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

TL;DR

This paper introduces a stochastic root finding method using the Homotopy Analysis Method for Dyson-Schwinger equations, demonstrating improved convergence and applicability to high-dimensional, non-perturbative problems in quantum field theory.

Contribution

The paper develops a novel stochastic summation algorithm based on rooted-tree diagrams, enhancing convergence and efficiency over existing diagrammatic Monte Carlo methods for solving Dyson-Schwinger equations.

Findings

Superior convergence compared to bold diagrammatic Monte Carlo

Applicable to high-dimensional integral equations

Accesses non-perturbative regimes in quantum field theories

Abstract

We present the construction and stochastic summation of rooted-tree diagrams, based on the expansion of a root finding algorithm applied to the Dyson-Schwinger equations (DSEs). The mathematical formulation shows superior convergence properties compared to the bold diagrammatic Monte Carlo approach and the developed algorithm allows one to tackle generic high-dimensional integral equations, to avoid the curse of dealing explicitly with high-dimensional objects and to access non-perturbative regimes. The sign problem remains the limiting factor, but it is not found to be worse than in other approaches. We illustrate the method for $\phi^4$ theory but note that it applies in principle to any model.

Tables

Table 1. Table 1: Comparing the tree expansion sampled with diagMC to the exact answer where every u m ​ ( x ) subscript 𝑢 𝑚 𝑥 u_{m}(x) is calculated and stored on an external x ∈ [ 0 , 1 ] 𝑥 0 1 x\in[0,1] grid. The last column indicates the probability p 𝑝 p of reaching the expansion order m 𝑚 m .

m	$u_m\left(x=0.5\right)$	$u_m\left(x=0.5\right)$ (diagMC)	p
1	-0.0339665
2	-0.0112517	-0.01127 $\pm$ 0.00004	0.587
3	0.00292815	0.00293 $\pm$ 0.00001	0.247
4	0.000486724	0.000488 $\pm$ 5e-06	0.096
5	-0.000359147	-0.000359 $\pm$ 5e-06	0.041
6	-2.77262e-06	-2.6 e-06 $\pm$ 3.6e-06	0.019
7	4.24149e-05	4.2e-05 $\pm$ 3e-06	0.009

Equations

$$S_E(\varphi )=\frac{1}{2}{\varphi }^2+\frac{\lambda }{4!}{\varphi }^4.$$

$$F(J)=\log ⟨e^{J\varphi }⟩=\log \int \text{d}\varphi e^{-S_E(\varphi )+J\varphi }.$$

$$\frac{\text{d}{S}_E}{\text{d}\varphi }\left[\frac{\text{d}}{\text{d}J}+\frac{\text{d}F(J)}{\text{d}J}\right]=J.$$

$$\begin{array}{rl}{\kappa }_2+\frac{\lambda }{6}{\kappa }_4+\frac{\lambda }{2}{\kappa }_2^2& =1\\ {\kappa }_4+\frac{\lambda }{6}{\kappa }_6+2\lambda {\kappa }_2{\kappa }_4+\lambda {\kappa }_2^3& =0.\end{array}$$

$$\begin{array}{rl}{\kappa }_2=& 1-\frac{\lambda }{2}{\kappa }_2^2-\frac{\lambda }{6}{F}^{\text{pert}}({\kappa }_2)\\ F^{\text{pert}}({\kappa }_2)=& -\lambda {\kappa }_2^3+2{\lambda }^2{\kappa }_2^4-\frac{17}{3}{\lambda }^3{\kappa }_2^5+\mathcal{O}({\lambda }^4).\end{array}$$

$${\kappa }_2^{(n+1)}$$

$${\kappa }_4^{(n+1)}$$

$${\kappa }_4^{(n)}$$

$${\kappa }_2^{(n+1/2)}$$

$${\kappa }_{4,i+1}^{(n)}={\kappa }_{4,i}^{(n)}+\frac{f({\kappa }_{4,i}^{(n)})}{f^{\prime }({\kappa }_{4,i}^{(n)})}.$$

$$\begin{array}{rl}{\kappa }_{4,i+1}^{(n)}& ={\kappa }_{4,i}^{(n)}+\frac{f({\kappa }_{4,i}^{(n)})}{f^{\prime }({\kappa }_{4,i}^{(n)})}\\ & ={\kappa }_{4,i-1}^{(n)}+\frac{f({\kappa }_{4,i-1}^{(n)})}{f^{\prime }({\kappa }_{4,i-1}^{(n)})}+\\ & +\frac{f({\kappa }_{4,i-1}^{(n)}+\frac{f({\kappa }_{4,i-1}^{(n)})}{f^{\prime }({\kappa }_{4,i-1}^{(n)})})}{f^{\prime }({\kappa }_{4,i-1}^{(n)}+\frac{f({\kappa }_{4,i-1}^{(n)})}{f^{\prime }({\kappa }_{4,i-1}^{(n)})})}\\ & =\cdots =F^{\text{tree}}({\kappa }_{4,0}^{(n)}).\end{array}$$

$${\kappa }_2^{(n+1/2)}+\frac{\lambda }{6}{F}^{\text{tree}}({\kappa }_{4,0}^{(n)})+\frac{\lambda }{2}({\kappa }_2^{(n+1/2)}{)}^2=1.$$

$$\mathcal{N}[f]=0.$$

$$(1-q)\mathcal{L}[\varphi (x,q)-f_0]+qh\mathcal{N}[\varphi (x,q)]=0.$$

$$\varphi (x,q)=\sum _m{\frac{1}{m!}\frac{{\text{d}}^m}{{\text{dq}}^m}\varphi (x,q)|}_{q=0}{q}^m.$$

$$\begin{array}{r}f(x)=\varphi (x,q=1)=\sum _m{u}_{f,m}\\ \text{with}{u}_{f,m}={\frac{1}{m!}\frac{{\text{d}}^m}{{\text{dq}}^m}\varphi (x,q)|}_{q=0}.\end{array}$$

$$f(x)=c(x)+{\int }_a^{b}K(x,t)n\left(f(t)\right)dt,$$

$$(1-q)\left[\varphi (x,q)-f_0(x)\right]+q\mathcal{N}[\varphi (x,q)]=0,$$

$$\mathcal{N}[\varphi (x,q)]=\varphi (x,q)-c(x)-{\int }_a^{b}K(x,t)n(\varphi (t,q))\text{d}t.$$

$$u_{f,m}(x)=\frac{1}{(m-1)!}{\frac{{\text{d}}^{m-1}}{{\text{dq}}^{m-1}}{\int }_a^{b}K(x,t)n(\varphi (t,q))\text{d}t|}_{q=0}.$$

$$u_{f,m}$$

$$=$$

$$\times {B_{m-1,k}\left(\varphi (t,q{)}^{\prime },\dots ,\varphi (t,q{)}^{(m-k-1)}\right)|}_{q=0}$$

$$=$$

$$\times B_{m-1,k}\left(u_{f,1}(t),\dots ,(m-k-1)!u_{f,m-k-1}(t)\right).$$

$$\begin{array}{rl}& B_{5,1}(x_1,x_2,x_3,x_4,x_5)=x_5\\ & B_{5,2}(x_1,x_2,x_3,x_4)=10x_2{x}_3+5x_1{x}_4,\end{array}$$

$$\begin{array}{rl}{u}_{f,6}=& \frac{1}{5!}{\int }_a^{b}K(x,t)\sum _{k=1}^2\sum _n^{n_{m-1,k}}{n}^{(k)}(u_0(t))\times \\ & b_n\left(u_{f,1}(t),\dots ,(m-k-1)!u_{f,m-k-1}(t)\right)dt.\end{array}$$

$$\frac{1}{5!}{\int }_a^{b}K(x,t)20\left(2!u_2(t)\right)\left(3!u_3(t)\right)dt$$

$$\begin{array}{rl}& \frac{1}{5!}{\int }_a^{b}K(x,t)20\times \\ & 2!\left({\int }_a^{b}K(t,t^{\prime })2u_0(t^{\prime })u_1(t^{\prime })u_1(t^{\prime })d{t}^{\prime }\right)\times \\ & 3!\left(\frac{1}{2!}{\int }_a^{b}K(t,t^{\prime \prime })2u_1(t^{\prime \prime })u_1(t^{\prime \prime })d{t}^{\prime \prime }\right)dt.\end{array}$$

$$u_1(x)={\int }_a^{b}K(x,t)u_0(t)u_0(t)dt.$$

$$\begin{array}{rl}\int & d{t}_{1}d{t}_{2}d{t}_{3}d{t}_{4}d{t}_{5}d{t}_6\frac{1}{5!}K(x,t_1)20\times \\ & 2!K(t_1,t_3)2u_0(t_3)K(t_3,t_4)2u_0(t_4)u_0(t_4)\times \\ & \frac{3!}{2!}K(t_1,t_2)2K(t_2,t_5)2u_0(t_5)u_0(t_5)\times \\ & K(t_2,t_6)2u_0(t_6)u_0(t_6).\end{array}$$

$$\begin{array}{rl}& d{t}_{1}d{t}_{2}d{t}_{3}d{t}_{4}d{t}_{5}d{t}_{6}32K(x,t_1)K(t_1,t_3)K(t_3,t_4)\times \\ & \times K(t_1,t_2)K(t_2,t_5)K(t_2,t_6)u_0(t_3)u_0(t_4{)}^2{u}_0(t_5{)}^2{u}_0(t_6{)}^2.\end{array}$$

$$f(x)=c(x)+{\int }_0^{1}K(x,t)n(f(t))dt.$$

Full text

I Introduction

Developing first-principles methods for strongly-interacting many-body systems remains an active field of research in theoretical physics. Quantum Monte Carlo algorithms are often the method of choice because of their versatility: Unprecedented insight has been gained with path integral Monte Carlo methods for (bosonic) cold atomic systems RecentDevelopmentsPIMC ; ComparisonMCExperiment ; HiggsBoseHubbard , superfluid and supersolid 4He FateVacancySuperSolid ; SFGrainBoundaries . Determinant Monte Carlo simulations remain indispensable for the nuclear shell model NuclearShell , lattice quantum chromodynamics HybridMonteCarlo , fermions at unitarity UnitaryFermionsI ; UnitaryFermionsII , the Hubbard model detQMCLatticeFermions , topological insulators KaneMeleHubbardAssad ; KaneMeleHubbardWu , and currently certain designer Hamiltonians with gauge fields GaugeFields have been added to this list. They are also used as impurity solvers in dynamical mean-field theory ImpuritySolvers . Diffusion Monte Carlo DiffusionMonteCarloElectronGroundState and full configuration interaction quantum Monte Carlo FICQMCcode ; FICQMCapplication are often used for electronic ground state properties and chemical molecules.

However, in the absence of a positive expansion scheme Monte Carlo algorithms scale exponentially. This is the case for interacting Fermi systems without special symmetry and frustrated spin-systems. Frustration can lead to bad sampling properties already for classical models. It is generally accepted that this exponential scaling is unlikely to be overcome in general SignProblem . One hence tries to develop methods such that the physics can be retrieved before the exponential scaling makes the calculations impossible.

Diagrammatic Monte Carlo (diagMC) is based on the stochastic evaluation of the perturbative Feynman series and recently gained attention: Its greatest appeal is that, for generic fermionic problems, the inevitable sign problem does not lead to an exponential scaling in the system volume but in the expansion order. This typically enables the evaluation of the diagrams up to order 5-10 depending on the problem under consideration, from which the converged answer can hopefully be extracted. When taking into account the range of problems, starting with the polaron problem FermiPolaron ; FermiPolaron3DMass ; FermiPolaron2DVlietinck over models of frustrated spins FrustratedSpins ; SpinIce and to strongly correlated electron systems Hubbard ; AnisotropicHubbard ; PseudoGapHubbard ; ColoumbElectronPhonon ; DiracLiquid , this shows that diagMC is a flexible and quite universal tool. Recently, it has also been used to systematically reintroduce non-local correlations in the dynamical mean-field theory framework DualFermionGull ; DualFermionGukelberger .

One of the biggest challenges for diagMC is the issue of a possible zero convergence radius for which a necessary condition is Dyson’s collapse argument: The Feynman series corresponds to a Taylor series in the coupling constant, and viewing the latter as a complex variable, the system is seen to be unstable against collapse with a zero convergence radius as a result. Dyson formulated the argument originally for quantum electrodynamics, but it applies generally to any bosonic system, in particular to $\phi^{4}$ theory which is studied in this paper. An alternative approach is provided by the skeleton technique known as Hedin’s equations in material science or the Dyson-Schwinger equations (DSEs) in physics. If one can prevent the perturbative expansion of the DSEs (which, in case of $\phi^{4}$ theory, would bring us back to Feynman diagrams subject to Dyson’s collapse Phi4Buividovich ) and instead solve these high-dimensional integral equations directly, progress is possible.

The main result of this paper is to present an algorithm which can be used to solve such high-dimensional integral equations stochastically. The main technique we develop in this paper is the expansion and stochastic summation of the Homotopy Analysis Method (HAM) HAMbook . The HAM is a numerical method to find solutions of non-linear differential and integral equations with a growing number of applications in science, finance, fluid dynamics and engineering HAMbook2 ; HAMbook3 . Therefore, series convergence problems, such as Dyson’s collapse, are avoided if this algorithm is applied to the DSEs as long as the HAM is able to find the solutions of the DSEs. A crucial aspect is that we can prevent the explicit storage and manipulation of $n-$leg vertices (as long as only integrals over such quantities are needed), which for $n\geq 4$ is all but impossible currently. This overcomes the main limitation of numerical methods dealing with DSEs or related skeleton techniques, e.g., functional renormalization group approaches fRGStatMech ; fRGFermi . We note in passing that the recently developed technique of Grassmannization GrassmanizationIsing also prevents Dyson’s collapse.

The paper is organized as follows, featuring an increasing degree of complexity ranging from a zero-dimensional problem to quantum field theory: In the next chapter (Sec. II) we illustrate the key ideas for a coupled set of algebraic equations and solve these issues by introducing 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, we apply the technique to the $\mathcal{Z}_{2}$ symmetric $\phi^{4}$ model in 1D (Sec. IV). Sec. III and IV work out the technical implementations of the ideas presented in Sec. II and can therefore be skipped by readers who are not interested in the detailed implementation. As a non-trivial example we apply the method to the $\mathcal{Z}_{2}$ symmetric $\phi^{4}$ model in 2D (Sec. V). We conclude and provide an outlook in Section VI. Appendix A contains further details of the toy model of Sec. II whereas Appendix B shows how the ideas developed in Sec. II can be formulated for the full DSEs, including the integro-differential equation for the vertex function, showing that the developed algorithm can, in principle, be applied to any model.

II Solution Strategy

The key idea can already be understood from the [math] space-time dimensional case of the $\phi^{4}$ model with action

[TABLE]

In this case the field is reduced to just a single variable $\phi$ and the connected n-point function $G^{(n)}_{c}(\mathbf{x}_{1},\dots,\mathbf{x}_{n})$ reduces to the cumulant $\kappa_{n}$ which is obtained from $\kappa_{n}=\left.\frac{\text{d}^{n}F(J)}{\text{d}J^{n}}\right|_{J=0}$ where

[TABLE]

The differential form of the DSEs can be written as

[TABLE]

From this form of the DSEs the first two non-zero cumulants can be derived by differentiating with respect to $J$,

[TABLE]

This builds an infinite tower of non-linear equations with an infinite number of unknown variables ($\kappa_{2}$ depends on $\kappa_{2}$ and $\kappa_{4}$; $\kappa_{4}$ depends on $\kappa_{2},\kappa_{4}$ and $\kappa_{6}$, etc.), known as the (integral) DSEs.

The infinite tower can be perturbatively expanded in the coupling constant $\lambda$ (cf. Appendix A),

[TABLE]

Such an expansion is in close analogy with existing bold diagrammatic Monte Carlo codes which presently rely on the Luttinger-Ward functional LuttingerWard (or a closely related functional). In these approaches the Luttinger-Ward functional is constructed as the sum of all possible closed diagrams built of (full) 2-point correlation functions and which do not fall apart when cutting any two 2-point correlation function lines. We shall refer to such an expansion as a skeleton series. These are, however, often asymptotic.

In contrast to such a perturbative expansion another strategy is to truncate the infinite tower, e.g. by setting $\kappa_{6}=0$ (this also assumes that 4-point vertices (cf. $\kappa_{4}$) can be dealt with appropriately). The system (4) then yields a closed set of two equations.

The usual procedure to solve such a system stochastically is by fixed point iterations: Starting from initial guesses $\kappa_{2}^{(0)}$ and $\kappa_{4}^{(0)}$ one obtains in iteration $n+1$ an approximate solution which depends on the solution in step $n$:

[TABLE]

Generically, such a fixed point iteration may not have a stable fixed point. It can be checked that already for $\lambda\sim 2$ the above fixed point iteration is diverging. To improve stability we consider each of the equations as non-linear (implicit) equations which are connected by the self-consistency loop shown in Fig. 1,

[TABLE]

The equations (7), (8) are solved by considering each as a root finding problem $f(\kappa^{(n)}_{4})=0$, $g(\kappa^{(n+1/2)}_{2})=0$. In Fig. 1 “Solver I” denotes a root finding algorithm which solves $f(\kappa^{(n)}_{4})=0$ for $\kappa^{(n)}_{4}$ with fixed $\kappa^{(n-1/2)}_{2}$ while “Solver II” solves $g(\kappa^{(n+1/2)}_{2})=0$ for $\kappa^{(n+1/2)}_{2}$ with fixed $\kappa^{(n)}_{4}$.

In order to find the root of $f(\kappa^{(n)}_{4})=0$ “Solver I” can employ the Newton-Raphson method which will be substituted by the Homotopy Analysis Method (HAM) in cases where the DSEs have the form of integral equations (cf. Sec. III). The Newton-Raphson method finds the root corresponding to the procedure,

[TABLE]

The index $i$ denotes the iteration number in the auxiliary 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 improves with the number of iterations $i$ but the result of each iteration step has to be stored in order to use it as a starting point for the next iteration step.

In case of a field theory, $\kappa_{4}$ should be thought of a hard-to-manipulate connected 4-point correlation function $G_{c}^{(4)}(\mathbf{x_{1}},\mathbf{x_{2}},\mathbf{x_{3}},\mathbf{x_{4}})$. To overcome the storage problem, we use a recursive expression for the $(i+1)$-th approximation to the root,

[TABLE]

The expanded result for $\kappa^{(n)}_{4,i+1}$ is a real-valued function $F^{\text{tree}}(\kappa^{(n)}_{4,0})$ depending on the initial guess $\kappa^{(n)}_{4,0}$ for the root finding. It also depends implicitly on $\kappa^{(n-1/2)}_{2}$ as this cumulant was taken fixed in the root finding of “Solver I” and can therefore be viewed as another parameter. In the following such an expansion is referred to as a tree expansion.

We can now proceed with “Solver II” using the tree expansion $F^{\text{tree}}$, which yields

[TABLE]

“Solver II” is used to solve for the unknown $\kappa_{2}^{(n+1/2)}$ after which iteration step $n$ is finished. Since the prime object of “Solver II” (i.e., the 2-point correlation function) is easy to store and manipulate, “Solver II” does not require a tree expansion. This notation makes clear that we can combine both solvers, shown as “Solver I+II” in Fig. 2, whenever the root finding of “Solver I” has been expanded in terms of the function $F^{\text{tree}}$. Compared to Fig. 1 the only change is that $\kappa_{4}$ is represented by a tree expansion and not by a single object.

By dropping the self-consistency index $n$ the final result of (11) can be compared with the skeleton series expansion of (5). The tree expansion inherits the properties of “Solver I” as it is constructed to exactly represent the root finding algorithm. Therefore, by using the tree expansion the often asymptotic skeleton series expansion of bold diagrammatic Monte Carlo is avoided and $F^{\text{tree}}$ is used instead of $F^{\text{pert}}$. It should also be noted that once the infinite tower is truncated by setting $\kappa_{n}=0$ for some fixed $n$ the perturbative expansion of the truncated tower (cf. Eq. (5)) has a finite convergence radius, e.g., for $\kappa_{6}=0$ the convergence radius $r$ is given by $|r|<2\lambda\kappa_{2}$ (cf. (41)).

We have glossed over the computation of the expansion $F^{\text{tree}}$ in Eq. (10), which requires keeping track and summing all the terms generated in the tree expansion to all orders in the recursion and which constitutes a formidable problem for a finite dimensional field theory. Here, we propose a stochastic approach to the tree expansion, in the same spirit as diagrammatic Monte Carlo samples Feynman diagrams. Consequently, the configuration space is given by a collection of diagrams where every diagram is in one to one correspondence to a single term in the tree expansion and a weight that is the product of the individual building blocks. How this Monte Carlo average is exactly implemented in practice is subject of the next section.

III Integral equations and Diagrammatic Monte Carlo

The aim of this section is to arrive at a practical scheme for the generalization of the tree expansion (10). The mathematical framework is provided by the Homotopy Analysis Method (Sec. III.1). A Monte Carlo updating scheme is presented in Sec. III.2, and finally, the method is illustrated for a one-dimensional integral equation where the answer can be compared with the exact one, see Sec. III.3.

III.1 Homotopy Analysis Method

We are interested in finding the roots of a non-linear equation $\mathcal{N}$,

[TABLE]

The idea of the Homotopy Analysis Method (HAM) HAMbook is to rewrite this problem with the help of an embedding parameter $q$ and a convergence control parameter $h$,

[TABLE]

Here, $\mathcal{L}$ is an arbitrary linear operator with $\mathcal{L}[0]=0$. This equation is called the [math]-th order deformation equation. By setting $q=0$ and $q=1$ one sees that the initial guess for the solution of the root finding problem $f_{0}$ is transformed to the full solution $f(x)=\phi(x,q=1)$. Under the assumption that this transformation is smooth and the Taylor expansion is well-defined, one can write

[TABLE]

Therefore,

[TABLE]

The Taylor coefficients $u_{f,m}$ can be obtained by differentiating the [math]-th order deformation equation $m$ times with respect to $q$ and setting $q=0$ afterwards. This can be done analytically yielding a set of deformation equations. It can also be represented graphically with the help of diagrams as will be shown in the following.

The non-linear equation under consideration typically has the form

[TABLE]

with given functions $c$, $n$ and the kernel of the integration $K$. Choosing the HAM convergence parameters $h=1$ and the linear operator $\mathcal{L}$ as the identity operator, the [math]-th order deformation equation is given by

[TABLE]

where

[TABLE]

Taking the derivative with respect to $q$ $m$ times in the [math]-th order deformation equation and taking $q=0$ afterwards gives an equation for $u_{m}=\frac{1}{m!}\frac{\text{d}^{m}\phi(x,q)}{\text{dq}^{m}}|_{q=0}$. For $m\geq 2$ one finds

[TABLE]

This is the starting point for the tree expansion of the HAM. The above equation can be written without specifying the non-linear function $n$,

[TABLE]

Here $B_{m-1,k}$ are the Bell polynomials of the second kind BellPolynomialRef encoding the combinatorial coefficients produced by the action of the derivative $\frac{\text{d}^{m-1}}{\text{d}q^{m-1}}$ on $n(\phi(t,q))$. It can be checked that for the initial guess of the root finding $f_{0}(x)=u_{f,0}(x)=c(x)$ Eq. (20) also holds for $m=1$.

When the root finding algorithm has converged in step $M$, all references to the previous iteration results $u_{f,m}$, $m<M$ should be eliminated from Eq. (20) by using Eq. (20) recursively for every $u_{f,m}$, $m\neq 0$, leading to the desired generalization of Eq. (10). This is again the general procedure leading to the tree expansion. To avoid confusion it should be noted that even though the solution of the non-linear equation is found in the form of $\sum_{m}u_{f,m}$ by the HAM this expansion does not constitute the tree expansion as in order to calculate $u_{f,j}$ all $u_{f,m}$ with $m<j$ have to be calculated and stored, cf. (20). In the sequel we restrict the discussion to a specific non-linear function $n(x)=x^{2}$, which reduces the sum over $k$ in Eq. (20) to $k_{\text{max}}=2$ as $n^{\prime\prime\prime}=0$.

Let us take $m=6$ to set the ideas. The expansion of $u_{f,6}$ involves the Bell polynomials $B_{5,1}$ and $B_{5,2}$ given by

[TABLE]

or $B_{m,k}=\sum_{n=1}^{n_{m,k}}b_{n}$, where $n_{m,k}$ is the number of monomials $b_{n}$ for the polynomial $B_{m,k}$. In this notation Eq. (20) can be written as

[TABLE]

In this form it is evident that all terms in the expansion can be obtained by writing down all possible configurations of $(k,n)$. For example, choosing the configuration $(k=2,n=1)$ yields the term

[TABLE]

by considering the first monomial of $B_{5,2}$ and $n^{\prime\prime}=2$. We have dropped the subscript $f$ because no ambiguity is possible. What we achieved so far is just the first step in obtaining a term in the tree expansion since $u_{2}(t)$ and $u_{3}(t)$ have to be expanded further. Writing down Eq. (22) for $u_{2}(t)$ and $u_{3}(t)$ and choosing for each a new configuration $(k,n)$ eliminates $u_{2}$ and $u_{3}$ from (23). To be concrete, choosing $(k=1,n=1)$ for $u_{2}(t)$ and $(k=2,n=1)$ for $u_{3}(t)$ yields

[TABLE]

At this point only $u_{1}$ remains to be eliminated, for which there is only one possibility,

[TABLE]

Graphically, this elimination procedure is depicted by rooted trees where the basic elements, see Fig. 3, are the roots and leafs of the tree which are connected by branches.

One term in the tree expansion corresponds to a fully grown rooted tree. A random term in the tree expansion is picked by growing a random rooted tree in the following way:

Select a random integer $m$ for the root. 2. 2.

Grow a branch from the root according to some random integer $k$, which fixes the Bell polynomial $B_{m-1,k}$. 3. 3.

Grow leafs from this branch according to some random integer $n$, corresponding to the monomials of the Bell polynomials $B_{m-1,k}$. 4. 4.

Regard every leaf of the branch of $B_{m-1,k}$ as a new root and go back to step 2 if $m>1$ or finish the recursion by using Eq. (25) if $m=1$.

Applied to the example discussed above (see also Fig. 4), the root is $m=6$ which has two different branch types $k$ where $k=2$ is picked. By selecting $n=2$ two leafs are grown on this branch (cf. Eq. 24), and the decomposition ends by invoking Eq. (25). In order to associate to each fully grown rooted tree a term in the tree expansion each element in the rooted tree must correspond to an element in expression (24) according to the following rules:

For each branch from a root with given $m$

1. (a)

there is a factor $\frac{1}{(m-1)!}$. 2. (b)

there is a factor $2u_{0}(t)$ if $k=1$ or $2$ if $k=2$. 3. (c)

there is the prefactor from the randomly picked monomial of $B_{m-1,k}$. 4. (d)

there is an integration over a new variable $t$ and a factor $K(x,t)$. 2. 2.

For each new leaf with label $m\neq 0$ there is a factor $m!$. 3. 3.

For each new leaf with label $m=0$ there is a factor $u_{0}(t)$.

Fig. 5 shows the fully grown, labelled tree, from which the integral can be read off,

[TABLE]

It corresponds to a single term in the tree expansion of the HAM.

The approximation of the root to the $M$-th order is given by $\sum_{m=0}^{M}u_{f,m}$ and therefore the tree expansion consists of generating and evaluating all integral expressions corresponding to all possible rooted trees. As the final answer is given by the sum over all deformations $u_{f,m}$ trees with variable height have to be considered.

In principle this expansion can be done explicitly by drawing and calculating each diagram. But already for moderate $M$ the exponential growth in the number of diagrams renders this approach impractical. Therefore only a stochastic evaluation of the the sum $\sum_{m=0}^{M}u_{f,m}$ is feasible.

III.2 Update structure for the Diagrammatic Monte Carlo sampling of rooted trees

This section introduces the Monte Carlo updates to stochastically sum the tree expansion of the HAM for the above example where $n(x)=x^{2}$. As in a generic diagMC sampling scheme the core of the algorithm is a Markov Chain which changes the topologies and integration variables of the diagrams according to their respective weights. The weight of a diagram is given by the integration kernel of the integral expression corresponding to the diagram. The rooted tree in Fig. 5 corresponds to the integral expression (26) and therefore its weight is

[TABLE]

Furthermore it has been shown that the integral expression can be read off from the rooted tree diagram by using a set of rules which associate to each elementary diagram element, Fig. 3, an analytic expression. With this set of rules changes in the topology and integration variables of a diagram by introducing or removing elementary diagram elements can be related to changes of the weight of the diagram. In the following a set of updates is introduced which generates a Markov Chain Monte Carlo in the space of all possible rooted trees and therefore stochastically sums the tree expansion of the HAM. The updates to calculate $\sum_{m=0}^{M}u_{f,m}(x)$ for a fixed external variable $x$ are:

change integration variable: This update performs a shift of the internal integration variables. It picks one of the internal integration variables and updates the variable according to standard detailed balance rules. As this update is balanced with itself and does not change the diagram structure only the weights of the diagrams have to be taken into account. For the above example the update is schematically depicted in Fig. 6. 2. 2.

shrink-tree/expand-tree: This update is changing the height of the rooted tree.

1. (a)

expand-tree: A new root with $m+1$ will be introduced. From this new root a branch with $k=1$ will be grown and the leaf of this branch is the old root. Furthermore a new integration variable $t_{new}$ will be seeded according to some probability density function $P(t_{new})$ and assigned to the old root. 2. (b)

shrink-tree: This update can only be performed if the branch growing from the root has $k=1$. The root of the rooted tree is deleted and the leaf of the root is assigned to be the new root. The acceptance ratio is the inverse of the expand-tree update. 3. 3.

shrink-tree-cluster/expand-tree-cluster: This is almost the same update pair as the update shrink-tree/expand-tree. Instead of introducing a new root with a new branch in the $k=1$ configuration a branch in the $k=2$ configuration is grown. As can be seen in Fig. 8 this can be done by first growing a new random rooted tree with random height and afterwards glue this tree to the root of the current random tree. These two roots can be regarded as the leafs of a $k=2$ branch grown form the new root of the combined rooted tree. 4. 4.

change-subtree: A random leaf of the rooted tree is chosen from which, regarded as a sub-root, a new subtree is grown. When accepted, the old subtree is deleted and replaced by the new subtree. This update is shown in Fig. 9.

III.3 Illustration

We illustrate the above algorithm for a one-dimensional integral equation,

[TABLE]

The kernel of the integration is $K(x,t)=(x-t)$ and $n(x)=x^{2}$. The function $c(x)$ is picked in such a way that the solution of Eq. (27) is given by $f(x)=\log(x+1)$. It can be checked that $c(x)=\log(x+1)+2\log(2)-2x(\log(2)-1)(\log(2)-1)-\frac{5}{4}$ satisfies this condition. The $m$-th order deformation equation is given by

[TABLE]

The initial approximation of the root is $u_{0}(x)=c(x)$ and $h=H(x)=1$. As is shown in Fig. 10 the result $f(x)\approx f_{M}(x)=\sum_{m=0}^{M}u_{m}(x)$ quickly converges to the correct answer as a function of maximum deformation order $M$. The results of the diagMC root finding calculations for $x=0.5$ are given in Table 1 and compared with the exact answer. As the expansion of $u_{1}$ into $u_{0}$ is immediate it has not been included into the diagMC procedure.

The 4-th column in Table 1 shows the probability $p$ to reach deformation order $m$, indicating the convergence of the diagMC root finding algorithm as $p\to 0$ for larger $m$.

In the next section we extend the proof of principle for the HAM method to the more challenging case of the DSEs for $\phi^{4}$ field theory.

IV Solving Dyson-Schwinger Equations for the One-Dimensional $\phi^{4}$ Theory

Since it is not a priori clear if (and how) the HAM equations and the stochastic sampling of rooted trees work in practice for the coupled set of DSEs, we analyze the HAM method for the DSEs of the 1D $\phi^{4}$ theory in two steps: First, in Sec. IV.2, we use the HAM for “Solver I” and “Solver II” separately, cf. Fig. 1. This is possible because the storage of the vertex function poses no problem in 1D. We demonstrate superior convergence properties compared to fixed point iterations. Second, in Sec. IV.3, we check the stochastic evaluation of the 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 the stochastic solution of the coupled set of DSEs in 2D including both 2-point and 4-point functions and which requires the combination of both steps, and which will be presented in Sec. V.

IV.1 Model and notation

The truncated set of DSEs for the $\phi^{4}$ theory PelsterGlaum is in the thermodynamic limit given by

[TABLE]

For later purposes (cf. Sec. V) we leave the dimensionality $D$ of the lattice, the bare mass $m^{2}$, and the bare coupling $\lambda$ as free parameters. The lattice constant is fixed, $a=1$ and the system size infinite unless otherwise specified. In the above equations the sum runs over all possible lattice points $\mathbf{r_{i}}$ and is denoted by $\int_{i}$. $G_{1,2}$ denotes the 2-point correlation function $G(\mathbf{r_{1}},\mathbf{r_{2}})$ and $\Gamma_{1,2,3,4}=\Gamma(\mathbf{r_{1}},\mathbf{r_{2}},\mathbf{r_{3}},\mathbf{r_{4}})$ denotes the 4-point vertex function. Terms of order $\mathcal{O}(\Gamma^{6})$ are neglected. $G_{0;1,2}$ is the 2-point correlation function of the Gaussian model ($\lambda$=0) given by

[TABLE]

The truncated DSEs, diagrammatically shown in Fig. 11, are regarded as root finding problems and solved by the HAM. The $m$-th order deformations of $G$ and $\Gamma$ are denoted by $u_{G,m}$ and $u_{\Gamma,m}$, respectively, and are given by

[TABLE]

The linear operator $\mathcal{L}$ (cf. Eq. (13)) is chosen to be the identity operator.

The notation is as follows:

The external coordinates are denoted as $\mathbf{x}=(\mathbf{s},\mathbf{t},\mathbf{u})$, where due to translational invariance $(\mathbf{s},\mathbf{t},\mathbf{u})=(\mathbf{r}_{1}-\mathbf{r}_{2},\mathbf{r}_{1}-\mathbf{r}_{3},\mathbf{r}_{1}-\mathbf{r}_{4})$. They are $D$-dimensional vectors. 2. 2.

The sum $\sum_{c}$ runs over the three possible symmetry channels $c=(\mathbf{s},\mathbf{t},\mathbf{u})$ which corresponds to the possible permutations of the external legs, cf. Eq. (29). 3. 3.

The kernel functions $K_{c}$ denote the contribution of the 2-point correlation functions to the vertex function. For example, $K_{c=s}(\mathbf{x},5,6)=G(\mathbf{r}_{1}-\mathbf{r}_{5})G(\mathbf{r}_{1}-\mathbf{r}_{6})$ and $K_{c=t}(\mathbf{x},5,6,7,8,9)=G(\mathbf{r}_{1}-\mathbf{r}_{5})G(\mathbf{r}_{1}-\mathbf{r}_{6})G(\mathbf{r}_{7}-\mathbf{r}_{8})G(\mathbf{r}_{1}-\mathbf{r}_{9})$. 4. 4.

The functions $f_{c}$ denote the linear combinations of positions at which the vertex is evaluated, e.g. $f_{c=s}(\mathbf{x},5,6)=(\mathbf{r}_{5}-\mathbf{r}_{6},\mathbf{r}_{5}-\mathbf{r}_{3},\mathbf{r}_{5}-\mathbf{r}_{4})$ and $f_{c=t}(\mathbf{x},5,6,7)=(\mathbf{r}_{5}-\mathbf{r}_{3},\mathbf{r}_{5}-\mathbf{r}_{7},\mathbf{r}_{5}-\mathbf{r}_{6})$. 5. 5.

$\chi_{m}=1$ for $m>1$ and 0 for $m=1$ and $\tilde{\chi}_{m}=1-\chi_{m}$. 6. 6.

The parameter $h$ is the convergence control parameter of the HAM (cf. Eq. (13)).

With these two equations it is possible to apply the self-consistency loop of Fig. 1. In the case of the 1D $\phi^{4}$ model the vertex depends only on three variables and therefore the vertex can still be stored and the HAM root finding of “Solver I” and “Solver II” in Fig. 1 can be implemented separately.

IV.2 Root finding without tree expansion

All results presented in this section are for the $D=1$ case where the bare mass and bare coupling are fixed to $m^{2}=1$ and $\lambda=10$, respectively. The HAM control parameter is taken as $h=\frac{1}{\lambda}$. The self-consistency loop is applied in the following way (cf. Fig. 1):

In step $n=0$ we start from the guess $G^{(n=0)}=G_{0}$. 2. 2.

With $G^{(n-1/2)}$ “Solver I” is searching for the root of the DSE for the vertex by using Eq. (30) with starting value the second order perturbative result for the vertex. Let the result be $\Gamma^{(n)}$. 3. 3.

With $\Gamma^{(n)}$ “Solver II” calculates the solution of the first DSE by using Eq. (31) with starting value $u_{G,0}=G_{0}$. Let the result be $G^{(n+1/2)}$.

Let us now discuss the results. Fig. 12 and Fig. 13 show the convergence of the HAM root finding of “Solver I” and “Solver II” for the first loop of the self-consistency, $n=1$. The convergence of “Solver I” is demonstrated at a single point $\Gamma^{(n=1)}(0,0,0)$ but holds for any point $(s,t,u)$. The inset shows the divergence of a fixed point iteration which tries to find the root by brute force iteration. Fig. 14 shows the convergence of the self-consistency loop in $G$ as a function of the loop index $n$. As can be seen in more detail in the inset, the self-consistency already converges at $n=3$. Convergence of the 2-point correlation function $G$ in the self-consistency loop implies convergence of $\Gamma$.

Finally, Fig. 15 shows the comparison of the 2-point correlation function with the one obtained by the classical worm algorithm ClassicalWorm , proving that the self-consistency loop converges to the correct answer.

IV.3 Root finding with tree expansion

The goal of this section is to show that the vertex function $\Gamma$ can be obtained by the stochastic evaluation of the tree expansion, i.e. only the DSE for the vertex function is solved with the 2-point correlation function fixed to be the result of the classical worm algorithm. We leave the discussion of the solution for the coupled set to the next section where the DSEs are considered for the 2D setup.

In the following the DSE for the vertex function is considered in momentum space where the $m$-th order deformation equation for $\Gamma(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})$ is given by the Fourier transform of (30). This is again the starting point for the tree expansion which eliminates all references to deformations $u_{\Gamma,i}$, $i<m$ from the $m$-th order deformation equation. Statistics for the tree expansion $\Gamma(0,0,0)=\sum_{m}u_{\Gamma,m}(0,0,0)$ is obtained from the estimator

[TABLE]

where $m(\nu_{\text{tree}})$ is the current height of the rooted tree diagram. The configuration space $\{\nu_{\text{tree}}\}$ of the Monte Carlo (MC) sampling includes rooted tree diagrams with a fixed external momentum configuration, $\mathbf{p}_{1}=\mathbf{p}_{2}=\mathbf{p}_{3}=0$. $\nu_{N_{\text{tree}}}$ is a normalization diagram whose numerical value $N_{\text{tree}}$ is calculated with a deterministic integration algorithm (in practice we take one of the leading terms in the tree expansion).

The tree expansion of (30) suffers from an alternating sign originating from the different signs of the linear and quadratic term with respect to $\Gamma$ in the DSE for the vertex function, cf. Eq (29). For this reason it is impossible to sample the tree expansion to arbitrary large orders. We will discuss the numerical sign problem for our algorithm in more detail in the next section.

The results from the sampling of the tree expansion are compared with the straightforward implementation of the HAM algorithm in Sec. IV.2 where $\Gamma=\sum_{m}u_{\Gamma,m}$ was calculated in real space by storing an manipulating the high-dimensional objects $u_{\Gamma,m}$. Fig. 16 compares the results from the tree expansion at $(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})=(0,0,0)$ with the Fourier transforms of the explicitly stored $u_{\Gamma,m}$. It clearly shows that it is possible to calculate the tree expansion for the case of the 1D $\phi^{4}$ model up to 8 iteration steps with high accuracy.

The tree expansion is a convergent expansion as long as the root finding algorithm is powerful enough to find the solution of the DSE. For the problem under consideration the skeleton series expansion of the Luttinger-Ward functional is asymptotic and therefore breaks down as soon as $\lambda\sim\mathcal{O}(1)$. Fig. 16 shows that the tree expansion considerably increases this parameter regime and seems to be only limited by the sign problem and/or a phase transition (cf. the next section).

V Solving Dyson-Schwinger Equations for the Two-Dimensional $\phi^{4}$ Theory

V.1 Problem Setup

In this section we will move on to the more challenging case of $\phi^{4}$ theory in 2D. For $m^{2}<0$ the theory undergoes a second order phase transition at some non-trivial coupling $\lambda_{c}(m^{2})$ from an ordered phase $0\leq\lambda<\lambda_{c}(m^{2})$ to an unordered phase $\lambda>\lambda_{c}(m^{2})$. The goal of this section is to study this phase transition for $m^{2}=-0.5$.

The bare 2-point correlation function in momentum space $G^{0}(\mathbf{p})$ has poles at $\sum_{i=1}^{D}\sin(\frac{p_{i}}{2})^{2}=\frac{|m^{2}|}{4}$. In order to avoid these poles it is convenient to solve the first DSE for the 2-point correlation function, Eq. (29), in the form $G^{-1}=G_{0}^{-1}-\Sigma$ leading to the root finding problem:

[TABLE]

where the DSE for the vertex function $\Gamma$ is still given by the Fourier transform of the second equation in (29).

In contrast to Sec. IV.2 it is impossible to apply the self-consistency of Fig. 1 straightforwardly as the high-dimensional object $\Gamma$ can no longer be stored and manipulated to high accuracy. In order to solve the coupled equations the stochastic evaluation of the tree expansion is used in the self-consistency of Fig. 2. How the algorithm discussed for the 1D example in Sec. IV.3 is used in the 2D case will be explained in the following.

V.2 Diagrammatic Monte Carlo setup

In the $n$-th self-consistency step, cf. Fig. 2, the solution of Eq. (33) is denoted as $\Gamma^{(2),(n+1/2)}$. This solution is found by using the HAM on (33) for fixed $\Gamma$ which leads to the $m$-th order deformation equation (the loop index $n$ is omitted in $u_{\Gamma^{(2),(n+1/2)},m}$):

[TABLE]

where $u_{\Gamma^{(2)},m}(\mathbf{p})=\frac{1}{m!}\frac{\text{d}^{m}}{\text{d}\text{q}^{m}}\left.\phi(\mathbf{p},q)\right|_{q}=0$, and

[TABLE]

In the function $S[u_{\Gamma^{(2)}},m](\mathbf{p})$ the vertex function $\Gamma$ is given by applying the HAM on the second DSE in (29). According to the self-consistency loop in Fig. 2 the tree expansion is constructed for $G^{-1}=\Gamma^{(2),(n-1/2)}$ in order to calculate $S[u_{\Gamma^{(2)}},m](\mathbf{p})$ in the $n$-th self-consistency step.

In practice the function $S(\mathbf{p})$ is calculated by discretizing the external momentum $\mathbf{p}$ on a grid and applying the diagMC algorithm of IV.3 to sample all possible rooted tree diagrams with variable heights and now variable external momentum. The external momentum of the rooted tree diagrams is updated by importance sampling of the variables $\mathbf{p}$, $\mathbf{k}$ and $\mathbf{l}$ with respect to the integral weight $W[u_{\Gamma^{(2)}},m]\,(\mathbf{p},\mathbf{k},\mathbf{l})$ and histograms are taken for the discrete external momentum points. These histograms are normalized and stored giving a discretized function $\hat{S}(\mathbf{p})$ which is used to compute $\hat{u}_{\Gamma^{(2)},m}(\mathbf{p})$ on the same momentum grid. For calculating the $j$-th order deformation $u_{\Gamma^{(2)},j}$ the deformations $u_{\Gamma^{(2)},m}$ with $m<j$ are needed which can be retrieved by bilinear interpolation of the stored results $\hat{u}_{\Gamma^{(2)},m}$. $H$ is calculated by an independent deterministic numerical integration algorithm which can be considered exact as the numerical errors are subleading to the stochastic diagMC errors.

V.3 Further approximations

Before discussing results for the full solution of the coupled DSEs (29) we introduce various truncations of $\Gamma$ which we will compare against in Sec V.5. Setting $\Gamma=0$ in Eq. (33) transforms the non-linear integral equation into the non-linear algebraic equation

[TABLE]

This equation can be easily solved by tabulating the integral for different $m_{R}$ and using standard numerical methods for solving algebraic non-linear equations.

Another simple truncation is to take $\Gamma=\text{const.}\neq 0$ in momentum space. In order to find a non-trivial fixed point the dimensionless renormalized coupling constant in 2D $\tilde{\lambda}_{R}=\Gamma(0,0,0)m_{R}^{-2}$ has to approach a non-trivial value if the system is tuned close to the phase transition. The next to leading order truncation to (29) satisfying this condition is

[TABLE]

The latter equation corresponds to a resummed ladder expansion at zero external momentum. This set of equations can be solved either by approximating $\Gamma^{(2)}(\mathbf{p})=4\sum_{i=1}^{D}\sin(\frac{p_{i}}{2})^{2}+m_{R}^{2}$, yielding a non-linear algebraic set of equations for $m_{R}$ and $\lambda_{R}$, or without further approximations by using the self-consistency loop, Fig. 1, with “Solver I” just calculating $\lambda_{R}$ and “Solver II” the HAM for $\Gamma^{(2)}(\mathbf{p})$.

V.4 Numerical sign problem

In order to obtain controlled diagMC simulation results to the coupled integral equations of the truncated tower of DSEs the following issues have to be taken into account. The best choice of the free simulation parameters is as follows: The starting value of $\Gamma^{(2)}$, $\Gamma^{(2),(0)}$ in the self-consistency loop, Fig.2, is taken to be the inverse of the 2-point correlation function obtained by the classical worm algorithm. The initial guess for the HAM root finding (34) in the $n$-th self-consistency loop is taken to be $u_{\Gamma^{(2),(n+1/2)},0}=\Gamma^{(2),(n-1/2)}$ and the initial guess for the tree expansion of the HAM is $u_{\Gamma,0}=\lambda_{R}[\Gamma^{(2),(n-1/2)}]$ a constant in momentum space. To determine the maximal order of the tree expansion $\sum_{m=0}^{M}u_{\Gamma,m}$ we calculate $u_{\Gamma^{(2),(1)},1}(\mathbf{p}=0)$, cf. Eq. (34), for different expansion orders $M$ and various convergence control parameters $h$, cf. (30). Examples of these calculations are shown in Fig. 17. The convergence control parameter must be chosen such that the convergence is as fast as possible, i.e. the expansion order is as low as possible. For $\lambda=7$ the choice $h=0.4$ is optimal while a higher value $h=0.6$ leads to an oscillating solution with high error bars. If $h$ is set to too small values (like for the $\lambda=2.5$, $h=0.1$ case) systematic errors in the calculation of the deformations of $\Gamma^{(2)}$ are introduced through the omission of higher order deformations.

The average sign $\langle s\rangle$ of the diagMC integration also depends on the convergence control parameter. While for $\lambda=7$, $h=0.3$, $\langle s\rangle=-0.18$, it goes down to $\langle s\rangle=-0.09$ for $h=0.4$ and drops further to $-0.03$ for $h=0.6$ in the case of an oscillating convergence. As already noted in Sec. IV.3 the sign problem originates from the different signs in front of the linear and quadratic term with respect to $\Gamma$ in the DSE for the vertex function, cf. Eq (29), and therefore all rooted tree diagrams can be classified by having either a positive or negative contribution to $S(\mathbf{p})=S^{+}(\mathbf{p})-\left|S^{-}(\mathbf{p})\right|$ in Eq. (34). Note that $\Gamma^{(2)}$ in the tree expansion is positive definite. The number of rooted tree diagrams in order $m$ grows exponentially with $m$. Therefore, without the weights of the rooted tree diagrams cancelling this exponential growth,

[TABLE]

will grow exponentially. Here, $u^{\pm}_{\Gamma,m}$ denotes the sum of all rooted tree diagrams contributing with a positive (negative) sign in the tree expansion of $u_{\Gamma,m}$. Assuming convergence in the form $S_{m}(\mathbf{p})\rightarrow 0$ for $m\rightarrow\infty$ while $S_{m}(\mathbf{p})$ having a definite sign for all $m$, which is actually the case for a carefully chosen $h$ as can be deduced from Fig. 17, both big numbers $S^{\pm}_{m}(\mathbf{p})$ have to be computed to very high precision in order for their difference to have a high precision. This also explains the fact why $\langle s\rangle$ is smaller for $h=0.4$ compared to $h=0.35$ as convergence sets in faster for $h=0.4$. Thus, with carefully chosen simulation parameters it is possible to extract meaningful results before the statistical errors become too large due to the numerical sign problem. If one has oscillating convergence the sign problem makes it, as expected, impossible to extract controlled results in a realistic simulation time.

Apart from the cancellation of terms in the calculation of $S(\mathbf{p})$ more cancellations occur in the calculation of the deformations $u_{\Gamma^{(2)},m}$. Assuming convergence of the self-consistency loop, i.e. $\Gamma^{(2)(n+1/2)}\approx\Gamma^{(2)(n-1/2)}$, it follows that $u_{\Gamma^{(2)(n+1/2)},1}\approx 0$. But $u_{\Gamma^{(2)(n+1/2)},1}\approx 0$ is obtained by an almost perfect cancellation of the individual contributions in (34) and clearly the absolute error $\Delta u_{\Gamma^{(2)(n+1)},1}$ is given by the statistical error of the diagMC calculation of $S$. As this statistical error is almost constant with respect to the self-consistency loop index $n$ the relative error $\delta u_{\Gamma^{(2)(n+1/2)},1}=\frac{\Delta u_{\Gamma^{(2)(n+1/2)},1}}{u_{\Gamma^{(2)(n+1/2)},1}}$ diverges. It is hence the sign problem that poses the strongest limitation on approaching the critical point more closely.

V.5 Results

In order to compare the quality of the truncations of the infinite tower of DSEs the susceptibility $\chi=\frac{1}{\Gamma^{(2)}(0)}=G(0)$ is calculated in each setup and compared to the results from the classical worm algorithm ClassicalWorm . As the worm algorithm works only on a finite lattice one must ensure that the system size is much larger than the correlation length. The results are shown in Fig. 18. The quality of the truncation increases smoothly with the number of terms taken exactly into account. Deviations with the classical worm remain nevertheless visible, which is attributed to the the omission of higher order vertices such as $\Gamma^{(6)}$, but we checked that the data show consistency (albeit within the large error bars) with the critical exponent $\gamma=7/4$ for the susceptibility.

VI Conclusion and Outlook

In this paper we have treated the Dyson-Schwinger equations (DSEs) as integral equations posing a root finding problem, which we solved by introducing a rooted tree expansion in the Homotopy Analysis Method (HAM) framework, which has better convergence properties than fixed point iterations. We introduced a Monte Carlo sampling procedure to deal with the proliferation of branches and leafs in the tree expansion. Storing high-dimensional objects such as 4-point correlation functions, which naturally appear in the DSEs, can be avoided as long as one is not directly interested in the full knowledge of these quantities. We compared this tree expansion with the skeleton series expansion of the Luttinger-Ward functional in bold diagrammatic Monte Carlo and showed that the convergence properties of this new expansion are superior for the case of $\phi^{4}$ theory. We could go up to correlation length $5$ for the 2D model where further increase is limited by the sign problem.

In future work, the following two major questions must be addressed: The first one is about the quality of the truncation of the DSEs, which has partly already been addressed in previous works on DSEs DSEQCD and within the functional renormalization group community fRGDelamotte ; fRGStatMech ; fRGFermi . It has been shown that truncating the average effective action in the fRG approach by using a low order derivative and field expansion yields already good results for the critical exponents of $\phi^{4}$ models fRGStatMech . We expect this to hold for the DSEs as well since both methods have to yield exactly the same results in the absence of approximations EllwangerDSEvs.fRG .

The second major question concerns the sampling and manipulation of higher order vertices in case one decides to keep them in the expansion. In principle the ideas developed in this paper can be generalized to construct a tree expansion of higher order vertices: There will be another tree each time a higher order vertex appears. If the sign problem remains manageable, the answer can in principle be obtained, but this seems questionable when fluctuations dominate.

Interestingly, the functional form of the integro-differential formulation of the DSEs (cf. PelsterGlaum ) is such that another approach is feasible: $n$-point vertices with $n>4$ can be written as functional derivative terms of the 4-point vertex with respect to the full 2-point correlation function. In the rooted tree expansion of the HAM for the 4-point vertex function the functional derivative terms can be included naturally, as is shown in Appendix B. Hence, high-dimensional objects do not appear and this may be an interesting avenue for future work.

Finally, we note that the DSE formulation can equally well be applied to fermionic systems. It may in particular be interesting to study the Hubbard model near half filling, where the breakdown of the Luttinger-Ward functional and the convergence of the bold diagrammatic Monte Carlo approach to an unphysical solution have been reported Kozik . Our approach should not suffer from these problems.

VII Acknowledgement

This work was supported by FP7/ERC Starting Grant No. 306897 (“QUSIMGAS”). We are grateful to P. Kroiß for useful discussions. TP acknowledges the help on technical questions from P. Kroiß. LP acknowledges support and hospitality from the Erwin Schrödinger Institute, Vienna.

Appendix A Perturbative expansion in $\lambda$ for the DSEs

In this section we derive a perturbative expansion in $\lambda$ for the infinite tower of DSEs. The DSEs for the first three non-zero cumulants are given by

[TABLE]

The general idea can easily be illustrated by first setting $\kappa_{6}=0$. In this case the solution for the second equation of (39) is

[TABLE]

This solution can also be found by recursively plugging the second equation of $(\ref{tower})$ into itself,

[TABLE]

With this result the perturbative solution of (39) up to order $\mathcal{O}(\lambda^{2})$ is given by

[TABLE]

In order to find all terms in such a perturbative expansion one has to know at which order in $\lambda$ the cumulant $\kappa_{2n}$ starts to contribute. From (40) it can be seen that $\kappa_{4}$ is at least of $\mathcal{O}(\lambda)$. The leading order term for $\kappa_{6}$ (cf. 39) is of $\mathcal{O}(\kappa_{4}\lambda)=\mathcal{O}(\lambda^{2})$. The general structure of the DSEs is such that the leading contribution of $\kappa_{2n}$ is $\mathcal{O}(\kappa_{2(n-1)}\lambda)=\mathcal{O}(\lambda^{n-1})$.

As an example, the perturbative expansion up to $\mathcal{O}(\lambda^{4})$ can be constructed by considering the leading order term of $\kappa_{6}$.

[TABLE]

Consequently, $\kappa_{4}$ is given by

[TABLE]

Plugging this result into the first equation of (39) yields the perturbative expansion in $\lambda$ up to order $\mathcal{O}(\lambda^{4})$,

[TABLE]

Appendix B Integro-differential formulation of DSEs

In the functional integro-differential formulation of the infinite tower of DSEs PelsterGlaum the equation for the 4-point vertex function can be written as

[TABLE]

where the kernel functions are given by

[TABLE]

The Homotopy Analysis Method can be applied to this equation leading to the $m$-th order deformation equation

[TABLE]

This equation is again the starting point for the tree expansion where all references to $u_{j}$, $j<m$ will be eliminated in the definition of the $m$-th order deformation equation. To see how the functional derivatives are handled consider for example the following term in the tree expansion of $u_{3}$:

[TABLE]

The functional derivatives can be evaluated easily by using $\frac{\delta u_{0}}{\delta G}=0$ and $\frac{\delta G_{1,2}}{\delta G_{3,4}}=\frac{1}{2}(\delta_{1,3}\delta_{2,4}+\delta_{1,4}\delta_{2,3})$. Due to the product rule for functional derivatives three different terms are generated.

[TABLE]

These terms can be accounted for in the stochastic evaluation of the tree expansion by extending the configuration space with all possible rooted tree diagrams generated by the functional derivative terms.