Diagrammatic Monte Carlo algorithm for the resonant Fermi gas
K. Van Houcke, F. Werner, T. Ohgoe, N. Prokof'ev, B. Svistunov
DDiagrammatic Monte Carlo algorithm for the resonant Fermi gas
K. Van Houcke,
1, 2
F. Werner,
3, 2
T. Ohgoe, N. V. Prokof’ev,
2, 5 and B. V. Svistunov
2, 5, 6 Laboratoire de Physique de l’Ecole Normale Sup´erieure, ENS - PSL, CNRS,Sorbonne Universit´e, Universit´e Paris-Diderot, Sorbonne Paris Cit´e, 75005 Paris, France Department of Physics, University of Massachusetts, Amherst, MA 01003, USA Laboratoire Kastler Brossel, Ecole Normale Sup´erieure - PSL, Sorbonne Universit´e,Coll`ege de France, CNRS, 24 rue Lhomond, 75005 Paris, France Department of Applied Physics, University of Tokyo,7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Russian Research Center “Kurchatov Institute”, 123182 Moscow, Russia Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute,Shanghai Jiao Tong University, Shanghai 200240, China (Dated: January 23, 2019)We provide a description of a diagrammatic Monte Carlo algorithm for the resonant Fermi gasin the normal phase. Details are given on diagrammatic framework, Monte Carlo moves, andincorporation of ultraviolet asymptotics. Apart from the self-consistent bold scheme, we also describea non-self-consistent scheme, for which the ultraviolet treatment is more involved.
I. INTRODUCTION
A major long-standing challenge is to find a methodfor solving a generic fermionic many-body problem inthe thermodynamic limit with controlled accuracy. Thediagrammatic technique is the most versatile quantum-field-theoretical tool allowing one to express the answersas series of integrals of a special structure. Each termin the series can be visualized with graphs—Feynmandiagrams—built using simple rules. In the absence ofsmall parameters, there is little hope to sum the dia-grammatic series analytically, and one commonly resortsto uncontrollable truncations. In contrast, the goal ofthe Diagrammatic Monte Carlo (DiagMC) approach isto sum up all Feynman diagrams in a systematic wayup to a controlled accuracy. Using an efficient MonteCarlo algorithm to evaluate all diagrams up to a highenough order N max , convergence as a function of N max can be observed, as first demonstrated for the Hubbardmodel . The thermodynamic limit is taken from theoutset since one works only with connected diagrams.Furthermore one can build diagrams with fully-dressedpropagators; this self-consistent formulation, called BoldDiagrammatic Monte Carlo (BDMC), was first demon-strated for the Fermi-polaron and the resonant Fermigas .The resonant Fermi gas is a three-dimensionalcontinuous-space model of great interdisciplinary inter-est. It features a smooth crossover between fermionicand bosonic superfluidity, as argued in the context ofcondensed matter physics and later observed exper-imentally in ultracold atomic Fermi gases near Fesh-bach resonances. The model is also relevant to neutronmatter and high-energy physics, in particular in theunitary limit where the scattering length is infinite.For the unitary Fermi gas in the normal unpolar-ized phase, first BDMC results for the equation of statewere reported in Ref. 4. Very recently, these resultswere confirmed using a much more advanced resumma- tion method, which was found to be necessary for con-trollability, due to the fact that the series has zero ra-dius of convergence . Contact and momentum distri-bution were also computed using the new resummationmethod . In the meantime, the DiagMC approach wasalso developed further and applied to frustrated quantummagnetism and various lattice models of interactingfermions including models with electron-phonon in-teraction and topological phase transitions .In this paper, we describe the numerical methodused for the equilibrium normal resonant Fermi gas inRefs. 4, 12, and 13, in particular how to evaluate theterms of the diagrammatic series to high orders (typicallyup to order 9) using a diagrammatic Monte Carlo algo-rithm, and how to properly incorporate large-momentumasymptotics coming from the contact interactions. An-other crucial ingredient is the proper resummation of thedivergent diagrammatic series which will be detailedelsewhere .We mostly use a bold diagrammatic scheme, where di-agrams are built with fully dressed single-particle propa-gators and pair propagators. We present a set of elemen-tary Monte Carlo updates to sample this diagrammaticspace. While some features of the updating scheme areanalogous to the ones introduced for the bare series of theHubbard model in Ref. 1, an important difference is thatonly fully irreducible skeleton diagrams are sampled, sothat ergodicity has to be carefully verified. Furthermore,resonant fermions feature specific ultraviolet singulari-ties governed by an observable called contact. Thisphysics manifests itself in a natural way within our skele-ton diagrammatic framework, and is readily incorporatedinto our BDMC scheme. For cross-validation, we use notonly the self-consistent bold scheme, but also a non self-consistent “ladder scheme”, in which case the ultravioletphysics governed by the contact can also be incorporatedsemi-analytically, using a more elaborate procedure.The paper is organized as follows. In Section II,the diagrammatic framework is constructed, arriving at a r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n the skeleton series for the single-particle and pair self-energies. Section III describes the diagrammatic MonteCarlo algorithm: The diagrammatic expansion is ex-pressed as a Monte Carlo average in subsection III A,precise descriptions of configuration space, probabilitydensity and measurement procedure are given in sub-sections III B, III C and III D, the update scheme is de-scribed in subsection III E, reducibility and ergodicity is-sues are discussed in subsection III F, the self-consistentiteration procedure is described in subsection III G, andresummation is briefly mentioned in Sec. III H. Section IVdescribes ultraviolet analytics and its incorporation intoBDMC. The ladder scheme is treated in Sec. V. II. DIAGRAMMATIC FRAMEWORKA. The resonant Fermi gas model
In the zero-range model, also known as the resonantgas model, the interaction is characterized by the s -wavescattering length a . The zero-range model is a universallimit of finite-range models. More precisely, a genericinteraction of range b can be replaced by the zero-rangemodel in the limit where b becomes much smaller thanother typical lengthscales of the problem, such as theinterparticle distance, the thermal wavelength, and | a | .For an atomic alkali Fermi gas near a broad Feshbachresonance, the range is set by the van der Waals length,and most current experiments are well within the zero-range limit, with finite-range corrections in the percentor sub-percent range. Even though our Monte Carlo scheme works directlywith the zero-range interaction in continuous space, it isconvenient to start with a lattice model, thereby elimi-nating ultraviolet divergences at the initial steps of con-structing the formalism. The Hamiltonian readsˆ H (cid:48) = ˆ H − (cid:88) σ = ↑ , ↓ µ σ ˆ N σ = (cid:88) σ = ↑ , ↓ (cid:88) k ∈B ( (cid:15) k − µ σ ) ˆ c † k ,σ ˆ c k ,σ + g b (cid:88) r ˆ n ↑ ( r )ˆ n ↓ ( r ) , (1)where the spin index σ takes on the values ↑ and ↓ , the op-erator ˆ c k ,σ annihilates a spin- σ fermion of momentum k ,ˆ ψ σ ( r ) is the corresponding position-space annihilationoperator, ˆ n σ ( r ) = ˆ ψ † σ ( r ) ˆ ψ σ ( r ) is the number-density op-erator, µ σ is the spin-dependent chemical potential, ˆ N σ is the number operator for spin- σ fermions, r is a positionvector whose components are integer multiples of the lat-tice spacing b (this b can also be viewed as the interactionrange since the interaction is on-site), B =] − π/b, π/b ]is the first Brillouin zone, and the dispersion relation is (cid:15) k = k / The barecoupling constant g is adjusted to have the desired scat-tering length a for two particles on the lattice in free space, namely 1 g = 14 πa − (cid:90) B d k (2 π ) k . (2)The zero-range limit corresponds to the continuum limit b →
0, with a fixed. One can note that g → − in thislimit. B. Single-particle propagator, self-energy, andladder summation
In the standard diagrammatic formalism for the many-body problem at finite temperature, the central ob-ject is the single-particle propagator G σ ( p , τ ) = − (cid:10) T ˆ c p ,σ ( τ )ˆ c † p ,σ (0) (cid:11) , (3)where τ is the imaginary time and T[ . . . ] is the time-ordered product. This Green’s function gives access tothe momentum distribution n σ ( p ) = G σ ( p , τ = 0 − ), andto the number density n σ = G σ ( r = , τ = 0 − ) . (4)In the series expansion of G in powers of the bare couplingconstant g , each term can be represented by a Feynmangraph: (5)where the bare interaction vertex • denotes g , the thinlines denote an ideal gas propagator G (0) , and the boldline denotes the fully dressed (i.e. exact) propagator G .The first natural step to organize the higher-orderterms is to introduce the self-energy Σ, which is relatedto G by the Dyson equation, given diagrammatically by(6)i.e. 1 G σ ( p , ω n ) = 1 G (0) σ ( p , ω n ) − Σ σ ( p , ω n ) (7)for any fixed momentum p and Matsubara frequency ω n . To avoid double counting, reducible diagrams areexcluded from Σ, so that (8)Another standard step is to perform summation of lad-der diagrams: (9)Physically, such a ladder summation is natural since invacuum it would correspond to the two-body scatteringamplitude or T -matrix. This allows one to take the zero-range limit and work directly with zero-range interac-tions in continuous space. Γ (0) is an approximate pairpropagator, which can also be viewed as a renormalizedinteraction vertex; eventually, Γ (0) will be replaced bya fully dressed pair propagator in our BDMC scheme.Summation of the geometric series in Eq. (9) gives1Γ (0) ( P , Ω n ) = 1 g − Π (0) ( P , Ω n ) (10)where Π (0) is the ( G (0) G (0) ) bubble given byΠ (0) ( P , Ω n ) = − β − (cid:88) m (cid:90) B d k (2 π ) G (0) ↑ ( P / k , ω m ) · G (0) ↓ ( P / − k , Ω n − ω m )= (cid:90) B d k (2 π ) − n (0) ↑ ( P / k ) − n (0) ↓ ( P / − k ) i Ω n + 2 µ − (cid:15) P / k − (cid:15) P / − k , with the Fermi factor n (0) σ ( k ) = [1 + e β ( k / − µ σ ) ] − . Theintegral over k is finite thanks to the restriction to thefirst Brillouin zone B . Here β is the inverse temperatureand µ = ( µ ↑ + µ ↓ ) / g in Eq. (10)in favor of the scattering length a —using relation (2)—finally yields1Γ (0) ( P , Ω n ) = 14 πa − (cid:90) d k (2 π ) (cid:34) k + 1 − n (0) ↑ ( P / k ) − n (0) ↓ ( P / − k ) i Ω n + 2 µ − P / − k (cid:35) , (11)where the integration domain for k is now taken to be R instead of B , i.e. the continuum limit is taken. Thediagrammatic expansion of the self-energy can then bewritten in terms of the vertex Γ (0) instead of g ; to avoid double counting one simply has to forbid diagrams con-taining ( G (0) G (0) ) bubbles: (12)Here each G (0) line is meant to have a fixed spin label,which is not shown for simplicity. We thus arrive atthe exact diagrammatic representation of the zero-rangecontinuous-space model to be used in what follows.Many diagrammatic studies of the BEC-BCS crossoverproblem are based on the bare T -matrix, Γ (0) , and thelowest-order diagram for Σ in terms of G (0) and Γ (0) , see,e.g., Refs. 7, 43, and 44. For example, this approximationis sufficient for obtaining the exponential scaling of thecritical temperature T c ∝ e − π/ (2 k F | a | ) in the BCS limit. C. Bold pair propagator
While the diagrammatic elements introduced in theprevious section are completely standard, a more originalaspect of our diagrammatic framework is the use of afully dressed (bold) pair propagator Γ. In the case of thepolaron problem, this was done in Refs. 3 and 45. Thepropagator Γ is defined byΓ( p , τ ) = g δ ( τ ) + g · P ( p , τ ) (13)with P ( r , τ ) ≡ − (cid:68) T ( ˆΨ ↓ ˆΨ ↑ )( r , τ )( ˆΨ †↑ ˆΨ †↓ )( , (cid:69) , (14)or, diagramatically, (15)One can note that the first term in Eqs. (13,15) goes tozero in the continuum limit.Similarly to the Dyson equation that expresses the boldsingle-particle propagator G in terms of the irreduciblesingle-particle self-energy Σ [Eq. (6)], we can write aDyson equation for the bold pair propagator Γ in termsof an irreducible pair self-energy Π: (16)i.e. 1Γ( p, Ω n ) = 1Γ (0) ( p, Ω n ) − Π( p, Ω n ) . (17) D. Feynman rules for the skeleton diagrams
Bold diagrammatic Monte Carlo works with skeletondiagrams built with fully dressed (bold) lines. For theunitary Fermi gas, we use diagrams built from the boldsingle-particle propagator G σ and the bold pair propa-gator Γ defined above. The first diagrams expressing thesingle-particle self-energy Σ in terms of G and Γ are(18)while the first diagrams for the pair self-energy Π are(19)In summary, the propagators G and Γ are expressed interms of the self-energies Σ and Π through the Dysonequations (7,17), and the self-energies are themselves ex-pressed in terms of the propagators through the diagram-matic expansions (18,19).Since the Feynman rules for these diagrammatic ex-pansions are the ones which our algorithm has to obey,we describe them in some detail. The goal is to expressthe sum Σ ( N ) σ or Π ( N ) of all order- N skeleton diagrams.We define the order N of a skeleton diagram through thenumber of Γ-lines: a Σ-diagram contains N such linesand a Π-diagram contains N − Q to denote either Σ ( N ) σ or Π ( N ) , and let S Q be the set of all skeleton diagram topologies for Q , mean-ing that all these diagrams are irreducible with respectto cutting any two internal lines of the same type ( i.e. the diagram should remain connected if one cuts two G σ lines or two Γ lines). We shall use the shorthand nota-tion Y = ( p ; τ , τ ) for the external diagram variables.Clearly, Q ( Y ) = Q ( p, τ − τ ).For a given topology, we can label each internal lineby an index l for a G -line (resp. λ for a Γ-line), and de-note the corresponding internal momentum by k l (resp. κ λ ), and the spins of G -lines by σ l . Similarly, the time-differences between the end and origin points of the lines are denoted by ∆ τ l (resp. ∆ τ (cid:48) λ ). It can be shown that forany topology T in a diagram of order N , one can alwaysfind N ‘loop momenta’ q , . . . , q N that, together with theexternal momentum, uniquely determine all the internalmomenta. More precisely, some of the internal momentaare equal to a loop momentum, while the others are linearcombinations of loop momenta such that momentum isconserved at each vertex. For our Feynman diagrams,the internal variables X can thus be parameterized by q , . . . , q N , as well as by the internal times τ , . . . , τ N which belong to [0 , β ] (these times are assigned to three-point vertices which connect a Γ-line with two G -lines).With these notations, the N -th order of the diagram-matic expansion simply reads Q ( Y ) = (cid:88) T ∈S Q (cid:90) dX D ( T , X, Y ) (20)with the differential measure dX ≡ d q . . . d q N dτ . . . dτ N (21)and D ( T , X, Y ) = ( − N ( − N loop (2 π ) N × (cid:34)(cid:89) l G σ l ( k l , ∆ τ l ) (cid:35) × (cid:34)(cid:89) λ Γ( κ λ , ∆ τ (cid:48) λ ) (cid:35) , (22)with N loop the number of closed fermion loops in thediagram of topology T . There is one exception: To avoiddouble counting, in the first-order diagram for Π, we haveto compensate for the fact that all ( G (0) G (0) ) bubbles arealready contained in Γ (0) :Π (1) ( p, τ ) = − π ) (cid:90) d q D (23)with D = G ↑ ( q , τ ) G ↓ ( p − q , τ ) − G (0) ↑ ( q , τ ) G (0) ↓ ( p − q , τ ) . (24)Note also that a diagram topology T is defined here bya graph with fixed spin labels.Note that if we restrict to the lowest order diagram inEq. (18) and (19), our framework becomes equivalent tothe approach introduced in Refs. 47 and 48. This ap-proach is called self-consistent T -matrix approximation,because Γ is then given by the ladder diagrams built with G . III. DIAGRAMMATIC MONTE CARLOALGORITHMA. From diagrams to Monte Carlo
In this section, we explain how the diagrammatic ex-pansion of the previous section can be formally rewrittenas a stochastic average. As in Refs. 1, 3, and 49, thegeneral idea is that the integral over internal variables X and the sum over topologies T will be evaluated stochas-tically, for all values of the external variables Y , througha single Monte Carlo process. Specifically, in order todetermine the function Q ( Y ), where Q stands as abovefor Σ ( N ) σ or Π ( N ) , we shall compute overlaps of the form A Q ,g ≡ (cid:90) dY Q ( Y ) g ( Y ) (25)for a set of functions g given below. Expanding Q ( Y ) interms of Feynman diagrams as in Eq. (20) yields A Q ,g = (cid:88) T ∈S Q (cid:90) dXdY D ( T , X, Y ) g ( Y ) . (26)Defining a configuration by C = ( T , X, Y ) , (27)i.e. by a given topology and given values of internal andexternal variables, the expression (26) can be rewrittenas a weighted average over configurations, A Q ,g = (cid:90) d C |D ( C ) | · sgn[ D ( C )] · g ( Y ) · T ∈S Q . (28)Here we introduced the indicator function1 T ∈S Q = (cid:40) , if T ∈ S Q , , otherwise , so that the integral over C can be extended to topolo-gies outside of S Q . Our choice of the extended space ofconfigurations will be discussed below.In order to evaluate (28) by Monte Carlo, it should berewritten in the form A Q ,g = (cid:90) d C w ( C ) A Q ,g ( C ) , (29)where w ( C ) ≥ Z ≡ (cid:90) d C w ( C ) (30)is finite so that w ( C ) / Z is a normalized probability dis-tribution. In practice we take w ( C ) = |D ( C ) | R ( C ) , (31)where R ( C ) is an arbitrary (non-negative) reweighingfunction. It is then clear that Eq. (28) can indeed berewritten as Eq. (29) provided we set A Q ,g ( C ) = sgn[ D ( C )] · g ( Y ) · T ∈S Q R ( C ) . (32)The Monte Carlo update scheme (described in Sec. III E)will generate a Markov chain of random configurations FIG. 1. Examples of diagrammatic topologies from the dif-ferent sectors of the Monte Carlo configuration space: (a) Σ-sector, (b): Π-sector, (c): Worm sector, (d): Normalizationsector. The dashed black line is the measuring line, whichhas the structure of a one-body propagator in the Σ-sector,and of a pair propagator in the Π-sector. In the Worm sector,the worm ends are represented by red dots connected with anextra unphysical thread (dashed red line). C , C , . . . with the stationary probability distribution w ( C ) / Z . The average over n generated configurationsthen converges to the true expectation value in the large n limit, A Q ,g = Z × lim n →∞ n n (cid:88) i =1 A Q ,g ( C i ) . (33)It remains to estimate Z , which can be done easily inthe following way. The trick is to have a subset S N ofthe configuration space, which we call the normalization-sector, whose total weight (cid:90) S N d C w ( C ) =: Z N , (34)is easy to calculate analytically. In our case, we artifi-cially create this normalization sector by enlarging theconfiguration space, as we shall see below (in contrast, inthe case of the Hubbard model algorithm of Ref. 1, thenormalization sector was that of the first-order diagram).Defining the “norm” N as the number of times that thenormalization sector was visited, N ≡ n (cid:88) i =1 C i ∈S N , (35) Z can be evaluated thanks to Z N / Z = lim n →∞ N /n .Inserting this into Eq. (33) yields the final expression A Q ,g = Z N lim n →∞ (cid:80) ni =1 A Q ,g ( C i ) N . (36) B. Configuration space
To be more specific, the allowed diagram topologies, T , belong to one of the following sectors (see Fig. 1 forexamples): • Σ-sectors ( S Σ ( N ) σ ): A self-energy diagram of order N contains N pair propagators Γ, N − G σ , and N single-particlepropagators G − σ . The open ends of the diagramare formally closed with some extra unphysical linewhich has the structure of a single-particle propaga-tor of spin σ . We refer to this line as the measuringline. • Π-sectors ( S Π ( N ) ): A pair self-energy diagram oforder N contains 2 N single-particle propagators G , N − • Worm sector: In addition to the above physical di-agrams, we also consider unphysical diagrams con-taining two vertices where the momentum conser-vation is not fulfilled. We will refer to these verticesas Worms, named Ira ( I ) and Masha ( M ). The mo-mentum conservation at I and M is restored if weconsider that a momentum δ is flowing from Irato Masha along some extra unphysical thread. Inthis sector, T includes the location of I and M ,while the momentum δ is included in the internalvariables X . • Normalization-sector ( S N ): The topology and vari-ables of the normalization diagram are the ones ofa fully closed N = 1 diagram. The lines in this dia-gram are certain “designed” simple functions ratherthan G and Γ propagators. C. Probability density
In order to precisely define the probability density w ( C ) d C on the above configuration space, we first specifywhat we mean by d C . For any function f ( C ), we set (cid:90) d C f ( C ) ≡ (cid:88) T (cid:90) dXdY f ( T , X, Y ) , (37)where dY = d p dτ dτ and dX depends on the topology T : it is given by Eq. (21) if T has no Worms, and by thesame expression with an additional factor d δ if T has apair of Worms.Alternatively, one can discretize the configurationspace. We emphasize that this introduces arbitrarilysmall discretization steps, which is really fundamentallyequivalent to working with continuous variables. In thiscase, all momentum coordinates and imaginary time areinteger multiples of some arbitrarily small δp and δτ . Wecan write, for topologies without Worms, (cid:90) dX dY f ( T , X, Y ) ≡ (cid:88) τ ,...,τ N δτ N (cid:88) ( p ,..., p N ) δp N +1) f ( T , X, Y ) , (38) where the sum over the momenta ( p , . . . , p N ) of alllines (internal and measuring) is constrained by the mo-mentum conservation at each vertex. For topologies withWorms, (cid:90) dX dY f ( T , X, Y ) ≡ (cid:88) τ ,...,τ N δτ N (cid:88) ( p ,..., p N ) δp N +1) (cid:88) δ δp f ( T , X, Y ) . (39)We then have (cid:90) d C f ( C ) ≡ (cid:88) C δ C f ( C ) , (40)where δ C is the “volume” of one discrete “cell” of theconfiguration space around the considered point C . Moreprecisely, if C is an order- N diagram, δ C is given by δτ N δp N +1) , multiplied by the additional factor δp ifthe Worms are present. A nice feature of the formula-tion (38,39) is that there is no need to introduce the loopmomenta; instead, all momenta are treated on equal foot-ing, which is also the case in the diagrammatic MonteCarlo code.We define the weighting function w ( C ) in the followingway. For a physical configuration, we take w ( C ) = |D ( C ) | R ( C ) , (41)where D ( C ) is given by the Feynman rules (see subsec-tion II D), and R ( C ) is an arbitrary non-negative reweigh-ing function. We take R ( C ) = W Q meas ( p ) O N , (42)where Q is equal to Σ or Π depending on the sector, W Q meas ( p ) is the weight of the measuring line, and O N is an order-dependent reweighting factor. We choose W Σmeas ( p, τ ) to be ∝ /p for intermediate momenta (tocompensate for the Jacobian), constant for small mo-menta (to avoid having rare configurations with a largeweight), and ∝ /p for large momenta. This is just oneof the many possible choices, subject to the conditionthat sampling of diagrams with large p has to be sup-pressed in order to have a normalizable distribution (i.e.the total weight Z has to be finite). For the Π-sector,the simplest option is W Πmeas = φ Π W Σmeas where φ Π is anoptimization factor controlling the relative weights of theΣ and Π sectors.The weight of unphysical configurations (belonging tothe Worm sector or to the normalization sector) is de-fined as follows. Formally, the weight of configurationscontaining Worms is arbitrary, since they do not con-tribute to the self-energy. These diagrams are auxiliaryand are only employed for obtaining an efficient updat-ing scheme. In order to have a good acceptance ratiowhen moving between the physical and Worm sectors wechoose the weights according to the Feynman rules for allpropagator lines, with the extra rule that the unphysicalthread contributes to w ( C ) a factor C ( δ ), i.e. w ( C ) = |D ( C ) | R ( C ) C ( δ ) . (43)The C ( δ ) function should be chosen to decay fast enoughat large δ to ensure that Z is finite and includes a con-stant prefactor to optimize the relative statistics of sam-pled diagrams with and without the Worms.In the normalization sector, w ( C ) is a simple expressionsuch that one can easily calculate analytically the totalweight of the normalization sector S N , Z N = (cid:90) S N d C w ( C ) . (44)We take w ( C ) = G N ( p ↑ , − τ ) G N ( p ↓ , − τ ) Γ N ( p ↑ + p ↓ , τ )with G N ( p, τ ) = exp( − p σ N ) and Γ N ( p, τ ) = φ N . Theparameters σ N and φ N can be freely chosen and opti-mized. D. Measuring
We recall that we determine the function Q = Σ ( N ) σ or Π ( N ) by computing its overlaps with a set of func-tions g . We now describe our specific choices of func-tions g ( p, τ ). We divide the space of all ( p, τ = τ − τ )into bins B = B p × B τ ⊂ [0 , p max ] × [0 , β ]. In practice, τ = τ − τ lies in the interval [ − β, β ], but thanks to the β -(anti-)periodicity of Π (Σ) we only need to consider τ ∈ [0 , β ]. In each bin B we define the ortho-normal setsof basis functions u k ( p ) and v l ( τ ) satisfying (cid:90) B p d p w( p ) u k ( p ) u k (cid:48) ( p ) = δ k,k (cid:48) (45) (cid:90) B τ dτ v l ( τ ) v l (cid:48) ( τ ) = δ l,l (cid:48) (46)where w( p ) >
0. Then the to-be-determined function Q can be expanded in the bin B as Q ( p, τ ) = (cid:88) k,l Q k,l u k ( p ) v l ( τ ) . (47)Setting g ( Y ) = 1 ( p,τ ) ∈B w( p ) u k ( p ) v l ( τ ), we obtain theexpansion coefficients, Q k,l = (cid:82) dY Q ( Y ) g ( Y ), by MonteCarlo as explained in Subsec. III A.We take the w( p ) function in the inner product to bew( p ) = 1 / (4 πp ) (except for the lowest bin, see below)so that the u k ’s and v l ’s can be chosen in the form ofLegendre polynomials up to the order 2. The procedurebecomes exact only in the limit of vanishing bin-size, butone can afford relatively large bins compared to the casewhen the function Q is approximated by a constant ineach bin (which would correspond to restricting to thepolynomial of order 0). In the lowest momentum-bin,we chose w( p ) = 1 / (4 π ). The reason for this choice is to avoid having a factor 1 /p in the right-hand side ofEq. (32), which would lead to huge contributions fromrare configurations with small p . The corresponding basisset of two functions is built from a constant and p .This choice ensures that for each consid-ered g , not only the mean value A Q ,g ( C ) =lim n →∞ n (cid:80) ni =1 A Q ,g ( C i ) is finite, but the corre-sponding variance [ A Q ,g ( C )] − [ A Q ,g ( C )] is also finite,where [ A Q ,g ( C )] = lim n →∞ n (cid:80) ni =1 [ A Q ,g ( C i )] . Thisfollows from the fact that for each bin B (including thespecial case of the lowest momentum-bin), [ A Q ,g ( C )] isbounded, because 1 /W Q meas ( p ) and g ( Y ) are bounded. E. Updates
Our updating scheme shares a number of features withthe one introduced for the Hubbard model in Ref. 1. Tosample the space of configurations with a variable num-ber of continuous variables, we use a Metropolis algo-rithm, with pairs of complementary updates. In ad-dition to the complementary pairs, a number of self-complementary updates are used. While not changingthe number of continuous variables, self-complementaryupdates allow us to efficiently sample diagram topologies.In this Section we present details of our specific imple-mentation of all Monte Carlo moves including expressionsfor their acceptance ratios. The updates presented inIII E 1, III E 2 and III E 3 suffice to perform the integra-tion over internal momenta and times while keeping theorder and topology of the diagram fixed. The updatesof III E 4 and III E 5 change the topology of the diagramwithout changing the order. The updates presented inIII E 6 and III E 7 allow one to change the diagram order.Finally the update of III E 8 allows to enter and leave thenormalization sector.
1. Create–Delete
In the complementary pair of updates
Create–Delete (see Fig. 2), a pair of Worms is created or deleted in thecurrent diagram. These updates are called with constantprobabilities p crt and p dlt , respectively. Delete (resp.
Create ) can only be called when the Worms are present(resp. absent). In
Create , a line is first chosen at random(i.e. with probability 1 / (3 N ) with N the order of thediagram). The chosen line can be a G σ , Γ or measur-ing line. Next, we choose with equal probability eitherIra or Masha to be located at the origin of the chosenline. An unphysical thread running from Ira to Mashais introduced and carries a momentum δ , chosen withprobability density W ( δ ). Note that, to optimize the ac-ceptance ratio, W is in principle allowed to depend on theimaginary-time difference between the ends of the chosenline (or any other configuration parameter).When the Worms are deleted, we first check whetherthere is at least one line connecting Ira and Masha. If FIG. 2. (color online) Graphical representation of the com-plementary pair of updates
Create-Delete . Only relevant (i.e.updated) parts of the Feynman diagrams are drawn. An un-physical thread (red dashed line) running from Ira to Mashaand carrying momentum δ is added to the graph in Create ,and is removed in the inverse update. Momentum conserva-tion is maintained when taking into account both the physicalpropagators and the unphysical thread. so, the Worms can be deleted with the inverse
Create ac-ceptance ratio. In case there is more than one line con-necting Ira and Masha, we choose one of these lines withequal probability 1 /N links where N links is the number ofconnections. When a G σ -line is chosen at the beginningof the update, the total acceptance ratio becomes q create = p dlt p crt NN links W ( δ ) | G σ ( p ± δ , τ ) | C ( δ ) | G σ ( p , τ ) | . (48)When the chosen line is of the Γ- or measuring-type,the acceptance ratios are constructed similarly. Recallthat C ( δ ) is an extra factor assigned to the diagram withWorms [see Eq. (43)]. The new momentum of the chosenline, p + δ or p − δ , depends on whether Ira is created atthe end or the origin of the line. Finally, we are left withthe choice for the probability density W ( δ ). We simplytake W ( δ ) ∝ C ( δ ).
2. Move In Move (see Fig. 3), one of the Worms is moved fromone three-point vertex to another along a single line. Thisline is chosen at random (i.e., with probability 1 /
3) andone has to ensure that Ira and Masha will not be placedon the same vertex (note that for this reason,
Move is im-possible for N = 1). Move is called with constant prob-ability p move , whenever Worms are present. If a Wormhappens to move along a G σ -line, the acceptance ratio is q move = | G σ ( p ± δ , τ ) || G σ ( p , τ ) | . (49) FIG. 3. Graphical representation of the self-complementaryupdate
Move . In this particular case, Ira is moved from onevertex to another along a G σ -propagator line. To preserve mo-mentum conservation, the momentum of this line is changed. The sign depends on the direction in which Ira or Mashais moved, and is chosen such that momentum conser-vation is preserved. As a result of
Move updates, Iraand Masha perform a random work within the Feynmangraph, while constantly updating line momenta.
3. Shift in time
This update can be called in every sector with con-stant probability p shift . It shifts the time variable of therandomly selected three-point vertex from τ to τ (cid:48) . Thenew variable τ (cid:48) is drawn from a distribution W ( τ (cid:48) ) on the(0 , β ) interval. The acceptance ratio is (here it is givenfor the case when all three lines attached to the shiftedvertex are physical propagators) q shift = W ( τ ) W ( τ (cid:48) ) · | G ↑ ( q ↑ , τ (cid:48) − τ ) G ↓ ( q ↓ , τ (cid:48) − τ )Γ( q , τ − τ (cid:48) ) || G ↑ ( q ↑ , τ − τ ) G ↓ ( q ↓ , τ − τ )Γ( q , τ − τ ) | , (50)where we have assumed that the propagators G σ are in-coming (see Fig. 4). We choose a seeding function W tak-ing into account the short-time behavior Γ( k, τ ) ∝ / √ τ [see Eqs. (62,64) below]: W ( τ (cid:48) ) = 1 / (2 √ β ∆ τ ) where∆ τ ≡ τ − τ (cid:48) for τ > τ (cid:48) , and ∆ τ ≡ τ − τ (cid:48) + β for τ < τ (cid:48) .
4. Reconnect
This self-complementary update changes the topol-ogy of the diagram without changing the order. It iscalled with constant probability p rec whenever Worms FIG. 4. Graphical representation of the self-complementaryupdate
Shift . The time τ of a three-point vertex is shifted to τ (cid:48) . are present and N >
1. The basic idea is that two spin- σ single-particle propagators that both leave from (or botharrive at) the Worm-vertices are reconnected, i.e. theirend-points are exchanged. This does not cause a prob-lem with momentum conservation: Only the unphysicalmomentum δ running from Ira to Masha changes. Reconnect is constructed as follows (see Fig. 5). First,we choose with equal probability whether to reconnectthe spin-up lines or the spin-down lines. These propaga-tors should be both arriving at or both leaving from theWorms, otherwise the update is rejected. In case theyboth arrive at the Worms, the acceptance ratio is q reconnect = | G σ ( q , τ M − τ ) G σ ( p , τ I − τ ) || G σ ( q , τ I − τ ) G σ ( p , τ M − τ ) |· C ( δ + p − q ) C ( δ ) . (51)
5. Swap measuring line
This update converts the measuring line into a realpropagator, while some other line becomes the new mea-suring line (see Fig. 6). Although very simple, this up-date changes the diagram topology and the values of in-ternal and external variables. The update is only calledin the Σ and Π sectors, since it is not useful in the Wormor normalization-sector. The update starts with choosingone of the lines at random (it should not be the measuringline). This line is proposed to become the new measuringline. The acceptance ratio is given by q swap = | Γ( q , τ ) | W Πmeas ( q , τ ) · W Σmeas ( q σ , τ (cid:48) ) | G σ ( q σ , τ (cid:48) ) | , (52)for the particular case which converts Π-sector to Σ-sector. For other cases, acceptance ratios are constructed FIG. 5. Graphical representation of the self-complementaryupdate
Reconnect . In this particular example, the two prop-agators incoming to the Ira and Masha-vertices are inter-changed. The momentum carried by the unphysical threadis changed from δ to δ + p − q .FIG. 6. Graphical representation of the Swap measuring line update. similarly.
6. Add–Remove
To add a pair-propagator line, the Worms should bepresent, and we should not be dealing with the normaliza-tion diagram. In this case, the update
Add is called withconstant probability p add . First, we choose the spin-upor spin-down line attached to the Ira-vertex. Let this linecorrespond to G σ . Next, we consider the opposite spinpropagator attached to the Masha-vertex, G − σ . Thesetwo propagators will be cut, and a new pair propagatorwill be inserted; see Fig. 7 for an illustration. The fi-nal diagram does not contain the Worms, which leavesus no freedom in choosing the momenta in the final di-agram. We propose initial and final times τ o and τ d forthe new pair-propagator line, from a probability density W ( τ o , τ d ), which in the current implementation is simply0the uniform distribution.For Remove , we need
N > p rm . Thepair-propagator line to be removed is chosen at random.If the topology of the diagram is such that the chosen Γ-line has the same G -propagator attached to its both ends,the update is immediately rejected since such a G -loopcannot be created through Add . An update trying to re-move a measuring Γ-line is also forbidden. Next, chooseone of the four lines attached to the pair-propagator lineat random. This will be the future G σ -line and the vertexit is connected to will become Ira. One of the remain-ing G − σ is also selected at random and the vertex it isconnected to will become Masha. If the same vertex ischosen for Ira and Masha, the move is rejected.The acceptance ratio for Add is q add = p rm p add O N +1 Γ( τ d − τ o )32 π ( N + 1) W ( τ o , τ d ) C ( δ ) O N G σ ( τ o − τ o ) G − σ ( τ o − τ o ) G σ ( τ d − τ d ) G − σ ( τ d − τ d ) G σ ( τ d − τ o ) G − σ ( τ d − τ o ) . (53)The momenta are omitted here for simplicity. There areseveral possibilities depending on the particular choice of G σ and G − σ and the positions of Ira and Masha. In allcases, however, the new momenta are completely deter-mined by the conservation laws. Fig. 7 shows a particularexample. If in Add the chosen propagator G σ (or G − σ )happens to be the measuring line, then a new measuringline will be chosen with equal probability among the twospin- σ propagators connected to Γ in the final diagram.The reverse is done in Remove .
7. Add–Remove loop
These updates are called in the Σ- and Π-sectors only.
Add loop (resp.
Remove loop ) is called with the probabil-ity p al (resp. p rl ). In Add loop a G σ propagator is chosenat random, and converted into the sequence G σ Σ (1) σ G σ where Σ (1) σ is the first-order self-energy diagram (Γ closedwith G − σ ), see Fig. 8 where we illustrate the setup for σ = ↓ . The initial and final times τ o and τ d for the newpair-propagator line are chosen from the probability den-sity W ( τ o , τ d ), and the momentum q ↑ for G ↑ is chosenfrom another distribution W ( q ↑ | τ o − τ d ). In Remove loop a pair-propagator line is first chosen at random. If thispropagator has the same G -line attached to its ends, thenit can possibly be removed by the update. If either the Γor G ↑ -line is the measuring line, the update is rejected.The acceptance ratio for Add loop is given by q add loop = p rl p al (2 N − O N +1 (2 π ) ( N + 1) W ( τ o , τ d ) W ( q ↑ | τ o − τ d ) O N · | G ↓ ( q ↓ , τ − τ d ) G ↓ ( q ↓ , τ o − τ ) || G ↓ ( q ↓ , τ − τ ) |· | Γ( q ↓ + q ↑ , τ d − τ o ) G ↑ ( q ↑ , τ o − τ d ) | , (54) FIG. 7. Graphical representation of the complementary pairof moves
Add-Remove . In this particular example, the G σ propagator in Eq. (53) corresponds to the spin- ↓ line carryingmomentum q ↓ , and G − σ is the spin- ↑ line with momentum q ↑ .FIG. 8. Graphical representation of the updates Add–Removeloop . with N the order of the diagram in Add loop . For W ( q ↑ | τ o − τ d ) we take a Gaussian distribution with vari-ance 1 /τ where τ = τ o − τ d for τ o > τ d and τ = τ o − τ d + β for τ o < τ d . This corresponds to the behavior of the vac-uum propagator G v .1
8. Swap to the normalization diagram
For normalization purposes, we introduce an unphysi-cal diagram, for which all integrals can be evaluated an-alytically (see subsections III B and III C). If the currentdiagram is the one-body self-energy diagram of order one,the
Norm update proposes to swap to the normalizationdiagram. The acceptance ratio is q Norm = G N ( p ↑ , − τ ) G N ( p ↓ , − τ )Γ N ( p ↑ + p ↓ , τ ) | W Σmeas ( p σ , − τ ) G − σ ( p − σ , − τ )Γ( p ↑ + p ↓ , τ ) | . (55)When the current diagram is the normalization dia-gram, Norm proposes to swap back to the physical self-energy diagram with the probability given by the inverseof Eq. (55).
F. Reducibility and ergodicity
The goal of our Monte Carlo setup is to sample thespace of one-body and two-body self-energy skeleton di-agrams in an ergodic way. These diagrams are connected,irreducible with respect to cutting a single G -propagatoror Γ-propagator, and irreducible with respect to cuttingany two G σ -propagators or any two Γ-propagators. Theset of updates presented in Section III E suffices to gener-ate this class of diagrams. In principle, the scheme couldbe used to generate a bigger class of diagrams (e.g., allconnected diagrams), but we focus the discussion here onsampling the skeleton diagrams only.Some of the updates of Section III E can propose to gofrom a skeleton diagram to a non-skeleton one. One pos-sibility is that all such proposals are simply rejected. Thisimmediately creates a problem with ergodicity: Sincethere is no skeleton diagram at order 2 and the diagramorder can only be changed by one, the simulation wouldnever leave the first-order diagram. Allowing some non-skeleton diagrams at orders 2 and 3 solves the problemand is sufficient for ergodicity (obviously, non-skeletondiagrams are excluded from the measurements). Beyondorder 3, we restrict sampling to skeleton diagrams onlywithout violating the ergodicity requirement.Explicitly checking the topology of high-order dia-grams at each update would be very time-consuming.Instead, our connectivity and reducibility checks rely onmomentum conservation. Let us start with discussing theconnectedness of the generated diagram. It is easy to seethat, by construction, the only moves that can possiblygenerate disconnected pieces are Reconnect and
Remove .The latter update, however, can only create a discon-nected piece if the initial diagram is not a skeleton dia-gram (since this diagram falls apart when cutting two G σ lines connected to the Γ-line that is removed). Reconnect ,on the other hand, can generate two disconnected piecesin the Worm sector starting from a skeleton diagram. Wesimply reject the update when this happens, which canbe straightforwardly done in the following way. When
N Σ σ [ G, Γ] Π[ G, Γ] Σ σ [ G (0) , Γ (0) ] Π[ G (0) , Γ (0) ]1 1 1 1 02 0 0 1 23 1 1 5 64 4 4 25 305 23 23 161 1866 168 168 1201 13627 1384 1384 10181 113828 12948 12948 96265 106446TABLE I. Number of diagram topologies contributing to theone-body self-energy Σ σ and two-body self-energy Π. In ad-dition to the number of skeleton diagrams built with G andΓ (first and second column), we also give for comparison thenumber of diagrams built with G (0) and Γ (0) (third and fourthcolumn). two disconnected pieces are generated by Reconnect , theWorms will be located on two three-point vertices whichare part of these two pieces. Due to momentum conser-vation, δ = 0. In this case, we reject the update.To test the topology of the diagram, we keep momentaof all lines in a hash table. The key point is that a di-agram has an irreducible skeleton topology if and onlyif no pair of lines (irrespective of the type of line: G , Γ,or measuring line) can have exactly the same momentum(or momenta which differ by ± δ in a Worm sector) withfinite probability. Indeed, such a pair of lines can only ex-ist if the two lines are of the same type and if the diagramfalls apart when cutting these two lines. The hash tableallows to find equal momenta in just a few operationsfor a sufficiently fine mesh in the hash table. Whenevera momentum of a line is changed, the hash table is up-dated. In each update, we ensure that the final diagramwill be of the skeleton type. Note that many of the up-dates cannot, by construction, result in disconnected ornon-skeleton topology, and we only check these proper-ties when there is a possibility that such a topology willbe created. For example, when adding a pair-propagatorin Add , there is only one way in which the diagram canbecome non-skeleton: when the final diagram falls apartby cutting the added pair-propagator and another line.This means that the added pair-propagator will be hav-ing the same momentum than another line.We have checked ergodicity explicitly using a dedicatedprogram which enumerates all topologies. In practice, weran these checks up to order 8. As a byproduct we get thenumber of topologies at each order, given in Table I. As mentioned earlier, we allow some non-skeleton dia-grams at order 2 and 3 to ensure ergodicity, namely theone-particle irreducible diagrams without ladders (theirnumber is also given in Table I). For this reason we haveintroduced the moves
Add Loop and
Remove Loop thatadd and remove loops. These updates should not becalled if the final (initial) order is bigger than 3.2
G. Bold diagrammatic Monte Carlo iterativescheme
The self-consistent nature of BDMC implies that thecalculation is performed iteratively. Starting from thepropagators G and Γ (for the first iteration, they are justsome initial guess), the self-energies Σ and Π are cal-culated by diagrammatic Monte Carlo. They are usednext in the Dyson equations to compute new values ofthe propagators, and the simulation continues with up-dated propagator lines. After a large enough number ofiterations, the process converges.A useful trick to accelerate this convergence is to per-form a weighted average over different iterations. Moreprecisely, the self-energy Σ j that we plug into the Dysonequation after iteration j is a weighted average of theMonte Carlo result of iteration j , and of Σ j − . The cor-responding weighting coefficients can be optimized to ob-tain small statistical errors as well as fast j -dependentconvergence.We have used the following weighting coefficients. TheMonte Carlo result of iteration j is an unnormalized his-togram ¯Σ ( h ) j and a norm ¯ N j ; instead of estimating theself-energy as Σ j = ¯Σ ( h ) j / ¯ N j , we use Σ j = Σ ( h ) j / N j withΣ ( h ) j = ¯Σ ( h ) j +(1 − f j )Σ ( h ) j − and N j = ¯ N j +(1 − f j ) N j − . Using a non-zero value of f j suppresses the contributionof older iterations, leading to a faster convergence. Inthe long-time limit j → ∞ , the statistical error still tendsto zero provided f j →
0. We used f j ∝ /j . The sameprocedure was applied for Π.The final error for each observable (density, contact,etc.) was estimated conservatively from its fluctuationsas a function of the iteration number j . More precisely,for a total number J of iterations, we estimated the erroras the maximal deviation between the final result afteriteration J and all intermediate results after the itera-tions j ∈ [ J/ , J ]. This automatically takes into accountthe combined effects of the statistical errors and the errordue to the finite number of self-consistency loops. H. Resummation
If the diagrammatic series was convergent, we wouldsimply have Q = lim N max →∞ N max (cid:88) N =1 Q ( N ) , (56)where Q stands either for the single-particle self-energyΣ or for the pair self-energy Π, Q ( N ) is the total con-tribution of the N -th order diagrams, and where it isimplicit that we consider arbitrary fixed values of the ex-ternal variables ( p, τ ). However, the series is divergent,as shown analytically in Refs. 12 and 31. To overcomethis difficulty, we employ a divergent-series-resummation method, of the form Q = lim N max →∞ N max (cid:88) N =1 R ( N max ) N Q ( N ) (57)where the R ( N max ) N are appropriate coefficients, cor-responding to a conformal-Borel transformation, seeRefs. 12 and 31.In practice, a full BDMC calculation must be per-formed for each value of N max , and the result is extrap-olated to N max → ∞ . This implies that the Q ( N ) arethemselves N max -dependent, and are assumed to tend tothe exact Q ( N ) when the N max → ∞ limit is taken. IV. ULTRAVIOLET PHYSICS
Zero-range interactions lead to a characteristic ultra-violet asymptotic behavior governed by the so-calledcontact.
This physics is expressed in a naturalway within the bold-line diagrammatic framework, aswe explain in subsection IV A (related discussions withinthe T -matrix approximation can be found in Refs. 55–58). Analytical understanding of the ultra-violet behav-ior is readily incorporated into our BDMC scheme, as de-scribed in subsection IV B. A short description of thesepoints was given in Ref. 13. A. Large-momentum analytics
1. The contact
The momentum distribution of the resonant gas hasthe power-law tail n σ ( k ) ∼ k →∞ C k . (58)In practice, this behavior holds for k much larger thanthe typical momentum k typ of the particles in the gas.(In the balanced unitary case, k typ is the maximum ofthe Fermi momentum and the thermal momentum.)In position space, the density-density correlation func-tion diverges at short distance as (cid:104) ˆ n ↑ ( r ) ˆ n ↓ ( ) (cid:105) ∼ r → C (4 π r ) . (59)An immediate consequence of the last equation is thatif one measures all the particle positions in a unit vol-ume, the number of pairs of particles whose interparticledistance is smaller than s is C s/ (4 π ) when s →
0; inthis sense, C can be viewed as a density of short-distancepairs. Furthermore, the contact can be directly expressed interms of the bold pair propagator C = − Γ( r = , τ = 0 − ) . (60)3This expression is analogous to the expression Eq. (4)of the single-particle density n in terms of the single-particle propagator G , which shows again that C controlsthe density of short-distance pairs. While Eq. (60) wasfirst obtained within the T -matrix approximations, it is actually an exact relation in terms of the fully dressedΓ.
2. Bold propagators at large momentum
At large momentum, the bold propagators can in somesense be replaced by vacuum propagators. More pre-cisely, when k → ∞ , G ( k, τ ) and Γ( k, τ ) become smallfor any τ in the interval ]0; β [, except in the narrow region0 < τ (cid:46) /k where G ( k, τ ) (cid:39) G v ( k, τ ) (61)Γ( k, τ ) (cid:39) Γ v ( k, τ ) (62)with G v ( k, τ ) ≡ − e − ( k / τ (63)Γ v ( k, τ ) ≡ − (cid:114) πτ e − ( k / τ . (64)This can be justified as follows. We first note that G (0) ( k, τ ) (cid:39) G v ( k, τ ) at large k , where we extend G v to negative times by β -antiperiodicity. To justify (61),we write ( G − G (0) )( k, τ ) = ( G (0) Σ G (0) )( k, τ ) + . . . = (cid:82) β dτ (cid:82) β dτ G (0) ( k, τ − τ )Σ( k, τ − τ ) G (0) ( k, τ ) + . . . .When k → ∞ , G (0) ( k, ∆ τ ) (cid:39) G v ( k, ∆ τ ) becomes a nar-row function of ∆ τ , so that the integrals over the internaltimes τ i are effectively restricted to narrow intervals ofwidth ∼ /k . This implies that G ( k, τ ) − G (0) ( k, τ )tends to zero uniformly in τ when k → ∞ .To derive (62), we first note that Γ (0) ( k, τ ) (cid:39) Γ v ( k, τ ) at large k and τ (cid:46) /k , as shown inAppendix A. Equation (62) then follows by writ-ing (Γ − Γ (0) )( k, τ ) = (Γ (0) Π Γ (0) )( k, τ ) + . . . = (cid:82) β dτ (cid:82) β dτ Γ (0) ( k, τ − τ )Π( k, τ − τ )Γ (0) ( k, τ ) + . . . .Again, when k → ∞ , the integrals over the internal times τ i are effectively restricted to narrow intervals, so thatΓ( k, τ ) − Γ (0) ( k, τ ) tends to zero uniformly in τ .We have also derived analytical expressions for G − G (0) at large momentum or short distance, which naturallydepend on the contact. These expressions are given inAppendix D and used in Appendix B.
3. Self-energy at large momentum
When k → ∞ , Σ σ ( k, τ ) becomes small for any τ in theinterval ]0; β [, except • for τ → + with 0 < τ (cid:46) /k , whereΣ σ ( k, τ ) (cid:39) Σ (+) σ ( k, τ ) (65) with Σ (+) σ ( k, τ ) ≡ − (cid:114) πτ n − σ e − ( k / τ (66) • for τ → β − with 0 < β − τ (cid:46) /k , whereΣ σ ( k, τ ) (cid:39) Σ ( − ) ( k, τ ) (67)with Σ ( − ) ( k, τ ) ≡ −C e − ( k / β − τ ) . (68)Furthermore, this behavior comes entirely from thelowest-order bold diagram Σ (1) σ . To justify these statements, let us first considerthe higher-order bold diagrams for Σ σ ( k, τ ). Theircontributions vanish uniformly in τ for k → ∞ . Indeed,they contain internal vertices, and at some of theseinternal vertices, a large momentum goes through andhence the integration over the internal time variableis restricted to a narrow range (because G and Γ arenarrow functions of imaginary time at large momentum,cf. Sec. IV A 2). We thus only need to consider thelowest-order bold self-energy diagram, represented inFig. 9. The momenta q and p of the G and Γ lines arerelated by momentum conservation, p = q + k . Thus,when k (cid:29) k typ , at least one of the momenta p and q hasto be (cid:29) k typ . Case 1: p (cid:29) k typ . Choosing q as the integrationvariable, we have Σ (1) σ ( k, τ ) = (cid:82) G − σ ( q , − τ ) Γ( p = k + q , τ ) d q/ (2 π ) . As discussed in Sec. IV A 2, Γ( p, τ )is small except in the relevant time-region 0 < τ (cid:46) /p where it can be replaced with Γ v . We further observethat the relevant values of q in the integral are (cid:46) k typ ,an assumption that will be justified a posteriori . Thisimplies that p (cid:39) k , and thus the relevant time-regionis 0 < τ (cid:46) /k . Therefore we can replace Γ( p, τ )with Γ v ( k, τ ) and G − σ ( q, τ ) with G − σ ( q, − ) = n − σ ( q ).Since the remaining integral over q gives us the particledensity n − σ , we arrive at the result (65,66). Finally, therelevant momenta in the integral for particle density are q (cid:46) k typ , which justifies the above assumption. Case 2: q (cid:29) k typ . We now choose p as theintegation variable, and write Σ (1) σ ( k, τ ) = − (cid:82) Γ( p , τ ) G − σ ( q = − k + p , β − τ ) d p/ (2 π ) . Ac-cording to Sec. IV A 2, G − σ ( q, β − τ ) is small exceptin the relevant time-region 0 < β − τ (cid:46) /q where itcan be replaced with G v ( q, β − τ ). We observe that therelevant values of p in the integral are (cid:46) k typ , whichimplies that q (cid:39) − k . Thus the relevant time-regionis 0 < β − τ (cid:46) /k , and we can replace G v ( q, β − τ )with − e − ( k / β − τ ) and Γ( p , τ ) with Γ( p , β − ). Theremaining integral over p gives us the contact, see (60),and we readily arrive at the result (67,68).4 FIG. 9. Lowest-order bold self-energy diagram, expressingΣ (1) in terms of G and Γ. This diagram contains the dominantcontributions to the self-energy at large momentum.FIG. 10. Leading diagrammatic contribution to the mo-mentum distribution n σ ( k ) at large k . The imaginary time isrunning from right to left. The single-particle lines propagateforward in time and can be replaced with the vacuum prop-agators. The pair propagator runs backwards in time and isfully dressed.
4. Tail of the momentum distribution
In short, the tail of the momentum distribution comesfrom the diagram depicted in Fig. 10, which can be in-terpreted physically as the simultaneous propagation oftwo opposite-spin particles of large and nearly oppositemomenta and of a missing pair with lower momentum.More precisely, for k → ∞ the Dyson equation simplifies: n σ ( k ) = G σ ( k , − ) = G (0) σ ( k , − )+ (cid:90) β − β dτ (cid:90) β − β dτ G (0) σ ( k , − τ )Σ σ ( k , τ − τ ) G σ ( k , τ ) (cid:39) (cid:90) β − β dτ (cid:90) β − β dτ G v ( k , − τ )Σ σ ( k , τ − τ ) G v ( k , τ ) . (69)Indeed, the ideal-gas momentum distribution decays ex-ponentially at large k so that we can neglect the term G (0) σ ( k , − ), and in the remaining term we can replace G with G v according to Subsec. IV A 2. Note that wetook the integration domain for the internal times τ and τ to be ] − β/ β/
2[ instead of the usual ]0; β [,which is allowed since the integrand is a periodic func-tion of τ and τ . As a result, the time-arguments ofthe G v factors never approach − β , and thus the G v canbe replaced by the retarded vacuum propagators, i.e. we have G v ( k , ∆ τ ) (cid:39) − θ (∆ τ ) e − ( k / τ for k → ∞ and∆ τ ∈ ] − β/ β/ θ ( . ) is the Heaviside function.Hence the integral is dominated by τ → + and τ → − ,and the imaginary time argument τ − τ of the self-energytends to 0 − . The asymptotic expression of Σ σ ( k , τ ) for k → ∞ , τ → − is known analytically, cf. Eq. (66). Aftersubstitution of this expression into the asymptotic Dysonequation given above, the large-momentum tail, Eq. (58),is recovered. B. Incorporating ultraviolet analytics into BDMC
A hallmark of BDMC is its unique capability to incor-porate analytical knowledge. The analytical considera-tions of the previous subsection have the following impli-cations for our BDMC calculation. Firstly, the contactcan be evaluated accurately from the bold pair propa-gator thanks to the relation Eq. (60), as was done inRef. 13. Furthermore, since the C /k tail of the mo-mentum distribution comes exclusively from the lowest-order self-energy diagram, this tail is automatically builtinto our self-consistent BDMC scheme provided this di-agram is evaluated with high precision. We achieve thisby using numerical Fourier transformations (rather thanMonte Carlo) and analytical treatments of leading-ordersingularities, in the spirit of Ref. 48, see Appendix Bfor details. As a result, in the BDMC data for themomentum distribution, the C /k tail is automaticallypresent and free of k -dependent noise. Note that here, C comes from the fully dressed pair propagator Γ, givenby the BDMC self-consistency which includes higher-order contributions; hence C differs from the one of theself-consistent T -matrix approximation of Refs. 48 and61. On the technical side, we mention that treating thelowest-order self-energy diagram separately (without us-ing Monte Carlo) has another advantage: the steep func-tions of τ in Eqs. (66,68) would be hard to capture byMonte Carlo sampling. V. LADDER SCHEME
As an alternative to the bold scheme discussed above,we also employ a partially dressed scheme, in which dia-grams are built from the bare single-particle propagator G (0) and the partially dressed pair-propagator Γ (0) , de-fined as the sum of ladder diagrams built with G (0) , seeEq. (9). For simplicity we will refer to this as the “lad-der scheme” (the ladder summation being the minimaldressing procedure allowing to work with zero-range in-teractions in continuous space). While the first diagramsof the ladder series for the single-particle self-energy Σare given by Eq. (12) above, the ones for the pair self-5energy Π are (70)A drawback of the ladder scheme is that it can onlybe used for temperatures above the approximate criticaltemperature T (0) c at which Γ (0) ( P = 0 , Ω n = 0) diverges.In the region T (0) c > T > T c , Γ (0) has a pole at finitemomentum so that the ladder scheme cannot be used. For the ladder scheme, no self-consistent iterations areneeded, which implies several advantages over the boldscheme: The ladder scheme is more practical for nu-merical computations; the justification of the conformal-Borel resummation method is more solid for the ladderscheme ; in particular, the ladder scheme is not sub-ject to the misleading-convergence problems that maypotentially affect the bold scheme (we also note thatmisleading convergence was observed in Ref. 63 only forfillings near one atom per lattice site, which is a regimevery different from the zero-filling limit corresponding tothe present continuous-space model). A. Dyson equations
In the ladder scheme, it is useful to consider the dia-grammatic series not only for the self-energies Σ and Π,but also for the propagators G and Γ. In this Section, letus denote by Σ ( N ) , Π ( N ) , G ( N ) and Γ ( N ) the sum of allorder- N diagrams in the ladder scheme for Σ, Π, G and Γ(the number of Γ (0) -lines in such diagrams is respectively N , N − N and N +1; accordingly the number of G (0) -lines is respectively 2 N −
1, 2 N , 2 N +1 and 2 N ). Notethat Π (1) = 0 (since all ( G (0) G (0) ) bubbles are alreadycontained in Γ (0) ).From the Dyson equations G σ ( p, ω n ) = ( G (0) σ + G (0) σ Σ σ G σ )( p, ω n ) (71)Γ( p, Ω n ) = (Γ (0) + Γ (0) Π Γ)( p, Ω n ) , (72) we have the order-by-order Dyson equations G ( N ) = N (cid:88) M =1 G (0) Σ ( M ) G ( N − M ) (73)Γ ( N ) = N (cid:88) M =1 Γ (0) Π ( M ) Γ ( N − M ) (74)for 1 ≤ N ≤ N max .As in the bold case, we need to apply a resummationprocedure to extract a result from the divergent diagram-matic series. The first way to do so is to proceed exactlyas in the bold case (Subsec. III H above): apply the re-summation procedure to Σ and Π, and plug the resultinto the Dyson equations (71,72) to get G and Γ. An-other way is to apply the resummation procedure to theseries (cid:80) N G ( N ) and (cid:80) N Γ ( N ) , i.e., to use Q = Q (0) + lim N max →∞ N max (cid:88) N =1 R ( N max ) N Q ( N ) (75)with Q = G or Γ. B. Ultraviolet physics
In the ladder scheme, the accurate incorporation ofultraviolet physics is more involved than in the bold case.Recall that Σ( p, τ ) and Π( p, τ ) are narrow functionsof τ when p is large. This would be difficult to captureby Monte Carlo. Our solution for the bold code wasvery simple: Given that this singular behavor is com-pletely contained (at leading order) in the lowest-orderbold diagrams, we compute these diagrams by Fouriertransformation rather than by Monte Carlo.For the ladder scheme, we have to do some extra workin order to achieve the same goal. Let us denote (in thepresent subsection) the lowest-order bold diagrams byΣ ,bold and Π ,bold . The problem is that these bold di-agrams contain an infinite number of ladder-scheme di-agrams. Our solution is as follows: During the MonteCarlo process, we do not measure the (ladder-scheme)diagrams that contribute to Σ ,bold and Π ,bold . Instead,we compute them by combining Fourier transformationwith order-by-order Dyson equations.More precisely, sinceΣ ,bold ; σ ( r, τ ) = Γ( r, τ ) G − σ ( r, − τ ) (76)Π ,bold ( r, τ ) = − G ↑ ( r, τ ) G ↓ ( r, τ ) , (77)we haveΣ ( N )1 ,bold ; σ ( r, τ ) = N (cid:88) M =1 Γ ( M − ( r, τ ) G ( N − M ) − σ ( r, − τ ) (78)for 1 ≤ N ≤ N max , andΠ ( N )1 ,bold ( r, τ ) = − N − (cid:88) M =0 G ( M ) ↑ ( r, τ ) G ( N − − M ) ↓ ( r, τ ) (79)6for 2 ≤ N ≤ N max . Here, Σ ( N )1 ,bold ; σ and Π ( N )1 ,bold denotethe sum of all ladder-scheme diagrams of order N thatare part of the lowest-order-bold diagram. The diagramscontributing to Σ ( N )1 ,bold ; σ up to N = 3 are the ones inEq. (12), except for the last diagram in Eq. (12) which isnot part of Σ (3)1 ,bold ; σ . Similarly, the diagrams contributingto Π ( N )1 ,bold up to N = 3 are the ones in Eq. (70), except forthe last diagram in Eq. (70) which is not part of Π (3)1 ,bold .We can thus perform the computations recursively inthe following order:( G (0) , Γ (0) ) −→ . . . −→ ( G ( N − , Γ ( N − ) −→ (Σ ( N ) , Π ( N ) ) −→ ( G ( N ) , Γ ( N ) ) −→ . . . −→ ( G ( N max ) , Γ ( N max ) ) (80)where at each order, the self-energies are obtained byadding up the Monte Carlo contribution with the (1 , bold )contribution. C. Monte Carlo
The diagrammatic Monte Carlo algorithm for sam-pling the ladder series is similar to the bold case de-scribed above in Sec. III, with the following differences.The iterative procedure (Subsec. III G) is not requiredany more. The topologies which are reducible with re-spect to cutting two internal G (0) σ lines, or two inter-nal Γ (0) lines, are sampled and measured. Accordingly,we perform the momentum-comparison checks describedin Subsec. III F only between one internal line and themeasuring line to omit one-particle reducible diagrams.Lastly, the (1 , bold ) diagrams are not measured; they areidentified in a way similar to detecting whether a dia-gram is non-skeleton, except now we only check whetherthe diagram falls apart if we cut two specific lines (theinternal G (0) -lines which are connected to the externalthree-point vertices). VI. CONCLUSION AND OUTLOOK
For spin-1 / .While the first numerical results presented in Refs. 4,12, and 13 are restricted to the unpolarized unitary gas,we expect the approach to be direcly applicable to the po-larized gas throughout the BEC-BCS crossover, as wellas to the mass-imbalanced case. Extension to two dimen-sions also seems feasible, as already demonstrated for thepolaron problem . A similar scheme may be used tostudy the leading finite-range correction.Another direction is the development of new algo-rithms to perform the summation over diagrams. Ratherthan sampling stochastically topologies, one may sum ex-actly over all topologies at each Monte Carlo update.With the efficient summation strategy that was recentlyintroduced for the Hubbard model , one obtains abetter computational complexity than for the originalDiagMC . It can also be advantageous to perform thisexact summation by brute-force enumeration providedthe momentum and time variables are chosen appropri-ately, as sucessfully demonstrated very recently for theelectron gas . A radically different approach would beto work with Schwinger-Dyson equations, for which newalgorithms were introduced and applied to bosonic mod-els . ACKNOWLEDGMENTS
To the memory of Joe Babcock, who has played a cru-cial role at the UMass cluster over many years.
We thank R. Rossi and E. Kozik for a fruitful collab-oration on the closely related Refs. 4, 12, 13, and 31.We are grateful to R. Haussmann, who provided us, forcomparison, with unpublished propagator data obtainedfrom the lowest-order bold diagrams as in Refs. 48 and75. We thank G. Bertsch, A. Bulgac, E. Mueller, G.Shlyapnikov and S. Tan for comments.This work was supported by the Research Founda-tion - Flanders FWO (K.V.H.), ERC grants Thermo-dynamix and Critisup2 (F.W.), a PICS from CNRS(F.W., N.P. and B.S.), National Science Foundation un-der grant DMR-1720465, MURI Program “Advancedquantum materials – a new frontier for ultracold atoms”from AFOSR and the Simons Collaboration on the ManyElectron Problem (N.P. and B.S.). T.O. was supportedby the MEXT HPCI Strategic Programs for InnovativeResearch (SPIRE), the Computational Materials ScienceInitiative (CMSI) and Creation of New Functional De-vices and High-Performance Materials to Support NextGeneration Industries (CDMSI), and by a Grant-in-Aid for Scientific Research (No. 22104010, 22340090,16H06345 and 18K13477) from MEXT, Japan. We ac-knowledge the hospitality of the Institute for NuclearTheory, Seattle (INT-10-1, INT-11-1).7
Appendix A: Ladder diagrams
In this appendix, we give some useful analytical prop-erties of the pair propagator Γ (0) defined by the sum ofladder diagrams [Eq. (9)] and describe its numerical cal-culation in frequency domain.The expression of Γ (0) ( P , Ω n ) was given in Eq. (11).For Ω n (cid:54) = 0 or P / − µ > (0) ( P , Ω n ) = 1˜Γ ( P , Ω n )+ (cid:90) d k (2 π ) n (0) ↑ ( P / k ) + n (0) ↓ ( P / − k ) i Ω n + 2 µ − P / − k (A1)where 1˜Γ (0) ( P , Ω n ) = 14 π (cid:18) a − (cid:112) P / − µ − i Ω n (cid:19) (A2)and we take the convention that the real part of thesquare root is positive.In time-domain, we get (after transforming the summa-tion over Matsubara frequencies into a contour integralusing the residue theorem)˜Γ (0) ( P , τ ) = − √ τ e − ( P / − µ ) τ (cid:90) ∞ dx e − x − e − β ( P / − µ ) − ( β/τ ) x , (A3)where, for simplicity, we restricted the analysis to theunitary case a = ∞ , and assumed that P / − µ > P → ∞ , τ → + , P τ (cid:46)
1, we haveΓ (0) ( P, Ω n ) (cid:39) ˜Γ (0) ( P, Ω n ) (cid:39) Γ v ( P, τ ) . (A4)Indeed, in this limit, in the integrand in Eq. (A3),the denominator tends to 1, which yields ˜Γ (0) ( P, τ ) (cid:39) Γ v ( P, τ ), where Γ v is defined in Eq. (64); moreover, inthis same limit, we have Γ (0) ( P, τ ) (cid:39) ˜Γ (0) ( P, τ ), becausein the large-momentum large-frequency limit, we haveΓ (0) ( P, Ω n ) (cid:39) ˜Γ (0) ( P, Ω n ) by neglecting the Fermi fac-tors compared to unity in Eq. (11).In practice we numerically compute and tabulateΓ (0) ( P , Ω n ). We distinguish between Ω n = 0 and Ω n (cid:54) =0. For Ω n (cid:54) = 0 we can use the expression (A1,A2).The angular integration is done analytically, and oneis left with a one-dimensional integral which is evalu-ated numerically. When Ω n = 0, we have to use thefull expression Eq. (11), whose integrand does not di-verge, because 2 µ − P / − k = 0 implies that also1 − n (0) ↑ ( P / k ) − n (0) ↓ ( P / − k ) = 0. The angularintegration is again done analytically.For the ladder scheme, we also need Γ (0) ( P , τ ), whichwe obtain from Γ (0) ( P , Ω n ) using the procedure describedfor Γ at the end of App. C. Appendix B: First order diagrams
The lowest-order diagram for the one-body and two-body self-energy is evaluated separately (without MonteCarlo), in order to accurately capture the singular be-havior coming from the zero-range interaction. Our pro-cedure, described in detail in the following, is similar tothe one of Ref. 48, in that it uses Fourier transformationbetween momentum and position space, with analyticaltreatment of singular pieces.In position space, we simply haveΣ (1) ( r, τ ) = Γ( r, τ ) G ( r, − τ ) . (B1)To Fourier transform the propagators G and Γ frommomentum space to position space, we write them as G = G v + δG and Γ = Γ v + δ Γ, where G v and Γ v capturethe leading-order large-momentum short-time behaviorof G and Γ, see Eqs. (61,62,63,64); the Fourier transformof G v and Γ v is then done analytically while δG and δ Γare Fourier transformed numerically. Furthermore, to en-sure that the Fourier transformation δG ( k → r ) is doneaccurately, we have derived analytical expressions for theleading-order ultraviolet behavior of δG both in momen-tum and position space, see Appendix D.Finally, Σ (1) has to be Fourier transformed back fromposition to momentum space. We again single out sin-gular parts which we transform analytically. We rewriteEq. (B1) asΣ (1) ( r, τ ) = Γ v ( r, τ ) G v ( r, − τ ) + Γ v ( r, τ ) δG ( r, − τ )+ δ Γ( r, τ ) G v ( r, − τ ) + δ Γ( r, τ ) δG ( r, − τ ) . (B2)The Fourier transform to momentum space is done ana-lytically for the first term, and numerically for the lastterm. For the cross-terms (second and third term), wesingle out a singular piece whose Fourier transform tomomentum space is done analytically:Γ v ( r, τ ) δG ( r, − τ ) = Γ v ( r, τ ) δG ( r = 0 , − τ )+ Γ v ( r, τ )[ δG ( r, − τ ) − δG ( r = 0 , − τ )] . (B3)The first term in Eq. (B3) is indeed singular for τ → + and r →
0, where Γ v ( r, τ ) becomes a sharply peakedfunction of r . Its Fourier transform simply gives the con-tribution Γ v ( p, τ ) δG ( r = 0 , − τ ) to Σ (1) ( p, τ ). The sec-ond term in Eq. (B3) is Fourier transformed numerically.The second cross-term in Eq. (B2) is treated similarly,by writing it as δ Γ( r, τ ) G v ( r, − τ ) = δ Γ( r = 0 , τ ) G v ( r, − τ )+ [ δ Γ( r, τ ) − δ Γ( r = 0 , τ )] G v ( r, − τ ) . (B4)We note that one could think of the following al-ternative procedure: subtract the analytical singularpieces Σ (+) ( r, τ ) + Σ ( − ) ( r, τ ) from Σ (1) ( r, τ ), do theFourier transform to momentum space, and then addback Σ (+) ( p, τ ) + Σ ( − ) ( p, τ ). Actually, this alternative8procedure would be essentially equivalent to the previ-ous one, since we haveΣ (+) ( r, τ ) = Γ v ( r, τ ) G ( r = 0 , − ) (B5)Σ ( − ) ( r, τ ) = Γ( r = 0 , β − ) G v ( r, − τ ) . (B6)The first-order pair self-energy Π (1) is computed simi-larly, by going to position space, the singular pieces beingtreated analytically.Finally, we note that it is important to use an ap-propriate numerical treatment of the functions and theirultraviolet singularities (even when the leading singu-larities are subtracted and treated analytically). Simi-larly to Ref. 48, we used non-linear grids to tabulate thefunctions, and we computed the Fourier transforms us-ing spline-interpolation and analytical evaluation of theresulting integrals. Appendix C: Dyson equations
To calculate the propagator G ( q, τ ) from Σ( q, τ ), wefirst Fourier transform Σ( q, τ ) to the frequency represen-tation. When doing so, we single out the singular partsΣ (+) ( q, τ ) and Σ ( − ) ( q, τ ) given in Eqs. (66,68), whoseFourier transforms are done analytically:Σ (+) ( q, ω n ) = − π n − σ erf (cid:0) √ β (cid:112) q / − iω n (cid:1)(cid:112) q / − iω n , (C1)Σ ( − ) ( q, ω n ) = C e − βq / iω n + q / . (C2)This way, we take care not only of the high-momentumleading behavior of Σ, but also of the short-time behaviorof Σ at any momentum, which is given byΣ σ ( q, τ ) (cid:39) τ → + − n − σ (cid:114) πτ (C3)see Eqs. (65,66).The propagator G is then given in frequency repre-sentation by the Dyson equation Eq. (7). When Fouriertransforming this back to time representation, we treatanalytically the singular piece given by G (0) .To calculate the dressed pair propagator Γ, we firstfourier transform Π( p , τ ) to the Matsubara frequencyrepresentation, Π( p , Ω n ), and insert this into the Dysonequation Eq. (17) to obtain Γ( p , Ω n ).Finally we need to take the Fourier transform to thetime domain to get Γ( p , τ ). In order to suppress nu-merical errors in the form of oscillations in Γ( P , τ ) asa function of τ , we treat the large-frequency short-timeand large-momentum singular part analytically. Moreprecisely, we write Γ = ˜Γ v + ˜ δ Γ in the momentum-time domain, where ˜Γ v is a simple function capturing the ultraviolet behavior of Γ whose Fourier transform tomomentum-frequency domain is done analytically, while˜ δ Γ is Fourier transformed numerically. We take˜Γ v ( P, τ ) = − (cid:114) πτ e − ( P / − µ ) τ − (cid:114) πβ e − β ( P / − µ ) (cid:20) e β ¯ E ( P ) − (cid:21) e − ¯ E ( P ) τ (C4)where ¯ E ( P ) = Max( p / − µ, k / v ( P, Ω n ) = − π erf (cid:16)(cid:112) ( p / − µ − i Ω n ) β (cid:17)(cid:112) p / − µ − i Ω n + 4 (cid:114) πβ e − β ( p / − µ ) i Ω n − ¯ E ( P ) . (C5)In this way, we take care of leading and higher-order sin-gular parts of Γ at short time and large momentum. Appendix D: Ultraviolet asymptotics for G − G (0) In this Appendix, we give large-momentum and short-distance asymptotic expressions for G − G (0) . The deriva-tions being rather long, we only present the final re-sults, which we obtained from the diagram G (0) [Σ (+) +Σ ( − ) ] G (0) where Σ ( ± ) are the analytical large-momentumexpressions given in Eqs. (66,68).
1. Momentum space
At large momentum, we already know that G ( q, τ = β − ) (cid:39) −C /q . The generalization to τ ∈ ]0; β [ isgiven by the following expression, valid when τ or β − τ are (cid:46) /q : ( G σ − G (0) σ )( q, τ ) (cid:39) q →∞ δG a ( q, τ ) (D1)where δG a ( q, τ ) = [ δG ( − ) + δG (+) A + δG (+) B + δG (+) C ]( q, τ ) (D2)with δG ( − ) ( q, τ ) = − C q e − q ( β − τ ) (D3) δG (+) A ( q, τ ) = 16 √ π n − σ e − q τ q (cid:34) q √ τ + i √ π (cid:18) q τ + 1 (cid:19) erf (cid:18) i q √ τ (cid:19) e − q τ (cid:35) (D4)9 δG (+) B ( q, τ ) = C τ e − q τ q (D5) δG (+) C ( q, τ ) = C q e − q τ . (D6)
2. Position space
The large-momentum behavior δG a ( q, τ ) of G ( q, τ ) ob-tained above gives rise to a short-distance singular behav-ior of G ( r, τ ). In order to obtain analytical expressionsfor this position-space behavior, one essentially needs totake the Fourier transform of δG a ( q, τ ) from momentumto position space. However, this would lead to infrareddivergences. To avoid this problem, we introduce a func-tion ˜ δG ( q, τ ) which has the same large- q behavior than δG a ( q, τ ) and is properly regularized at low q . More pre-cisely, we define˜ δG ( q, τ ) = [ ˜ δG ( − ) + δG (+) A + ˜ δG (+) B + ˜ δG (+) C ]( q, τ ) (D7)with ˜ δG ( − ) ( q, τ ) ≡ δG ( − ) ( q, τ ) (cid:104) − e − ( q/q m ) (cid:105) (D8)˜ δG (+) B ( q, τ ) ≡ δG (+) B ( q, τ ) (cid:104) − e − ( q/q m ) (cid:105) (D9)˜ δG (+) C ( q, τ ) ≡ δG (+) C ( q, τ ) (cid:104) − e − ( q/q m ) (cid:105) (D10)where q m is a lower momentum cutoff whose precise valueis arbitrary ( e.g. , one can take q m = k typ ). These four terms have the following expressions in po-sition space:˜ δG ( − ) ( r, τ ) = C π r [ F ( X ) − F ( Y ) + F ( Z )] (D11)where F ( x ) = I ( x ) (cid:18) x (cid:19) + (cid:114) π e − x / x , I ( X ) = π (cid:18) X √ (cid:19) ,X ≡ r √ ∆ τ , Y ≡ r (cid:112) ∆ τ + 2 /q m , Z ≡ r (cid:112) ∆ τ + 4 /q m , and ∆ τ ≡ β − τ ; δG (+) A ( r, τ ) = 4 n − σ π (cid:90) ∞ dx e − x / (2 X (cid:48) ) i erf (cid:16) i x X (cid:48) (cid:17) × (cid:18) xx − sin xX (cid:48) (cid:19) (D12)with X (cid:48) ≡ r/ √ τ ,˜ δG (+) B ( r, τ ) = C τ π r I ( X (cid:48) ) − I (cid:113) X (cid:48) + q m r ) , (D13)˜ δG (+) C ( r, τ ) = − ˜ δG ( − ) ( r, β − τ ) . (D14) K. Van Houcke, E. Kozik, N. Prokof’ev, and B. Svis-tunov,
Diagrammatic Monte Carlo , in Computer Simula-tion Studies in Condensed Matter Physics XXI. CSP-2008.Eds. D.P. Landau, S.P. Lewis, and H.B. Sch¨uttler, PhysicsProcedia , 95 (2010). E. Kozik, K. Van Houcke, E. Gull, L. Pollet, N. Prokof’ev,B. Svistunov, and M. Troyer, EPL , 10004 (2010). N. Prokof’ev and B. Svistunov, Phys. Rev. B , 125101(2008). K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svis-tunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schi-rotzek, and M. W. Zwierlein, Nature Phys. , 366 (2012). A. J. Leggett, in: A. Pekalski, J. Przystawa (eds.),
ModernTrends in the Theory of Condensed Matter , p. 13. Springer,New York (1980). A. J. Leggett, J. Phys. (Paris) , C7 (1980). P. Nozi`eres and S. Schmitt-Rink, J. Low Temp. Phys. ,195 (1985). A. J. Leggett and S. Zhang, Lecture Notes in Physics ,33 (2012), in . The BCS-BEC Crossover and the Unitary Fermi Gas , Lec-ture Notes in Physics , W. Zwerger ed. (Springer, Hei-delberg, 2012). J. Carlson, S. Gandolfi, and A. Gezerlis, Prog. Theor. Exp.Phys. , 01A209 (2012). New J. Phys (2011), Focus on Strongly CorrelatedQuantum Fluids: from Ultracold Quantum Gases to QCDPlasmas , A. Adams, L. D. Carr, T. Schaefer, P. Steinberg,J. E. Thomas (eds.). R. Rossi, T. Ohgoe, K. Van Houcke, and F. Werner, Phys.Rev. Lett , 130405 (2018). R. Rossi, T. Ohgoe, E. Kozik, N. Prokof’ev, B. Svistunov,K. Van Houcke, and F. Werner, Phys. Rev. Lett ,130406 (2018). S. Kulagin, N. Prokof’ev, O. Starykh, B. Svistunov, andC. Varney, Phys. Rev. B , 024407 (2013). S. Kulagin, N. Prokof’ev, O. Starykh, B. Svistunov, andC. Varney, Phys. Rev. Lett. , 070601 (2013). Y. Huang, K. Chen, Y. Deng, N. Prokof’ev, and B. Svis-tunov, Phys. Rev. Lett. , 177203 (2016). J. Gukelberger, E. Kozik, L. Pollet, N. Prokof’ev,M. Sigrist, B. Svistunov, and M. Troyer, Phys. Rev. Lett. , 195301 (2014). Y. Deng, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov,EPL , 57001 (2015). J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bu-lik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero,T. M. Henderson, C. A. Jim´enez-Hoyos, E. Kozik, X.-W.Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria,H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R.White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Si-mons Collaboration on the Many-Electron Problem), Phys.Rev. X , 041041 (2015). J. Gukelberger, S. Lienert, E. Kozik, L. Pollet, andM. Troyer, Phys. Rev. B , 075157 (2016). I. Tupitsyn and N. Prokof’ev, Phys. Rev. Lett. , 026403(2017). W. Wu, M. Ferrero, A. Georges, and E. Kozik, Phys. Rev.B , 041105 (2017). J. Gukelberger, L. Wang, and L. Pollet, Phys. Rev. B ,205121 (2017). J. Carlstr¨om, Phys. Rev. B , 075119 (2018). S. Iskakov, A. E. Antipov, and E. Gull, Phys. Rev. B ,035102 (2016). J. Gukelberger, E. Kozik, and H. Hafermann, Phys. Rev.B , 035152 (2017). M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez,E. Gull, S. Guo, C. A. Jim´enez-Hoyos, T. N. Lan, J. Li,F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuse-ria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn,S. R. White, D. Zgid, and S. Zhang (Simons Collaborationon the Many-Electron Problem), Phys. Rev. X , 031059(2017). A. S. Mishchenko, N. Nagaosa, and N. Prokof’ev, Phys.Rev. Lett. , 166402 (2014). I. S. Tupitsyn, A. S. Mishchenko, N. Nagaosa, andN. Prokof’ev, Phys. Rev. B , 155145 (2016). I. Tupitsyn and N. Prokof’ev, “Phase diagram topology ofthe haldane-hubbard-coulomb model,” arXiv:1809.01258. R. Rossi, T. Ohgoe, K. Van Houcke, and F. Werner, inpreparation. S. Tan, Ann. Phys. , 2952 (2008). S. Tan, Ann. Phys. , 2971 (2008). E. Braaten, Lecture Notes in Physics , 193 (2012), in . Y. Castin and F. Werner, Lecture Notes in Physics ,127 (2012). F. Werner and Y. Castin, Phys. Rev. A , 013626 (2012). As usual, it is implicit that, in the sum over momentum inEq. (1), formally considering a finite system before even-tually taking the thermodynamic limit, the coordinates of k are integer multiples of 2 π/L where L is the length ofthe cubic box with periodic boundary conditions. Also, forthe lattice model, (cid:15) k = k / k ∈ B , while (cid:15) k is extended outside of B by periodicity. A. L. Fetter and J. D. Walecka,
Quantum Theory of Many-Particle Systems (Dover, Mineola, New York, 2003). L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii,
Sta-tistical Physics part 2 (Butterworth-Heinemann, Oxford,2000). A. Abrikosov, L. Gor’kov, and I. Y. Dzyaloshinskii,
Meth-ods of Quantum Field Theory in Statistical Physics (Dover,New York, 1975). We use the following standard convention for the Fouriertransformation between position and momentum space: f ( r ) = (cid:82) f ( p ) e i k · r d p / (2 π ) . It is implicit that the inte-grals over momenta run over the first Brillouin zone B whenworking with the lattice model, and over the entire space R when working with the zero-range model in continuousspace. We use the following standard notations for the Fouriertransformation between imaginary time and Matsubarafrequencies: for a β -antiperiodic function f ( τ ), f ( ω n ) = (cid:82) β dτ f ( τ ) e iω n τ where ω n = (2 n + 1) π/β are the fermionicMatsubara frequencies; while for a β -periodic function f ( τ ), f (Ω n ) = (cid:82) β dτ f ( τ ) e i Ω n τ where Ω n = 2 nπ/β arethe bosonic Matsubara frequencies. G. C. Strinati, Lecture Notes in Physics , 99 (2012). R. Combescot, X. Leyronas, and M. Y. Kagan, Phys. Rev.A , 023618 (2006). J. Vlietinck, J. Ryckebusch, and K. Van Houcke, Phys.Rev. B , 115133 (2013). The well-known fact that one can always choose N inde-pendent momenta in this way can be proven rigorously byusing a covering tree of the diagram (J. Magnen, privatecommunication ). R. Haussmann, Z. Phys. B , 291 (1993). R. Haussmann, Phys. Rev. B , 12975 (1994). N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. Lett. ,2514 (1998). N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn,JETP , 310 (1998). These numbers were used and partially checked in Ref. 76. N. Prokof’ev and B. Svistunov, Phys. Rev. Lett. , 250201(2007). The relation with the notations of Sec. III A is: ¯ N j = N = (cid:80) ni =1 C i ∈S N and ¯Σ ( h ) j = Z N (cid:80) ni =1 A Σ ( N ) σ ,g ( C i ), with n the number of Monte Carlo steps per iteration. This was demonstrated to us by E. Kozik. R. Combescot, F. Alzetto, and X. Leyronas, Phys. Rev.A , 053640 (2009). P. Pieri, A. Perali, and G. C. Strinati, Nature Physics ,736 (2009). R. Haussmann, M. Punk, and W. Zwerger, Phys. Rev. A , 063612 (2009). H. Hu, X.-J. Liu, and P. D. Drummond, New J. Phys. ,035007 (2011). E. Braaten and L. Platter, Phys. Rev. Lett. , 205301(2008). Starting from Eq. (3.36) of Ref. 77, our Eqs. (66,68) canbe rederived [Y. Nishida, private communication]. T. Enss, R. Haussmann, and W. Zwerger, Ann. Phys. ,770 (2011). For the construction of the skeleton formalism in Sec. II,this is only a formal issue that does not cause any problems:Equation (18) can be directly derived from Eq. (8) withoutusing Eq. (12), and only the inverse of Γ (0) appears in thefinal expressions Eqs. (6,17,18,19). E. Kozik, M. Ferrero, and A. Georges, Phys. Rev. Lett. , 156402 (2015). J. Vlietinck, J. Ryckebusch, and K. Van Houcke, Phys.Rev. B , 085119 (2014). P. Kroiss and L. Pollet, Phys. Rev. B , 104510 (2014). R. Rossi, Phys. Rev. Lett. , 045701 (2017). F. Simkovic and E. Kozik, arXiv:1712.10001. A. Moutenet, W. Wu, and M. Ferrero, Phys. Rev. B ,085117 (2018). R. Rossi, N. Prokof’ev, B. Svistunov, K. Van Houcke, and F. Werner, EPL , 10004 (2017). K. Chen and K. Haule, “Feynmann’s solution ofthe quintessential problem in solid state physics,”arXiv:1809.04651. P. V. Buividovich, Nucl. Phys. B , 688 (2011). P. V. Buividovich and A. Davody, Phys. Rev. D , 114512 (2017). T. Pfeffer and L. Pollet, New J. Phys. , 043005 (2017). T. Pfeffer and L. Pollet, Phys. Rev. B , 195104 (2018). R. Haussmann, W. Rantner, S. Cerrito, and W. Zwerger,Phys. Rev. A , 023610 (2007). F. B. Kugler, Phys. Rev. E , 023303 (2018). Y. Nishida, Phys. Rev. A85