Bold diagrammatic Monte Carlo: A generic technique for polaron (and many-body?) problems

0.33 MB·English

by  Nikolay Prokof'ev

Bold diagrammatic Monte Carlo:
A generic technique for polaron (and many-body?) problems
N.V. Prokof’ev1,2,3 and B.V. Svistunov1,3
1Department of Physics, University of Massachusetts, Amherst, MA 01003, USA
2Theoretische Physik, ETH Zu¨rich, CH-8093 Zu¨rich, Switzerland
3Russian Research Center “Kurchatov Institute”, 123182 Moscow, Russia
We develop a Monte Carlo scheme for sampling series of Feynman diagrams for the proper self-
energy which are self-consistently expressed in terms of renormalized particle propagators. This
approach is used to solve the problem of a single spin-down fermion resonantly interacting with
8 the Fermi gas of spin-up particles. Though the original series based on bare propagators are sign-
0 alternating and divergent one can still determine the answer behind them by using two strategies
0 (separatelyortogether): (i)usingproperseriesre-summationtechniques,and(ii)introducingrenor-
2 malized propagators whicharedefinedintermsofthesimulated properself-energy, i.e. makingthe
entire scheme self-consistent. Our solution is important for understanding the phase diagram and
n
properties of theBCS-BEC crossover in the strongly imbalanced regime. On thetechnical side, we
a
developagenericsign-problemtolerantmethodforexactnumericalsolutionofpolaron-typemodels,
J
and,possibly, of theinteracting many-body Hamiltonians.
7

PACSnumbers: 05.30.Fk,05.10.Ln,02.70.Ss
]
l
e
-
r I. INTRODUCTION coupling between the particle and bosonic environment,
t but it can not be generalized to fermionic environment,
s
t. Modern science has radically changed our view of the sign-alternating coupling (as in the t-J model [11]), nor
a physical vacuum. Instead of na¨ıve “empty space” we is it suitable for continuous-space models.
m
have to deal with a complex groundstate of an interact- In this article (which follows a short communication
- ingsystem,and,strictlyspeaking,thereisnofundamen- [12]), we develop a Monte Carlo technique which sim-
d
n tal difference between the outer space, helium, or any ulates series of Feynman diagrams for the proper self-
o other condensed matter system. With this point of view energy. ThediagrammaticMonteCarlo(Diag-MC)tech-
c comesunderstandingthatthenotionofa“bare”particle nique was used previously to study electron-phonon po-
[ is somewhat abstract since its coupling to vacuum fluc- larons [13, 14]. The essence of Diag-MC is in interpret-

1 tuations, orenvironment,maystrongly(sometimes radi- ing the sum of all Feynman diagrams as an ensemble
v cally)changeparticlepropertiesatenergiesaddressedby averaging procedure over the corresponding configura-
1 the experimental probes. The polaron problem [1] is by tion space. It was consideredessentialfor the method to
1 now canonical across all of physics with the same ques- workthatthe seriesofdiagramsforthe Green’s function
9 tions about particle energy, effective mass, residue, etc., be convergent and sign-positive. Though the configura-
0
being askedfor differenttypes of particles,environments tion space of diagrams for polarons in Fermi systems is
.
1 and coupling between them [2]. In a broader context, more complex, similar methods of generating and sam-
0 particles are not necessarily external objects unrelated pling the correspondingconfigurationspace canbe used.
8
to a given vacuum, but quasiparticle excitations of the ThecrucialdifferenceisthatintheFermisystemwehave
0
very same ground state. Thus, the solution of the po- to deal with the sign-alternating and divergent (at least
:
v laron problem paves the way to the effective low energy for strong coupling) series. Under these conditions a di-
Xi theoryofagivensystem,and,tolargeextent,determines rectsummationofallrelevantFeynmandiagramsforthe
basic properties of all condensed matter systems at low Green’s functions is not possible,and one has to develop
r
a temperature. additional tools for (i) reducing the number of diagrams
Atthemoment,thereisnogenericanalyticornumeric bycalculatingself-energiesratherthanGreen’sfunctions,
technique to study quasiparticle properties for arbitrary (ii) employing the “bold-line” trick in the form of the
strongly interacting system. Analytic solutions are typi- Dyson equation which allows self-consistent summation
cally(withfewexceptions,see,e.g.,[3])basedonpertur- of infinite geometric series and further reduces the num-
bative corrections to certain limiting cases [1, 2, 4, 5, 6] berofself-energydiagrams,and,ifnecessary,(iii)extrap-
or variational treatment [7]. Several numeric schemes olatingDiag-MCresultsto theinfinite diagramorderfor
were suggested in the past, but all of them have limita- a divergentseries. At the moment we do notsee any ob-
tions either in the system size, system dimension, inter- viouslimitationsofthe newmethodsince evendivergent
action or environment type. Exact diagonalization and signalternatingseriescanbedealtwithtoobtainreliable
variationalmethodsinthelow-energysubspace[8, 9]are results. We believe that our findings are important in a
mostly restricted to one-dimensional models with short- muchbroadercontextsincetheDiag-MCapproachtothe
range interactions. The continuous time path integral many-body problem has essentially the same structure.
formulation [10] works for the lattice model with linear Asapracticalapplicationofthe methodweconsidera
2
particle coupled to the ideal Fermi sea via a short-range Monte Carlo simulations [19]. On the experimental side,
potential with divergent s-wave scattering length. This MIT experiments [26] are in good agreement with the
problemwasrecognizedasthekeyoneforunderstanding predictionsmadeinRef.[19],whileRiceexperiments[27]
thephasediagramofthepopulationimbalancedresonant are not. The origin of discrepancy between the two ex-
Fermigas [15, 16]. Inparticular,to constructthe energy periments is not understood. Our results for polaron
functional describing dilute solutions of minority (spin energies are in excellent agreement with Ref. [19].
down)fermions resonantlycoupled to the majority(spin It is worth noting that the model (1) (in general, the
up) fermions one has to know precisely the quasiparticle particle mass is different from that of the Fermi gas) is
parameters of spin-down fermions since they determine also known as the Anderson orthogonality problem with
the coefficientsin the energyexpansionin the spin-down recoil [28, 29]. It can be also considered as a specific
concentration x↓: the linear term is controlled by the example of a particle coupled to the Ohmic dissipative
polaron energy, and the x5/3 term is determined by the bath (see monograph [30] for numerous other examples
↓
and connections to realistic systems).
polaron mass [16].
The general Hamiltonian we deal with in this article The paper is organized as follows. In Sec. II we dis-
can be written as cuss the configuration space of Feynman diagrams for
self-energyinmomentum–imaginary-timerepresentation
H = H ∆ /2m + drV(r R)n(r), (1) (both in the particle and molecule channels), and ex-
F R
− − plain how polaron parameters can be obtained in this
Z
representation. In Sec. III we describe a Monte Carlo
where H is the Hamiltonian of the ideal spin-up Fermi
F
algorithm for generating and sampling the correspond-
gas with density n and Fermi momentum k , R is the
F
ing diagrammatic space. A small technical Section IV
particle coordinate, and V(r R) is the interaction po-
− deals with numerically evaluating the effective T-matrix
tential of finite range r between the particle and the
0
bybolddiagrammaticMonteCarlo. Wepresentanddis-
spin-up Fermi gas. In what follows we refer to (1) as the
cussresultsinSec.V. Inparticular,weshowthatonecan
Fermi-polaron problem. The specifics of the BCS-BEC
use re-summation techniques for divergent series of dia-
crossover physics in the strongly imbalanced regime is
grams based on bare propagators to determine the final
two-fold: one is that the particle and the Fermi gas have
answer. In Sec. VI we further advance the algorithm by
the same bare mass m (in what follows we will use units
employingbold-line approachinwhichthe entirescheme
such that m = 1/2 and k = 1), and the other is that
F
is self-consistently based on renormalized (“bold-line”)
one has to take explicitly the so-called zero-range reso-
propagators. We present our conclusions and perspec-
nant limit when r 0, but the s-scattering length a
0
→ tives for future work in Sec. VII.
remains finite, i.e. k a remains fixed for k r 0. In
F F 0
→
this limit, the nature of the interactionpotential is irrel-
evant, and the same results will be obtained e.g. for the
neutron matter and Cesium atoms. We note, however, II. CONFIGURATION SPACE OF
that the method we develop for the numeric solution of SELF-ENERGY DIAGRAMS
the resonant Fermi-polaron problem is absolutely gen-
eralandcanbe usedforanarbitrarymodeldescribedby As mentioned above, when coupling between spin-
Eq. (1). down and spin-up fermions is strong enough they from
It turns out that the structure of the phase diagram a composite boson, or molecule state. In what follows,
is very sensitive to polaronparameters. If the state with we will be using the term “polaron” in a narrow sense,
a dilute gas of spin-down fermions is stable at all val- i.e. only for the unbound fermionic spin-down excita-
ues of kFa then the solution of the single-particle prob- tion. For the composite boson we will be using the
lem would define the phase diagram in the vicinity of term “molecule”. Since our goal is to calculate parti-
themulticriticalpointdiscussedrecentlybySachdevand cle properties for arbitrary coupling strength we have to
Yang[17],wherefourdifferentphasesmeet. Atthispoint considerone-andtwo-particlechannelsonequalfooting.
the spin-down fermion forms a bound state with a spin- In the rest of this section we render standard diagram-
up fermion thus becoming a spin-zero composite boson matic rules for irreducible self-energy diagrams in both
(“molecule”); i.e. quasi-particles radically change their channels,withanemphasisonspecificsofworkinginthe
statistics. Themulticriticalpoint,however,maybether- imaginary-time representation.
modynamicallyunstableifthe effective scatteringlength
between the composite bosons and spin-up electrons is
large enough, and the analysis of Refs. [18, 19] based on
A. Polaron channel
thefixed-nodeMonteCarlosimulationsfindsevidencein
favor of this scenario. Phase separation was also found
incalculationsbasedonmean-fieldandnarrow-resonance We start from the definition of the single-particle
approaches(bothatfiniteandzerotemperature)seee.g. Green’s function (see e.g. [31])
Refs. [20, 21, 22, 23, 24, 25], though with strong quan-
titative deviations from results based on the fixed-node G(τ,r)= T ψ(τ,r)ψ¯(0) , (2)
τ
−h i
3
and its frequency-momentum representation
-G = -G(0) + -G(0) -G
G(ξ,p)= ei(ξτ−p·r)G(τ,r)drdτ . (3) -Σ
Z
FIG. 1: Dyson equation for the single-particle Green’s func-
Here ψ(τ,r) is the fermion annihilation operator at the
tion.
space-time point (τ,r). For the ideal spin-up Fermi gas
at T =0 we have
Equation (11) allows one to solve for E provided
1
G (ξ,p)= . (4) Σ(τ,p,µ) is known. All we have to do is to calculate
↑ iξ p2/2m+ǫ
F the the zero-frequency value of Σ for µ=E
−
The vacuum Green’s function for the spin-down po- ∞
E =p2/2m+ Σ(τ,p,µ)e(E−µ)τdτ . (13)
laron is
Z0
G(0)(τ,p,µ) = θ(τ)e−(p2/2m−µ)τ , (5) After E is found, the residue is obtained from Eq. (12)
↓ − using
where θ is the step function, and µ is a free parameter. ∞
From Dyson’s equation for the polaron, see Fig. 1, one A(p,E)= τΣ(τ,p,µ)e(E−µ)τ dτ . (14)
−
finds Z0
Note alsothatthe dependence onµ dropsoutfromboth
1
G (ξ,p,µ) = , (6) (13) and (14).
↓ iξ p2/2m+µ Σ(ξ,p,µ) A comment is in order here. Strictly speaking, po-
− −
laron and molecule poles exist only for p = 0 because
where self-energyΣ is givenby the sumofallirreducible
the fermionic bath they couple to is gapless. However,
diagrams (i.e. diagrams which can not be made discon-
the spectrum E(p) is well-defined in the limit p 0, as
nected but cutting through the G(0) line) taken with the decay width vanishes faster than [E(p) E(→0)]. To
↓
the negative sign. Taking into account that in the τ- have stable quasi-particles, one can use a tr−ick of intro-
representationbothG andΣdependonµonlythrough ducing a gap ∆ in the environment spectrum, e.g. by
↓
exponential factors exp(µτ), one obtains G↓(ξ,p,µ) redefining the dispersion relation for spin-up fermions:
≡
G↓(ξ−iµ,p) and Σ(ξ,p,µ)≡Σ(ξ−iµ,p). εk → max(εk,∆). In the p → 0 limit, the system-
If polaron is a well-defined quasi-particle, then its en- atic error vanishes faster than [E(p) E(0)], provided
ergy E(p) and residue Z(p) can be extracted from the ∆ [E(p) E(0)]. It should be also−possible to work
∼ −
asymptotic decay with∆’sessentiallylargerthan[E(p) E(0)]andextrap-
−
olate to ∆ 0. In particular, such an extrapolation is
G↓(τ,p,µ) Ze−(E−µ)τ , τ . (7) possible(an→disimplicitlyimplied)attheanalyticallevel
→− →∞
in the relation for the effective mass, which we consider
This asymptotic behavior immediately implies that the
below.
function G (ξ iµ,p) has a pole singularity
↓ − Onewaytodeterminetheeffectivemassistocalculate
the quasi-particle energy as a function of momentum for
Z(p)
G (ξ iµ,p)= +regular part. (8) a set of small but finite values of p and fit it with the
↓ − iξ+µ E(p)
parabola. It is, however, possible to skip this procedure
−
and to write a direct estimator for the effective mass in
Now setting µ = E(p) in (8) and comparing the result
terms of the calculated self-energy. Acting with the op-
to (6), we conclude that
erator 2 onboth sides of Eq.(11) and taking the limit
∇P
iξ/Z =iξ p2/2m+E Σ(0,p,E)+iξA(p,E), (9) p→0 we get
− −
1+A 1
0
= +B , (15)
where (we take into account that ∂Σ/∂ξ =i∂Σ/∂µ) m m 0
∗
A(p,E)= ∂Σ(0,p,µ) . (10) B = 1 ∞dτe(E0−µ)τ [ 2 Σ(τ,p,µ)] , (16)
− ∂µ (cid:12)(cid:12)µ=E 0 3Z0 ∇P p=0
This yields two important relations (s(cid:12)(cid:12)ee also [31]): where A0 A(p=0) and E0 E(p=0).(cid:12)(cid:12)
≡ ≡
E =p2/2m+Σ(0,p,E), (11)
B. Molecule channel
and
In this case we start with the two-particle propagator
1
Z = . (12)
1+A(p,E) K(τ,p)=−hTτΦp(τ)Φ†p(0)i, (17)
4
dq/(2π)3
-K Π¯(η,p)= . (24)
-Q q2/2m+(p q)2/2m η
= + Zq≤kF − −
Forfinitedensityofspin-upfermionsconvertingEq.(23)
~ to the imaginary time domain has to be done numeri-
-K
-Q -Γ cally. One possibility is to use the inverseLaplace trans-
= +
form. We employedthe bold diagrammatic Monte Carlo
technology [32] to achieve this goal and further details
are provided in Sec. IV. The two-dimensional function
FIG.2: DefiningfunctionsQandK˜ inthetwo-particlechan- Γ(τ,p) is tabulated prior to the polaron simulation.
nel. In Fig.2 we define function Q thatcan be viewedas a
renormalized pair propagator related to Γ by the Dyson
equation. Intheupperpanel,weshowthediagrammatic
where
structure for K, which includes dotted lines represent-
dq
Φp = ϕqψ↑(p q)ψ↓(q), (18) ingexternalfunctionsϕq,greysquaresrepresentingsums
2π − of all Γ-irreducible diagrams, and the renormalized pair
Z
propagatorQ. By Γ-irreduciblediagramswe understand
andϕ isthepairwavefunction(inmomentumrepresen-
q
diagramswhichcannotbemadedisconnectedbycutting
tation) that localizes the relative distance between two
them through a single Γ-line. All Γ-reducible diagrams
particles. If there is a bound state (molecule), then
are absorbed in the Q function which is shown in the
K(τ,p,µ) Z e−(Emol−µ)τ , τ , (19) lower panel in Fig. 2. The grey circle has nearly the
mol
→− →∞
samestructureasthegreysquare(the zerothorderterm
and the pair propagator in the frequency representation is present in the crossed square, but not in the crossed
has a pole: circle): since Γ is defined as the sum of ladder diagrams,
all terms which include ladder-type structures based on
Z (p)
K(ξ iµ,p)= mol +regular part. (20) freeone-particlepropagatorshavetobeexcludedfromQ
− iξ+µ E (p)
− mol and K. With the replacements G Q, G(0) Γ, and
↓ → ↓ →
Now we introduce yet another function, that features Σ K˜ we find an exact analogy between the one- and
thesamemolecularpole,buthasasimplerdiagrammatic →
two-particle propagators.
structure. The specifics of the resonant zero-range limit
The analogycan be carriedout further by noting that
is that the sum of all ladder diagrams for the interac-
thestructureofdiagramsinFig.2impliesthatQhasthe
tion potential V(r) has to be considered as a separate
same pole as K, while the rest of the functions simply
diagrammatic element. We denote this sum by Γ(τ,p)
change the value of the quasi-particle residue. Thus, if
and consider it as a “bare pair propagator”. Of course,
molecule is a well-defined excitation we expect that
the same approach can be taken in a general case to re-
placethebareinteractionpotentialwiththescatteringT-
matrix,but inthe zero-rangelimit we really do not have Z˜ (p)
Q(ξ iµ,p)= mol +regular part. (25)
any other alternative for the ultra-violet-divergence-free − iξ+µ E (p)
mol
formulation. Thesumofladderdiagramstakestheultra- −
violet physics into account exactly and allows to express
This explains why introducing the function Q is conve-
Γ(τ,p) in terms of the s-scattering length a. The ladder
nient: nowEqs.(12)-(16)areimmediatelygeneralizedto
structure of diagrams absorbed in Γ(τ,p) also explains
themoleculecase(uptoreplacementsmentionedabove).
why we treat it as a “pair propagator” (we will depict
it with a double-line, see Fig 2). The exact expression
for Γ is readily obtained from the geometrical series (in
frequency domain)
III. WORM ALGORITHM FOR FEYNMAN
Γ = U + ( U)2Π + ... , (21) DIAGRAMS
− − −
where the polarization operator is defined by
InthisSectionwedescribehowtheconfigurationspace
dq/(2π)3 of Feynman diagrams for Σ and K˜ is parameterizedand
Π(η,p)= , (22)
q2/2m+(p q)2/2m η updatedusingDiag-MCrules. Ouralgorithmisdesigned
Zq>kF − − to treatpolaronand molecule channels on equalfooting.
and η =ω+εF +µ. Using suitable ultra-violet regular- We achieve this goal by introducing auxiliary diagrams
ization,onecancastthesameexpressionintheuniversal whichcontaintwo“loose”endscalled“worms”—thiswas
form which depends only on the s-scattering length: proven to be an efficient strategy for reducing the MC
m m autocorrelation time when simulations are performed in
Γ−1(η,p)= p2 4mη Π¯(η,p), (23) the configuration space with complex topology [33, 34].
4πa − 8π − −
p
5
τ τ τ τ τ τ τ τ τ τ τ τ
a b c d e f c d e f a b
FIG. 3: The backboneof thecyclic diagram. FIG.5: Backward connection. Thepairof Γ-endsbeingcon-
nected is precisely the same as in Fig. 4, but the direction is
opposite, and thearc represents −G↑(τ =−τf −τa−τb).
τ τ τ τ τ τ
a b c d e f I M
FIG. 4: Forward connection. The arc represents −G↑(τ =
τc+τd+τe).
FIG. 6: A diagram with two worms, I and M.
A. Cyclic diagrams
manipulations with and . As it will become clear
I M
soon, of special importance for normalization of MC re-
We start with the introduction of cyclic diagrams.
sultsisthefirst-orderdiagramwiththeworm,seeFig.7.
Though we work in the imaginary time representa-
Itsweightconsistsofjusttwofactors,G(0)(τ )andΓ(τ ).
tion at T = 0 when τ [0, ), it is still conve- ↓ a b
∈ ∞
nient not to specify the time origin and to consider di-
agrams on the imaginary time circle. The backbone
of each cyclic diagram is a pre-diagram illustrated in C. Parametrization of diagrams
Fig. 3. It consists of a cyclic chain of the struc-
ture G(0)(τ )Γ(τ )G(0)(τ )Γ(τ )G(0)(τ )Γ(τ ) ... (all Apartfrom the diagramorderand its topology,we se-
↓ a b ↓ c d ↓ e f
thetimesarepositive). Wedonotexplicitlyshow“direc- lecttime intervals ofΓ’s andG(0)’s, andmomenta ofthe
↓
tions”ofthe propagators,sincethese areunambiguously spin-up propagators as independent variables. The mo-
fixed by the global direction of all the backbone lines, menta of Γ’s and G(0)’s are then unambiguously defined
↓
which we select—without loss of generality—to be from by the momentum conservation law, while the time in-
right to left. With this convention, the left (free) spin- terval of a spin-up propagator is obtained by summing
up end of any Γ-line is outgoing, while the right end is up the time intervals of Γ’s and G(0)’s covered by this
↓
incoming. Aphysicaldiagramisobtainedbypairwisere- propagator. Technically,wefinditconvenienttoworkex-
placingfreespin-upendswithpropagatorsG↑. Thereare plicitlyintheparticle-holerepresentationforthespin-up
twowaystoconnectincomingandoutgoinglines: (i)for- propagatorswhenbackwardspin-uppropagatorisunder-
ward, i.e., in the direction of the backbone propagators, stoodasaforwardholepropagatorwiththeoppositemo-
and(ii) backward(opposite to forward). Forward(back- mentum. This is achievedby introducing a non-negative
ward)connectionsresultinpropagatorsG↑ withpositive function,
(negative)times,seeFigs.4and5forillustrations. They
represent particle (hole) excitations in the fermionic en- G (τ,p), p p ,
vironment. It is important to emphasize that in cyclic G˜(τ,p)= G− (↑ τ, p), p≥<Fp , (26)
↑ F
diagrams the only time-variables are the positive time- (cid:26) − −
lengthsofG(0)’sandΓ’s. Thereisnoabsolutetime,and, whichisassignedtoallspin-uplines(theglobalfermionic
↓
correspondingly, all moments in time are equivalent. sign of the diagram is defined separately, by standard
diagrammatic rules). All momenta assigned to the spin-
uplinesareunderstoodasmomentaofthecorresponding
B. Worms G˜-propagators;i.e. they are either momenta of particles
(for forward propagators they are non-zero only for p
≥
To have a MC scheme which simulates diagrams in pF),ormomentaofholes(forbackwardpropagatorsthey
one- and two-particle channels on equal footing we ex- are non-zero only for p<pF). An explicit formula for G˜
tendthespaceofphysicaldiagramsbyallowingdiagrams (subject to the above conditions) is
with two special end-points, or “worms”. They will be
denoted and andstandfor the unconnectedincom- G˜(τ,p)=θ(τ)e−|p2/2m−εF|τ . (27)
I M
ing(outgoing)spin-upends,seeFig.6foranillustration.
Correspondingly,the entire updating scheme is based on
6
as a parameter. The acceptance ratio for this update is
M I Γ(τ ,p)G(0)(τ ,p)
τ τ P =NC 1 ↓ 2 , (28)
ins N+1
b a W (τ )W (τ )
Γ 1 ↓ 2
where C is an artificial weighing factor ascribed to all
N
worm diagrams of order N (it can be used for optimiza-
FIG.7: “Normalization” diagram. Itisthesimplest diagram
with theworm; itsweight isaproduct of G(↓0)(τa)and Γ(τb). tion purposes and depend on p too). A natural choice
for W-functions is to make them proportional to BBL,
i.e.
To simplify the description of updates below we will
genericallyrefer to G(0)- andΓ-propagatorsas backbone W (τ)= Γ(τ,p) , W (τ)= G(↓0)(τ,p) .
lines (BBLs) and den↓ote them as D. The diagram order Γ Γ(τ′,p)dτ′ ↓ G(↓0)(τ′,p)dτ′
N isgivenbythetotalnumberofspin-downpropagators. R R (29)
Wealsousespecialcounterstocharacterizethe topology Then, to have an acceptance ratio of order unity and
of the diagram. For each BBL we define a cover num- independent of p we choose
ber, N , equal to the total number of G˜-lines covering a
c
1
gcoivveenreBd.BLF.inAalblya,ckpbhoynsiecallindeiawgirtahmNs—c t=he0oinsecsalwleidthuount- CN = NΛ, Λ= Γ(τ1,p)dτ1 G(↓0)(τ2,p)dτ2. (30)
Z Z
worms—are divided into polaron (0) and molecule (1)
In the rest of the paper, we will refer to this choice of
sectors; the diagram sector is defined by the difference
W andC as “optimized”,though we do notmean that
between the number of particle and hole spin-up propa- N
it is the best one possible for the entire scheme. For the
gators covering any of the G(0)-lines (the same result is
↓ “optimized” choice
obtained by analyzing propagatorscovering Γ-lines after
adding unity for the spin-upparticle participatingin the
P =N/(N +1). (31)
ins
ladder diagrams).
In the molecule sector, we essentially repeat all steps,
up to minor modifications. Now the propagator being
D. Updates
selectedisΓ (once again,the update isrejectedif the se-
lected propagatoris covered.) Then, a pair of new prop-
The cyclic structure of diagrams in combination with agators, G(0)(τ ,p) and Γ(τ ,p), is inserted in front of
the possibility of considering non-physical diagrams al- ↓ 1 2
the selected uncovered propagator. The new propagator
lowsonetoconstructaverysimpleergodicsetofupdates. Γ(τ ,p)inheritstheoutgoingspin-uplinepreviouslycon-
2
Theminimalsetconsistsoftwocomplementarypairs,In-
nected to the selected Γ; the latter gets instead the -
sert/DeleteandOpen/Close,andoneself-complementary M
end of the worm while the -end is attached to the new
update Reconnect. The description that follows aims at I
Γ. The acceptance ratio is identical to (28) [for the opti-
the most transparent and straightforward realization of
mizedchoiceitisN/(N+1)]. The polaronandmolecule
updates. Standard performance enhancement tricks and
sectors are mutually exclusive due to particle conserva-
optimization protocols are not mentioned. In particular,
tion,thusonlyonetypeoftheInsertupdateisapplicable
webaseourconsiderationsontheupdatingschemewhich
to a given diagram.
may propose a change leading to a forbidden configura- Delete: Thisupdateconvertswormdiagramstophys-
tion. Such proposalsare rejectedas if they resultin zero
ical ones while reducing the diagram order. It applies
acceptance ratio.
only to diagrams of order N > 1 with worms being sep-
Insert: This update applies only to physical—no
arated by one uncovered BBL. Otherwise, the update is
worms—diagrams and is rejected otherwise. First, con-
rejected. If worms are separated by an uncovered BBL,
siderthepolaronsector. SelectatrandomoneoftheG(0) thentheleftandrightneighborBBLsarealsouncovered,
↓
propagators;if it is covered,the update is automatically and an update opposite to Insert is possible. In Delete
rejected. If the selected propagator is uncovered, insert we simply remove two consecutive BBL and worms from
a pairofnew propagators,Γ(τ ,p)andG(0)(τ ,p), right the diagram. Its acceptance probability is the inverse of
1 ↓ 2
after the selected one. The new Γ(τ ,p)-propagator is Eq. (28),
1
supposed to contain and at its ends. Worms radi-
I M 1 W (τ )W (τ )
cally simplify this diagram-orderincreasingupdate since Γ 1 ↓ 2
P = . (32)
due to conservation laws the momenta of new BBLs are del (N −1)CN Γ(τ1,p)G(↓0)(τ2,p)
equal to the global momentum of the diagram p. The
times τ1 and τ2 are drawn from normalized probabil- With Eqs. (29), (30), we have
ity distributions W (τ ) and W (τ ) (arbitrary at this
Γ 1 ↓ 2
point). Note that W (τ ) and W (τ ) can depend on p P =N/(N 1). (optimized) (33)
Γ 1 ↓ 2 del
−
7
Close: The update applies only to diagrams with the the selected spin-up propagator. Primes indicate new
worms. The proposal is to connect and with a line values of the BBL momenta:
G˜(τ,q)andeliminatewormsfromthIediagMram. Themo-
mentum variable q is drawn from the probability distri- p′ν =pν +q. (40)
butionW (q),whileτ isthetimeintervalbetween and
↑
I In the optimized version, we have
to be covered by the new propagator. There are two
M
ways of connecting and , forwards and backwards.
The ambiguity is auItomaticMally resolved by the absolute P = (2π)3 Dν(τν,p′ν) . (41)
op 2ΛΩ(τ) D (τ ,p )
value of the momentum variable q: if q ≥ pF (q < pF), Yν ν ν ν
the propagator is supposed to go forwards (backwards).
Reconnect: The purpose of this update is to change
Inpractice,wefirstselectthedirection(withequalprob-
thetopologyofdiagramswiththeworm. Theproposalis
abilities), and then generate the momentum variable q
toselectatrandomoneoftheG˜-propagatorsandswapits
accordingly: either q p or q <p .
F F
≥ outgoing end with ; the momentum of the propagator
The acceptance ratio for this update is M
remains the same,only its time variable changesfrom τ
0
P = 2 G˜(τ,q) Dν(τν,p′ν) , (34) to τ0′. The acceptance ratio is given by:
cl (2π)3NCN W↑(q) Yν Dν(τν,pν) P = G˜(τ0′,q) Dν(τν,p′ν) . (42)
where the subscript ν runs through all BBLs to be cov- rec G˜(τ0,q) ν Dν(τν,pν)
Y
eredbythenewpropagator(clearly,τ = τ ). Primes
ν ν The subscript ν runs through all BBLs that will change
indicate new values of the corresponding momenta:
P their momenta (and cover numbers Nc’s) as a result of
p′ =p q. (35) the update. Topologically, there are two different situa-
ν ν − tions (the two a complementary to each other in terms
Asusual,the distributionfunctionW (q)candependon of the update): (i) is covered by the propagator in
↑
M
τ andthedirectionofthepropagator. Thenaturalchoice question, and (ii) is not covered by the propagator.
M
for this function would be Correspondingly,thenewvaluesofthediagramvariables
are calculated as
G˜(τ,q)
W (q)= , (36)
↑ Ω(τ) (τ′,p′)= (τ0−τ, pν +q) (i), (43)
0 ν (τ +τ, p q) (ii),
0 ν
(cid:26) −
Ω(τ)= q≥pF GG˜˜↑((ττ,,qq))ddqq ((bfoarcwkawradr)d,), (37) witThhτe =abovνeτsνe.t of updates is ergodic. It can be sup-
( Rq<pF ↑ plementedPby additional updates that may improve the
R algorithm performance by more efficient sampling of the
leading to the optimized acceptance ratio
diagram variables and topologies. Over-complete sets of
2ΛΩ(τ) D (τ ,p′) updates are also useful for meaningful tests of the de-
P = ν ν ν . (38)
cl (2π)3 D (τ ,p ) tailedbalance. The possibilities are endless,andhere we
ν ν ν
Yν simply mention two updates we have been using.
Time shift: Weproposenewtime variables,τ τ′,
In this Section we deal with diagrams based on bare ν → ν
for all uncovered BBLs (labeled here with the subscript
propagators. To avoid double-counting, we have to ex-
ν). The acceptance probability is
clude all cases which can be reduced to ladders already
summed in Γ’s. When and are on the nearest-
neighbor BBL the propoMsal to coInnect them with the P = Wν(τν,p)Dν(τν′,p) . (44)
sh W (τ′,p)D (τ ,p)
spin-up particle propagator is rejected. The last rule to ν ν ν ν ν
Y
be monitored it to restrict all physical diagrams to be
All uncovered propagators have the same momentum p.
either in the polaronor molecule sectors,i.e. sectors dif-
Withthe optimizedchoiceforW (τ ,p) D (τ ,p)the
ferent from 0 and 1 are not allowed. ν ν ν ν
∝
acceptance ratio becomes unity.
Open: The update applies only to physical diagrams
Redirect: Here we propose to select at random one
and proposes to create a worm by selecting at random
of the G˜-propagatorsand change its direction to the op-
and removing one of the spin-up propagators. The ac-
posite. Simultaneously, we change the momentum of the
ceptance ratio is given by the inverse of (34), (38)
selected propagator (resulting in new momenta for all
(2π)3NC W (q) D (τ ,p′) BBLitcoversorwillcoverasaresultoftheupdate). Let
Pop = 2 N G˜↑(↑τ,q) Yν Dνν(τνν,pνν) , (39) Gt˜h(eτ′s,eqle′c)twedithprτop=agatνoτrνb,eanG˜d(ττ,′q=), anλdτλth,ewnheewreoνnreubnes
where the subscript ν runs through all BBLs covered by throughallBBLscoveredbythepropagatorG˜ (τ,q)and
P P ↑
the propagator, τ = τ , and q is the momentum of λrunsthroughallBBLstobecoveredbythepropagator
ν ν
P
8
τ τ τ τ τ τ τ τ τ τ
a b c d e f b c d e
FIG. 8: A GΣ-diagram contributing to the polaron self- FIG. 9: A ΓK˜-diagram contributing to the molecule self-
energy. It factorizes into a product of G(↓0)(τa) and Σ(τ = energy. It factorizes into a product of Γ(τb) and K˜(τ =
τb+τc+τd+τe+τf). τc +τd + τe). The vertical dashed lines cut (for the sake
of better visual perception) thesame Γ(τ )-line.
b
G˜ (τ′,q′). We assume that q′ is drawn from the distri-
↑
either filtered out at the time when the contribution to
butionW introducedabove. Inthiscase,theacceptance
↑
the self-energy histogram is made, or, are not produced
ratio is given by
inthesimulationatall.] Theyfactorizeintoaproductof
W (τ,q)G˜(τ′,q′) G(0) andsomediagramcontributingtotheself-energyΣ.
P = ↑ ↓
rdr W (τ′,q′)G˜(τ,q) × The utility of cyclic representationis that the uncovered
↑
propagator can be anywhere on the backbone. In view
Dν(τν,p′ν) Dλ(τλ,p′λ) , (45) of factorization, it is easy to write the MC estimator for
D (τ ,p ) D (τ ,p ) the integral (to simplify notations we omit irrelevant to
" ν ν ν #" λ λ λ #
ν λ
Y Y the discussion momentum p)
p′ =p +q, p′ =p q′ . (46) ∞
ν ν λ λ− I = f(τ)Σ(τ)dτ , (48)
Z0
For the optimized choice of W , see Eq. (36),
↑
where f(τ) is some function [see, e.g., Eqs. (13), (14)].
W (τ,q)G˜(τ′,q′) Ω(τ′) Indeed, consider the estimator
↑
. (47)
W (τ′,q′)G˜(τ,q) → Ω(τ)
↑ = δ f(τ). (49)
I GΣ
E
Diagram sign: The sign of a diagram with worms
which counts all instances of GΣ-diagrams with an ad-
is somewhat arbitrary since it is not physical. The am-
ditional weight f(τ). Here δ is unity for each GΣ-
biguity is resolved by assuming that is always con- GΣ
M diagram and zero otherwise, and τ is the total duration
nected to in the backward direction by an auxiliary
I in time of the Σ-part of the GΣ-diagram. The Monte
unity propagator. Then, to comply with the diagram-
Carlo average of this estimator is
matic rules, we change the configuration sign each time
anyofthefollowingupdatesareaccepted: (i)Reconnect, ∞ ∞
(ii) Open/Close updates dealing with spin-up propaga- hEIiMC ∝ G(↓0)(τ′)dτ′ f(τ)Σ(τ)dτ . (50)
tors in the forward direction, (iii) Insert/Delete in the Z0 Z0
molecule sector, and (iv) Redirect. For Open/Close up- Similarly,withinthesamescheme,wecollectstatisticsof
dates dealing with spin-up propagators in the backward all“normalization”diagrams,seeFig.7,usinganestima-
direction the sign remains the same, because here the tor projecting to the first-order diagram with the worm,
sign coming from changing the number of closedspin-up δ . Then,
norm
loopsiscompensatedbythe signinEq.(26); the sameis
∞ ∞
also true for Insert/Delete updates in the polaron sector δ C G(0)(τ′)dτ′ Γ(τ)dτ . (51)
(due to our choice of the auxiliary propagatorsign). For h normiMC ∝ 1 ↓
Z0 Z0
precisely the same sign compensation reason, the Redi-
rect update does change the sign despite the fact that it The proportionalitycoefficientcancelsinthe ratioofthe
preserves the number of loops. two averages,leading to
∞
I MC
I = C hE i Γ(τ)dτ . (52)
1
E. Estimators hδnormiMC Z0
In particular, for the optimal choice of C , we have
OnlyphysicaldiagramswithoneuncoveredG(0) prop- N
↓
agator contribute to the polaron self-energy. We will ∞ −1
call them GΣ-diagrams. An example is shown in Fig. 8. I = hEIiMC G(0)(τ)dτ . (53)
δ ↓
[Depending on restrictions imposed on the configuration h normiMC (cid:20)Z0 (cid:21)
space (easy to implement in any scheme) physical dia- Imaginary-time integrals for the product of Σ(τ) and
gramswithmorethanoneuncoveredG(0) propagatorare exponentials, see Eqs. (13) and (14), is all we need to
↓
9
determine the polaronenergy and residue. For the bold-
-Γ -Γ
line generalization of the scheme described in the next = 0 -
Section we need the entire dependence of self-energy on
time and momentum. This is achieved by differentiating
partial contributions of GΣ-diagrams. For example, FIG. 10: Diagrammatic equation for the T-matrix Γ(τ,p).
Thearcisthevacuumspin-uppropagatorwiththeconstraint
EΣ,i = δGΣ δτ∈bini , (54) qno<thkinFgobnuittscomrroemcteinntgumtheq.vaTchueummeraensiunlgt oΓf tbhyisseuqbutartaicotninigs
0
is an estimator counting contributions with τ within the contributions from the spin-up fermions with momenta q <
i-th imaginary-time bin of width ∆τ centered at the kF.
i
point τ . Due to linear relation between I and Σ(τ) we
i
immediately realize that (for optimized choice)
Equation(58) is one of the simplest examples ofprob-
∞ −1 lems solvable by bold diagrammatic Monte Carlo [32].
Σ(τ ) = hEΣ,iiMC ∆τ G(0)(τ)dτ . (55) We refer the reader to Ref. [32], where the algorithm of
i δ i ↓
h normiMC (cid:20) Z0 (cid:21) solvingsuchequationsisdescribedindetail. Herewejust
mention some specific details. We find it helpful to start
Incompleteanalogywiththepolaroncase,weconsider
from a good trial function for obtaining high-accuracy
ΓK˜-diagrams that contain one, and only one, uncovered
results in a short simulation time. We achieve this goal
Γ-propagator, see Fig. 9, and use them to collect statis-
using the following protocol. We start the simulation by
tics for the molecule self-energy. Up to straightforward
restricting imaginary time to be smaller than τ and
replacements G(0) Γ, Σ K˜ all relations of this sub- max
↓ ↔ ↔ select relatively short τmax. When the result is accurate
section hold true.
enough, we extrapolate it to longer times, increase τ ,
max
and restart the simulation with the extrapolated func-
tion, Γ serving as the trial function, i.e. we substitute
ext
IV. SIMULATING T-MATRIX Γ(τ,p) BY BOLD Γ=Γ +δΓ to Eq. (58) and solve for δΓ. If necessary,
ext
DIAGRAMMATIC MONTE CARLO
this procedure can be repeated several times.
Despite relatively simple form of Eqs. (23) and (24),
tabulatingthetwo-dimensionalfunctionΓ(τ,p)withhigh
V. SIMULATION RESULTS AND SERIES
accuracy using the inverse Laplace transform of Γ(ω,p) CONVERGENCE PROPERTIES
turns out to be a time consuming job. In this work we
have used an alternative route based on the bold dia-
Nearlyallresults inthis paper wereobtainedby simu-
grammatic Monte Carlo technology introduced recently
latingdiagramsbuiltonbareone-andtwo-particleprop-
inRef.[32]. ThecrucialobservationisthattheT-matrix
agators G(0) and Γ. We observed that the correspond-
Γ(τ,p) can be diagrammatically related to its vacuum ↓
ing series are likely to be divergent. This, however, does
counterpartΓ (τ,p),seeFig.10,withthelatterisknown
0
notmeanthattheentireideaofcalculatingcontributions
analytically:
fromdiagramsofhigherandhigherorderandextrapolat-
Γ (τ,p) = 4π e(ǫF+µ−p2/4m)τg (τ), (56) ing results to the infinite order is useless and ill-defined.
0 m3/2 ∓ On the contrary, it was recognized long ago that appro-
priate re-summation techniques allow one to determine
were
reliablythefunctionstandingbehindthedivergentseries.
1 Moreover, all re-summation techniques (formally, there
g∓(τ) = √E eEτerfc( √Eτ), (57) are infinitely many!), if applicable, have to agree with
−√πτ ± ±
each other on the final result. This important observa-
tionvastlyincreasestheutilityoftheDiag-MCtechnique
fornegative/positivevaluesofthescatteringlength,E =
1/ma2, and erfc(x) is the error-function. [The Fermi- we are developing here. In the next Section we demon-
strate that making the series for Σ self-consistent with
energyandthe chemicalpotentialinEq.(56) come from
the use of Dyson equation—bold-line technique—is an-
the globalenergyshiftnecessaryforcompliancewiththe
other way to improve series convergence properties.
Dyson equation shown in Fig. 10.]
Withtheexplicitexpressionfortheproductoftwovac- For the resonant Fermi-polaron considered here the
uumpropagatorstherelationshowninFig.10reads(the Ces`aro-Riesz summation method solves the convergence
momentum argumentof Γ(τ,p) is suppressedfor clarity) problem. In general,for any quantity of interest—in our
case they are polaron or molecule self-energy—one con-
τ τ structs partial sums
Γ(τ) = Γ (τ) + ds ds′Γ (s)Γ(τ s′)
0 0
− − −
Z0 Zs N∗
dq e−[(p−q)2+q2](s′−s)/2m . (58) Σ(N∗) = DNFN(N∗) , (59)
× (2π)3 −
Zq<kF (cid:16) (cid:17) NX=1
10
-2.44 E
Em
-2.46 a- 1 = 0
-2.48 kFa = 1.0
-0.608
-2.50
-0.61
-2.52
-2.54 -0.612
-2.56
-0.614
-2.58
-0.616
-2.60
-2.62 -0.618
-2.64
0 0.1 0.2 0.3 1/N
*
0.0 0.1 0.2 0.3 0.4
1/N*
FIG. 11: The molecule energy (at kFa = 1) as a function FIG. 12: (Color online) The polaron energy (at the uni-
of the maximum diagram order N∗ for different summation tarity point a−1 = 0) as a function of the maximum dia-
techniques: Ces`aro (open squares), Riesz δ=2 (filled circles, gram order N∗ using Eq. (61). The data are fitted using
fittedwiththeparabolay=−2.6164+0.28013x+0.01638x2), linear −0.618+0.033/N∗ (red) and exponential −0.6151+
Riesz δ = 4 (open circles, fitted with the parabola y = 0.026e−0.39N∗ (black) functions to have an estimate of sys-
−2.6190+0.61635x −0.3515x2), and Eq. (61) (stars fitted tematic errors introduced by theextrapolation procedure.
with thehorizontal dashed line). Reproduced from Ref. [12].
in Fig. (11)]:
defined as sums of all terms up to order N with the N-
∗
th order terms being multiplied by the factor FN(N∗). In F(N∗) = C(N∗) N∗ exp (N∗+1)2 , (61)
thelimitoflargeN andN N themultiplicationfac- N −m(N m+1)
∗ ≪ ∗ m=N (cid:20) ∗− (cid:21)
tors F approach unity while for N N they suppress X
∗
→
higher-ordercontributions in sucha waythat Σ(N∗) has where C(N∗) is such that F(N∗) = 1. The most impor-
1
a well-defined N∗ limit. There are infinitely many tant conclusion we draw from Fig. 11 is that in our case
→∞
ways to construct multiplication factors satisfying these the series are subject to re-summation methods and the
conditions. This immediately leadstoanimportantcon- result of extrapolation is method independent. We con-
sistency check: Final results have to be independent of sider small variations in the final answer due to different
the choice of F. In the Ces`aro-Rieszsummation method re-summation techniques and extrapolation methods as
we have our systematic errors. An example is shown in Fig. 12.
In the next section we will present evidence that the ac-
F(N∗) =[(N N +1)/N ]δ , (Cesa`ro-Riesz). (60)
N ∗− ∗ tual answer is closer to the upper bound of 0.615. We
−
seethatinthe absenceofadditionalinformationone has
Here δ >0 is an arbitrary parameter (δ =1 corresponds
to allow for different ways of extrapolating the answer.
to the Ces`aro method). The freedom of choosing the
Apart from consistency checks, one can test numeri-
value of Riesz’s exponent δ can be used to optimize the
cal results against an analytic prediction for the strong
convergence properties of Σ(N ).
∗
coupling limit k a 0 corresponding to a compact
We proceed as follows. For the series truncated at or- F
→
molecule scattering off majority spins. In this limit the
derN wefirstdeterminethepolaronandmoleculeener-
∗
molecule energy is given by the expression
giesandthenstudytheirdependenceonN asN .
∗ ∗
→∞
In Fig. 11 we show results for the molecule energy at
1 2πa˜
kF a = 1. Without re-summation factors the data are Em =−ma2 −εF + (2/3)mn↑ (kFa →0), (62)
oscillating so strongly that any extrapolationto the infi-
nite diagramorderwouldbeimpossible;weconsiderthis wherethefirsttermisthemoleculebindingenergyinvac-
asanindicationthattheoriginalseriesaredivergent. Os- uum,thesecondtermreflectsfinitechemicalpotentialof
cillations remain pronounced for δ = 1, but are strongly spin-up fermions, and the last term comes from the in-
suppressedforlargervaluesofδ,sothatforδ =4wewere teractionbetweenthecompositemoleculewiththeFermi
notable to resolveodd-evenoscillationsanymore. How- gas. The molecule-fermion s-scattering length a˜ 1.18a
ever, the smoothness of the curve for large δ = 4 comes [35,36]isobtainedfromthenon-perturbativesol≈utionof
at the expense of increasedcurvature, which renders the the three-body problem. Agreement with Eq. (62) pro-
extrapolation to the 1/N∗ 0 limit more vulnerable to vides a robust test for the entire numerical procedure of
→
systematic errors. Empirically, we constructed a factor sampling asymptotic diagrammatic series. Our data are
F(N∗)whichleadstoafasterconvergence[seeanexample in a perfect agreement with the a˜ 1.18a result within
N ≈