Hypergeometric continuation of divergent perturbation series. II. Comparison with Shanks transformation and Pad\'e approximation

We explore in detail how analytic continuation of divergent perturbation series by generalized hypergeometric functions is achieved in practice. Using the example of strong-coupling perturbation series provided by the two-dimensional Bose-Hubbard model, we compare hypergeometric continuation to Shanks and Pad\'e techniques, and demonstrate that the former yields a powerful, efficient and reliable alternative for computing the phase diagram of the Mott insulator-to-superfluid transition. In contrast to Shanks transformations and Pad\'e approximations, hypergeometric continuation also allows us to determine the exponents which characterize the divergence of correlation functions at the transition points. Therefore, hypergeometric continuation constitutes a promising tool for the study of quantum phase transitions.


Introduction
Quantum mechanical perturbation theory usually is applied to a Hamiltonian of the form where the "unperturbed system" H 0 can be diagonalized exactly, and the "perturbation" V is supposed to be small in some suitable sense [1]. The dimensionless parameter λ, which is set equal to one at the end of the calculation, connects the perturbed system H to the unperturbed one when varying between zero and unity. Employing the customary Rayleigh-Schrödinger perturbation series, one may then compute, e.g., the perturbed energy eigenvalues and eigenstates as a formal power series in λ [2][3][4]. In practice, the evaluation of the perturbation series may pose insurmountable technical difficulties in higher orders, so that one often restricts oneself to the lowest few terms, hoping that the series converges sufficiently fast for such a truncation to be meaningful. If, however, the perturbation series has only a finite radius of convergence, the formal series may still bear significance even beyond that radius. As an illustrative example, consider the one-dimensional linear harmonic oscillator [4] H 0 = p 2 2m where q is the position operator, p its conjugate momentum operator, m denotes the mass of the oscillator particle, and ω 0 is the (positive) oscillator angular frequency, so that the unperturbed energy eigenvalues are given by E n (λ = 0) = ω 0 (n + 1/2) with quantum numbers n = 0, 1, 2, . . . . Suppose further that this system is perturbed by another oscillator potential with (positive) frequency ω 1 , Evidently, the exact perturbed energy eigenvalues then are given by E n (λ = 1) = Ω(n + 1/2) with Ω = ω 2 0 + ω 2 1 .
If, on the other hand, one works out the perturbation series to fourth order, one finds E n (λ = 1) = ω 0 (n + 1/2) f ((ω 1 /ω 0 ) 2 ) , where the function f is given by Now one faces two closely related tasks. Firstly, one needs to deduce that this function f , completed to all orders of x, takes the form This series converges for |x| ≤ 1, but becomes a diverging, asymptotic series for |x| > 1 [5,6]. Therefore, the full perturbation series E n (λ = 1) = ω 0 (n + 1/2) ∞ ν=0 1/2 ν ω 1 ω 0 2ν (10) converges only for ω 1 ≤ ω 0 . Secondly, bearing in mind that one has to realize that the expression on the right-hand side of this equation also constitutes the analytic continuation of the left-hand side for x > 1. Therefore, for both ω 1 ≤ ω 0 and ω 1 > ω 0 the summed perturbation series yields E n (λ = 1) = ω 0 (n + 1/2) 1 + which, of course, equals the above expression (5). This pedagogical example sets the stage for the current work. The plan of the present paper is to subject a recently suggested powerful technique for the analytic continuation of divergent perturbation series based on the use of generalized hypergeometric functions [7][8][9] -dubbed hypergeometric continuation for short -to a comprehensive test, thereby demonstrating its outstanding value for practical calculations. Instead of aiming for general theorems, here we consider a definite system of particular significance, the two-dimensional Bose-Hubbard model. This model, which describes interacting Bose particles on an infinite two-dimensional lattice, shows a quantum phase transition from a Mott insulator to a superfluid at zero temperature when the relative strength of the interparticle interaction is reduced, or the strength of the tunneling contact between neighboring lattice sites is increased [10,11]. This Mott insulator-tosuperfluid transition reflects itself in a divergence of the strong-coupling perturbation series of certain correlation functions of the Bose-Hubbard model [9]. Therefore, in order to obtain information on system properties in the superfluid phase, these series need to be analytically continued beyond the point of divergence into the superfluid regime. Hence, we again face the two tasks exemplified by equations (9) and (11): Starting from a series of which only the first few terms are at our disposal, we need to guess the underlying systematics to all orders, and to deduce the corresponding analytic continuation. In the Bose-Hubbard case, however, the leading contributions to the perturbation series are given only numerically, and their extension to infinite orders is not obvious. For tackling these tasks, we proceed as follows: In Sec. 2 we briefly review the formulation of the Bose-Hubbard model, and give a formal definition of its k-particle correlation functions c 2k in terms of strong-coupling perturbation series. In Sec. 3 we then explain the basics of their analytic continuation by means of the familiar Shanks transformation and Padé approximation methods [12][13][14][15], and through the novel hypergeometric technique. All three schemes are applied in Sec. 4 to the one-particle correlation function c 2 of the Bose-Hubbard model, in order to deduce its phase diagram.
Since there exist fairly precise previous computations of this diagram, comparison of such reference data with our results allows us to gauge the accuracy of the respective method, and thereby to confirm that hypergeometric continuation is at least competitive with its well-established rivals. In Sec. 5 we then address a subject which highlights a particular strength of hypergeometric continuation, and which is not amenable to any of the other two methods: Using hypergeometric continuation, we determine the exponents with which the one-particle correlation function c 2 and the two-particle correlation function c 4 diverge at the phase boundary. Finally, some conclusions are drawn in Sec. 6. Thus, the present case study serves a twofold purpose. On a fairly general level, we wish to establish hypergeometric continuation as a reliable and highly versatile tool for extracting observable quantities from divergent perturbation series, corroborating the pioneering work by Mera et al [7,8]. More specifically, the results obtained in Sec. 5 also serve as input data for our preceding paper [16], in which we have determined the critical exponent of the order parameter of the Mott transition on a two-dimensional lattice. While the present investigation is essentially independent of that preceding paper [16], it provides the technical background information required for assuring the correctness of the data employed therein.

Perturbative evaluation of correlation functions
Consider a d-dimensional lattice of arbitrary geometry with sites labeled by an index i, and let b † i and b i be the Fock space operators which create and annihilate, respectively, a Bose particle at the ith site, obeying the commutation relation [b i , b † j ] = δ ij , so that n i = b † i b i is the operator which counts the number of particles placed on site i. Assume further that neighboring sites are connected by a tunneling contact with hopping matrix element J, and that particles occupying a common site repel each other, with each pair of particles contributing an amount U to the total potential energy, while interaction between particles sitting on different sites is neglected, as sketched in Fig. 1. This system is subjected to a chemical potential µ, which enables one to control its total particle content. After scaling with respect to U , the dimensionless Hamiltonian of this Bose-Hubbard model is cast into the form where the first part is site-diagonal, accounting for the total repulsion energy and the coupling to the chemical potential, while the second term incorporates the tunneling links between neighboring sites i and j, as indicated by the symbol i, j . This fairly minimalistic model provides an excellent description of ultracold atoms in deep optical lattice potentials [17,18], since the effective range of the two-particle scattering potential for neutral atoms is much shorter than the optical lattice constant, as given by half the wavelength of the lattice-generating laser radiation [19]. Assuming integer filling with g particles per site at zero temperature, the ground state of the system in the absence of the tunneling contact, that is, for J/U = 0, is given by the product state with µ/U < g < µ/U + 1 [10], and |vac denoting the empty-lattice state. Now there are two alternative approaches to account for finite tunneling strength by means of perturbation theory in J/U . The first strategy is to compare the energy of the ground state which emerges from the product state (16) with the energy of two defect states which contain either an additional particle or an additional hole moving coherently over the lattice; the boundary between the incompressible Mott phase and the superfluid phase is reached when the energy difference between the Mott state and a defect state vanishes. This strategy has been adopted and implemented to third order already in Ref. [20]. The technical inconvenience of having to deal with the defect states is avoided if instead one probes the response of the basic system (13) to external sources and drains, thus breaking its particle number conservation. Taking these sources and drains to be spatially uniform with strength η, we therefore consider the extended system where and determine its intensive ground state energy E in the thermodynamic limit, where M is the number of lattice sites. Observing that the sought-for response then is quantified by for arbitrary i, where the expectation value again is taken with respect to the ground state of the extended system (17), the responseless Mott phase is characterized by b i = 0, whereas one finds b i = 0 in the superfluid phase, reflecting spontaneous symmetry breaking. Expanding the ground state energy (19) within the Mott phase, which has to be an even function of η, in the form the quantity e 0 (µ/U, J/U ) merely represents the intensive ground state energy of the basic model (13), whereas the desired information is contained in the k-particle correlation functions c 2k (µ/U, J/U ). When these correlation functions, in their turn, are expanded for given µ/U in powers of the scaled hopping strength J/U , writing each term in such a series corresponds to various chains of k creation processes, ν tunneling events, and k annihilation processes, as visualized in Fig. 2 for k = 1 and k = 2, respectively. We pay particular attention to these two correlation functions c 2 and c 4 : When varying J/U at fixed µ/U close to the respective transition point (J/U ) c , they exhibit power-law behaviors and with positive divergence exponents 2 (µ/U ) and 4 (µ/U ) which enable one, for d = 2, to determine the critical exponent of the order parameter [16]. From the viewpoint of systematic many-body perturbation theory, this approach means that the Hamiltonian (17) is split according to . Graphical representation of the lowest terms of the perturbation series (22) for the one-particle correlation function c 2 (above), and for the two-particle correlation function c 4 (below), with weight factors pertaining to a d-dimensional hypercubic lattice. Open circles symbolize a creation process, crosses an annihilation process, and arrows indicate a tunneling event which connects neighboring lattice sites.
where the "unperturbed system" H 0 is given by the site-diagonal part (14), and the "perturbation" consists of both the tunneling contacts (15) and the sources and drains (18). The perturbative evaluation of the series (22), starting from a ground state (16) with appropriate filling factor g for J/U = η = 0, then proceeds by means of the process chain approach as devised by Eckardt [21], which constitutes an adaption of Kato's non-recursive formulation of the general perturbation series [22] to many-body lattice models. This enables one to represent physical observables by appropriate sequences of processes as exemplified above in Fig. 2; the evaluation of such diagrams is sketched briefly in Appendix A. The nth order of perturbation theory thus comprises all connected diagrams which consist of n processes of any kind; the computational bottleneck is caused by the necessity to evaluate all processes of a given chain in all possible permutations. The details of this procedure have been communicated elsewhere [23] and do not need to concern us here; suffice it to state that we are able to compute the series (22) up to ν max,2 = 10 for k = 1 and up to ν max,4 = 7 for k = 2 (corresponding, respectively, to orders n = 12 and n = 11 of perturbation theory). Thus, in the present study we fully account for all process chains contributing to the one-particle correlation function c 2 with up to ten tunneling events, and for all chains contributing to c 4 with up to seven tunneling events.
As an example, Fig. 3 shows the coefficients −α  7,8,9). Observing the logarithmic scale of the ordinate, one deduces from the exponential growth of the coefficients that simply terminating the series (22) at the respective ν max,2k will not yield reasonable approximations. Hence, we need some kind of extrapolation scheme for estimating the coefficients α (ν) 2k (µ/U ) for ν > ν max,2k , and some means of analytic continuation for giving meaning to the series beyond their radius of convergence.

Analytic continuation
We start by inspecting two standard techniques for handling slowly converging or diverging series: the Shanks transformation and the Padé approximation, before turning to analytic continuation utilizing hypergeometric functions. We discuss the three methods using c 2 as example; c 4 can be treated analogously.

Shanks transformation
We first consider the Shanks transformation [12,13], which transforms a sequence (s n ) n∈N into the sequence (S n ) n∈N given by  The Shanks transformation is constructed in such a way that for s n = A + αq n with geometrically decaying transient (|q| < 1) the Shanks transform (S n ) n∈N is a constant sequence S n = A which equals the limit of the original sequence (s n ) n∈N . More generally, in many cases the Shanks transform has better convergence properties than the original sequence. Applied to the sequence of partial sums s n = n ν=0 b ν of some sequence (b ν ) ν∈N , the transformation formula (27) takes the form We now apply this general prescription to the truncation s n = n ν=0 α (ν) 2 (µ/U ) (J/U ) ν of the perturbation series for c 2 . Since we always monitor the phase transition such that µ/U is kept fixed while J/U is varied, we omit the argument µ/U in the following for the sake of clear notation. Thus we obtain The underlying hypothesis now is that the limit of this Shanks transform (29) is c 2 , the object we are interested in. Then approximants to the radius of convergence of c 2 are given by the respective zero of the denominator of the second term, so that we have the easy-to-calculate approximation for the phase boundary (J/U ) c = lim ν→∞ (J/U ) c,ν . Exemplarily we have displayed these ratios for the tip of the lowest Mott lobe in Fig. 4. The approximants (J/U ) c,ν seem to  converge towards a limit consistent with the corresponding phase boundary (J/U ) c,QMC obtained from Quantum Monte Carlo (QMC) calculations [24]. Since the limit is not reached in finite order, we have linearly extrapolated the inverse of these finite-order approximants over the reciprocal perturbation order 1/ν to the limit 1/ν → 0, as displayed in Fig. 5. This latter figure also shows the ratios α corresponding to the two-particle correlation function c 4 . Observe that the fit then yields a notably different numerical value for (J/U ) c , although we expect the same (J/U ) c for both c 2 and c 4 [16]. This indicates that the present method is not optimal.
Although the Shanks transformation thus allows us to compute the phase boundary, it does not provide access to the divergence exponents 2k , as introduced in the relations (23) and (24), because the transform (29) has a pole of order 1. Therefore, the Shanks transform yields only trivial, integer-valued estimates of the divergence exponents 2k , independent of the perturbation coefficients α (ν) 2k (µ/U ). We remark that in some cases it is possible to further increase the convergence properties of the Shanks transform S n by iterating the Shanks transformation. While this may lead to a better approximation of c 2 in general, its radius of convergence and consequently our estimate (J/U ) c,ν for the phase boundary are not affected by such an iteration. This is a consequence of the fact that any singularity of s n is as well a singularity of the Shanks transform S n , due to the numerator's quadratic and the denominator's linear dependence on s n , so that the Shanks transform of S n shares the singularity of S n itself at (J/U ) c,ν = α To support this theoretical discussion, we have sketched in Appendix A the analytical calculation of the first three coefficients α 2 , and α 2 of c 2 , in order to evaluate S 1 . In the limit of high dimensionality the previously derived formula (30) then reproduces the mean-field phase boundary, which is exact in this limit, thus serving as a showcase for the quality of the Shanks transformation.

Padé approximation
Next, we address the Padé approximation [13][14][15]. The underlying hypothesis here is that the coefficients α (ν) 2k of the perturbation series correspond to coefficients of the series expansion of a rational function Given the degrees L and M of the polynomials in the denominator and numerator of A L/M , the coefficients p 0 , p 1 , ..., p L and q 1 , ..., q M of the polynomials can be calculated by solving the system of equations The degrees L and M are restricted through the set of available coefficients α To obtain an approximation to the phase boundary, as given by the radius of convergence of c 2 , we observe that the rational function A L/M has up to M isolated poles of integer order, corresponding to the zeros of the denominator 1 + M k=1 q k · (J/U ) k . Thus, if we assume that c 2 is given by the Padé approximation A L/M , the smallest positive pole of A L/M corresponds to the phase boundary. Again, the integer order of the poles of A L/M implies that the Padé approach only yields integer estimates for the divergence exponents. Shanks transformation and Padé approximation are directly related to each other. It can be shown that the Shanks transform S n corresponds to the Padé approximation A n/1 , while Padé approximations A L/M with M ≤ L equal so-called generalized Shanks transforms [13].

Hypergeometric continuation
Summing up, we stress that neither Shanks transformation nor Padé approximation are able to yield non-trivial estimates for the divergence exponents 2k . The idea to overcome this deficiency is to replace the rational functions used in the Padé approximation by functions with an essential singularity. In a recent series of papers [7,8], Mera et al have suggested to employ hypergeometric functions for the analytic continuation of typical perturbation series in quantum mechanics. As we will discuss below, these hypergeometric functions possess a tunable singularity, which is exactly what is needed for our purposes. Moreover, we have already given a heuristic argument why hypergeometric functions are particularly well suited for the analytic continuation of strong coupling perturbation series for the Bose-Hubbard model [9]. For instance, we may assume that the coefficients c 2k are given by Gaussian hypergeometric functions with parameters a, b, c and (J/U ) c still to be determined, where (a) ν = a(a + 1) · · · (a + ν − 1) is the usual Pochhammer symbol [25]. Beyond this Gaussian hypergeometric function 2 F 1 with 4 degrees of freedom, we also consider generalized hypergeometric functions p F q a 1 , a 2 , ..., a p , b 1 , ..., b q , It is known that such generalized hypergeometric functions converge for all arguments if p < q + 1 and diverge for any value of J/U = 0 if p > q + 1 [26]. As a finite radius of convergence is precisely what we are interested in, we employ only generalized hypergeometric functions with p = q + 1, as natural generalizations of 2 F 1 . As defined in equation (35), these functions have (J/U ) c as their radius of convergence. In the following we will describe our approach for the Gaussian hypergeometric function 2 F 1 ; the adjustments necessary for generalized hypergeometric functions q+1 F q are straightforward. For calculating the parameters a, b, c and (J/U ) c we utilize all perturbatively available coefficients α (ν) 2k in a least-squares fit as displayed in Fig. 6. As the coefficients vary over at least thirteen orders of magnitude (cf. Fig. 3), an unweighted fit yields huge relative deviations for the smaller coefficients α (ν) 2k , and has a correspondingly bad performance. To overcome this deficiency we have experimented with various weights. We obtained good results with least-squares fits of the relative deviations, thus weighting the absolute deviations by the respective coefficient 1/α (ν) 2k . However, we achieved the best overall performance when fitting the ratios α of subsequent coefficients by corresponding ratios of terms of the hypergeometric function, These ratios are of the same order of magnitude, so that we do not have to introduce artificial weights when fitting the data. In Appendix B we have compiled a comparison between both procedures, employing the ratios (36) on the one hand, and relative deviations on the other, and show that both methods lead to more or less the same results.
Depending on the details of the fitting algorithm used, unwanted local minima of the sum of least squares can cause problems which worsen with increasing number of degrees of freedom. To deal with this we calculate every fit multiple times with randomized initial values, and choose the fit with the smallest squared deviations. In Fig. 6 we have displayed the ratios α of the coefficients near the tip of the first Mott lobe together with the fits we obtained for the Gaussian hypergeometric function 2 F 1 , and for the generalized hypergeometric function 1 F 0 , which is the well-known binomial series. We now turn to the evaluation of the divergence exponents 2 and 4 of the correlation functions (23) and (24). To this end we apply the linear transformation formula (see, e.g., [25]) which is valid for any z ∈ C with | arg(1 − z)| < π. The hypergeometric functions 2 F 1 on the right-hand side tend to 1 as their arguments 1 − z tend to 0. Therefore, the first term is asymptotically constant for z → 1. In contrast, for a + b > c the factor (1 − z) c−a−b diverges, and the asymptotics are given by 2 F 1 a, b; c; Therefore, within the scope of hypergeometric analytic continuation by means of 2 F 1 , the divergence exponent of c 2k at (J/U ) c is estimated by When employing the generalized hypergeometric functions q+1 F q as defined in equation (35), this estimate becomes In order to ensure the validity of the transformation formula (37) for arguments z = J/U (J/U ) c > 1 we have to add to J/U an arbitrarily small imaginary part. Moreover, beyond the radius of convergence of its series representation, the hypergeometric functions q+1 F q adopt complex values. If real-valued quantities are required, we may utilize the replacement beyond the phase transition. This construction resolves both issues. In this work, however, we will not study the correlation functions in the superfluid phase.

Phase diagram
We now report the results we have obtained for the phase diagram by evaluating the radius of convergence of the one-particle correlation function c 2 according to the three different schemes. This phase diagram has previously been calculated in various ways [9,20,24,[27][28][29] and, therefore, allows us to assess the quality of the respective scheme. While the process-chain approach easily allows us to calculate the phase diagram for arbitrarily high filling factors, this proves difficult with other methods, so that we focus our comparison on the first Mott lobe (g = 1).

Shanks transformation
When resorting to the Shanks transformation, we calculate the phase boundary by means of formula (30) and extrapolate the results linearly in 1/ν, as indicated in Fig. 5. We have depicted the highest numerically accessible approximation S 9 together with the extrapolation in Fig. 7. Our results agree qualitatively very well with the benchmark expression for the phase boundary stated by Freericks et al [29]. Generally, the finite-order Shanks transform S 9 seems to significantly underestimate the boundary (J/U ) c around the tip of the Mott lobe, while the 1/ν-extrapolation seems to slightly overestimate the boundary away from the tip. Quantitatively, we find the tip of the first Mott lobe at µ/U = 0.3740 compared to µ/U = 0.3724 as stated by Freericks et al , and there is a small deviation between both estimates of the phase boundary with a maximum difference of ∆(J/U ) c = 0.00105, which corresponds to a relative difference of 1.9 %, half-way between µ/U = 0 and the tip of the first Mott lobe. From the trend in the finite orders available to our perturbational treatment we expect that the phase boundary derived from the Shanks transform S ν by applying formula (30) converges for ν → ∞ to a phase boundary very similar to the one obtained by Freericks et al . This is supported by the very good agreement emphasized in the inset of Fig. 7.

Padé approximation
When employing the Padé approximation, we obtain the phase boundary by calculating the zeros of the denominator 1 + M k=1 q k · (J/U ) k appearing in equation (

Hypergeometric continuation
For each value of µ/U separately, we have fitted the ratios of subsequent coefficients α (ν) 2 (µ/U ) α (ν−1) 2 (µ/U ) to the generalization of expression (36) in order to obtain the parameters of the generalized hypergeometric functions q+1 F q . The phase boundary (J/U ) c is then simply given by one of the fit parameters. In this manner we obtain the phase diagram depicted in Fig. 9 for the first four Mott lobes, and magnified in Fig. 10 for the first Mott lobe only. Diagrams for 1 F 0 , 2 F 1 and 3 F 2 agree very well with each other, and with the benchmark expression. Exemplarily, we have calculated the largest difference between the phase boundary obtained by continuation based on 2 F 1 and that benchmark. This gives ∆(J/U ) c = 0.00151, corresponding to a relative deviation of 2.8 %. When inspecting the different phase diagrams obtained with an increasing number of degrees of freedom using generalized hypergeometric functions q+1 F q , the boundary (J/U ) c decreases and the difference between our findings and the benchmark becomes smaller. Similar to the Shanks transformation and the Padé approximation the finite-order estimates seem to converge. In view of the minor differences between the phase diagrams obtained by 3 F 2 and 4 F 3 , as shown in Fig. 10, it appears unlikely that the limit for q → ∞ agrees exactly with the benchmark; however, the relative deviation is less than 2.0 %. Especially at the tip of the first Mott lobe the results hardly change when going from 3 F 2 to 4 F 3 , with less than 0.25 % relative deviation, and to 5 F 4 . The hopping strength J U chemical potential µ U 1 F 0 continuation 2 F 1 continuation 3 F 2 continuation benchmark expression [29]  precise values of (J/U ) c at the tip provided by the various approximations are collected in Tab. 2. Thus, it appears that with 4 F 3 we are very close to the limit q → ∞, so that it is not necessary to calculate fits beyond 5 F 4 .

Comparison
Comparing the estimates for the phase boundary obtained from all three analytic continuation schemes, we conclude that they agree qualitatively very well with each other, yielding the well-known Mott lobes. As illustrated in Fig. 11, the finite-order  the computation of the phase boundary. This finding instills great confidence that hypergeometric analytic continuation also provides a reliable basis for the computation of divergence exponents.

Divergence exponents of correlation functions
The key advantage of hypergeometric continuation, as compared to the Shanks or Padé approach, lies in the fact that the hypergeometric functions themselves possess a nontrivial divergence exponent at the border of their convergence regime. Therefore, by fitting these functions to the perturbatively calculated coefficients α (ν) 2k , the divergence exponents 2k of the correlation functions c 2k are estimated directly from the parameters of the fit functions via equation (39) for Gaussian hypergeometric functions, or equation (40) for generalized hypergeometric functions. In this manner, the divergence exponents are obtained from the fits already presented in Sec. 4.3, without the need for additional work. In Fig. 12 we show the divergence exponent 2 for scaled chemical potentials ranging from µ/U = 0 to µ/U = 4, as resulting from fits with 1 F 0 , 2 F 1 , and 3 F 2 . Evidently, 3 F 2 . To obtain these fits each transition point (J/U ) c has been set to the value calculated before from c 2 , instead of treating it as an independent fit parameter. This is justified by the knowledge that c 2 and c 4 share the same radius of convergence [16], but can also be checked numerically: If this assumption is not made, but (J/U ) c is determined independently from c 2 and c 4 , the results agree to within less than 3 %.
To get an impression of how strongly our results depend on the parameter q of the fitting functions q+1 F q , Tab. 3 lists the divergence exponents at the tip of the first Mott lobe, as estimated for the q numerically accessible to us. While the results for 2 are quite stable, similar stability for 4 is achieved only with q = 1 and q = 2, whereas the value for q = 0 appears to be somewhat off, indicating that 1 F 0 does not offer a sufficient number of degrees of freedom. From the connection of the divergence exponents to the critical exponents of the Mott insulator-to-superfluid transition [16], we know that both 2 and 4 each adopt a common value at all lobe tips. This expectation is confirmed by Tab. 4, which has been compiled using 2 F 1 , on the sub-percent accuracy level.    Table 3. Comparison of the divergence exponents 2 and 4 at the tip of the first Mott lobe, as obtained by hypergeometric analytic continuation with various q+1 F q .

Conclusions
The present study serves both a technical and a conceptual purpose. On the technical level, we have employed the strong-coupling perturbation series for the correlation functions (22) of the two-dimensional Bose-Hubbard model at zero temperature for demonstrating how analytic continuation of divergent perturbation series by means of hypergeometric functions [7][8][9] is achieved in practice. To this end, we have evaluated the strong-coupling series numerically to high orders, and fitted their coefficients to Mott lobe g generalized hypergeometric functions q+1 F q , providing 2q + 2 degrees of freedom. This necessarily enforces a compromise between the quest for high flexibility, as increasing with q, and the task to provide a correspondingly large number of data points, i.e., to compute high-order coefficients of the perturbation series. Fortunately, it turns out that already the Gaussian hypergeometric function 2 F 1 with its four degrees of freedom is suitable for most purposes. By comparing the phase diagram of the Mott insulator-to-superfluid transition computed in this manner with the results provided by the well known Shanks and Padé continuation techniques, we have confirmed that hypergeometric continuation constitutes an accurate and reliable tool. On the conceptual level, it is of profound theoretical interest to explore just how the strong-coupling perturbation series diverge at the phase boundary. Namely, the Mott insulator-to-superfluid transition shown by the two-dimensional Bose-Hubbard model generally is of the mean field type, with the exception of the tips of the Mott lobes, which constitute multicritical points [10]. Thus, the question arises how the switch from "mean field-like" to "multicritical" in response to a variation of the chemical potential affects the perturbation series. Hypergeometric continuation is ideally suited to address this question in detail, since generalized hypergeometric functions provide a tunable singularity which can be employed to determine the exponents which characterize the divergence of the correlation functions at the phase boundary. Our numerical results for the divergence exponents 2 and 4 of the one-and two-particle correlation function, respectively, furnish the starting point of a novel strategy for determining critical exponents [16]. Thus, the idea to utilize hypergeometric functions for analytic continuation of divergent perturbation series, originally put forward by Mera et al [7], not only provides an efficient alternative to the tools already at hand, such as the Shanks and Padé techniques, but it also provides insights not obtainable with these older methods. Similar to the present study, the solution of a quantum mechanical single-particle or many-body problem can often be separated into two different steps: (i) The evaluation of a systemspecific perturbation series, and (ii) the analytic continuation of that series. While step (i) may strongly depend on the physical nature of the respective system under study, hypergeometric continuation offers a powerful universal tool to perform step (ii).
of the perturbation series for c 2 , and utilize the Shanks transformation to obtain an estimate for the phase boundary. As detailed in Sec. 2, we have to evaluate the correlation function  2 (µ/U ), which corresponds to the energy shift proportional to η 2 and (J/U ) 0 , is interpreted to originate from a process chain creating a Bose particle at an arbitrary site i, and annihilating it at the same site. According to the Kato perturbation theory we have to account for both permutations: Starting from the ground state (16), either a particle is first created and then a particle is annihilated. In this case the sum over all intermediate states |m reduces to a single term, involving a state |. . . , g, g + 1, g, . . . with one additional particle at site i, giving . . , g, site i g + 1, g, . . . . . . , g, Or a particle is annihilated first and then a particle is created, . . , g, Analogously, we compute the first-order contribution: With i and j labeling neighboring sites, this requires the calculation of together with the five other permutations of the processes. As there are 2d directions for the tunneling process from a fixed site i to a neighboring site j on a d-dimensional hypercubic lattice, we have to multiply their sum by 2d, and obtain α (1) To second order in J/U we have to consider the two diagrams depicted in Fig. 2: Either the second tunneling process ends at the origin of the first one or it goes to a different lattice site. For both diagrams we have to count the number of corresponding paths on our d-dimensional hypercubic lattice: While the first diagram has a multiplicity of 2d, the multiplicity for the second diagram is 2d · (2d − 1). For general spatial dimension d both diagrams have to be taken into account, and we have done so in our numerical evaluation for the two-dimensional square lattice. Here, however, we focus on the limit d → ∞, which allows us to neglect the first diagram. Evaluating and summing over all 4! = 24 permutations, we find α 2 (µ/U ) = −2d(2d − 1) (µ/U + 1) 3 With these first three coefficients (A.4), (A.6) and (A.7) we are now equipped to calculate an approximation for the phase boundary by means of the Shanks transformation. According to equation (30) this approximation is given by 2 (µ/U ) α divergence exponent of c 2 chemical potential µ U relative deviations fit ratio fit Figure B1. (Color online) Zero-temperature phase diagram of the 2d Bose-Hubbard model (left) and divergence exponents 2 for the Mott insulator-to-superfluid transition (right), as obtained from the hypergeometric functions 2 F 1 with a least-squares fit minimizing the relative deviations, and with a least-squares fit to the ratios α 2,ratio = 1.392 and 2,relative = 1.372 deviate from each other by less than 1.4 %. This remarkably good agreement achieved with the two different fitting strategies further increases our confidence in the viability of the hypergeometric function approach for analytic continuation: Regardless of the particular fitting strategy one obtains stable estimates both for the phase boundary and even the divergence exponent.