Diagrammatic Monte Carlo study of mass-imbalanced Fermi-polaron system

We apply the diagrammatic Monte Carlo approach to three-dimensional Fermi-polaron systems with mass-imbalance, where an impurity interacts resonantly with a noninteracting Fermi sea whose atoms have a different mass. This method allows to go beyond frequently used variational techniques by stochastically summing all relevant impurity Feynman diagrams up to a maximum expansion order limited by the sign problem. Polaron energy and quasiparticle residue can be accurately determined over a broad range of impurity masses. Furthermore, the spectral function of an imbalanced polaron demonstrates the stability of the quasiparticle and allows to locate in addition also the repulsive polaron as an excited state. The quantitative exactness of two-particle-hole wave-functions is investigated, resulting in a relative lowering of polaronic energies in the mass-imbalance phase diagram. Tan's contact coefficient for the mass-balanced polaron system is found in good agreement with variational methods. Mass-imbalanced systems can be studied experimentally by ultracold atom mixtures like $^6$Li-$^{40}$K.


I. INTRODUCTION
One of the most general and successful concepts in physics is the separation of a physical system into a simpler, controlled subsystem that is interacting with a perturbing subsystem. A specific example of this method is given by a basic impurity problem, consisting of a noninteracting homogeneous medium and one particle disturbing it. In the case of a noninteracting Fermi gas, this is called Fermi-polaron problem 1 . This theoretical model can help to map out the phase diagram of a strongly population-imbalanced Fermi gas 2 , where the quasiparticle energy and effective mass serve as input parameters for Landau-Pomeranchuk Hamiltonians 3,4 helping to quantify zero temperature phase separation and the ground state energy of different phases. Moreover, the N +1 Fermi-polaron system was shown to undergo a transition of its own, featuring as possible ground states 5,6 a polaronic spin-1/2 quasiparticle and the composite spin-0 molecule, consisting of the impurity and a single bath atom. However, this transition does not generalize easily to population-imbalanced Fermi gases as it might be preempted by phase separation.
One of the most common approaches to the Fermipolaron problem is the use of a variational ansatz 2,7-12 which can be translated into diagrammatic language. This method was able to map the transition between polaronic and molecular states and helped to evaluate key quasiparticle properties, e.g., effective mass or polaron residue. Mathy et al. extended these ideas to the case of a mass-imbalanced Fermi-polaron problem and depicted the ground state phase diagram with respect to polaronic, molecular, trimer and tetramer states 4 , where a trimer (tetramer) is the bound state of two (three) bath particles with the impurity. Concerning other techniques, Functional Renormalization Group 13 and fixed-node diffusion Monte Carlo 14 have been successfully applied. Experi-mental works include Refs. 15-17. Another approach was presented by Prokof'ev and Svistunov in 2008, diagrammatic Monte Carlo 6 (di-agMC). Its key ingredient is the sampling of Feynman diagram integrals by a set of ergodic updates linking all topologies and internal variables. By reducing the diagrammatic space, they managed to reach sufficiently high expansion orders allowing them to extract energies and effective masses in very good agreement with other techniques despite the fermionic sign problem. Recently, the method was used [18][19][20] for the extraction of polaron quasiparticle residues and two-dimensional geometries.
Up to now, these diagMC implementations have only been applied to the special case of equal masses of impurity and bath atoms. This is important as such a system can be created by different spin states of a homogeneous atomic gas. However, also mixtures like 6 Li-40 K are experimentally realizable in ultracold atom systems. In our work, we extend diagMC to the case of arbitrary massimbalance and present the dependence of polaron energy and residue on the imbalance ratio. Determining the polaronic spectral function for mass-imbalance helps understanding the stability of quasiparticles close to the limit of a heavy impurity. It features the repulsive polaron, an excited state with finite lifetime. We show that twoparticle-hole wave-functions remain essentially exact in three dimensions and demonstrate the implications for the mass-imbalanced phase diagram. We also present results for Tan's contact parameter for a mass-balanced polaron system. This paper is structured as follows: Section II presents the basic Fermi-polaron model summarizing the diagrammatic ingredients for imbalanced masses. In Sec. III, some changes of the diagrammatic routine are proposed in order to increase performance, while Sec. IV introduces the enhancements of bold diagMC we use in our code. Section V will explain a newly developed regrouping technique helping to increase extrapolation speed of the series. Finally, Sec. VI exhibits our results. A brief conclusion is given in Sec. VII, while the appendices give a detailed derivation of the mass-imbalanced T matrix and explain our extrapolation procedure.

II. MODEL
The system referred to as a Fermi-polaron problem consists of two main ingredients: a noninteracting Fermi bath and an impurity resonantly interacting with it. This can be realized at ultracold temperatures because s-wave scattering between any bath particles is forbidden due to Pauli's principle and p-wave scattering is energetically suppressed. The simplest Hamiltonian compatible with these requirements is 2 c k,σ labels a field operator for a particle in state σ and with momentum k, g is the bare coupling constant and k,σ = k 2 2mσ incorporates the energy dispersion. For convenience, we label the impurity by ↓ and bath particles by ↑, although this does not necessarily refer to pure spin states. k F is the bath particle Fermi momentum, E F its Fermi energy. We are working in units assuring = 1. For the diagrammatic series, we employ the imaginary time representation and establish the following particlehole convention: Bath lines with momentum |k| > k F are denoted particles possessing positive time while lines with |k| < k F are called holes and propagate with negative time. The positive time direction is defined to be from left to right. The inter-particle interaction can be modeled by a pseudo-potential in case of a zero-range interaction. It is important, though, that this potential assures the correct two-particle scattering length 21,22 a. We introduce the T matrix Γ(τ, p) in conventional form 6 and tabulate it prior to the Monte Carlo run. This is tantamount to replacing the bare coupling g by the experimentally accessible scattering length a. We refer to appendix A for details. The Green's functions are given by where µ 0 ↓ is a tuning parameter used for convergence reasons. Reduced mass m r and total mass M are introduced canonically FIG. 1. The first-order diagram is used for normalization purposes. The (local) appearance of this diagram topology as part of the whole diagram will be used to identify reducible diagrams in partially bold diagMC in Sec. IV.
FIG. 2. This (unphysical) second-order diagram connects the first-order fake diagram with higher-order diagrams.

III. DIAGRAMMATIC FRAMEWORK
The set of updates we use is different from the approach of Prokof'ev-Svistunov 6 . Rather than linking different orders by worm diagrams, we prefer to implement these transitions by direct updates. We propose the following update pairs: • Go-to-fake and Go-from-fake, • Insert-mushroom and Remove-mushroom, • Insert-T matrix and Remove-T matrix, • Reconnect (self-inverse).
A fake diagram is used for normalization purposes and is graphically identical to the first-order diagram, but with analytically easy weights, cf. Fig. 1. Its internal variables are updated by the update Change-Fake. The updates Go-to-fake and Go-from-fake connect this fake diagram with the lowest-order diagram we sample, presented in Fig. 2. Note that this diagram is unphysical if no boldification scheme is used. The first-order diagram is not included in our simulation because its contribution experiences a 1 √ τ -behavior for τ → 0, thus forcing the program to spend a lot of time on small times. It is straightforward to include the first-order self-energy by a numerical tabulation in ω space 7 .
There are four updates linking different orders: Insertmushroom, Insert-T matrix and their inverse updates Remove-mushroom and Remove-T matrix. Insert-T matrix chooses any T matrix of the current diagram and splits it into two linked T matrices. The resulting (unphysical) diagram has good overlap with the previous configuration if G 1 ↓ −G 0 ↓ is artificially attributed as weight of the underlying impurity propagator. Here, G 1 ↓ denotes the impurity Green's function evaluated by plugging the first-order self-energy contribution into Dyson's equation.
Last, an update called Reconnect ensures that all different topologies of a certain order are sampled.
This set of updates is ergodic and avoids sampling of first-order contributions. The last remaining unphysical diagrams connect two adjacent T matrices -however, this is important for partially bold sampling (cf. section III). If no boldification is used, sampling of relevant diagrams can be enforced by assigning an additional penalty weight to those diagrams. We will present the updates Insert-mushroom, Remove-mushroom and Reconnect in the next subsections. All other updates were designed in the same spirit.

A. Insert-mushroom
This update is available for impurity propagator lines. It attempts to insert the diagrammatic structure of Fig. 1 (called mushroom) on one of these lines. If the current diagram order is denoted by N , there are N − 1 propagators available for insertion. Having selected one of those with imaginary time τ and momentum p, internal time slices τ 1 and τ 2 are uniformly seeded (probabilities: dτ1 τ and dτ2 τ −τ1 ), as well as a bath propagator momentum q with |q| < k F (probability: . This fixes the time variable of the last piece to τ 3 = τ − τ 1 − τ 2 . The whole process is illustrated in Fig. 3. The acceptance factor P IM is

(4)
The factor (2π) 3 in the denominator is part of the diagrammatic weight of the new configuration. p IM and p RM are the probabilities of selecting the updates Insertmushroom or Remove-mushroom, respectively.

B. Remove-mushroom
Remove-mushroom is the inverse update for Insertmushroom. If the current diagram order is denoted by N , there are N T matrices which could be removed together with the corresponding impurity propagator. However, the first T matrix can not be removed, because it is never constructed by Insert-Mushroom. The same is true for the last T matrix. That leaves N − 2 possible T matrices and balances the selection factors in the Metropolis algorithm. Being the inverse update of Insert-Mushroom, the acceptance factor of Remove-Mushroom is given by Reconnect is the key update of our procedure. It randomly selects one of the T matrices, except the last one. In diagram order N , there are N − 1 possible choices. Suppose that a T matrix with parameters (τ t , p t ) and an impurity neighbor adjacent to the right with parameters (τ p , p ↓ ) is selected. The update then proposes to swap the incoming bath propagator with the incoming bath propagator of its right neighbor. There is a unique way of swapping as we are not working in cyclical representation. Accordingly, the former bath propagator times τ 1 and τ 2 are exactly mapped on new times τ 1 and τ 2 . Index 1 labels the bath propagator linked to the selected T matrix. Note that the mapping of the bath propagator momenta to the corresponding new momenta is not clear at this instant, as the shape of the current topology has to be reflected. In the moment of linking the new propagator configuration, the underlying momenta have to be adjusted in a manner described below. Last, the resulting diagram has to be checked for one-particle-irreducibility -the update has to be rejected if any impurity propagator line turns out uncovered. Subsequent application of Reconnect updates allows to reach every bath propagator configuration and guarantees ergodicity.
More precisely, the update separates into three different cases depending on the current diagram configura-tion. The first diagram topology (cf. Fig. 4) is identified by having a mushroom-structure on the selected T matrix -its incoming bath line is connected with its outgoing bath line. Since swapping will transfer a hole into a bath particle, it is necessary to create its new particle momentum q. This is done by uniform seeding on the interval [−k max , k max ] for each component of q, where k max introduces the momentum cutoff of our procedure. The update is rejected if |q| > k max or if |q| < k F . Concerning underlying momenta, the selected T matrix is assigned the momentum p of the right neighboring T matrix, while its right impurity neighbor obtains p − q. It is easy to compute final times The acceptance factor is The second topology is identified by a link between the outgoing end of the selected T matrix and its right neighbor. Being the inverse of the latter update, only one more step is necessary. Instead of seeding new particle momentum, now the hole momentum has to be created on the selected T matrix, thus explaining the factor of 2k F in Eq. 6. This equation yields the inverse acceptance factor for the second topology, (P 2 R ) −1 . All other cases are included in the third topology (cf. Fig. 5), defined by neither connecting the selected T matrix with its right neighbor nor with itself. Such cases are self-inverse. No seeding is necessary, all bath momenta are purely swapped or added. Determining final times is straightforward Note that holes are defined to have negative times. Concerning momenta, only the selected T matrix and its right impurity neighbor have to be considered, yielding new momenta P t for T matrix and P ↓ for impurity line: This results in the acceptance factor . (9)

IV. PARTIALLY BOLD DIAGRAMMATIC MONTE CARLO
In addition to the modifications of the bare code, we propose some changes only affecting the bold diagram-matic Monte Carlo routine. Using full Green's functions as self-consistent input has the disadvantage of increasing sampling space drastically by forcing the tabulation of the full Green's function in both imaginary time and momentum. It is therefore beneficial to use the fact that first-order contributions dominate the fully bold propagator and construct a partially bold diagrammatic series out of quantities that are easily tabulated. It is possible to put the analytically known first-order self-energy into Dyson's equation in Matsubara frequency space to obtain the first-order Green's function. A Fourier transform yields the basic propagator G 1 ↓ (τ, p) without stochastic errors. Only a few modifications are necessary to adjust the basic Monte Carlo routine: First, every diagram containing at least one first-order self-energy diagram (cf. Fig. 1) has to be excluded from measurements. Second, it is no longer forbidden to connect adjacent T matrices, just as in fully bold code 6 -the only difference is that the associated impurity weight is now G 1 ↓ − G 0 ↓ . We also extended the molecule code to include partially bold propagators. This extension is almost identical to the bold polaron code. Linking T matrices is readmitted with the same (now non-artificial) weight G 1 ↓ − G 0 ↓ for impurity propagators. The only new feature concerns the first-order molecule diagram that is now sampled and has to be calculated with impurity weight G 1 ↓ − G 0 ↓ . By this means, the ultraviolet divergence is cured and a second, independent molecule series is constructed which helps confirming robustness and reproducibility of the results. Resummation is still needed for this new series.
As a last comment on bold sampling, we would like to draw attention to the flaws of bold diagMC. First, if the Dyson series is not absolutely convergent, the rearrangement of this series into the fully or partially bold series is potentially harmful, as it constitutes a (though physically motivated) regrouping of terms which can yield any result for non-convergent series 23 . Second, many series in quantum field theory are asymptotic expansions, implying that results begin to get better and better with increasing expansion order until some maximum expansion order N max is reached, after which the factorial growth of the number of diagrams leads to huge oscillations. In such a case, using a bare series is a common procedure, whereas the bold approach captures diagrams of higher order from the start. Summarizing this line of argument, a bold diagrammatic approach seems only reasonable if the underlying series is convergent.

V. DIAGRAM REGROUPING AND RESUMMATION
The resummation procedure requires a discussion in more detail. Typically we find the molecular energies to be stable, but the polaron energies in the Bose-Einsteincondensate limit are harder. Sharp resummations are potentially dangerous if the maximum sampling order is not high enough as can be seen as follows: With a strong resummation method, the produced curve is almost flat for low expansion orders and then bends down sharply for higher orders. Weak resummation methods on the other hand have more curvature for low expansion orders and flatten off if the order of divergence of the series is weak enough. There is thus a risk with strong resummation methods if only low expansion orders can be reached in the sense that a possibly strong curvature is missed, resulting in an apparently converging but wrong extrapolated result in close vicinity to the first-order result. This effect is demonstrated in Fig. 6. Note that the bare series is mainly decreasing with increasing expansion order. At unitarity, it is even monotonously decreasing which makes a high maximum resummation order nec-essary to extract the correct answer. Summing it up: the more sign-alternating the bare series is, the better resummation works. A typical resummation method is the Riesz resummation method. These resummations will act upon selfenergy series which can be written as S = N S (N ) , where S (N ) contains all contributions of diagrams of order N . The order of a self-energy diagram is defined as the number of interaction T matrices. The resummed self-energy series S for some given maximum order L is defined 6 as where F is given by the Riesz coefficients δ fixes the strength of the resummation: For δ = 0, no resummation is performed at all, while only the first-order contribution is maintained in the limit δ → +∞. If this method is used on molecular series, it might be beneficial to set F (L) 2 = 1 as this ensures that the first contributing diagram (which is of second order for molecules) always contributes with full weight.
We introduce a regrouping technique which seems to saturate much faster. Provided the series is absolutely convergent, this is always allowed. It is based on a regrouping of terms in the bare series in such a way that sign-alternation is maximized. The technique consists in splitting S (N ) , the self-energy contributions of order N , into two parts: where S r (N ) collects all diagrams containing at least one T matrix linked by a hole to itself, cf. Fig. 1, and S ir (N ) gathers the rest. We propose a new series S = N S (N ) which is sign-alternating by construction. The coefficients in the resummation procedure depend on the expansion order for k ∈ N as These coefficients are in principle arbitrary -our choice was designed to show fast saturation as can be seen in same answer can be found but only after a stronger resummation. Note that although the sum of diagrams of a specific order of the unitary polaron series is vanishing within error bars, the different terms are not small. As the regrouped series agrees with the results of the standard bare series, this might indicate that the polaron Dyson series is convergent or that the maximum expansion order N max of its asymptotic expansion is essentially infinity at unitarity. Nevertheless, Fig. 7 of Ref. 18 suggests that the series might be asymptotic in fact, as a doubly bold regrouping shows clear signs of growing fluctuations for increasing expansion order after an initial improvement of results.

VI. RESULTS
In this section, our diagrammatic Monte Carlo results for a mass-imbalanced polaron are presented. here labeled as first order, captures the whole curve qualitatively, its quantitative accuracy gets less precise for low r, i.e., a light polaron. Note that this energy curve reproduces the correct infinite mass limit 7 E pol = −0.5 for an imbalance ratio as low as r = 2. For the case of an immobile impurity, the polaron is subject to Anderson's orthogonality catastrophe 24  The quasiparticle residue Z can be extracted from the asymptotic decay of the full propagator 6

A. Polaron energy and residue at unitarity
where E is the ground state energy of the quasiparticle. This form implies that a suitable Monte Carlo estimator for the residue is S is the sampled self-energy.
Our results for Z are depicted in Fig. 9. In this case, the first-order ansatz is both quantitatively and qualitatively different, predicting a different position of the mass-imbalance ratio of maximum residue. Therefore the quasiparticle with maximum residue can be found at higher r than estimated by first order. At low r, higher orders affect Z stronger and stronger, down to a ratio as low as r = 0.125 which is sufficient for most mixtures, e.g., 6 Li-40 K. The measured residue is strictly lower than the first order variational result. This is remarkable as the Functional Renormalization Group analysis of Ref. 13 predicts a higher residue for unitarity at r = 1. Further investigation is needed to understand this discrepancy.
No resummation was used to extract the quasiparticle residues. For r 0.5, the series seemed to saturate within our maximum expansion order. The extrapolation error was approximated to be twice the fluctuation of the saturating points. For r 0.5, the series changed and the saturation was not visible anymore. These points are therefore only valid if a linear extrapolation to infinite expansion order is appropriate. This extrapolation error was approximated by the method explained in appendix B.  10. (Color online) Different contributions to the full polaron energy are compared for three maximum expansion orders N . En denotes contributions including up to n-particlehole diagrams. Two of the curves have been offset by ±0.03EF for clarity. The points were measured for a mass-imbalance ratio r = 0.125 at unitarity.

B. Two-particle-hole channel
For quasi-two-dimensional geometries, a remarkable quantitative precision of two-particle-hole wave-functions was found for polaron energies 20,25 . A natural question is whether this approach remains valid in three dimensions. Fig. 10 compares different particle-hole channels for three maximum expansion orders. For the selected measurement point (unitarity with imbalance r = 0.125), the three-particle-hole channel contributes with slight quantitative differences, whereas four-particle-hole and fiveparticle-hole diagrams vanish within error bars. This confirms the observation that the two-particle-hole result is essentially correct and can be used for quantitative calculations. The classification of diagrams into particle-hole classes breaks down for the fully bold approach, because each bold diagram captures bare diagrams of different particle-hole order. For the partially bold scheme, this does not apply since the diagrammatic structure of the partially bold Green's functions ensures that the propagation will start with zero holes and is guaranteed to switch to the 1-ph sector at least once.

C. Spectral function
It is possible to extract the polaronic spectral function a mass-imbalance ratio r = 2. Although the energy corresponds to the infinite mass limit 7 , the polaron remains a stable quasiparticle. The dispersion follows a parabola with positive effective mass, whereas at higher energies, the repulsive polaron (a metastable eigenstate of the Hamiltonian) can clearly be seen 11,13 .

D. Quantitative exactness of variational energies
In this subsection, the extraction of polaron and molecule energies by resummation is compared to variational one-particle-hole wave-functions. We choose the point k F a = 0.5 and r = 0.25 of the mass-imbalanced phase diagram 4 as it stays away from the peculiarities of unitarity and the trimer threshold. Molecular energies (Fig. 12) yield perfect agreement with the variational ansatz and motivate the quantitative correctness of variational wave-functions for the Fermi-polaron problem. For the polaron (Fig. 13), using a one-particle-hole wave-function underestimates its energy slightly. Hence, it would be beneficial to use at least two-particle-hole precision for the polaron sector for a precise mapping of the phase diagram. Altogether, the phase diagram of Ref. 4 can be expected to be nearly quantitatively exact. Nevertheless, as the polaron phase is underestimated, it will be shifted into the molecular sector. As this will reduce the small size of the nonzero-momentum molecular phase (FFLO) further, it remains open whether this phase really exists. We suggest to compute the phase boundary with high precision within a variational 2-ph polaron approach.

E. Contact coefficient
Tan's contact density C for the strongly populationimbalanced Fermi gas is linked to the dimensionless contact density s by Here, k F,↓ is the Fermi momentum of the minority species which is finite for the strongly imbalanced Fermi gas. The dimensionless contact density s can be accessed easily by calculating the derivative of the polaron energy with respect to the dimensionless coupling 9 (k F,↑ a) −1 . The re-sulting contact curve agrees with first order calculations within error bars.

VII. CONCLUSION
Our work extends the diagrammatic Monte Carlo polaron routines to the more general case of a massimbalanced polaron. This is an important limiting case of population imbalanced Fermi gases and allows to estimate key properties of its phase diagram. The diagrammatic space can be drastically reduced by sampling the self-energy instead of the full Green's functions and by using the ladder approximation as basic interaction element. We presented a critical analysis and alternative check of diagrammatic Monte Carlo as well as a partially bold approach, thus broadening the toolbox of the diagrammatic Monte Carlo method. An additional regrouping technique was presented to speed up extrapolation to infinite diagram order for absolutely convergent series. While the first-order variational ansatz could give qualitative and quantitative good results for the polaron energy at different polaron masses, discrepancies are more pronounced for the polaron residue. For this quantity, higher orders have to be included in order to capture the whole physics. Concerning Tan's contact coefficient, an excellent agreement was found with the Chevy variational wave-function. The polaronic spectral function was extracted from imaginary time representation of diagrammatic Monte Carlo data by means of analytic continuation. It demonstrates a clean parabolic dispersion as well as the existence of the repulsive polaron. The two-particle-hole wave-function ansatz, which was shown to be essentially exact in quasitwo-dimensional geometries 20 , provides an equally good description of quasi-particle energies in three dimensions. Therefore, using one-particle-hole trial functions will lead to a phase diagram which overestimates molecular contributions and might lead to a weakening of the FFLO state found in Ref (A2) The integral is calculated by residue calculus; as only the pole of G 0 ↓ is in the upper plane, it follows 1 The problem of the latter expression is an ultraviolet divergence in q, consequently it has to be regularized.
A similar series appears in the standard scattering Lippmann-Schwinger equation and can be used for regularization. Starting point is the Lippmann-Schwinger form of the relative motion Schrödinger equation ImposingV on the left and introducing the two-body T matrixT 2B viaT 2B |k ≡V |ψ k , it followŝ In the next step, k | is multiplied from the left with |k| = |k | and a complete set of eigenfunctions ofĤ 0 is inserted The matrix elements of the potential are trivial for a pseudo-potential and the newly inserted h| is an eigenstate ofĤ 0 Looking at h|T 2B |k more closely (A8) Remarkably, this expression does not depend on the first element. As there is no h-dependence left, the T matrix can be extracted from the integral (A9) The two-body T matrix is related to the scattering length 29 if the effective range can be set to zero k is related to the energy of relative motion and can be generalized to off-shell behavior Note that the total energy ω is measured with respect to E F and µ 0 ↓ . p is the total momentum of the system. Inserting this into equation A9 yields Finally, a suitable expression for the T matrix can be found by inserting equation A12 into equation A3 This can be cast into its final form by a shift q → q +

Appendix B: Extrapolation of resummed data
In this appendix, we explain the details of our resummation procedure. Its use is most delicate for cases where the maximum diagram order is small. Therefore, it is illustrative to use a quasi-two-dimensional Fermipolaron series (characterized by the dimensionless parameter η = k 2D F a 2D , where a 2D is the two-dimensional scattering length and k 2D F is the two-dimensional Fermi momentum) to explain this technique, since the maximum expansion order is approximately 8. For these systems, it is additionally necessary to deal with large binding energies, hence aggressive resummation has to be applied  to the bare series in order to be able to extrapolate to infinite expansion order. However, this tends to conceal the curvature of the series in the first points, leading to an initially flat curve.
Our extrapolation procedure is the following: For the upper value of the error bar, we apply linear extrapolation on the Riesz-resummed data with Riesz exponent δ. In this linear extrapolation, only the two points corresponding to highest and second highest expansion order are taken into account. For the lower value of the error bar, we assume a worst-case scenario with large curvature of the extrapolated curve, according to the following fit function f of parameters a and b: N denotes the maximum expansion order. We emphasize that the curvature of this function is empirically set by us. This curve includes only the highest and second highest expansion points. An important feature of f is the dependence on the Riesz exponent. This ensures that a stronger resummation results in a bigger error bar due to extrapolation errors. The fit f can also be used for bare data (δ = 0). In this case, we replace δ by -1 in Eq. B1. Note that the error bars represent a variability of results due to systematic origins corresponding to the one-σ-interval. The result of this technique is shown in Fig. 14. For a maximum expansion order of 7, the two ways of extrapolating are shown in comparison. If the maximum expansion order was 6, then this error bar would increase, just as one would expect regarding the loss of information. Fig. 15 uses a sharper resummation on the same data, demonstrating that the error bar increases with δ.
Therefore, it becomes clear that the weakest possible resummation procedure (among the ones resulting in a monotonously decaying series) should be applied. As a final example, we show our resummation for a polaron point inside the quasi-two-dimensional transition region in Fig. 16. Here, resummation with δ = 2 is too weak, so δ = 4 has to be used, resulting in a stronger curvature. It is important to stress that these two extrapolations represent assumed worst-case scenarios. Finally, as always for diagMC, the extrapolation result has to be checked with available experimental or theoretical results, thus justifying its application in retrospect. In our case, these results would be variational 2-ph results which we expect to be quantitatively exact. As an example for a system in which extrapolated error bars were underestimated for a similar system, consult Fig. 22 of Ref. 30.
For molecular energies, resummation is more straightforward. As this series is typically alternating, resummed curves can often be extrapolated linearly. This is demonstrated in Fig. 17: The last four points are well fit by a straight line. However, as this resummation involves the same dangers as described above, we try to vary both the fitting (e.g., fitting three of the four last points) and the resummation technique to test the variability of this result.