Diagrammatic Monte Carlo study of mass-imbalanced Fermi-polaron system
DDiagrammatic Monte Carlo study of mass-imbalanced Fermi-polaron system
Peter Kroiss and Lode Pollet
Department of Physics, Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience,University of Munich, Theresienstrasse 37, 80333 Munich, Germany (Dated: October 11, 2018)We apply the diagrammatic Monte Carlo approach to three-dimensional Fermi-polaron systemswith mass-imbalance, where an impurity interacts resonantly with a noninteracting Fermi sea whoseatoms have a different mass. This method allows to go beyond frequently used variational techniquesby stochastically summing all relevant impurity Feynman diagrams up to a maximum expansionorder limited by the sign problem. Polaron energy and quasiparticle residue can be accurately deter-mined over a broad range of impurity masses. Furthermore, the spectral function of an imbalancedpolaron demonstrates the stability of the quasiparticle and allows to locate in addition also the re-pulsive polaron as an excited state. The quantitative exactness of two-particle-hole wave-functionsis investigated, resulting in a relative lowering of polaronic energies in the mass-imbalance phasediagram. Tan’s contact coefficient for the mass-balanced polaron system is found in good agreementwith variational methods. Mass-imbalanced systems can be studied experimentally by ultracoldatom mixtures like Li– K. PACS numbers: 02.70.Ss, 05.10.Ln, 05.30.Fk
I. INTRODUCTION
One of the most general and successful concepts inphysics is the separation of a physical system into a sim-pler, controlled subsystem that is interacting with a per-turbing subsystem. A specific example of this method isgiven by a basic impurity problem, consisting of a nonin-teracting homogeneous medium and one particle disturb-ing it. In the case of a noninteracting Fermi gas, this iscalled Fermi-polaron problem . This theoretical modelcan help to map out the phase diagram of a stronglypopulation-imbalanced Fermi gas , where the quasipar-ticle energy and effective mass serve as input param-eters for Landau-Pomeranchuk Hamiltonians helpingto quantify zero temperature phase separation and theground state energy of different phases. Moreover, the N +1 Fermi-polaron system was shown to undergo a tran-sition of its own, featuring as possible ground states apolaronic spin-1/2 quasiparticle and the composite spin-0 molecule, consisting of the impurity and a single bathatom. However, this transition does not generalize eas-ily to population-imbalanced Fermi gases as it might bepreempted by phase separation.Two of the most common approaches to the Fermi-polaron problem are variational ans¨atze and the use ofan approximate diagrammatic technique . Thesemethods are able to calculate key quasiparticle proper-ties, e.g., effective mass or polaron residue and helpedto map the transition between polaronic and molecularstates. Mathy et al. extended these ideas to the case ofa mass-imbalanced Fermi-polaron problem and depictedthe ground state phase diagram with respect to pola-ronic, molecular, trimer and tetramer states, wherea trimer (tetramer) is the bound state of two (three) bathparticles with the impurity. Concerning other techniques,Functional Renormalization Group and fixed-node dif-fusion Monte Carlo have been successfully applied. Ex- perimental works include Refs. 20–22.Another approach was presented by Prokof’ev andSvistunov in 2008, diagrammatic Monte Carlo (di-agMC). Its key ingredient is the sampling of Feynmandiagram integrals by a set of ergodic updates linking alltopologies and internal variables. By reducing the dia-grammatic space, they managed to reach sufficiently highexpansion orders allowing them to extract energies andeffective masses in very good agreement with other tech-niques despite the fermionic sign problem. Recently, themethod was used for the extraction of polaron quasi-particle residues and two-dimensional geometries.Up to now, these diagMC implementations have onlybeen applied to the special case of equal masses of impu-rity and bath atoms. This is important as such a systemcan be created by different spin states of a homogeneousatomic gas. However, also mixtures like Li– K are ex-perimentally realizable in ultracold atom systems. In ourwork, we extend diagMC to the case of arbitrary mass-imbalance and present the dependence of polaron energyand residue on the imbalance ratio. Determining the po-laronic spectral function for mass-imbalance helps under-standing the stability of quasiparticles close to the limitof a heavy impurity. It features the repulsive polaron ,an excited state with finite lifetime. We show that two-particle-hole wave-functions remain essentially exact inthree dimensions and demonstrate the implications forthe mass-imbalanced phase diagram. We also presentresults for Tan’s contact parameter for a mass-balancedpolaron system.This paper is structured as follows: Section II presentsthe basic Fermi-polaron model summarizing the diagram-matic ingredients for imbalanced masses. In Sec. III,some changes of the diagrammatic routine are proposedin order to increase performance, while Sec. IV intro-duces the enhancements of bold diagMC we use in ourcode. Section V will explain a newly developed regroup-ing technique helping to increase extrapolation speed of a r X i v : . [ c ond - m a t . qu a n t - g a s ] A p r the series. Finally, Sec. VI exhibits our results. A briefconclusion is given in Sec. VII, while the appendices givea detailed derivation of the mass-imbalanced T matrixand explain our extrapolation procedure. II. MODEL
The system referred to as a Fermi-polaron problemconsists of two main ingredients: a noninteracting Fermibath and an impurity resonantly interacting with it. Thiscan be realized at ultracold temperatures because s-wavescattering between any bath particles is forbidden dueto Pauli’s principle and p-wave scattering is energeticallysuppressed. The simplest Hamiltonian compatible withthese requirements is ˆ H = (cid:88) k ,σ (cid:15) k,σ ˆ c † k ,σ ˆ c k ,σ + g (cid:88) k , k (cid:48) , q ˆ c † k + q , ↑ ˆ c † k (cid:48) − q , ↓ ˆ c k (cid:48) , ↓ ˆ c k , ↑ . (1)ˆ c k ,σ labels a field operator for a particle in state σ andwith momentum k , g is the bare coupling constant and (cid:15) k,σ = k m σ incorporates the energy dispersion. For con-venience, we label the impurity by ↓ and bath particlesby ↑ , although this does not necessarily refer to pure spinstates. k F is the bath particle Fermi momentum, E F itsFermi energy. We are working in units assuring (cid:126) = 1.For the diagrammatic series, we employ the imaginarytime representation and establish the following particle-hole convention: Bath lines with momentum | k | > k F are denoted particles possessing positive time while lineswith | k | < k F are called holes and propagate with neg-ative time. The positive time direction is defined to befrom left to right. The inter-particle interaction can bemodeled by a pseudo-potential in case of a zero-rangeinteraction. It is important, though, that this potentialassures the correct two-particle scattering length a .We introduce the T matrix Γ( τ, p ) in conventional form and tabulate it prior to the Monte Carlo run. This istantamount to replacing the bare coupling g by the ex-perimentally accessible scattering length a . We refer toappendix A for details. The Green’s functions are givenby G ↑ ( τ, k ) = − θ ( τ ) θ ( k − k F ) e − ( (cid:15) k, ↑ − E F ) τ + θ ( − τ ) θ ( k F − k ) e − ( (cid:15) k, ↑ − E F ) τ G ↓ ( τ, k ) = − θ ( τ ) e − ( (cid:15) k, ↓ − µ ↓ ) τ , (2)where µ ↓ is a tuning parameter used for convergence rea-sons. Reduced mass m r and total mass M are introducedcanonically M = m ↓ + m ↑ m r = m ↓ m ↑ M . (3)
FIG. 1. The first-order diagram is used for normalizationpurposes. The (local) appearance of this diagram topology aspart of the whole diagram will be used to identify reduciblediagrams in partially bold diagMC in Sec. IV.FIG. 2. This (unphysical) second-order diagram connectsthe first-order fake diagram with higher-order diagrams.
III. DIAGRAMMATIC FRAMEWORK
The set of updates we use is different from the ap-proach of Prokof’ev-Svistunov . Rather than linking dif-ferent orders by worm diagrams, we prefer to implementthese transitions by direct updates. We propose the fol-lowing update pairs: • First-to-fake and Fake-to-first, • Change-fake (self-inverse), • Insert-mushroom and Remove-mushroom, • Insert- T matrix and Remove- T matrix, • Reconnect (self-inverse).A fake diagram is used for normalization purposes andis graphically identical to the first-order diagram, butwith analytically easy weights, cf. Fig. 1. Its internalvariables are updated by the update Change-Fake. Theupdates First-to-fake and Fake-to-first connect this fakediagram with the lowest-order diagram we sample, pre-sented in Fig. 2. Note that this diagram is unphysical ifno self-consistent bold scheme is used. The first-orderdiagram is not included in our simulation because its con-tribution experiences a √ τ -behavior for τ →
0, thus forc-ing the program to spend a lot of time on small times.It is straightforward to include the first-order self-energyby a numerical tabulation in ω space .There are four updates linking different orders: Insert-mushroom, Insert- T matrix and their inverse updatesRemove-mushroom and Remove- T matrix. Insert- T ma-trix chooses any T matrix of the current diagram andsplits it into two linked T matrices. The resulting (un-physical) diagram has good overlap with the previousconfiguration if G ↓ − G ↓ is artificially attributed as weightof the underlying impurity propagator. Here, G ↓ denotesthe impurity Green’s function evaluated by plugging thefirst-order self-energy contribution into Dyson’s equation.Last, an update called Reconnect ensures that all dif-ferent topologies of a certain order are sampled.This set of updates is ergodic and avoids sampling offirst-order contributions. The last remaining unphysicaldiagrams connect two adjacent T matrices – however,this is important for partially bold sampling (cf. sectionIII). If no self-consistent bold scheme is used, sampling ofrelevant diagrams can be enforced by assigning an addi-tional penalty weight to those diagrams. We will presentthe updates Insert-mushroom, Remove-mushroom andReconnect in the next subsections. All other updateswere 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 currentdiagram order is denoted by N , there are N − τ and momentum p , internal timeslices τ and τ are uniformly seeded (probabilities: d τ τ and d τ τ − τ ), as well as a bath propagator momentum q with | q | < k F (probability: d q (2 k F ) ). This fixes the timevariable of the last piece to τ = τ − τ − τ . The wholeprocess is illustrated in Fig. 3. The Metropolis accep-tance ratio P IM ismin (cid:32) , p RM p IM G ↓ ( τ , p )Γ( τ , p + q ) G ↓ ( τ , p ) G ↑ ( τ , q )(2 π ) G ↓ ( τ, p ) · τ τ − τ k F ) (cid:33) . (4)The factor (2 π ) in the denominator is part of the di-agrammatic weight of the new configuration. p IM and p RM are the probabilities of selecting the updates Insert-mushroom or Remove-mushroom, respectively. B. Remove-mushroom
Remove-mushroom is the inverse update for Insert-mushroom. If the current diagram order is denoted by N ,there are N T matrices which could be removed togetherwith the corresponding impurity propagator. However,
FIG. 3. Illustration of inverse updates Insert-mushroom andRemove-mushroom. FIG. 4. Illustration of the first and second case of updateReconnect. The dotted vertical bath propagator line is con-nected to an arbitrary T matrix in the diagram.FIG. 5. Illustration of the third case of update Reconnect.The dotted and full vertical bath propagator lines are con-nected to arbitrary T matrices in the diagram (as long asneither the full line of the left figure nor the dotted line of theright figure is connected to itself. the first T matrix can not be removed, because it is neverconstructed by Insert-Mushroom. The same is true forthe last T matrix. That leaves N − T matricesand balances the selection factors in the Metropolis al-gorithm. Being the inverse update of Insert-Mushroom,the acceptance ratio of Remove-Mushroom is given bymin (cid:32) , p IM p RM (2 π ) G ↓ ( τ, p ) · τ τ − τ k F ) G ↓ ( τ , p )Γ( τ , p + q ) G ↓ ( τ , p ) G ↑ ( τ , q ) (cid:33) . (5) C. Reconnect
Reconnect is the key update of our procedure. It ran-domly selects one of the T matrices, except the last one.In diagram order N , there are N − T matrix with parameters ( τ t , p t ) and animpurity neighbor adjacent to the right with parameters( τ p , p ↓ ) is selected. The update then proposes to swapthe incoming bath propagator with the incoming bathpropagator of its right neighbor. There is a unique wayof swapping as we are not working in cyclical represen-tation. Accordingly, the former bath propagator times τ and τ are exactly mapped on new times τ (cid:48) and τ (cid:48) .Index 1 labels the bath propagator linked to the selected T matrix. Note that the mapping of the bath propagatormomenta to the corresponding new momenta is not clearat this instant, as the shape of the current topology hasto be reflected. In the moment of linking the new propa-gator configuration, the underlying momenta have to beadjusted in a manner described below. Last, the resultingdiagram has to be checked for one-particle-irreducibility– the update has to be rejected if any impurity propaga-tor line turns out uncovered. Subsequent application ofReconnect updates allows to reach every bath propagatorconfiguration and guarantees ergodicity.More precisely, the update separates into three differ-ent cases depending on the current diagram configura-tion. The first diagram topology (cf. Fig. 4) is identifiedby having a mushroom-structure on the selected T matrix– its incoming bath line is connected with its outgoingbath line. Since swapping will transfer a hole into a bathparticle, it is necessary to create its new particle momen-tum q . This is done by uniform seeding on the interval[ − k max , k max ] for each component of q , where k max in-troduces the momentum cutoff of our procedure. Theupdate is rejected if | q | > k max or if | q | < k F . Con-cerning underlying momenta, the selected T matrix isassigned the momentum p of the right neighboring T matrix, while its right impurity neighbor obtains p − q .It is easy to compute final times τ (cid:48) = τ − τ t − τ p τ (cid:48) = τ p . (6)The acceptance ratio ismin (cid:32) , Γ( τ t , p ) G ↓ ( τ p , p − q ) G ↑ ( τ (cid:48) , p ) G ↑ ( τ (cid:48) , q ) k Γ( τ t , p t ) G ↓ ( τ p , p ↓ ) G ↑ ( τ , p ) G ↑ ( τ , p ) k F (cid:33) . (7)The second topology is identified by a link betweenthe outgoing end of the selected T matrix and its rightneighbor. Being the inverse of the latter update, only onemore step is necessary. Instead of seeding new particlemomentum, now the hole momentum has to be createdon the selected T matrix, thus explaining the factor of k F in Eq. 7. The acceptance ratio for the second topology ismin (cid:32) , Γ( τ t , p t ) G ↓ ( τ p , p ↓ ) G ↑ ( τ , p ) G ↑ ( τ , p ) k F Γ( τ t , p ) G ↓ ( τ p , p − q ) G ↑ ( τ (cid:48) , p ) G ↑ ( τ (cid:48) , q ) k (cid:33) . (8)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. Suchcases are self-inverse. No seeding is necessary, all bathmomenta are purely swapped or added. Determining fi-nal times is straightforward τ (cid:48) = τ − τ t − τ p τ (cid:48) = τ + τ t + τ p . (9)Note that holes are defined to have negative times. Con-cerning momenta, only the selected T matrix and its rightimpurity neighbor have to be considered, yielding newmomenta P t for T matrix and P ↓ for impurity line: P t = p t + p − p P ↓ = p ↓ + p − p p (cid:48) = p p (cid:48) = p . (10) This results in the acceptance ratiomin (cid:32) , Γ( τ t , P t ) G ↓ ( τ p , P ↓ ) G ↑ ( τ (cid:48) , p ) G ↑ ( τ (cid:48) , p )Γ( τ t , p t ) G ↓ ( τ p , p ↓ ) G ↑ ( τ , p ) G ↑ ( τ , p ) (cid:33) . (11) IV. PARTIALLY BOLD DIAGRAMMATICMONTE CARLO
In addition to the modifications of the bare code, wepropose some changes only affecting the bold diagram-matic Monte Carlo routine. Using full Green’s functionsas self-consistent input has the disadvantage of increas-ing sampling space drastically by forcing the tabulationof the full Green’s function in both imaginary time andmomentum. It is therefore beneficial to use the fact thatfirst-order contributions dominate the fully bold propa-gator and construct a partially bold diagrammatic seriesout of quantities that are easily tabulated. It is possibleto put the analytically known first-order self-energy intoDyson’s equation in Matsubara frequency space to ob-tain the first-order Green’s function. A Fourier transformyields the basic propagator G ↓ ( τ, p ) without stochasticerrors. Only a few modifications are necessary to ad-just the basic Monte Carlo routine: First, every dia-gram containing at least one first-order self-energy dia-gram (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 – the only differenceis that the associated impurity weight is now G ↓ − G ↓ .We also extended the molecule code to include partiallybold propagators. This extension is almost identical tothe bold polaron code. Linking T matrices is readmittedwith the same (now non-artificial) weight G ↓ − G ↓ forimpurity propagators. The only new feature concerns thefirst-order molecule diagram that is now sampled and hasto be calculated with impurity weight G ↓ − G ↓ . By thismeans, the ultraviolet divergence is cured and a second,independent molecule series is constructed which helpsconfirming robustness and reproducibility of the results.Resummation is still needed for this new series.As a last comment on bold sampling, we would like todraw attention to the flaws of bold diagMC. First, if theDyson series is not absolutely convergent, the rearrange-ment of this series into the fully or partially bold seriesis potentially harmful, as it constitutes a (though physi-cally motivated) regrouping of terms which can yield anyresult for non-convergent series . Second, many seriesin quantum field theory are asymptotic expansions, im-plying that results begin to get better and better with in-creasing expansion order until some maximum expansionorder N max is reached, after which the factorial growthof the number of diagrams leads to huge oscillations. Insuch a case, using a bare series is a common procedure,whereas the bold approach captures diagrams of higherorder from the start. Summarizing this line of argument,a bold diagrammatic approach seems only reasonable if N − E p o l / E F δ = 1 δ = 2 δ = 4 FIG. 6. (Color online) Resummation methods of type Rieszwith different exponents δ are compared for unitarity in athree-dimensional setup. The plot shows polaron energies de-pending on maximum sampling order. N − E p o l / E F δ = 1 δ = 2 δ = 4 FIG. 7. (Color online) Resummation methods of type Rieszwith different exponents δ applied to the modified bare series(see Eq. 15) are compared at unitarity in a three-dimensionalsetup. The plot shows polaron energies depending on maxi-mum sampling order. the underlying series is convergent. V. DIAGRAM REGROUPING ANDRESUMMATION
As diagrammatic expansions are in general not abso-lutely convergent, an important tool to study the un-derlying series is resummation . This resummation pro-cedure requires a discussion in more detail. Typicallywe find the molecular energies to be stable, but thepolaron energies in the Bose-Einstein-condensate limitare harder. Sharp resummations are potentially danger-ous if the maximum sampling order is not high enough as can be seen as follows: With a strong resummationmethod, the produced curve is almost flat for low ex-pansion orders and then bends down sharply for higherorders. Weak resummation methods on the other handhave more curvature for low expansion orders and flat-ten off if the order of divergence of the series is weakenough. There is thus a risk with strong resummationmethods if only low expansion orders can be reached inthe sense that a possibly strong curvature is missed, re-sulting in an apparently converging but wrongly extrap-olated result in close vicinity to the first-order result.This effect is demonstrated in Fig. 6 at unitarity. Notethat the bare series is monotonously decreasing whichmakes a high maximum resummation order necessary toextract the correct answer. Summing it up: the moresign-alternating the bare series is, the better resumma-tion works.A typical resummation method is the Riesz resumma-tion method. These resummations will act upon self-energy series which can be written as S = (cid:80) N S ( N ) ,where S ( N ) contains all contributions of diagrams of or-der N . The order of a self-energy diagram is defined asthe number of interaction T matrices. The resummedself-energy series S (cid:48) for some given maximum order L isdefined as S (cid:48) ( L ) = L (cid:88) N =1 S ( N ) F ( L ) N , (12)where F is given by the Riesz coefficients F ( L ) N = (cid:18) L − N + 1 L (cid:19) δ . (13) δ fixes the strength of the resummation: For δ = 0, no re-summation is performed at all, while only the first-ordercontribution is maintained in the limit δ → + ∞ . If thismethod is used on molecular series, it might be beneficialto set F ( L )2 = 1 as this ensures that the first contributingdiagram (which is of second order for molecules) alwayscontributes with full weight.We introduce a regrouping technique which seems tosaturate much faster. Provided the series is absolutelyconvergent, this is always allowed. It is based on a re-grouping of terms in the bare series in such a way thatsign-alternation is maximized. The technique consists insplitting S ( N ) , the self-energy contributions of order N ,into two parts: S ( N ) = S r ( N ) + S ir ( N ) , (14)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 (cid:48) = (cid:80) N S (cid:48) ( N ) which aims to maximize sign-alternation in S (cid:48) ( N ) . Thecoefficients in the resummation procedure depend on theexpansion order for k ∈ N as S (cid:48) (1) = S (1) , S (cid:48) ( N =2 k +1) = S ir ( N ) + 12 S ir ( N − + 12 S r ( N ) , S (cid:48) ( N =2 k ) = S r ( N ) + 12 S r ( N − + 12 S ir ( N ) . (15)These coefficients are in principle arbitrary – our choicewas designed to show fast saturation as can be seen inFig. 7 for the case of the Fermi-polaron at unitarity:The application of conventional Riesz resummations onthe reordered series illustrates that the manually imple-mented oscillations make it possible to read off polaronenergies reliably and allow for a clear statement whetherthe expansion order is high enough or not. When revert-ing the roles of reducible and irreducible in Eq. 15, thesame answer can be found but only after a stronger re-summation. Note that although the sum of diagrams ofa specific order of the unitary polaron series is vanishingwithin error bars, the different terms are not small.As the regrouped series agrees with the results of thestandard bare series, this might indicate that the polaronDyson series is convergent or that the maximum expan-sion order N max of its asymptotic expansion is essentiallyinfinity at unitarity. Nevertheless, Fig. 7 of Ref. 23 sug-gests that the series might be asymptotic in fact, as adoubly bold regrouping shows clear signs of growing fluc-tuations for increasing expansion order after an initialimprovement of results. VI. RESULTS
In this section, our diagrammatic Monte Carlo resultsfor a mass-imbalanced polaron are presented.
A. Polaron energy and residue at unitarity
Fig. 8 plots the polaron energy at unitarity for differentmass ratios r = m ↓ m ↑ . While the variational Chevy ansatz,here labeled as first order, captures the whole curve qual-itatively, its quantitative accuracy gets less precise forlow r , i.e., a light polaron. Note that this energy curvereproduces the correct infinite mass limit E pol = − . r = 2. For the case ofan immobile impurity, the polaron is subject to Ander-son’s orthogonality catastrophe and the quasiparticledescription is no longer appropriate. For a light impu-rity, the polaron energy decreases rapidly as the effectiveinteraction is stronger for smaller reduced mass (see Ap-pendix A) at unitarity. The next subsection will showthat this effect is weakened for finite scattering length.Eventually, a very light impurity will be subject to rel-ativistic effects so that our description will no longer beappropriate. The extraction of error bars was based onconservative extrapolations of light Riesz resummationwith δ = 1. For details, consult appendix B. r E p o l [ E F ] r E p o l [ E F ] Extrapolated DiagMCFirst order DiagMC
FIG. 8. (Color online) Polaron energy at unitarity for differ-ent mass ratios r . The inset shows the flat part of the figure.Many-body effects get more pronounced for a lighter polaron. r Z Extrapolated DiagMCFirst order DiagMC r Z FIG. 9. (Color online) Polaron residue at unitarity for dif-ferent mass ratios r . The whole range of masses shows a cleardifference between the first order result and the full diagram-matic answer. The inset shows the curve for additional valuesof r . The quasiparticle residue Z can be extracted from theasymptotic decay of the full propagator G ↓ ( τ, k ) −→ τ →∞ − Ze − ( E − µ ↓ ) τ , where E is the ground state energy of thequasiparticle. This form implies that a suitable MonteCarlo estimator for the residue is Z = (1 + A ( E, p )) − with A ( E, p ) = − (cid:90) ∞ τ S ( τ, p ) e ( E − µ ↓ ) τ d τ. (16) 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 quali-tatively different, predicting a different position of themass-imbalance ratio of maximum residue. This mightindicate that the variational wave-function descriptionworks particularly well for energy based quantities, whileit might take further particle-hole terms to capture theresidue equally well. Therefore the quasiparticle withmaximum residue can be found at higher r than esti-mated by first order. At low r , higher orders affect Z stronger and stronger, down to a ratio as low as r = 0 . Li– K. Inthis regime, the differences between the diagrammatic an-swer and the first order result are most pronounced. Forhigh r , the diagMC solution yields a roughly constantshift to the first order answer.As a trimer state acquires more and more strength withrespect to the polaron state for decreasing r and as themolecular state is strengthened for increasing r , it seemsnatural that the polaron residue takes on a maximumvalue in between – once it is no longer the ground state,its residue will decay quickly (although it will not be zerobecause there is no simple decay channel ). The mea-sured residue is strictly lower than the first order varia-tional result. This is remarkable as the Functional Renor-malization Group analysis of Ref. 18 predicts a higherresidue for unitarity at r = 1. Further investigation isneeded to understand this discrepancy.No resummation was used to extract the quasiparti-cle residues. For r (cid:38) .
5, the series seemed to saturatewithin our maximum expansion order. The extrapolationerror was approximated to be twice the fluctuation of thesaturating points. For r (cid:46) .
5, the series changed andthe saturation was not visible anymore. These points aretherefore only valid if a linear extrapolation to infiniteexpansion order is appropriate. This extrapolation errorwas approximated by the method explained in appendixB.
B. Polaron energy beyond unitarity
Fig. 10 shows polaron energies at varying couplingstrength in the BEC-regime for two different mass-imbalance ratios r . Both curves experience a peak ofmaximum dressing around ( k F a ) − = 0 .
4. Decreasingthe coupling further, this relative energy is decreased forboth masses, although the light impurity is affected morestrongly. Eventually, the heavy impurity has a higher ef-fective dressing (relative to the binding energy) than thelight impurity. This is a consequence of the m r / (2 πa )term in the denominator of the T matrix (see AppendixA) that strengthens the effective interaction between im-purity and bath atoms for increasing r at a given in-teraction strength. Note that these curves extend intothe molecular sector where the polaron ceases to be theground state. Concerning the residue, no qualitative dif-ference could be seen between the r = 0 . r = 2curves. ( k F a ) − ( E p o l − E B ) / E F r = 2r = 0.5 FIG. 10. (Color online) The polaron energy of two valuesof r is plotted for various interaction parameters k F a . Notethat the binding energy E B = (2 m r a ) − is subtracted. Forstrong interactions, the light impurity acquires a higher ef-fective dressing relative to the binding energy than the heavyimpurity. This is eventually reversed in the BEC-regime. N ( E − E )[ E F ]( E − E )[ E F ]( E − E )[ E F ] FIG. 11. (Color online) Different contributions to the full po-laron energy are compared for three maximum expansion or-ders N . E n denotes contributions including up to n-particle-hole diagrams. Two of the curves have been offset by ± . E F for clarity. The points were measured for a mass-imbalanceratio r = 0 .
125 at unitarity.
C. Two-particle-hole channel
For quasi-two-dimensional geometries, a remarkablequantitative precision of two-particle-hole wave-functionswas found for polaron energies . A natural question iswhether this approach remains valid in three dimensions.Fig. 11 compares different particle-hole channels for threemaximum expansion orders. For the selected measure-ment point (unitarity with imbalance r = 0 . p [ k F ] ω [ E F ] FIG. 12. (Color online) The polaron spectral function isplotted for a mass-imbalance r = 2 at unitarity. Note thequadratic dispersion and the repulsive polaron at positive en-ergies. particle-hole diagrams vanish within error bars. Thisconfirms the observation that the two-particle-hole re-sult is essentially correct and can be used for quanti-tative calculations. The classification of diagrams intoparticle-hole classes breaks down for the fully bold ap-proach, because each bold diagram captures bare dia-grams of different particle-hole order. For the partiallybold scheme, this does not apply since the diagrammaticstructure of the partially bold Green’s functions ensuresthat the propagation will start with zero holes and isguaranteed to switch to the 1-ph sector at least once. D. Spectral function
It is possible to extract the polaronic spectral function A ( ω, p ) = − (cid:16) ω + i − (cid:15) p + µ ↓ − S ( ω, p ) (cid:17) − fromthe sampled self-energy S . The full Green’s functioncan be calculated by Fourier transform, subsequent ap-plication of Dyson’s equation and another Fourier trans-form. Then, analytic continuation to real quantities isapplied . Fig. 12 presents the spectral function fora mass-imbalance ratio r = 2. Although the energycorresponds to the infinite mass limit , the polaron re-mains a stable quasiparticle. The dispersion follows aparabola with positive effective mass, whereas at higherenergies, the repulsive polaron (a metastable eigenstateof the Hamiltonian) can clearly be seen . E. Quantitative exactness of variational energies
In this subsection, the extraction of polaron andmolecule energies by resummation is compared to vari-ational one-particle-hole wave-functions. We choose the E m o l [ E F ] FIG. 13. (Color online) Molecule energy extrapolation for k F a = 0 . r = 0 .
25. The maximum expansion order isdenoted by N . The solid line corresponds to a linear fit of thelast five points. The crosses mark the one-particle-hole resultof Ref. 4. Riesz resummation with exponent δ = 6 was used. E p o l [ E F ] FIG. 14. (Color online) Polaron energy extrapolation for k F a = 0 . r = 0 .
25. The maximum expansion order isdenoted by N . The solid line corresponds to a linear fit of thelast five points, the dashed line is a special fitting function ex-plained in Appendix B. The one-particle-hole result of Ref. 4is identically with the point N = 2. Riesz resummation withexponent δ = 4 was used. point k F a = 0 . r = 0 .
25 of the mass-imbalancedphase diagram as it stays away from the peculiaritiesof unitarity and the trimer threshold. Molecular ener-gies (Fig. 13) yield perfect agreement with the varia-tional ansatz and motivate the quantitative correctnessof variational wave-functions for the Fermi-polaron prob-lem. For the polaron (Fig. 14), using a one-particle-holewave-function underestimates its energy slightly. Hence,it would be beneficial to use at least two-particle-hole pre-cision for the polaron sector for a precise mapping of thephase diagram. Altogether, the phase diagram of Ref. 4can be expected to be nearly quantitatively exact. Never-theless, as the polaron phase is underestimated, it will beshifted into the molecular sector. As this will reduce thesmall size of the nonzero-momentum molecular phase (la-beled as Fulde-Ferrell-Larkin-Ovchinnikov phase, FFLO)further, it remains open whether this phase really exists.We suggest to compute the phase boundary with highprecision within a variational 2-ph polaron approach. F. Dimensionless contact coefficient
Tan’s contact coefficient C for a strongly population-imbalanced Fermi gas is linked to the dimensionless con-tact coefficient s by C = s · k F, ↓ k F, ↑ . (17)Here, k F, ↓ is the Fermi momentum of the minority specieswhich is finite for the strongly imbalanced Fermi gas. Thedimensionless contact coefficient s can be accessed easilyby calculating the derivative of the polaron energy withrespect to the dimensionless coupling ( k F, ↑ a ) − . The re-sulting contact curve agrees with first order calculationswithin error bars. VII. CONCLUSION
Our work extends the diagrammatic Monte Carlo po-laron routines to the more general case of a mass-imbalanced polaron. This is an important limiting caseof population imbalanced Fermi gases and allows to esti-mate key properties of its phase diagram. The diagram-matic space can be drastically reduced by sampling theself-energy instead of the full Green’s functions and byusing the ladder approximation as basic interaction el-ement. We presented a critical analysis and alternativecheck of diagrammatic Monte Carlo as well as a par-tially bold approach, thus broadening the toolbox of thediagrammatic Monte Carlo method. An additional re-grouping technique was presented to speed up extrap-olation to infinite diagram order for absolutely conver-gent series. While the first-order variational ansatz couldgive qualitative and quantitative good results for thepolaron energy at different polaron masses, discrepan-cies are more pronounced for the polaron residue. Forthis quantity, higher orders have to be included in or-der to capture the whole physics. Concerning Tan’scontact coefficient, an excellent agreement was foundwith the Chevy variational wave-function. The pola-ronic spectral function was extracted from imaginarytime representation of diagrammatic Monte Carlo databy means of analytic continuation. It demonstrates aclean parabolic dispersion as well as the existence of therepulsive polaron. The two-particle-hole wave-functionansatz, which was shown to be essentially exact in quasi-two-dimensional geometries , provides an equally good description of quasi-particle energies in three dimensions.Therefore, using one-particle-hole trial functions will leadto a phase diagram which overestimates molecular con-tributions and might lead to a weakening of the FFLOstate found in Ref. 4.We are grateful to M. Bauer, T. Enss, J. Levinsen,M. Parish, N. Prokof’ev, R. Schmidt, B. Svistunov, K.Van Houcke and W. Zwerger for valuable discussions.We would also like to thank M. Parish for sharing herdata with us. This work is supported by the ExcellenceCluster NIM, FP7/Marie-Curie Grant No. 321918 (FDI-AGMC), FP7/ERC Starting Grant No. 306897 (QUSIM-GAS) and by a grant from the Army Research Office withfunding from DARPA. Appendix A: Many-body T matrix for unequalmasses It is possible to write the T matrix in presence of theFermi sea in a self-consistent way − i Γ( ω, p ) = − ig + (cid:90) q>k F d q (2 π ) (cid:90) dq π ( − ig ) iG ↑ ( q , q ) iG ↓ ( ω − q , p − q )( − i Γ( ω, p )) . (A1)Since T is replacing an interaction line in first approx-imation it has to follow the same sign convention as g .Also note that for a point interaction, g ( k ) is constantand independent of momentum. Rewriting yields1Γ( ω, p ) =1 g − i (cid:90) q>k F d q (2 π ) (cid:90) dq π G ↓ ( ω − q , p − q ) G ↑ ( q , q ) . (A2)The integral is calculated by residue calculus; as only thepole of G ↓ is in the upper plane, it follows1Γ( ω, p ) =1 g − (cid:90) q>k F d q (2 π ) q m ↑ + ( p − q ) m ↓ − E F − µ ↓ − ω . (A3)The problem of the latter expression is an ultraviolet di-vergence in q , consequently it has to be regularized.A similar series appears in the standard scatteringLippmann-Schwinger equation and can be used for reg-ularization. Starting point is the Lippmann-Schwingerform of the relative motion Schr¨odinger equation | ψ k (cid:105) = | k (cid:105) + 1 E − ˆ H + i V | ψ k (cid:105) . (A4)Imposing ˆ V on the left and introducing the two-body T matrix ˆ T B via ˆ T B | k (cid:105) ≡ ˆ V | ψ k (cid:105) , it followsˆ T B | k (cid:105) = ˆ V | k (cid:105) + ˆ V E − ˆ H + i T B | k (cid:105) . (A5)0In the next step, (cid:104) k (cid:48) | is multiplied from the left with | k | = | k (cid:48) | and a complete set of eigenfunctions of ˆ H isinserted (cid:104) k (cid:48) | ˆ T B | k (cid:105) = (cid:104) k (cid:48) | ˆ V | k (cid:105) + (cid:90) d h (2 π ) (cid:104) k (cid:48) | ˆ V | h (cid:105) (cid:104) h | E − ˆ H + i T B | k (cid:105) . (A6)The matrix elements of the potential are trivial for apseudo-potential and the newly inserted (cid:104) h | is an eigen-state of ˆ H (cid:104) k (cid:48) | ˆ T B | k (cid:105) = g + g (cid:90) d h (2 π ) E − h m r + i (cid:104) h | ˆ T B | k (cid:105) . (A7)Looking at (cid:104) h | ˆ T B | k (cid:105) more closely (cid:104) h | ˆ T B | k (cid:105) = (cid:68) h (cid:12)(cid:12)(cid:12) ˆ V (cid:12)(cid:12)(cid:12) ψ k (cid:69) = (cid:90) d k (cid:48)(cid:48) (2 π ) (cid:104) h | ˆ V | k (cid:48)(cid:48) (cid:105) (cid:104) k (cid:48)(cid:48) | ψ k (cid:105) = (cid:90) d k (cid:48)(cid:48) (2 π ) g (cid:104) k (cid:48)(cid:48) | ψ k (cid:105) = (cid:90) d k (cid:48)(cid:48) (2 π ) (cid:104) k (cid:48) | ˆ V | k (cid:48)(cid:48) (cid:105) (cid:104) k (cid:48)(cid:48) | ψ k (cid:105) = (cid:104) k (cid:48) | ˆ T B | k (cid:105) . (A8)Remarkably, this expression does not depend on the firstelement. As there is no h -dependence left, the T matrixcan be extracted from the integral1 (cid:104) k (cid:48) | ˆ T B | k (cid:105) = 1 g − (cid:90) d h (2 π ) E − h m r . (A9)The two-body T matrix is related to the scatteringlength if the effective range can be set to zero1 (cid:104) k (cid:48) | ˆ T B | k (cid:105) = m r πa (1 + iak ) . (A10) k is related to the energy of relative motion and can begeneralized to off-shell behavior k m r = E = E tot − E com = ω + E F + µ ↓ − p M . (A11)Note that the total energy ω is measured with respectto E F and µ ↓ . p is the total momentum of the system.Inserting this into equation A9 yields1 g = m r πa − m r π (cid:114) m r M p − m r ( ω + E F + µ ↓ )+ (cid:90) d h (2 π ) p M + h m r − ω − E F − µ ↓ . (A12)Finally, a suitable expression for the T matrix can be E p o l [ E F ] FIG. 15. (color online) The extrapolation procedure onresummed data is demonstrated for η = − .
248 and δ = 2.The dotted curve shows the linear part, whereas the dashedline is a fit with f . found by inserting equation A12 into equation A3Γ( ω, k ) − = m r πa − m r π (cid:114) m r M k − m r ( E F + µ ↓ + ω )+ (cid:90) d q (2 π ) k M + ( q − k ) m r − E F − µ ↓ − ω − (cid:90) q>k F d q (2 π ) q m ↑ + ( k − q ) m ↓ − E F − µ ↓ − ω . (A13)This can be cast into its final form by a shift q → q + k − m ↑ k M in the first integralΓ( ω, k ) − = m r πa − m r π (cid:114) m r M k − m r ( E F + µ ↓ + ω )+ (cid:90) q In this appendix, we explain the details of our re-summation procedure. Its use is most delicate for caseswhere the maximum diagram order is small. Therefore,it is illustrative to use a quasi-two-dimensional Fermi-polaron series (characterized by the dimensionless param-eter η = ln ( k F a ), where a is the two-dimensionalscattering length and k F is the two-dimensional Fermimomentum) to explain this technique, since the maxi-mum expansion order is approximately 8. For these sys-tems, it is additionally necessary to deal with large bind-ing energies, hence aggressive resummation has to be ap-1 E p o l [ E F ] FIG. 16. (color online) The extrapolation procedure onresummed data is demonstrated for η = − . 248 and δ = 4.The dotted curve shows the linear part, whereas the dashedline is a fit with f . Note that the resulting error bar exceedsthe corresponding δ = 2 error bar. E p o l [ E F ] FIG. 17. (color online) The extrapolation procedure onresummed data is demonstrated for η = − . 02 and δ = 4.The dotted curve shows the linear part, whereas the dashedline is a fit with f . plied to the bare series in order to be able to extrapolateto infinite expansion order. However, this tends to con-ceal the curvature of the series in the first points, leadingto an initially flat curve.Our extrapolation procedure is the following: For theupper value of the error bar, we apply linear extrapola-tion on the Riesz-resummed data with Riesz exponent δ .In this linear extrapolation, only the two points corre-sponding to highest and second highest expansion orderare taken into account. For the lower value of the errorbar, we assume a worst-case scenario with large curva-ture of the extrapolated curve, according to the following E m o l [ E F ] FIG. 18. (color online) The extrapolation procedure on re-summed data is demonstrated for a molecule with η = − . δ = 6. fit function f of parameters a and b : f ( N − , a, b ) = 4 δ (cid:18) N − − N − (cid:19) + aN − + b. (B1) N denotes the maximum expansion order. We empha-size that the curvature of this function is empirically setby us. This curve includes only the highest and secondhighest expansion points. An important feature of f isthe dependence on the Riesz exponent. This ensures thata stronger resummation results in a bigger error bar dueto extrapolation errors. The fit f can also be used forbare data ( δ = 0). In this case, we replace δ by -1 inEq. B1. Note that the error bars represent a variabil-ity of results due to systematic origins corresponding tothe one- σ -interval. The result of this technique is shownin Fig. 15. For a maximum expansion order of 7, thetwo ways of extrapolating are shown in comparison. Ifthe maximum expansion order was 6, then this error barwould increase, just as one would expect regarding theloss of information. Fig. 16 uses a sharper resummationon the same data, demonstrating that the error bar in-creases with δ .Therefore, it becomes clear that the weakest possibleresummation procedure (among the ones resulting in amonotonously decaying series) should be applied. As afinal example, we show our resummation for a polaronpoint inside the quasi-two-dimensional transition regionin Fig. 17. Here, resummation with δ = 2 is too weak, so δ = 4 has to be used, resulting in a stronger curvature. Itis important to stress that these two extrapolations rep-resent assumed worst-case scenarios. Finally, as alwaysfor diagMC, the extrapolation result has to be checkedwith available experimental or theoretical results, thusjustifying its application in retrospect. In our case, theseresults would be variational 2-ph results which we expectto be quantitatively exact. As an example for a systemin which extrapolated error bars were underestimated for2a similar system, consult Fig. 22 of Ref. 38.For molecular energies, resummation is more straight-forward. As this series is typically alternating, resummedcurves can often be extrapolated linearly. This is demon-strated in Fig. 18: The last four points are well fit by a straight line. However, as this resummation involves thesame dangers as described above, we try to vary boththe fitting (e.g., fitting three of the four last points) andthe resummation technique to test the variability of thisresult. P. Massignan, M. Zaccanti and G. M. Bruun, Rep. Prog.Phys. , 034401 (2014). F. Chevy, Phys. Rev. A , 063628 (2006). S. Pilati and S. Giorgini, Phys. Rev. Lett. , 030401(2008). C. J. M. Mathy, M. M. Parish and D. A. Huse, Phys. Rev.Lett. , 166404 (2011). N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. B ,020408(R) (2008). N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. B ,125101 (2008). R. Combescot, A. Recati, C. Lobo and F. Chevy, Phys.Rev. Lett. , 180402 (2007). R. Combescot and S. Giraud, Phys. Rev. Lett. , 050404(2008). M. Punk, P. T. Dumitrescu and W. Zwerger, Phys. Rev.A , 053605 (2009). C. Mora and F. Chevy, Phys. Rev. A , 033607 (2009). P. Massignan and G. M. Bruun, Eur. Phys. J. D , 83(2011). C. Trefzger and Y. Castin, Europhys. Lett. , 3006(2013). R. Combescot, S. Giraud and X. Leyronas, Europhys. Lett. , 60007 (2009). J. E. Baarsma, J. Armaitis, R. A. Duine, and H. T. C.Stoof, Phys. Rev. A , 033631 (2012). O. I. Kartavtsev and A. V. Malykh, J. Phys. B , 1429(2007). D. Blume, Phys. Rev. Lett. , 230404 (2012). J. Levinsen and M. M. Parish, Phys. Rev. Lett. ,055304 (2013). R. Schmidt and T. Enss, Phys. Rev. A , 063620 (2011). C. Lobo, A. Recati, S. Giorgini, S. Stringari, Phys. Rev.Lett. , 200403 (2006). A. Schirotzek, C. H. Wu, A. Sommer and M. W. Zwierlein,Phys. Rev. Lett. , 230402 (2009). S. Nascimb`ene, N. Navon, K. J. Jiang, L. Tarruell, M. Teichmann, J. McKeever, F. Chevy, and C. Salomon, Phys.Rev. Lett. , 170402 (2009). C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P.Massignan, G. M. Bruun, F. Schreck and R. Grimm, Na-ture (London) , 615 (2012). J. Vlietinck, J. Ryckebusch and K. Van Houcke, Phys. Rev.B , 115133 (2013). J. Vlietinck, J. Ryckebusch and K. Van Houcke, Phys. Rev.B , 085119 (2014). P. Kroiss and L. Pollet, Phys. Rev. B , 104510 (2014). X. Cui and H. Zhai, Phys. Rev. A , 041602(R) (2010). I. Bloch, J. Dalibard and W. Zwerger, Rev. Mod. Phys. , 885 (2008). M. Inguscio, W. Ketterle and C. Salomon Ultra-Cold FermiGases: Varenna 20-30 June 2006 , Ios PressInc, ISBN9781586038465 (2007). We define the order of a polaron (molecule) diagram bythe number of T matrices (impurity propagators). E. Kozik, M. Ferrero, A. Georges, arXiv:1407.5687. We use the same energy estimator as Ref. 6 throughoutthe paper. P. Anderson, Phys. Rev. Lett. M. M. Parish and J. Levinsen, Phys. Rev. A , 033616(2013). M. Jarrell and J. E. Gubernatis, Phys. Rep. , 133(1996). B. Bauer, L. Carr, H. G. Evertz, A. Feiguin, J. Freire, S.Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler etal. , J. Stat. Mech. (2011) P05001. This section is taken from P. Kroiss, diploma thesis (2013)and follows Refs. 6 and 37 in large parts. H. Stoof, K. Gubbels and D. Dickerscheid Ultracold quan-tum fields (1st Edition), Springer, ISBN 9781402087639(2009).38