A stochastic root finding approach: The Homotopy Analysis Method applied to Dyson-Schwinger Equations

2.9 MB·

by  Tobias Pfeffer

A stochastic root finding approach: The Homotopy Analysis Method applied to
Dyson-Schwinger Equations
Tobias Pfeffer and Lode Pollet
Department of Physics, Arnold Sommerfeld Center for Theoretical Physics,
University of Munich, Theresienstrasse 37, 80333 Munich, Germany
(Dated: January 11, 2017)
Wepresenttheconstructionandstochasticsummationofrooted-treediagrams,basedontheex-
pansion of a root finding algorithm applied to the Dyson-Schwinger equations (DSEs). The math-
ematical formulation shows superior convergence properties compared to the bold diagrammatic
7 Monte Carlo approach and the developed algorithm allows one to tackle generic high-dimensional
1 integral equations, to avoid the curse of dealing explicitly with high-dimensional objects and to ac-
0 cessnon-perturbative regimes. The signproblem remainsthelimiting factor, butitis notfound to
2 be worse than in other approaches. We illustrate the method for φ4 theory but note that it applies

n in principle to any model.
a
J
0 I. INTRODUCTION shows that diagMC is a flexible and quite universal
1 tool. Recently, it has also been used to systematically
reintroduce non-local correlations in the dynamical
Developing first-principles methods for strongly-
] mean-field theory framework [29, 30].
h interactingmany-bodysystemsremainsanactivefieldof
One of the biggest challenges for diagMC is the issue of
c research in theoretical physics. Quantum Monte Carlo
e algorithms are often the method of choice because of a possible zero convergence radius for which a necessary
m condition is Dyson’s collapse argument: The Feynman
their versatility: Unprecedented insight has been gained
- with path integral Monte Carlo methods for (bosonic) series corresponds to a Taylor series in the coupling
at cold atomic systems [1–3], superfluid and supersolid 4He constant, and viewing the latter as a complex variable,
t the system is seen to be unstable against collapse with
s [4, 5]. Determinant Monte Carlo simulations remain
. indispensable for the nuclear shell model [6], lattice a zero convergence radius as a result. Dyson formulated
t
a quantum chromodynamics [7], fermions at unitarity the argument originally for quantum electrodynamics,
m but it applies generally to any bosonic system, in
[8, 9], the Hubbard model [10], topological insulators
particular to φ4 theory which is studied in this paper.
- [11, 12], and currently certain designer Hamiltonians
d An alternative approach is provided by the skeleton
with gauge fields [13] have been added to this list. They
n
techniqueknownasHedin’sequationsinmaterialscience
are also used as impurity solvers in dynamical mean-
o
or the Dyson-Schwinger equations (DSEs) in physics. If
c field theory [14]. Diffusion Monte Carlo [15] and full
[ configuration interaction quantum Monte Carlo [16, 17] one can prevent the perturbative expansion of the DSEs
(which, in case of φ4 theory, would bring us back to
are often used for electronic ground state properties and
1
Feynman diagrams subject to Dyson’s collapse [34]) and
v chemical molecules.
instead solve these high-dimensional integral equations
0 However, in the absence of a positive expansion scheme
directly, progress is possible.
8 Monte Carlo algorithms scale exponentially. This is
6 the case for interacting Fermi systems without special The main result of this paper is to present an algorithm
2 symmetry and frustrated spin-systems. Frustration can whichcanbeusedtosolvesuchhigh-dimensionalintegral
0 equationsstochastically. Themaintechniquewedevelop
lead to bad sampling properties already for classical
1. models. It is generally accepted that this exponential in this paper is the expansion and stochastic summa-
0 scaling is unlikely to be overcome in general [18]. One tion of the Homotopy Analysis Method (HAM) [37].
7 The HAM is a numerical method to find solutions of
hence tries to develop methods such that the physics
1 non-linear differential and integral equations with a
can be retrieved before the exponential scaling makes
v: the calculations impossible. growing number of applications in science, finance, fluid
Xi Diagrammatic Monte Carlo (diagMC) is based on the dynamics and engineering [38, 39]. Therefore, series
convergence problems, such as Dyson’s collapse, are
stochastic evaluation of the perturbative Feynman series
r avoided if this algorithm is applied to the DSEs as long
a and recently gained attention: Its greatest appeal is
as the HAM is able to find the solutions of the DSEs.
that, for generic fermionic problems, the inevitable
A crucial aspect is that we can prevent the explicit
sign problem does not lead to an exponential scaling
storage and manipulation of n−leg vertices (as long as
in the system volume but in the expansion order.
only integrals over such quantities are needed), which
This typically enables the evaluation of the diagrams
for n≥4 is all but impossible currently. This overcomes
up to order 5-10 depending on the problem under
the main limitation of numerical methods dealing with
consideration, from which the converged answer can
DSEs or related skeleton techniques, e.g., functional
hopefully be extracted. When taking into account the
renormalization group approaches [44, 45]. We note
range of problems, starting with the polaron problem
in passing that the recently developed technique of
[19–21] over models of frustrated spins [22, 23] and
Grassmannization [33] also prevents Dyson’s collapse.
to strongly correlated electron systems [24–28], this
2
Thepaperisorganizedasfollows,featuringanincreasing
degree of complexity ranging from a zero-dimensional
problem to quantum field theory: In the next chapter
Solver II (n+1/2) (n) Solver I
(oSfeacl.geIIb)rawiceeiqlluuasttrioantes athnedkseoylveidtehaesseforissauecsoubpyleidntsreot- 2 ��(4n) 4 ��(2n�1/2)
� �
ducing the zero-dimensional equivalent of rooted trees.
We proceed with discussing the Homotopy Analysis
Method mathematically and illustrate it for the case of
a one-dimensional integral equation in Sec. III. Next,
FIG. 1. The self-consistency loop to solve the coupled set of
we apply the technique to the Z symmetric φ4 model
2 DSEs. The solution of the coupled equations, cf. Eqs. (7)
in 1D (Sec. IV). Sec. III and IV work out the technical
and(8),afterthen-thiterationintheself-consistencyloopis
implementationsoftheideaspresentedinSec.IIandcan denoted as κ(n), κ(n+1/2).
thereforebeskippedbyreaderswhoarenotinterestedin 4 2
the detailed implementation. As a non-trivial example
we apply the method to the Z2 symmetric φ4 model in coupling constant λ (cf. Appendix A),
2D (Sec. V). We conclude and provide an outlook in
Section VI. Appendix A contains further details of the λ λ
κ =1− κ2− Fpert(κ )
toy model of Sec. II whereas Appendix B shows how the 2 2 2 6 2
(5)
ideas developed in Sec. II can be formulated for the full 17
Fpert(κ )=−λκ3+2λ2κ4− λ3κ5+O(λ4).
DSEs, including the integro-differential equation for the 2 2 2 3 2
vertex function, showing that the developed algorithm
Such an expansion is in close analogy with existing bold
can, in principle, be applied to any model.
diagrammatic Monte Carlo codes which presently rely
on the Luttinger-Ward functional [36] (or a closely re-
lated functional). In these approaches the Luttinger-
II. SOLUTION STRATEGY
Ward functional is constructed as the sum of all possible
closed diagrams built of (full) 2-point correlation func-
The key idea can already be understood from the 0
tions and which do not fall apart when cutting any two
space-time dimensional case of the φ4 model with action
2-point correlation function lines. We shall refer to such
an expansion as a skeleton series. These are, however,
1 λ
S (φ)= φ2+ φ4. (1) often asymptotic.
E 2 4!
In contrast to such a perturbative expansion another
In this case the field is reduced to just a single variable strategy is to truncate the infinite tower, e.g. by setting
φ and the connected n-point function G(n)(x ,...,x ) κ6 = 0 (this also assumes that 4-point vertices (cf. κ4)
c 1 n
can be dealt with appropriately). The system(4) then
reduces to the cumulant κ which is obtained from κ =
dnF(J)(cid:12)(cid:12) where n n yieldsaclosedsetoftwoequations. Theusualprocedure
dJn (cid:12) to solve such a system stochastically is by fixed point it-
J=0
erations: Starting from initial guesses κ(0) and κ(0) one
(cid:90) 2 4
F(J)=log(cid:104)eJφ(cid:105)=log dφe−SE(φ)+Jφ. (2) obtains in iteration n+1 an approximate solution which
depends on the solution in step n:
The differential form of the DSEs can be written as λ λ
κ(n+1) =1− κ(n)− (κ(n))2
(cid:20) (cid:21) 2 6 4 2 2
dS d dF(J)
E + =J. (3) κ(n+1) =−2λκ(n)κ(n)−λ(κ(n))3. (6)
dφ dJ dJ 4 2 4 2
Generically, such a fixed point iteration may not have
From this form of the DSEs the first two non-zero cu-
a stable fixed point. It can be checked that already for
mulantscanbederivedbydifferentiatingwithrespectto
λ ∼ 2 the above fixed point iteration is diverging. To
J,
improve stability we consider each of the equations as
non-linear (implicit) equations which are connected by
λ λ
κ + κ + κ2 =1 the self-consistency loop shown in Fig.1,
2 6 4 2 2
(4)
κ + λκ +2λκ κ +λκ3 =0. κ(n) =−2λκ(n−1/2)κ(n)−λ(cid:16)κ(n−1/2)(cid:17)3 (7)
4 6 6 2 4 2 4 2 4 2
λ λ(cid:16) (cid:17)2
Thisbuildsaninfinitetowerofnon-linearequationswith κ(2n+1/2) =1− 6κ(4n)− 2 κ(2n+1/2) . (8)
an infinite number of unknown variables (κ depends on
2
κ and κ ; κ depends on κ ,κ and κ , etc.), known as The equations (7), (8) are solved by considering each as
2 4 4 2 4 6
the (integral) DSEs. a root finding problem f(κ(n)) = 0, g(κ(n+1/2)) = 0. In
4 2
Theinfinitetowercanbeperturbativelyexpandedinthe Fig.1 “Solver I” denotes a root finding algorithm which
3
solves f(κ(n)) = 0 for κ(n) with fixed κ(n−1/2) while
4 4 2
“Solver II” solves g(κ(n+1/2))=0 for κ(n+1/2) with fixed
2 2
κ(n).
4
In order to find the root of f(κ(n)) = 0 “Solver I” can
4
employ the Newton-Raphson method which will be sub- Solver I+II (n+1/2) Ftree((n))
2 4,0
stituted by the Homotopy Analysis Method (HAM) in
⇣ ⌘
cases where the DSEs have the form of integral equa-
tions (cf. Sec. III). The Newton-Raphson method finds
the root corresponding to the procedure,
f(κ(n))
κ(n) =κ(n)+ 4,i . (9)
4,i+1 4,i f(cid:48)(κ(n)) FIG.2. “SolverI”and“SolverII”fromFig1canbecombined
4,i into “Solver I+I” here. “Solver I+II” is in fact “Solver II”
where κ is represented by the tree expansion of “Solver I”,
Theindexidenotestheiterationnumberintheauxiliary 4
cf. Eq. (11).
Newton-Raphson process and is distinct from the index
n for the iterations in the main self-consistency problem
(i.e., the DSEs). The quality of the solution generically
Fig.2, whenever the root finding of “Solver I” has been
improves with the number of iterations i but the result expanded in terms of the function Ftree. Compared to
of each iteration step has to be stored in order to use it
Fig.1 the only change is that κ is represented by a tree
4
as a starting point for the next iteration step.
expansion and not by a single object.
In case of a field theory, κ should be thought of a
4 By dropping the self-consistency index n the final result
hard-to-manipulate connected 4-point correlation func-
of (11) can be compared with the skeleton series expan-
tion G(c4)(x1,x2,x3,x4). To overcome the storage prob- sion of (5). The tree expansion inherits the properties
lem, we use a recursive expression for the (i+1)-th ap- of “Solver I” as it is constructed to exactly represent the
proximation to the root, root finding algorithm. Therefore, by using the tree ex-
pansion the often asymptotic skeleton series expansion
f(κ(n)) of bold diagrammatic Monte Carlo is avoided and Ftree
κ(n) =κ(n)+ 4,i
4,i+1 4,i f(cid:48)(κ(n)) is used instead of Fpert. It should also be noted that
4,i once the infinite tower is truncated by setting κ =0 for
n
f(κ(n) ) somefixedntheperturbativeexpansionofthetruncated
=κ(n) + 4,i−1 +
4,i−1 f(cid:48)(κ(n) ) tower (cf. Eq. (5)) has a finite convergence radius, e.g.,
4,i−1 (10) forκ6 =0theconvergenceradiusrisgivenby|r|<2λκ2
f(κ(n) + f(κ(4n,i)−1)) (cf. (A3)).
4,i−1 f(cid:48)(κ(n) ) We have glossed over the computation of the expan-
+ 4,i−1
f(cid:48)(κ(n) + f(κ(4n,i)−1)) sion Ftree in Eq. (10), which requires keeping track and
4,i−1 f(cid:48)(κ(n) ) summing all the terms generated in the tree expansion
4,i−1
to all orders in the recursion and which constitutes a
=···=Ftree(κ(n)).
4,0 formidable problem for a finite dimensional field theory.
Here, we propose a stochastic approach to the tree ex-
The expanded result for κ(n) is a real-valued function pansion,inthesamespiritasdiagrammaticMonteCarlo
4,i+1
Ftree(κ(n))dependingontheinitialguessκ(n)fortheroot samples Feynman diagrams. Consequently, the configu-
4,0 4,0
ration space is given by a collection of diagrams where
finding. Italsodependsimplicitlyonκ(n−1/2) asthiscu-
2 everydiagramisinonetoonecorrespondencetoasingle
mulant was taken fixed in the root finding of “Solver I”
terminthetreeexpansionandaweightthatistheprod-
and can therefore be viewed as another parameter. In
uct of the individual building blocks. How this Monte
the following such an expansion is referred to as a tree
Carlo average is exactly implemented in practice is sub-
expansion.
ject of the next section.
We can now proceed with “Solver II” using the tree ex-
pansion Ftree, which yields
λ λ III. INTEGRAL EQUATIONS AND
κ(2n+1/2)+ 6Ftree(κ(4n,0))+ 2(κ(2n+1/2))2 =1. (11) DIAGRAMMATIC MONTE CARLO
“SolverII”isusedtosolvefortheunknownκ(n+1/2)after Theaimofthissectionistoarriveatapracticalscheme
2
whichiterationstepnisfinished. Sincetheprimeobject for the generalization of the tree expansion (10). The
of “Solver II” (i.e., the 2-point correlation function) is mathematical framework is provided by the Homotopy
easy to store and manipulate, “Solver II” does not re- Analysis Method (Sec. III.1). A Monte Carlo updat-
quire a tree expansion. This notation makes clear that ing scheme is presented in Sec. III.2, and finally, the
we can combine both solvers, shown as “Solver I+II” in methodisillustratedforaone-dimensionalintegralequa-
4
tion where the answer can be compared with the exact m≥2 one finds
one, see Sec. III.3.
(cid:12)
1 dm−1 (cid:90) b (cid:12)
u (x)= K(x,t)n(φ(t,q))dt(cid:12) .
f,m (m−1)! dqm−1 a (cid:12)(cid:12)
III.1. Homotopy Analysis Method q=0
(19)
This is the starting point for the tree expansion of the
We are interested in finding the roots of a non-linear
HAM. The above equation can be written without spec-
equation N,
ifying the non-linear function n,
N[f]=0. (12)
1 (cid:90) b dm−1n(φ(t,q))(cid:12)(cid:12)
The idea of the Homotopy Analysis Method (HAM) [37] uf,m(x)= (m−1)! dtK(x,t) dqm−1 (cid:12)(cid:12)
a q=0
is to rewrite this problem with the help of an embedding
parameter q and a convergence control parameter h, 1 (cid:90) b m(cid:88)−1
= dtK(x,t) n(k)(φ(t,q))
(m−1)!
(1−q)L[φ(x,q)−f ]+qhN[φ(x,q)]=0. (13) a k=1
0 (cid:16) (cid:17)(cid:12)
× B φ(t,q)(cid:48),...,φ(t,q)(m−k−1) (cid:12) (20)
Here, L is an arbitrary linear operator with L[0] = 0. m−1,k (cid:12)q=0
This equation is called the 0-th order deformation equa- 1 (cid:90) b m(cid:88)−1
tion. By setting q =0 and q =1 one sees that the initial = dtK(x,t) n(k)(u (t))
(m−1)! 0
guess for the solution of the root finding problem f0 is a k=1
transformed to the full solution f(x)=φ(x,q =1). Un- ×B (u (t),...,(m−k−1)!u (t)).
m−1,k f,1 f,m−k−1
der the assumption that this transformation is smooth
and the Taylor expansion is well-defined, one can write HereB aretheBellpolynomialsofthesecondkind
m−1,k
φ(x,q)=(cid:88) 1 dm φ(x,q)(cid:12)(cid:12)(cid:12) qm. (14) [t4h0e]aecntcioodninogf tthhee cdoemrivbaintiavteoridaml−c1oeoffincinen(φts(tp,rqo))d.ucIetdcabny
m!dqm (cid:12) dqm−1
m q=0 be checked that for the initial guess of the root finding
f (x)=u (x)=c(x) Eq.(20) also holds for m=1.
Therefore, 0 f,0
When the root finding algorithm has converged in step
(cid:88) M, all references to the previous iteration results u ,
f(x)=φ(x,q =1)= u f,m
f,m
m < M should be eliminated from Eq.(20) by using
m
1 dm (cid:12)(cid:12) (15) Eq.(20) recursively for every uf,m, m (cid:54)= 0, leading to
with uf,m = m!dqmφ(x,q)(cid:12)(cid:12) . the desired generalization of Eq. (10). This is again the
q=0 generalprocedureleadingtothetreeexpansion. Toavoid
confusion it should be noted that even though the solu-
TheTaylorcoefficientsu canbeobtainedbydifferen-
f,m tion of the non-linear equation is found in the form of
tiatingthe0-thorderdeformationequationmtimeswith
(cid:80)
u by the HAM this expansion does not consti-
respect to q and setting q = 0 afterwards. This can be m f,m
tute the tree expansion as in order to calculate u all
doneanalyticallyyieldingasetofdeformationequations. f,j
u with m < j have to be calculated and stored, cf.
Itcanalsoberepresentedgraphicallywiththehelpofdi- f,m
(20). Inthesequelwerestrictthediscussiontoaspecific
agrams as will be shown in the following.
non-linear function n(x) = x2, which reduces the sum
The non-linear equation under consideration typically
over k in Eq.(20) to k =2 as n(cid:48)(cid:48)(cid:48) =0.
has the form max
Letustakem=6tosettheideas. Theexpansionofu
f,6
(cid:90) b involves the Bell polynomials B and B given by
5,1 5,2
f(x)=c(x)+ K(x,t)n(f(t))dt, (16)
a
B (x ,x ,x ,x ,x )=x
5,1 1 2 3 4 5 5
(21)
withgivenfunctionsc,nandthekerneloftheintegration
B (x ,x ,x ,x )=10x x +5x x ,
5,2 1 2 3 4 2 3 1 4
K. Choosing the HAM convergence parameters h = 1
and the linear operator L as the identity operator, the orB =(cid:80)nm,kb , wheren isthenumberofmono-
0-th order deformation equation is given by m,k n=1 n m,k
mials b for the polynomial B . In this notation
n m,k
Eq.(20) can be written as
(1−q)[φ(x,q)−f (x)]+qN[φ(x,q)]=0, (17)
0
where 1 (cid:90) b (cid:88)2 nm(cid:88)−1,k
u = K(x,t) n(k)(u (t))×
(cid:90) b f,6 5! a k=1 n 0 (22)
N[φ(x,q)]=φ(x,q)−c(x)− K(x,t)n(φ(t,q))dt. (18)
b (u (t),...,(m−k−1)!u (t))dt.
a n f,1 f,m−k−1
Taking the derivative with respect to q m times in the In this form it is evident that all terms in the expansion
0-th order deformation equation and taking q = 0 after- can be obtained by writing down all possible configura-
wards gives an equation for u = 1 dmφ(x,q)| . For tions of (k,n). For example, choosing the configuration
m m! dqm q=0
5
1. Choose a random m for the root.
u
6
2. Grow a branch with random k.
a) b) c)
u
6
FIG. 3. The basic elements of the rooted trees. a) The roots k=2
andleafsaredrawnasopencircles. b)Branchesareconnect-
ingtherootswiththeleafs. c)Asetconsistingofaroot,leafs
3. Grow leafs according to a random monomial of B .
and branches constitutes a rooted tree. 6,2
u
6
(k =2,n=1) yields the term k=2
u u
2 3
1 (cid:90) b
5! K(x,t)20 (2!u2(t))(3!u3(t))dt (23) 4.Regard the leafs as new roots and start again
a
with Step 2 for each new root.
by considering the first monomial of B and n(cid:48)(cid:48) = 2.
5,2
We have dropped the subscript f because no ambiguity u
6
is possible. What we achieved so far is just the first k=2
stepinobtainingaterminthetreeexpansionsinceu2(t) u2 u3
and u (t) have to be expanded further. Writing down
3
Eq.(22) for u (t) and u (t) and choosing for each a new
2 3
configuration (k,n) eliminates u and u from (23). To
2 3
be concrete, choosing (k = 1,n = 1) for u (t) and (k =
2
2,n=1) for u (t) yields
3
1 (cid:90) b u
6
K(x,t)20×
5! k=2
a u u
(cid:32) (cid:33) 2 3
(cid:90) b
k=1 k=2
2! K(t,t(cid:48))2u0(t(cid:48)) u1(t(cid:48))u1(t(cid:48))dt(cid:48) × (24) u1 u1 u1
a k=1 k=1 k=1
u u u
(cid:32) (cid:33) 0 0 0
1 (cid:90) b
3! K(t,t(cid:48)(cid:48))2u (t(cid:48)(cid:48))u (t(cid:48)(cid:48))dt(cid:48)(cid:48) dt.
2! 1 1 FIG.4. Theconstructionofarootedtreewithheightm=6.
a
Afullygrownrootedtreecorrespondstoasingleterminthe
Atthispointonlyu remainstobeeliminated,forwhich tree expansion which is constructed by recursively applying
1
there is only one possibility, the definition of the root finding, cf. Eq. (20).
(cid:90) b
u (x)= K(x,t)u (t)u (t)dt. (25)
1 0 0 Applied to the example discussed above (see also Fig.4),
a
the root is m=6 which has two different branch types k
Graphically, this elimination procedure is depicted by
where k = 2 is picked. By selecting n = 2 two leafs are
rooted trees where the basic elements, see Fig.3, are
grownonthisbranch(cf. Eq.24),andthedecomposition
the roots and leafs of the tree which are connected by
ends by invoking Eq. (25). In order to associate to each
branches. One term in the tree expansion corresponds
fullygrownrootedtreeaterminthetreeexpansioneach
to a fully grown rooted tree. A random term in the tree
elementintherootedtreemustcorrespondtoanelement
expansion is picked by growing a random rooted tree in
in expression (24) according to the following rules:
the following way:
1. For each branch from a root with given m
1. Select a random integer m for the root.
2. Grow a branch from the root according to some (a) there is a factor 1 .
(m−1)!
random integer k, which fixes the Bell polynomial
(b) there is a factor 2u (t) if k =1 or 2 if k =2.
B . 0
m−1,k
(c) there is the prefactor from the randomly
3. Growleafsfromthisbranchaccordingtosomeran-
picked monomial of B .
m−1,k
dom integer n, corresponding to the monomials of
the Bell polynomials B . (d) there is an integration over a new variable t
m−1,k
and a factor K(x,t).
4. RegardeveryleafofthebranchofB asanew
m−1,k
root and go back to step 2 if m > 1 or finish the 2. Foreachnewleafwithlabelm(cid:54)=0thereisafactor
recursion by using Eq. (25) if m=1. m!.
6
u
6 x u6 x uu76 xtnew
u2 u3 t1 u2 u3 t1
u2 u3 t1 t3 u1 u1 u1 t2 t3 u1 u1 u1 t2
u u u t4 u0 u0 t5 u0 t6 t4 u0 u0 t5 u0 t6
t 1 1 1 t
3 2 FIG. 7. The update to change the height of a rooted tree by
a branch with k=1.
u u u
t t
FIG. 54. Each fully0grown ro0oted trte5e corres0ponds to an6inte- tt43 u0uu12uu01u6u3u1u0 xttt261 uuuu3201 tttt17980 ttt431 u0uu12u0uu1u63uuu1u1300t2 u2uu01 txttt17890
gralexpression,e.g. Eq.(26),inthetreeexpansionwhichcan t5 t5 t6
be read off from the labeled rooted tree.
FIG. 8. The update to change the height of a rooted tree by
u6 x u6 x a branch with k=2.
t3 uu12u1 u3u1 tt21 t3!t03 t03 uu12u1 u3u1 tt21 the above example where n(x) = x2. As in a generic
t4 u0 u0 u0 t6 t4 u0 u0 u0 t6 diagMC sampling scheme the core of the algorithm is a
t5 t5 MarkovChainwhichchangesthetopologiesandintegra-
tion variables of the diagrams according to their respec-
FIG. 6. The update to change an integration variable. The
acceptance ratio is only determined by the different weights tive weights. The weight of a diagram is given by the
of the rooted tree. integration kernel of the integral expression correspond-
ingtothediagram. TherootedtreeinFig.5corresponds
totheintegralexpression(26)andthereforeitsweightis
3. Foreachnewleafwithlabelm=0thereisafactor
u0(t). dt1dt2dt3dt4dt5dt632K(x,t1)K(t1,t3)K(t3,t4)×
Fig.5 shows the fully grown, labelled tree, from which ×K(t ,t )K(t ,t )K(t ,t )u (t )u (t )2u (t )2u (t )2.
1 2 2 5 2 6 0 3 0 4 0 5 0 6
the integral can be read off,
(cid:90) 1 Furthermore it has been shown that the integral expres-
dt dt dt dt dt dt K(x,t )20× sioncanbereadofffromtherootedtreediagrambyusing
1 2 3 4 5 6 5! 1
asetofruleswhichassociatetoeachelementarydiagram
2!K(t1,t3)2u0(t3)K(t3,t4)2u0(t4)u0(t4)× element, Fig.3, an analytic expression. With this set of
(26)
3! rules changes in the topology and integration variables
K(t ,t )2K(t ,t )2u (t )u (t )×
2! 1 2 2 5 0 5 0 5 of a diagram by introducing or removing elementary di-
agram elements can be related to changes of the weight
K(t ,t )2u (t )u (t ).
2 6 0 6 0 6
of the diagram. In the following a set of updates is in-
It corresponds to a single term in the tree expansion of troduced which generates a Markov Chain Monte Carlo
the HAM. in the space of all possible rooted trees and therefore
TheapproximationoftheroottotheM-thorderisgiven stochastically sums the tree expansion of the HAM. The
by (cid:80)Mm=0uf,m and therefore the tree expansion consists updates to calculate (cid:80)Mm=0uf,m(x) for a fixed external
of generating and evaluating all integral expressions cor- variable x are:
responding to all possible rooted trees. As the final an-
swerisgivenbythesumoveralldeformationsu trees 1. change integration variable: This update performs
f,m
with variable height have to be considered. ashiftoftheinternalintegrationvariables. Itpicks
In principle this expansion can be done explicitly by one of the internal integration variables and up-
drawing and calculating each diagram. But already for dates the variable according to standard detailed
moderate M the exponential growth in the number of balance rules. As this update is balanced with
diagrams renders this approach impractical. Therefore
oisnfleyaasibslteo.chastic evaluation of the the sum (cid:80)Mm=0uf,m u6 x uu32 tt021 u6 x
III.2.MoUnptedaCtaerslotrsuacmtuprleinfgorofthreooDteiadgtrraemesmatic tt43 u0uu12uu01 ut53u1u0 ttt261 uu01 tt0506 tt43 u0uu12 uuuu3201 tttt0205061
This section introduces the Monte Carlo updates to FIG. 9. The update to change a random subtree of a rooted
stochastically sum the tree expansion of the HAM for tree.
7
itself and does not change the diagram structure
0.7
only the weights of the diagrams have to be taken
0.6
into account. For the above example the update is
schematically depicted in Fig.6. 0.5
2. shrink-tree/expand-tree: This update is changing 0.4 fM=0
the height of the rooted tree. fM0.3 fM=1
fM=2
(a) expand-tree: A new root with m+1 will be 0.2 fM=3
introduced. Fromthisnewrootabranchwith 0.1 fM=4
k =1willbegrownandtheleafofthisbranch
0.0
istheoldroot. Furthermoreanewintegration
variablet willbeseededaccordingtosome 0.0 0.2 0.4 0.6 0.8 1.0
new
probability density function P(t ) and as- x
new
signed to the old root.
FIG. 10. The successive approximations f (x) =
M
(b) shrink-tree: This update can only be per- (cid:80)M u (x) to the solution of the example integral equa-
formed if the branch growing from the root tiomn=E0q.m(27). Theinitialguessu =f forthesolutionofthe
0 0
has k = 1. The root of the rooted tree is integralequationconvergestothecorrectsolutionalreadyfor
deleted and the leaf of the root is assigned to low M ≈ 4 iteration steps in the HAM.
be the new root. The acceptance ratio is the
inverse of the expand-tree update.
f (x) = (cid:80)M u (x) quickly converges to the correct
3. shrink-tree-cluster/expand-tree-cluster: This is al- M m=0 m
answer as a function of maximum deformation order M.
most the same update pair as the update shrink-
The results of the diagMC root finding calculations for
tree/expand-tree. Instead of introducing a new
x=0.5aregiveninTableIandcomparedwiththeexact
rootwithanewbranchinthek =1configurationa
answer. As the expansion of u into u is immediate it
1 0
branch in the k =2 configuration is grown. As can
has not been included into the diagMC procedure.
be seen in Fig.8 this can be done by first growing
The 4-th column in Table I shows the probability p to
a new random rooted tree with random height and
reach deformation order m, indicating the convergence
afterwards glue this tree to the root of the current
of the diagMC root finding algorithm as p→0 for larger
random tree. These two roots can be regarded as
m.
the leafs of a k = 2 branch grown form the new
In the next section we extend the proof of principle for
root of the combined rooted tree.
the HAM method to the more challenging case of the
DSEs for φ4 field theory.
4. change-subtree: Arandomleafoftherootedtreeis
chosen from which, regarded as a sub-root, a new
subtree is grown. When accepted, the old subtree
IV. SOLVING DYSON-SCHWINGER
is deleted and replaced by the new subtree. This
EQUATIONS FOR THE ONE-DIMENSIONAL φ4
update is shown in Fig.9.
THEORY
III.3. Illustration Since it is not a priori clear if (and how) the HAM
equations and the stochastic sampling of rooted trees
We illustrate the above algorithm for a one-
dimensional integral equation,
m u (x = 0.5) u (x = 0.5) (diagMC) p
m m
(cid:90) 1
1 -0.0339665
f(x)=c(x)+ K(x,t)n(f(t))dt. (27)
0 2 -0.0112517 -0.01127 ± 0.00004 0.587
The kernel of the integration is K(x,t) = (x−t) and 3 0.00292815 0.00293 ± 0.00001 0.247
n(x)=x2. Thefunctionc(x)ispickedinsuchawaythat
4 0.000486724 0.000488 ± 5e-06 0.096
the solution of Eq. (27) is given by f(x) = log(x+1).
5 -0.000359147 -0.000359 ± 5e-06 0.041
It can be checked that c(x) = log(x + 1) + 2log(2) −
2x(log(2)−1)(log(2)−1)−5 satisfiesthiscondition. The 6 -2.77262e-06 -2.6 e-06 ± 3.6e-06 0.019
4
m-th order deformation equation is given by
7 4.24149e-05 4.2e-05 ± 3e-06 0.009
(cid:90) 1 (cid:32)m−1 (cid:33)
(cid:88)
u (x)= dtK(x,t) u (t)u (t) . (28) TABLE I. Comparing the tree expansion sampled with di-
m k m−k−1
0 k=0 agMC to the exact answer where every um(x) is calculated
and stored on an external x ∈ [0,1] grid. The last column
Theinitialapproximationoftherootisu (x)=c(x)and indicates the probability p of reaching the expansion order
0
h = H(x) = 1. As is shown in Fig.10 the result f(x) ≈ m.
8
in two steps: First, in Sec. IV.2, we use the HAM for
= −12 +16 Γ4 “Solver I” and “Solver II” separately, cf. Fig.1. This is
possible because the storage of the vertex function poses
noproblemin1D.Wedemonstratesuperiorconvergence
properties compared to fixed point iterations. Second,
in Sec. IV.3, we check the stochastic evaluation of the
Γ4 = −123× Γ4 + tree expansion up to 8th order (given the exact 2-point
correlation function) against the result obtained in the
first step (Sec. IV.2). These are preliminary steps for
Γ4 the stochastic solution of the coupled set of DSEs in 2D
+613× Γ4 including both 2-point and 4-point functions and which
requires the combination of both steps, and which will
be presented in Sec. V.
FIG. 11. The diagrammatic representation of the truncated
towerofDSEfortheφ4 theory. Thepermutationofexternal
legs is assumed implicitly and indicated by 3×.
IV.1. Model and notation
work in practice for the coupled set of DSEs, we analyze The truncated set of DSEs for the φ4 theory [41] is in
the HAM method for the DSEs of the 1D φ4 theory the thermodynamic limit given by
(cid:90)
G =G + G Σ G
1,2 0;1,2 0;1,3 3,4 4,2
3,4
λ λ(cid:90)
where Σ =− δ(1−2)G + G G G Γ
1,2 2 1,1 6 1,3 1,4 1,5 3,4,5,2
3,4,5
(29)
λ(cid:90)
Γ =λδ(1−2)δ(1−3)δ(1−4)− G G [Γ +Γ +Γ ]+
1,2,3,4 2 1,5 1,6 5,6,3,4 5,6,2,4 5,6,3,2
5,6
λ(cid:90)
G G G G [Γ Γ +Γ Γ +Γ Γ ] .
6 1,5 1,6 7,8 1,9 5,2,7,6 9,8,3,4 5,3,7,6 9,8,2,4 5,4,7,6 9,8,3,2
5,6,7,8,9
For later purposes (cf. Sec. V) we leave the dimension- the 2-point correlation function of the Gaussian model
ality D of the lattice, the bare mass m2, and the bare (λ=0) given by
coupling λ as free parameters. The lattice constant is
fierxwedis,easp=ecifi1eda.ndInthtehesyasbtoevme seiqzueatiniofinnsitteheunsluemss routnhs- G (r)=(cid:90) dDp e−ipr(cid:32)4(cid:88)D sin2(pi)+m2(cid:33)−1 .
(cid:82) 0 (2π)D 2
over all possible lattice points ri and is denoted by i. i=1
G denotes the 2-point correlation function G(r ,r )
1,2 1 2
ThetruncatedDSEs,diagrammaticallyshowninFig.11,
and Γ =Γ(r ,r ,r ,r ) denotes the 4-point vertex
1,2,3,4 1 2 3 4
are regarded as root finding problems and solved by the
function. Terms of order O(Γ6) are neglected. G is
0;1,2
HAM. The m-th order deformations of G and Γ are de-
noted by u and u , respectively, and are given by
G,m Γ,m
9
(cid:34)
λ(cid:88)(cid:90)
u (x)=χ u (x)−h u (x)−λδ(x)χ˜ − K (x,5,6)u (f (x,5,6))+
Γ,m m Γ,m−1 Γ,m−1 m 2 c Γ,m−1 c
c 5,6
λ(cid:88)(cid:90) m(cid:88)−1 (cid:35)
K (x,5,6,7,8,9) u (f (x,8,9))u (f (x,5,6,7)) . (30)
6 c Γ,k c Γ,m−1−k c
c 5,6,7,8,9 k=0
u (1−2)=χ u (1−2)−h[u (1−2)+χ˜G (1−2)
G,m m G,m−1 G,m−1 0
λm(cid:88)−1(cid:90)
− G (1−3)u (0)u (3−2)
2 0 G,m G,m
k=0 3
λm(cid:88)−1(cid:88)k m(cid:88)−1−k(cid:90)
+ u (3−5)u (3−6)×
6 G,i G,k−i
k=0 i=0 j=0 3−7
u (3−7)u (4−2)G (1−3)Γ(7−6,7−5,7−4)]. (31)
G,j G,m−1−k−j 0
0.7 ● ● 5. χm = 1 for m > 1 and 0 for m = 1 and χ˜m =
1−χ .
0.6 1014 ● m
0,0,0) 00..45 ● 0,0,0Γ()|k 1111000015810 ● ● 6. eTtherepofartahmeeHtAerMhi(scft.hEeqco.n(v1e3r)g).encecontrolparam-
u(,mΓ 0.3 ● | 0.1●0 ● ●2 ● ●4 ● 6 8 With these two equations it is possible to apply the self-
consistencyloopofFig.1. Inthecaseofthe1Dφ4 model
0.2 ● k
the vertex depends only on three variables and therefore
●
0.1 ● thevertexcanstillbestoredandtheHAMrootfindingof
●
0.0 ● ● ● ● ● ● ● ● ● ● ● ● ● “Solver I” and “Solver II” in Fig.1 can be implemented
0 5 10 15 20 separately.
m
FIG. 12. The convergence of the deformations u at
Γ,m IV.2. Root finding without tree expansion
(s,t,u)=(0,0,0). The inset shows the divergence of a fixed
pointiterationforthesolutionoftheDSEofthevertexfunc-
tion. Here k is the number of fixed point iterations. All results presented in this section are for the D = 1
case where the bare mass and bare coupling are fixed to
m2 = 1 and λ = 10, respectively. The HAM control
The linear operator L (cf. Eq. (13)) is chosen to be the parameter is taken as h = 1. The self-consistency loop
λ
identity operator. is applied in the following way (cf. Fig.1):
The notation is as follows:
1. Instepn=0westartfromtheguessG(n=0) =G .
0
1. The external coordinates are denoted as x =
(s,t,u), where due to translational invariance
d(sim,te,nus)io=na(lrv1e−ctro2rs,.r1−r3,r1−r4). They are D- 0.00●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●
-0.01●● ● ● ●
2. Tcsihbhalenesnpueemlrsmc(cid:80)u=tcar(tusio,ntns,souov)fewrthhtehiceehxtthceorrernreaepslpolesogsnisbd,lsceft.soyEtmhqme.(pe2to9rs)y-. i()G,m -0.02●● ●● ● i0(=)m ---0000....00003210 ●●●●●●●●●●●●●●●●●●
u -0.03 uG, -0.04 ●
● ● -0.05●
3. The kernel functions K denote the contribution 5 10 15 20
c -0.04
of the 2-point correlation functions to the vertex m
function. For example, K (x,5,6) = G(r −
c=s 1 -0.05●
r )G(r − r ) and K (x,5,6,7,8,9) = G(r −
5 1 6 c=t 1 0 2 4 6 8 10 12
r5)G(r1−r6)G(r7−r8)G(r1−r9). i
4. The functions f denote the linear combinations FIG. 13. Contributions of the m-th order deformation
c
of positions at which the vertex is evaluated, e.g. uG,m(i = |r1 −r2|) for increasing order m, bottom (lowest
f (x,5,6) = (r − r ,r − r ,r − r ) and value of m) to top (highest value of m). The convergence at
c=s 5 6 5 3 5 4
f (x,5,6,7)=(r −r ,r −r ,r −r ). the origin i=0 is shown in the inset.
c=t 5 3 5 7 5 6
10
0.5 proving that the self-consistency loop converges to the
● ● G(n=0) correct answer.
0.4
0.280 ■ G(n=1)
0.3 0) 0.275 ◆▲ ◆ G(n=2)
Gi() ◆■▲ Gi(= 00..226750 ■ ▲ G(n=3) IV.3. Root finding with tree expansion
0.2 0.260
● -0.10-0.05 0.00 0.05 0.10 Thegoalofthissectionistoshowthatthevertexfunc-
i
tionΓcanbeobtainedbythestochasticevaluationofthe
0.1
◆■▲ ● tree expansion, i.e. only the DSE for the vertex function
0.0 ◆■▲ ●◆■▲ ◆●▲■ ◆●▲■ ●◆▲■ ●◆▲■ ●◆▲■ ●◆▲■ ●◆■▲ issolvedwiththe2-pointcorrelationfunctionfixedtobe
0 2 4 6 8 10 the result of the classical worm algorithm. We leave the
i discussion of the solution for the coupled set to the next
section where the DSEs are considered for the 2D setup.
FIG. 14. The approximations for the full 2-point correlation In the following the DSE for the vertex function is con-
function G(i = |r1 −r2|) after n = 0,...,3 self-consistency sidered in momentum space where the m-th order defor-
steps. The inset shows in detail the fast convergence of the mation equation for Γ(p ,p ,p ) is given by the Fourier
1 2 3
self-consistency loop at i=|r −r |.
1 2 transform of (30). This is again the starting point for
the tree expansion which eliminates all references to de-
0.500● formations u , i<m from the m-th order deformation
◆■ Γ,i
● equation. Statistics for the tree expansion Γ(0,0,0) =
00..015000 ◆■ ● (cid:80)muΓ,m(0,0,0) is obtained from the estimator
Gi() 0.010 ◆■ ● ● (cid:104)uΓ,m(0,0,0)(cid:105)MC = (cid:104)(cid:104)δδm,m(νtree)(cid:105)(cid:105)MCNtree, (32)
0.005 ● GaussianModelλ=0 ◆■ ● νtree,νNtree MC
0.001 ■ ClassicalWormλ=10 ◆■ wherem(νtree)isthecurrentheightoftherootedtreedia-
gram. Theconfigurationspace{ν }oftheMonteCarlo
5.×10-4 ◆ DSE+HAMλ=10 tree
◆■ (MC)samplingincludesrootedtreediagramswithafixed
external momentum configuration, p = p = p = 0.
0 1 2 3 4 5 1 2 3
i νNtree is a normalization diagram whose numerical value
N is calculated with a deterministic integration algo-
tree
FIG.15. Theresultforthe2-pointcorrelationfunctionG(i= rithm(inpracticewetakeoneoftheleadingtermsinthe
|r −r |) from the self-consistency loop (HAM+DSE) and tree expansion).
1 2
from the classical worm algorithm [47]. The tree expansion of (30) suffers from an alternating
sign originating from the different signs of the linear and
quadratic term with respect to Γ in the DSE for the ver-
2. With G(n−1/2) “Solver I” is searching for the root
tex function, cf. Eq (29). For this reason it is impossible
of the DSE for the vertex by using Eq. (30) with
to sample the tree expansion to arbitrary large orders.
starting value the second order perturbative result
We will discuss the numerical sign problem for our algo-
for the vertex. Let the result be Γ(n).
rithm in more detail in the next section.
The results from the sampling of the tree expansion are
3. WithΓ(n) “SolverII”calculatesthesolutionofthe
compared with the straightforward implementation of
first DSE by using Eq. (31) with starting value
(cid:80)
the HAM algorithm in Sec.IV.2 where Γ = u
uG,0 =G0. Let the result be G(n+1/2). m Γ,m
was calculated in real space by storing an manipulating
Let us now discuss the results. Fig.12 and Fig.13 show the high-dimensional objects u . Fig.16 compares the
Γ,m
the convergence of the HAM root finding of “Solver I” results from the tree expansion at (p ,p ,p )=(0,0,0)
1 2 3
and “Solver II” for the first loop of the self-consistency, withtheFouriertransformsoftheexplicitlystoredu .
Γ,m
n = 1. The convergence of “Solver I” is demonstrated It clearly shows that it is possible to calculate the tree
at a single point Γ(n=1)(0,0,0) but holds for any point expansion for the case of the 1D φ4 model up to 8 itera-
(s,t,u). The inset shows the divergence of a fixed tion steps with high accuracy.
point iteration which tries to find the root by brute The tree expansion is a convergent expansion as long as
force iteration. Fig.14 shows the convergence of the the root finding algorithm is powerful enough to find the
self-consistency loop in G as a function of the loop index solution of the DSE. For the problem under considera-
n. As can be seen in more detail in the inset, the self- tion the skeleton series expansion of the Luttinger-Ward
consistency already converges at n = 3. Convergence of functional is asymptotic and therefore breaks down as
the 2-point correlation function G in the self-consistency soon as λ∼O(1). Fig.16 shows that the tree expansion
loop implies convergence of Γ. Finally, Fig.15 shows considerably increases this parameter regime and seems
the comparison of the 2-point correlation function with to be only limited by the sign problem and/or a phase
the one obtained by the classical worm algorithm [47], transition (cf. the next section).