From weak to strong: constrained extrapolation of perturbation series with applications to dilute Fermi systems

We develop a method that uses truncation-order-dependent re-expansions constrained by generic strong-coupling information to extrapolate perturbation series to the nonperturbative regime. The method is first benchmarked against a zero-dimensional model field theory and then applied to the dilute Fermi gas in one and three dimensions. Overall, our method significantly outperforms Pad\'e and Borel extrapolations in these examples. The results for the ground-state energy of the three-dimensional Fermi gas are robust with respect to changes of the form of the re-expansion and compare well with quantum Monte Carlo simulations throughout the BCS regime and beyond.


Introduction.
A common situation in physics is that properties of a system can be computed analytically in a weakcoupling expansion, but only numerically at discrete points in the nonperturbative regime. The constrained extrapolation problem is to construct approximants that combine these two sources of information.
Consider an observable F(x) defined relative to the noninteracting system, e.g., the ground-state energy E/E 0 . Its perturbation series (denoted PT), truncated at order N in the coupling x, reads While Eq.
(1) provides precise information about the behavior of F(x) as x → 0, it generally fails to yield viable approximations away from weak coupling. Indeed, the PT is often a divergent asymptotic series, with large-order coefficients obeying, e.g., c k k→∞ ∼ k! [1,2]. Experiment or computational methods can give access to the behavior of F(x) at a specific point x 0 . Since x 0 can be mapped to infinity by a conformal transformation, we may take x 0 = −∞. Weak-to-strong-coupling extrapolants can then be defined as classes of functions F N (x) that reproduce both the PT to order N and the strong-coupling limit F(−∞) = ξ. F N (x) may also incorporate available information on the leading coefficient(s) d k in the strong-coupling expansion (SCE): The goal is then to find extrapolants F N (x) that converge rapidly and smoothly to the correct F(x) as N → ∞. As often only a few PT coefficients are known, from a practical perspective, F N (x) should be well converged at low orders.
A textbook example that has been the focus of many experimental and theoretical studies in the past two decades is the dilute Fermi gas [3,4]. Here x = k F a s , where k F is the Fermi momentum and a s is the s-wave scattering length. Due to its universal properties, the three-dimensional (3D) Fermi gas serves as an important benchmark for neutron matter [5] and neutron stars [6]. Its ground-state energy F = E/E 0 has been studied from weak attractive coupling through the BCS-BEC crossover with ultracold atoms [7] and sophisticated quantum Monte Carlo (QMC) simulations [8][9][10]. The Bertsch parameter ξ has been determined experimentally as ξ = 0.376(4) [7] and from QMC as ξ = 0.372(5) [9]. Moreover, recently the weak-coupling expansion has been calculated to N = 4 [11]. On the SCE side, one has viable estimates for d 1 and d 2 only, with d 1 known more precisely [12]. Such a situation is typical when only limited data are available in the nonperturbative region.
Padé approximants [13,14] are a standard approach to the extrapolation problem. However, when applied to the dilute Fermi gas, several Padés give flawed approximants with poles in the BCS region, see Fig. 1(a). Therefore, in this Letter we develop a new extrapolation method. Our method improves on the order-dependent mapping (ODM) approach introduced by Seznec and Zinn-Justin [15] (see also Ref. [16]) and builds in information on the leading strong-coupling coefficients, d 1 and d 2 . We refer to our method as order-dependent-mapping extrapolation (ODME).
For benchmarks we consider, in addition to the dilute Fermi gas in 3D, also its 1D variant (both for spin 1/2), and a 0D model problem that has long been a proving ground for extrapolation methods. We first summarize the PT and SCE information available for these problems, and then briefly discuss Padé extrapolants. This is followed by an explanation of the ODME and how it incorporates strong-coupling constraints. We then test the ODME for the 0D model and the 1D Fermi gas and find it outperforms Padés and produces very accurate approximants already at low orders. As a preview, this is shown in Figs. 1(b) and (c) and in Fig. 2. Our main results are for the 3D Fermi gas, where we find that the ODME, in contrast to Padés, leads to well-converged extrapolants that predict results consistent with QMC within uncertainties, see Fig. 1(a) and Fig. 3. The Supplemental Material (SM) provides additional details and shows that ODME also outperforms various Borel extrapolants. 3D Fermi gas exact ODME [4] ODME [2] PT [4] PT [2] 2pt-Padé [3,3] 2pt-Padé [ [10]. The errors of the QMC data include, in addition to the statistical uncertainty [10], an uncertainty based on QMC systematics [17,18]. The latter is taken to be ∆F QMC (x) = q(1 − F QMC (x)), with q = 0.038 obtained from the difference between ξ = 0.390 from Ref. [10] and the updated value ξ = 0.372(5) [9].
0D model. A well-known benchmark for resummation methods is the 0D field theory model [2,15,[19][20][21][22][23][24] Its weak-coupling coefficients are given by and the SCE involves fractional powers of g, The SCE of the 0D model has infinite radius of convergence, but Fig. 2 [and Fig. 1(c)] shows that for low truncation orders it is not very accurate, even for comparatively large values of g.
3D Fermi gas. For the 3D case, the PT coefficients have recently been calculated to fourth order [11]: Note that for spins higher than 1/2, logarithmic terms (and three-body parameters) appear in the PT [11]. On the SCE side, from QMC one finds d 1 ≈ −0.9 and d 2 ≈ −0.8 [12]. Comparing the gap between PT and SCE in Fig. 1(a) and (b) suggests that the PT and the SCE constrain F(x) less in the 3D case than in 1D. Padé approximants. The strong-coupling limit of the 3D and 1D Fermi gas can be described by "diagonal Padés" only (see also Refs. [11,33]), which are given by "Two-point Padés" are constructed by matching a k and b k to both the PT and the SCE up to specified orders. The twopoint Padé results shown in Fig. 1 are obtained by matching to ξ and d 1 and N = 2n − 2 PT coefficients. This is equivalent to matching a Padé[n − 1,n] to a rescaled version of F(x) that approaches 0 as x → −∞: In the 0D case we use square-roots of Padé[n,n + 1] functions such that successive orders in the SCE (g −1/4 , g −3/4 , etc.) are correctly reproduced [22,23].  A problem with Padé approximants is that they can have spurious poles in the region of interest [13,14]. In the case of the 3D Fermi gas, the two-point Padé [2,2] provides good results in the BCS region, but not beyond. However, matching a two-point Padé [3,3] to either (ξ, d 1 , c 1,2,3,4 ) or (ξ, d 1 , d 2 , c 1,2,3 ) produces poles at negative real coupling and hence a flawed extrapolant, see Fig. 1(a). In the 1D case many two-point Padés give good results at negative x (although some higherorder ones have poles there). But their continuation beyond the strong-coupling limit produces spurious poles at a small positive value of 1/x, see Fig. 1 Method of order-dependent mappings (ODM). We now consider extrapolants for f (x) of Eq. (9) of the form Here, w(x) satisfies w(0) = 0 and is chosen such that no poles occur on the negative real axis. Further, by Eq. (9), w(−∞) = 1 and the large-x behavior of w is taken consistent with the (rescaled) SCE (2). The coefficients h k are chosen to reproduce the first N terms in PT. For this, we rescale the PT (1) to Eq. (9), multiply by 1/(1−w), substitute x = x(w), and expand in powers of w to determine: . An approximant f N (x) is specified through the mapping w(x), which can contain control parameters {α i }. In the following, a single parameter α is used.
For the 0D model, approximants Z N (g) consistent with the SCE (5) are obtained if we use the mapping [15] and construct A possible choice for the mapping in Eq. (10) that builds in the SCE for the 1D and 3D Fermi gas, Eq. (2), is The method is called ODM because the parameter α is a function of N, i.e., it will be adjusted at each truncation order according to some criterion. This is similar to perturbation theory with an order-dependent reference point [15,[34][35][36][37]. The values of α(N) can be complex, in which case the ODM approximant is defined as the real part of f N (x) [or Z N (g)], see also Refs. [38][39][40].
Constrained extrapolation with order-dependent mappings (ODME). In the original ODM by Seznec and Zinn-Justin [15] the parameter α(N) is fixed by requiring that h N = 0, corresponding to the notion of "fastest apparent convergence" (FAC) [41]. This yields several possibilities for α(N), of which one is selected according to an additional criterion, e.g., smallest h N−1 . For the 0D model this approach converges to the exact solution as N → ∞ (the imaginary part of Z N (g) converges to zero for g ∈ [0, ∞]); the FAC criterion is not crucial for this, but suitable N dependence of α is [15,19,24,36,42,43].
In the ODME we fix α(N) by ensuring the SCE of f N (x) has a first-order coefficient equal to d 1 . This again yields several possibilities for α(N); we select the one that minimizes the difference between d 2 and the corresponding coefficient in f N (x). (A different approach to include SCE information in an order-dependent re-expansion was proposed in Ref. [44].) The set of approximants { f N (x)} corresponding to the uncertainty in the input (ξ, d 1 , d 2 ) and different choices for w(x) is then assessed according to the convergence behavior of f N (x), see below.
0D and 1D benchmarks. Figure 1(c) shows that in the 0D case, the ODME leads to high-precision approximants, already at low truncation orders. These significantly outperform two-point Padés. The ODME precision at higher N is examined in Fig. 2. While for large g and N the OMDE converges less rapidly and smoothly than the SCE, at low orders it is The second panel is the mapping used in Fig. 1(a). The bands for given N result from the uncertainties in the SCE coefficients d 1 = −0.90(5) and d 2 = −0.8(1) used to constrain the mappings.
more accurate even at relatively large coupling. At g = 10 ODME outperforms the SCE for N 4 (similarly at g = 100, see SM).
For the 1D Fermi gas, Figs. 1(b) and 2 show that also there ODME [with the mapping (14)] produces excellent approximants, again already at low orders. The ODME convergence is less pronounced compared to 0D, but Fig. 2 shows this can be improved by using mappings other than Eq. (14). More details are given in the SM.
Altogether, the study of the 0D model and the 1D Fermi gas suggests that ODME can produce accurate approximants already at low N, and is more broadly applicable than Padé (and Borel) extrapolation.
Next, we explore how robust the ODME predictions are with respect to the choice of mapping and what the sensitivity of ODME extrapolants is to the values chosen for (d 1 , d 2 ). To address this, we have investigated the class of mappings of the form w( consistent with the SCE of the 1D and 3D Fermi gas is, e.g., D(x; α) = κ 1 α − κ 2 x + (κ 3 α µ + (−x) ν ) 1/ν . In principle large sums of such terms are permitted. However, we find that to have well-converged results at low N, excessively complicated forms of D(x; α) are disfavored.
The ODME results for the 3D Fermi gas for different A more sophisticated algorithm would be to select sequences of ODME approximants according to their convergence for each (ξ, d 1 , d 2 ) input value. Fully implementing this requires uncertainties for (ξ, d 1 , d 2 ) based on a detailed analysis of the QMC data around strong coupling, in particular regarding error correlation. This requires further QMC input and is left to future work.
The ODME approximant sequences F N (x) with good convergence behavior produce results that are fully consistent with the QMC data, but generally lie below the central (variational) QMC values. This is most pronounced for the third mapping in Fig. 3 Assessing the variability in the ODME result over the four mappings of Fig. 3 and from uncertainty in the input (ξ, d 1 , d 2 ) specified above, we predict (N → ∞ extrapolated via Shanks transformation) for example: F ODME (−2) = 0.664(7) at x = −2.
Conclusions and future directions. We have developed the ODME method to provide powerful weak-to-strong coupling extrapolants constrained by limited data on strong-coupling behavior. For a 0D model and the 1D Fermi gas the ODME produces very accurate approximants already for low PT truncation orders. We then focused on the dilute Fermi gas in 3D. For this universal many-body system the weak-coupling PT is known to fourth order and limited strong-coupling data is available from experiment and QMC computations. With this input, the ODME yields robust approximant sequences with good convergence properties (in contrast to Padés). The predicted ground-state energies agree very well with the available QMC data over the entire range of intermediate couplings and even into the BEC side.
It is important to understand for which conditions the ODME works so well, especially in regard to the smoothness of the behavior from weak to strong coupling. This question is related to the way that nonperturbative features are encoded in weak-coupling asymptotics, and how different resummation methods capture them. For Padé [13,14] and Borel methods [2,[45][46][47][48][49] general results regarding such issues have been obtained. The comparison of our work to these methods shows that the ODME often performs better. It will be interesting to see if the ODME can be studied in similar mathematical generality as these standard resummation methods. In particular, the ODME can likely be improved by incorporating further analyticity constraints for the specific system of interest.
The ODME method is very flexible and broadly applicable. Interesting future applications include the unitary Fermi gas at finite temperature, where the virial expansion has recently been extended to fifth order [50], as well as hot and dense QCD matter from strong to weak coupling. Accurate methods that connect these regimes will also enable further progress for the nuclear equation of state in astrophysical simulations.
We thank R. J. Furnstahl for useful discussions and a detailed reading of the manuscript. We are also grateful to S. Gandolfi and A. Gezerlis for discussions regarding this work. DRP is grateful for the warm hospitality of the IKP Theoriezentrum Darmstadt. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project-ID 279384907 -SFB 1245, the US Department of Energy (contract DE-FG02-93ER40756), and by the ExtreMe Matter Institute.

SUPPLEMENTAL MATERIAL
Here, we provide more details on the performance of the ODME. We first compare ODME results against those from various Borel methods for the 0D model as well as the 1D (and 3D) Fermi gas. We then provide a more detailed assessment of the sensitivity of the ODME for the 3D and 1D Fermi gas to the choice of mapping and the input (ξ, d 1 , d 2 ).

Brief review of Borel extrapolation methods
Borel extrapolation is based on the Borel(-Leroy) transformed perturbation series which is constructed to have a finite convergence radius. That is, the large-order behavior c k k→∞ ∼ a k Γ(k + 1 + β), together with the choice of β 0 , determines the nature of the leading singularity of B(t) at t = 1/a [45,48,49]. From a given approximant B N (t) for B(t), constructed from the truncated-atorder-N Borel transformed perturbation series (see below), the corresponding approximant B N (x) for F(x) is obtained via the inverse Borel transform: If B N (tx) has poles on the positive real axis one can shift the integration path infinitesimally off the real axis. In this case the approximant for F(x) may be taken as the real part of B N (x); when applied to the (analytically continued) exact B(tx) this prescription often gives the correct result [2,15,30,46]. There are several methods to construct approximants B N (t) from incomplete perturbative information. The most straightforward is Padé-Borel, i.e., matching Padé approximants to the (truncated version of the) Borel series (A.1). If one has knowledge of the large-order behavior (specifically, if a is known), more sophisticated methods are available. In the "conformal-Borel" [2,48,49] approach one constructs B N (t) by re-expanding the (truncated) Borel series in terms of the conformal mapping that maps the cut Borel t plane to the interior of the unit disc [2,49]; i.e., Further, in the "Padé-conformal-Borel" [48] method one uses for B N (t) Padé approximants matched to Eq. (A.4). The Borel approximants discussed so far are "pure extrapolants" in that they include no strong-coupling constraints. One can also construct "constrained Borel extrapolants" where the strong-coupling limit F(−∞) = ξ is incorporated. "Constrained-conformal-Borel" ("CCB") approximants for f (x) = F(x)−ξ 1−ξ are obtained by re-expanding the (truncated) Borel series of f (x) as [45] and choosing η such that the known analytic structure at infinity is best reproduced. We choose η = 1/2 for the 0D model, and η = 1 otherwise. Finally, "Padé-constrained-conformal-Borel" ("PCCB") corresponds to matching Padé approximants to Eq. (A.5). The implementation of further SCE constraints is less straightfoward, and not considered here. A study of this problem can be found in Ref. [23], where it was found that two-point Padé-Borel approximants do not improve upon two-point Padés. Standard Borel resummation corresponds to β 0 = 0 in Eq. (A.1). The conformal transformation (A.3) yields a function that has a square-root branch point at t = 1/a. Based on this, a refinement of conformal Borel approximants corresponds to setting β 0 = β + 3/2, since then the exact Borel transform has the same feature [45]. With this, for the 0D model where a = −4 and β = 0, CCB (and PCCB) gives for all N 1 (resp. N 2) the exact result with K 1/4 (x) a modified Bessel function. Equation (A.6) matches the exact Z(g) given by Eq. (3) for g ∈ [0, ∞] and also provides its complex analytic continuation. However, apart from the case of (P)CCB extrapolants for the 0D model, we found that using β 0 = β + 3/2 does not yield much improvement compared to the standard choice β 0 = 0. Therefore the Borel results shown in Figs. A.1 and A.2 are obtained using β 0 = 0. g ∈ [0, ∞], provided the mapping parameter α scales appropriately with N: α(N) N→∞ ∼ 1/N γ , with 1 ( ≤ ) γ < 2 [36]. Indeed, Z N (g) then converges to the complex analytic continuation of Z(g), Eq. (A.6) [36], see also Ref. [20]. In the original ODM method [15], this is implemented by fixing α(N) via the FAC ("fastest apparent convergence") criterion h N = 0 (see also Refs. [19,24]). Another heuristic prescription is the "principle of minimal sensitivity" (PMS) [36,41], meaning that α(N) should be chosen such that Z N (g; α) is least sensitive to variations of α about its chosen value (see also Refs. [16,20,[38][39][40][42][43][44]).
Clearly, the optimal choice of α(N) is that which gives the most accurate results, with a smoothly converging sequence of approximants Z N (g). In Fig. A.1 we show that our ODME method-which fixes α(N) by matching to the SCE coefficients (d 1 , d 2 ), see main text-produces better approximants than the FAC criterion. We tried other prescriptions (e.g., PMS), which were similarly outperformed by ODME. Of course, this is not really surprising: ODME includes more information about the exact Z(g) than FAC and PMS.
In Fig. A.1, we also compare ODME against the various Borel extrapolants discussed above (using β 0 = 0). For Padé-Borel, Padé-conformal-Borel and PCCB we use Padé[n, m] functions with n = m − 1 = (N − 1)/2 and n = m = N/2, respectively, for odd and even truncation orders N. The Borel extrapolants all perform better than simple (i.e., non-Borel) one-point Padés, and exceptionally good results are obtained from the PCCB method. (By comparison, Padé-conformal-Borel does not improve much upon conformal-Borel.) For small coupling g 1 several Borel extrapolants are more accurate than ODME, but for g 1 ODME is outperformed only by PCCB (and the SCE) at large orders. However, for low PT truncation orders N 6 the ODME method gives the best approximants. This ability to produce accurate approximants at low N is a crucial asset for applications in realistic problems.

Results for the 1D (and 3D) Fermi gas
The exact ground-state energy density E(x) of the 1D Fermi gas can be computed with the Bethe ansatz, i.e., by solving numerically a Fredholm integral equation of the second kind [28,30,32]. From this, we compute the errors | f N (x) − f (x)| of approximants f N (x) to the exact solution for the rescaled energy f (x) given by Eq. (9).
Our results are shown in Fig. A.2. For ODME and FAC approximants we use the mapping (14). One sees that again the ODME leads to much better approximants than FAC, even for smaller x where one might expect that additional strongcoupling information does not improve the accuracy. For the conformal Borel extrapolants we use the recently determined large-order behavior, a = 1/π 2 (and β = −2) [29,30]. At small couplings |x| 1 the various Borel extrapolants are very precise, but their accuracy decreases with increasing coupling strengths; for |x| 10 they fail badly. (The Borel extrapolants with the correct strong-coupling limit often have local extrema at large x. Note also that here the conformal mapping technique does not improve upon Padé-Borel.) We have applied the various Borel extrapolants also to the 3D Fermi gas. [For the conformal Borel methods we have used, e.g., the conjectured large-order behavior a = −1/π (and β = 0) [30].] The results are similar to the 1D case: while accurate results are obtained for |x| 1, for larger couplings the various Borel extrapolants disperse strongly.
In summary, compared to the 0D model the Fermi gas in 1D (and, even more so, in 3D) represents a more difficult extrapolation problem. Nevertheless, although there the ODME is not as precise as in the 0D case (and the decrease of the errors with increasing N is diminished), it gives accurate results in the whole range x ∈ [0, −∞], in contrast to Borel methods. In addition, the ODME also reliably extrapolates the 1D (and 3D) Fermi gas to positive x, see Fig. 1.

Sensitivity to SCE input and mapping choice
Here, we study in more detail for the 3D and 1D Fermi gas the class of one-parameter mappings w(x) If the inverse mapping x(w) is not available in closed form, the coefficients γ n,m in Eq. (11) can be calculated iteratively starting from The iterations can be formulated in terms of polylogarithms, i.e., starting fromγ The γ n,m (x) are then obtained from theγ n,m (x) by substituting (A.10) In Fig. A.3 we compare the ODME results for the 1D Fermi gas for different D(x; α). One sees that many mappings perform better than our initial choice D(x; α) = α + (α 2 + x 2 ) 1/2 [see Eq. (14)]. The overall trend of the results is however similar for all D(x; α), i.e., the increase in precision with increasing N diminishes at larger couplings. We note that, while for x ∈ [0, −∞] some two-point Padés are more accurate than ODME with the mapping (14), the better mappings of In Fig. A.4, the mappings are ordered according to the convergence of the BCS results with increasing N, i.e., from smallest to largest deviations.
A weighted average of the deviation |F N (x) − F N−1 (x)| over orders N ∈ {2, 3, 4} together with an average over x ∈ [0, ∞] and the input (d 1 , d 2 ) is used for this purpose; see also the main text. While our focus here is on the BCS region, note that the ODME predictions for the BEC region may be improved by extending the convergence analysis to include values 1/x > 0.
The obtained ordering depends to some degree on the values of (d 1 , d 2 ) as well as the precise form of the quantitative convergence criterion. Qualitatively, the ordering in Fig. A.4 is as follows. For the two best-converged mappings, D(x; α) = α + (α + x 2 ) 1/2 and D(x; α) = α + (α 2 + x 2 ) 1/2 , the N = 2, 3, 4 results are very similar. For the third mapping, D(x; α) = 2α − x + (α + x 2 ) 1/2 , the deviations F N (x) − F N−1 decrease monotonically. The results of the fourth mapping are very similar for N = 3 and N = 4. The fifth, sixth, and seventh mappings appear about as well-converged as the third or fourth. On the other hand, the eighth mapping D(x; α) = α − x clearly has worse convergence behavior, see in particular the change from N = 3 to N = 4 in the plot with d 1 = −0.90(1) (last panel in Fig. A.5).
The sensitivity to mapping choice is more pronounced in the 3D case than in 1D. This reflects the fact the 1D extrapolation problem is more strongly constrained by the PT and SCE. The convergence behavior of different mappings deviates from the 1D case also in terms of which mappings perform better. In particular, for the 3D Fermi gas the simple mapping with D(x; α) = α − x gives approximant sequences with unfavorable convergence properties. (Note that this mapping has also the most irregular dependence on the ODME approximants F N (x) approach the QMC data with increasing N, and the extrapolated (N → ∞) values are well within the QMC errors.
We have examined several mappings other than those shown in Figs. A.4 and A.5. Among the ones not shown, those that have good convergence properties give results for the 3D Fermi gas similar to the results obtained from the first seven mappings of Fig. A.4. The input sensitivity of the ODME extrapolants is well-controlled for many mappings, such as the ones used in Figs. A.4 and A.5, for (ξ, d 1 , d 2 ) varied in ranges comparable to the ones specified there. Approximant sequences with poor convergence behavior can appear for these mappings if one allows larger input variations, but this can be dealt with by selecting sequences of ODME approximants according to their convergence for each (ξ, d 1 , d 2 ) input; see also the main text.