Full and unbiased solution of the Dyson-Schwinger equation in the functional integro-differential representation

We provide a full and unbiased solution to the Dyson-Schwinger equation illustrated for $\phi^4$ theory in 2D. It is based on an exact treatment of the functional derivative $\partial \Gamma / \partial G$ of the 4-point vertex function $\Gamma$ with respect to the 2-point correlation function $G$ within the framework of the homotopy analysis method (HAM) and the Monte Carlo sampling of rooted tree diagrams. The resulting series solution in deformations can be considered as an asymptotic series around $G=0$ in a HAM control parameter $c_0G$, or even a convergent one up to the phase transition point if shifts in $G$ can be performed (such as by summing up all ladder diagrams). These considerations are equally applicable to fermionic quantum field theories and offer a fresh approach to solving integro-differential equations.


I. INTRODUCTION
Despite decades of research there continues to be a need for developing novel methods for strongly correlated systems. The standard Monte Carlo approaches [1][2][3][4][5] are convergent but suffer from a prohibitive sign problem, scaling exponentially in the system volume [6]. Diagrammatic Monte Carlo simulations [7][8][9] were developed to prevent this, scaling exponentially only in the expansion order [10]. However, this happens at the expense of substantially worse series convergence properties [11][12][13]. After nearly ten years and despite recent and tremendous progress [14][15][16], one may well fear that the combination of an asymptotic/divergent series and a mild sign problem is as prohibitive as the standard approaches.
Recently [17], we therefore suggested to use the more flexible Dyson-Schwinger equation (DSE) instead of self-consistent Feynman diagrams [18] to provide a fully self-consistent scheme on the one and two particle level. Furthermore, we extended the homotopy analysis method (HAM) [19,20] to φ 4 field theory in two dimensions (2D) (providing us with more tools to enhance the convergence properties in a systematic way), and showed how the expansion in terms of rooted trees is amenable to a systematic Monte Carlo sampling. This expansion is furthermore convenient when dealing with multidimensional objects such as the four-point vertex function. Clearly, this is a radically different way at looking at interacting field theories. However, in our previous work [17], we truncated the DSE at the level of the six-point vertex. The infinite tower of equations for n-point correlation functions was not solved and differences with the full, exact answer could be seen when the correlation length increases.
In this work, we solve the full DSE by writing them as a closed set of integro-differential equations. Within the HAM theory there exists a semianalytic way to treat the functional derivatives without resorting to an infinite expansion of the successive n-point correlation functions, cf. Refs. [21,22] where the DSE has been used to generate new weak coupling expansions. The unbiased numerical solution of the DSE is largely unexplored as it was considered to be too complex to be solved even in the simplest cases [23][24][25]. Furthermore, taking into account the functional derivatives deteriorates the convergence properties of the field theory substantially compared to the truncated case considered in Ref. [17]. Here we show how the remaining theory can be brought under control within the HAM as an asymptotic expansion of the HAM deformations in terms of an auxiliary convergence control parameter c 0 (times the two-point correlation function G) around G = 0, or even as a convergent expansion in the HAM deformations when a shift of G is possible, e.g., by solving the ladder equations. Although the ideas are illustrated for a 1d integral and φ 4 theory in 2D the convergence considerations and the methodology are generically applicable as long as the two-point and four-point correlation functions are bounded, and may just as well be applied to the better known Hedin [26] and parquet equations [27] or to the functional renormalization group equation [28,29].
The paper is organized as follows. The main ideas and results are introduced in Secs. II-IV. In Secs. II and III, we first analyze the toy model (i.e., the 1d integral, cf. Ref. [30]) to illustrate the approach and then proceed in Sec. IV with the full solution of the DSE for the φ 4 model in 2D. We provide further technical details, derivations, and information about the developed algorithm in Appendixes A-C. We do not include these details in the main text in order to provide the reader with a clear introduction to our new approach. The appendixes should be consulted in order to judge the validity and broad applicability of the results presented in the main text.

II. FUNCTIONAL CLOSURE
The functional closure is most easily demonstrated in the case of a 0D field theory. Consider the 1D integral where Z[J ] is the generating functional of the n-point correlation functions G (n) , . (2) Although in this example the generating functional Z[J ] is a real-valued function and the n-point correlation functions are real numbers, we will, nevertheless, use the terminology of (quantum) field theory (FT) as there will be no ambiguities. Moreover, we use the shorthand notation G = G (2) for the two-point correlation function. The DSE [31] can be derived by introducing an infinitesimal shift δ in the integration variable, φ → φ + δ, and expanding the resulting expression in powers of δ. This yields The first derivative of the action S with respect to the field φ is promoted to a differential operator by the substitution φ = d dJ . The DSE (3) is a definition of the generating functional in terms of a differential equation equivalent to the definition of Z[J ] through (1). For a realistic FT (3) turns into a functional integro-differential equation, whereas (1) turns into a functional integral.
Instead of focusing on the solution of (1) we focus on (3). Differentiating (3) once with respect to J and setting J = 0 afterwards yields, after introducing the connected four-point correlation function G (4) c = G (4) − 3G 2 and the four-point vertex function = −G −4 G (4) c , This is the first equation in the expansion of the differential equation (3) into an infinite hierarchy of coupled equations for the correlation functions. We close the hierarchy by considering the generating functional to be a functional of the inverse noninteracting two-point correlation function k, Z = Z[J = 0, k]. Therefore we can write G (4) = −2 dG dk + G 2 and We use (4) to expand the derivative of the inverse of the twopoint correlation function with respect to k. This together with the chain rule where we denote the derivative with respect to G as , leads to the differential equation A solution to (4) requires knowledge of the universal functional [G], which can be obtained by solving the differential equation (7).
[G] is subsequently used in (4) and solved for G for a given inverse of the noninteracting two-point correlation function k. The combined system (4) and (7) is typically solved by a fixed point iteration. This gives the physical two-point correlation function G. The physical fourpoint vertex function follows by evaluating the universal functional [G] for the physical G. In Ref. [17], we have already shown how to solve a realistic FT when only taking the first three terms on the right-hand side of (7) into account, i.e., working with a truncated version of [G], which has a finite convergence radius. Next, we examine the deterioration in convergence properties when taking the derivatives in (7) into account.

III. VERTEX EQUATION
The exact solution of (7) obtained by the implicit Euler method is shown in Fig. 1 ]. The additional vertical (horizontal) grid line corresponds to the physical G( ) obtained for the model parameters k = 1, λ = 10. In case of a FT, (7) takes the form of a functional integro-differential equation and the implicit Euler method cannot be applied. It is therefore necessary to study semianalytic solutions to the differential equation (7).
The power series expansion [G] = λ + lim M→∞ M n=1 c n G n has zero convergence radius [32]: bringing (7) into the form of [G] = F [G, [G]], we find that F is not holomorphic at G = 0, = λ. An expansion around G = 0 corresponds to approaching the noninteracting limit by fixing λ and taking k → ∞, i.e., , the ladder approximation to the vertex equation (7). The functional [G] is evaluated at the physical G for k = 1. The inset shows the convergence of the main plot on a much finer scale. solution of the differential equation (7), [G, 1] = [G], at q = 1. u ,0 [G] is the initial guess for the solution of N [ [G]] = 0. The convergence control parameter c 0 controls the rate at which the deformation happens. The HAM attempts to find the solution of (9) through a Taylor series expansion in q, i.e., [G, q] = u ,0 [G] + m=1 u ,m [G] q m . The expansion coefficients are given by u ,m [G] = (m) [G, q = 0]/m! and can be obtained by the mth derivative of (9) with respect to q. Therefore the HAM gives a series solution of the differential equation (7) in terms of the deformation coefficients u ,m [G], We first use the easiest possible linear operator . This choice can be straightforwardly generalized to functional integro-differential equations. In this case, the mth-order deformation in the HAM series solution is given in simple powers of G and additionally depend on the auxiliary parameter c 0 . Picking the identity operator, L = id, is a too simplistic choice and does not lead to a convergent HAM series solution. Nevertheless, the homotopy introduces the convergence parameter c 0 which gives us additional freedom since the effective parameter c 0 G can always be made small as long as G remains finite and therefore the limitations of the standard power series approach can be overcome. As is illustrated in Figs. 1(a) and 1(c), the asymptotic nature can then be postponed to larger expansion orders for smaller c 0 but with larger deviations from the exact result at low expansion orders. This approach can be generalized to realistic field theories in which case c 0 ||G|| can be chosen to be small where ||G|| denotes the L p norm of the two-point correlation function. In the next section, we show numerically that this holds for a generic FT and that it is possible to get accurate results even close to a second-order phase transition at which ||G|| gets unbounded. Before extending this approach to a realistic FT we like to emphasise the great potential and novelty of the Dyson-Schwinger and homotopy approach to FT. Instead of constructing the HAM series expansion around the noninteracting limit, G = 0, we can use the linear operator L = L ladder to construct a homotopy with respect to the analytically solvable ladder approximation, In an FT, the term 3λ 2 G 2 [G] sums up all RPAladder diagrams. The resulting homotopy can be written as which should be compared to the closed DSE (7). Taking into account the equivalence between the functional integral and the functional integro-differential approach, an important question is if (11) can be related back to an auxiliary action such as S = S 0 (q ) + qS int , which has been considered in variational perturbation theory approaches [33,34] or shifted action approaches [35]. S 0 (q ) is some quadratic auxiliary action and S int is the interaction part such that S(q = 1) = S with S the physical action of the theory under consideration. A general feature of DSEs is that each term on the righthand side of (7) contributes exactly with one factor of the bare coupling constant. In contrary, for the homotopy (11), there are terms contributing with a modified bare coupling constantλ and therefore (11) can not be associated with the closed DSE of an auxiliary action. In conclusion, we see that the HAM series solution is effectively an expansion with respect to the truncated interacting ladder model which has no correspondence in a functional integral representation and can only be written down in the equations of motion approach. Figure 1(d) shows that the HAM series with respect to the ladder approximation yields a convergent series solution. Moreover, we find that the HAM series solution is globally convergent, i.e., it is not restricted by a finite, c 0 -dependent convergence interval. Therefore we also find a convergent solution for values of G, which correspond to G[k] with k < 0. In this parameter regime, the action (1) has degenerate minima and the standard perturbative series for the vertex function is not Boreal summable [36]. For further details see Appendix A.

IV. φ 4 THEORY IN 2D
Consider the Z 2 -symmetric φ 4 theory on a 2D lattice with action The inverse noninteracting two-point correlation function is given by G −1 0;i,j = − i,j + m 2 δ i,j , where i,j denotes the discretized Laplace operator in 2D. For m 2 < 0, this model undergoes a second order phase transition from a magnetically ordered to an unordered phase at a critical coupling constant λ c (m 2 ). The phase transition is signalled by the divergence of the magnetic susceptibility χ = G(p = 0), which leads to the same critical exponents as for the 2D Ising model previously studied diagrammatically with grassmannization techniques [37]. In our simulation, we use m 2 = −0.5, which corresponds to λ c (m 2 ) ∼ 2. The coupled set of equations for (12), generalization of (4) and (7) and derived in Appendix B, are graphically depicted in Fig. 2.
We use an extension of the Monte Carlo algorithm developed in Ref. [17] to sample the expansion of the HAM in rooted tree diagrams. The detailed algorithm, especially the correct implementation of the functional derivatives, is discussed in Appendix C.  (12) is closed through a functional integro-differential equation. Each diagram is in one-to-one correspondence with the terms in (4) and (7). The noninteracting two-point correlation function G 0;i,j is denoted by a thin line, the two-point correlation function G i,j by a bold line, and the bare vertex by a dot. The correct convolution of lattice indices can be obtained by standard diagrammatic rules. The terms without functional derivatives involve the permutation of external indices and is denoted by the factor 3. set of equations, Fig. 2, are defining the FT (12) directly in the thermodynamic limit, i.e., results obtained by our method are for infinite system sizes. The results are compared to the numerically exact simulation of model (12) with the classical Worm algorithm [38] on system sizes which are much larger than the correlation length. We find that it is possible to obtain controlled results with our current Monte Carlo algorithm up to χ ∼ 10, which corresponds to a correlation length of ξ ∼ 3. Our diagrammatic Monte Carlo sampling is based on a direct sampling of all topologies of rooted tree diagrams at a given order. This naive approach leads to sampling problems at higher deformation orders and restricts the order to 5-6. Due to the sign alternation between the linear and quadratic terms in the vertex equation the sampling suffers from a sign problem discussed in Ref. [17] for the case of the stochastic construction of a truncated functional. We find no qualitative difference if the functional derivative terms are added. We performed extensive Monte Carlo simulations such that the error bars are exclusively determined by the extrapolation of the lowest deformation orders as is shown in Fig. 4. It is only possible to extrapolate the inverse of the two-point correlation function and therefore the error introduced through the extrapolation is growing rapidly as the susceptibility diverges, cf. Fig. 3. It should be pointed out that we have no global excess to [G] but only to a stochastic, local evaluation of the universal functional through the diagrammatic Monte Carlo sampling of the HAM series solution in terms of rooted tree diagrams. Further details about the algorithm are presented in Appendix C. In order to obtain the results in  Fig. 2 at the numerically exact two-point correlation function. This reduces the computational time as we have to perform only a single fixed point iteration step. For more general starting points of the fixed point iteration we refer to Ref. [17]. Figures 3 and 4 show that, although there is no single physical small parameter close to the phase transition, we obtain controlled results as we can construct the homotopy (9) always with respect to a small enough c 0 . In the vicinity of the second-order phase transition, ||G|| → ∞ implies that c 0 → 0 but in order to get meaningful results, the number of required deformations increases rapidly. Therefore, due to the limited number of deformations, which can be computed with our algorithm, we can not go closer to the phase transition. For transitions where ||G|| remains finite (the divergence may then occur in the four-point vertex function), this argument is not applicable.

V. OUTLOOK
Although the approach to FT introduced in this work is illustrated for the simple though representative case of φ 4 theory, models with arbitrary two-body interactions in its symmetric phase can be tackled on the same footing. As shown in Appendix B, the transition from the functional integral representation for general models to functional integro-differential equations yields exactly the same set of equations as depicted in Fig. 2. The only difference is that the convolution of indices in the diagrams representing the functional integro-differential equation is with respect to a collective index i, which summarizes all possible field labels of the considered model. The statistics of fermionic fields translate into additional signs for the permutations of external indices in the diagrams of Fig. 2 and into sign alternating two-point correlation functions.

VI. CONCLUSION
In conclusion, we have introduced a general approach to tackle FT through full and unbiased solutions of functional integro-differential equations derived from the DSE. We showed for a toy model that by using the semianalytic HAM we can solve the differential equation for the vertex function. The HAM gives a convergent series solution if the homotopy is constructed with respect to the analytically solvable ladder approximation for the four-point vertex function or an asymptotic series solution, which can be controlled by the convergence control parameter c 0 if the homotopy is constructed around the noninteracting limit. The latter result found for the toy model can be readily applied to a simple FT where the asymptotic series solution can be controlled by the convergence control parameter c 0 even close to a second-order phase transition.
The open data for this project can be found at Ref. [39].

ACKNOWLEDGMENTS
The authors would like to thank D. Hügel, A. Toschi and the participants of the Diagrammatic Monte Carlo Workshop held at the Flatiron Institute, June 2017 for fruitful discussions and valuable input. This work is supported by FP7/ERC Consolidator Grant No. 771891, and the Nanosystems Initiative Munich (NIM).

APPENDIX A: ANALYTIC STRUCTURE OF THE UNIVERSAL FUNCTIONAL [G]: CONVERGENT HAM SERIES SOLUTION
In order to obtain the analytic structure of [G] near G = 0, we first consider the analytic solution to the integral The analytic solution to this integral for k ∈ C with fixed λ ∈ R + is given as Here, I n (z)/K n (z) are the modified Bessel functions of the first/second kind and is the Gamma function. The integral (A1) can be regarded as an integral representation of the function Z defined in (A2). The same holds for the two-point correlation function G[k], which we do not show here. The important point is that both Z[k] and G[k] are piecewise-defined functions in k, which can be represented by a single integral expression. There exists also a differential equation equivalent to the integral representation of the function this is the coupled system of differential equations: There are two independent solutions to this equation for G and, consequently, also for denoted as G ± [k] and ± [k], respectively. While (A1) is single-valued, the solution to the system of equations (A3) in principle leads to two independent solutions. They are shown for G ± [k], ± [k] on the real axis in Fig. 5. We will show in the following that the analytic structure of [G] near G = 0 can be understood from the existence of the two independent solutions to (A3). Formally, [G] can be obtained from ± [k] by inverting the relation . This leads to two branches for [G] where [G] must evaluate to (4) + [k 1 ] on the first branch and to (4) − [k 2 ] on the second branch. In order to obtain the analytic structure of [G] near |G| = 0, it is necessary to study the intersection of the images of G ± for elements which satisfy |G| → 0. The complete circle in the complex G plane |G|e iφ , φ ∈ [0, 2π ] for |G| → 0 in Fig. 6 is included in the image of G + [k] and can be obtained by the parametrization k = |k|e −iφ where |k| → ∞. On the other hand, the image of G − [k] includes two segments of the circle |G|e iφ , namely, with φ ∈ [ π 4 , 3π 4 ] and φ ∈ [ 5π 4 , 7π 4 ] and can also be obtained by the parametrization k = |k|e −iφ . Therefore there are branch cuts starting from G = 0 cutting the complex G plane in the ± π 4 direction. According to the above result a power series solution aroundG ∈ R + , should have convergence radius R = sin( π 4 )G. The convergence radius can also be determined numerically by considering the coefficients c n in the power series [G] = n c n (G −G) n . The power series can be obtained by using the analytically known result for (4) ] on the real axis, i.e., c 0 = [G]. As shown in Fig. 7, the numerically obtained convergence radius from the large order behavior of the coefficients c n agrees with the expected convergence radius.
Moreover, we numerically find a convergent HAM series solution if we are considering a linear operator L such that the HAM does not yield a series solution around the noninteracting limit G = 0. With the simple choice we obtain the result shown in Fig. 8(a). The linear operator is the ladder approximation to the functional [G] and for a field theory  sums up all RPA-ladder diagrams. Therefore the HAM series solution is an expansion with respect to the interacting ladder approximation. With the linear operator L ladder the HAM series solution gives a globally convergent solution. Therefore the series also converges for large values of G. This regime corresponds to the double well potential, i.e., k < 0. The convergence is demonstrated in Figs. 8(b) and 8(c). The HAM series solution is converging for arbitrary G without tuning the convergence control parameter c 0 .

APPENDIX B: FUNCTIONAL CLOSURE FOR FIELD THEORY
In this appendix, we derive a closed system of functional integro-differential equations for general interacting manybody models with arbitrary two-body interaction. These equations define the many-body theory under consideration in a functional integro-differential formulation.

Correlation functions and functional derivatives
To find a closed set of equations for correlation functions we write high-order correlation functions as functional derivatives of low order correlation functions. In this section, we show how we define the functional derivatives in a collective index notation.
In Sec. II, the algebraic 0D quantum field theory (FT) was considered, where Z[J ] is the generating functional of the n-point correlation functions We extend the definition of the generating functional Z[J ] with a one-point source term J to a two-point source term. The inverse noninteracting two-point correlation function k can be thought of as a two-point source term. Therefore the partition function Z = Z[J = 0] is considered to be a generating functional with respect to k, Z = Z[k]. The two-point correlation function G = G (2) can be obtained by differentiation with respect to the two-point source term k, In order to practically use (B3), Z[k] has to be computed for arbitrary k ∈ R and afterwards the derivative has to be evaluated at k = k phys . We first discuss the extension to a FT for the prototypical Z 2 symmetric φ 4 model on a lattice in arbitrary dimensions D, In this case, G −1 0 ∈ R L×L with L the number of lattice sites and G −1 where is the discretized Laplace operator in D dimensions. The measure of the functional integral is denoted by d(φ). The n-point correlation function carries now additional lattice indices and is defined as In analogy with the 0D case, we consider and write the two-point correlation function as Therefore Z[G −1 0 ] has to be computed for arbitrary G −1 0 ∈ R L×L and only after the derivative has been taken the expression is evaluated in the subspace of physical G −1 0 = G −1 0;i,j ;phys . In particularly, although G −1 0;i,j ;phys lies inside the subspace of translational invariant G −1 0 the calculation of Z[G −1 0 ] must not be restricted to this physical subspace. Translational invariance can only be restored after the functional derivative of FIG. 9. Feynman diagrams in momentum space contributing to different correlation functions for φ 4 theory with action (B4). The lines represent noninteracting two-point correlation functions G 0 whereas the interaction vertex is represented by the dot. Due to the local interaction, momentum is conserved at each vertex. (a) The second-order contribution to the self-energy. Due to the translational symmetry of the φ 4 model incoming and outgoing momenta have to be equal and G 0 depends only on a single momentum variable. (b) A second-order contribution to the four-point vertex function depending on the four external momentum variables k i . Due to translational symmetry, δ(k 1 + k 2 + k 3 + k 4 ), the four-point vertex function effectively depends on three momentum variables. (c) The second-order diagram contributing to the self-energy without assuming translational symmetry. Therefore incoming and outgoing momenta do not have to be equal and the two-point correlation function lines depend on two momentum variables.
The following argument [40] illustrates that this is indeed indispensable. The four-point correlation function Model (B4) is translational invariant and therefore G (4) depends on three relative distances, e.g., If translational invariance is already restored before taking the functional derivative in (B7), i.e., taking the functional derivative only in the restricted subspace of translational invariant G −1 0 , the right-hand side depends only on two relative distances i 1 − i 2 , i 3 − i 4 . In this case, the functional derivative does not provide the full information about all values of the four-point correlation function. This shows that it is crucial to consider G −1 0 ∈ R L×L and not perform the calculation in a restricted subspace. We can also consider an argument based on diagrammatic reasonings which is closer to the actual computational techniques introduced in this chapter. The above result can be obtained by taking into account the following considerations. Formally, a FT can also be defined through its series expansion of correlation functions in terms of Feynman diagrams. The functional derivative with respect to G 0 (related to the one with respect to G −1 0 through the chain rule of functional differentiation) can be considered for each diagram. The action of this derivative on a single Feynman diagram is diagrammatically represented by removing one G 0 line in all possible ways. The result of cutting a single G 0 line is a new diagram with two additional external legs, which can be a Feynman diagram in the series expansion of a higher order correlation function. The precise relation will be derived in the following chapter. We consider this process for a specific example where translational invariance is already restored before taking the functional derivative. The functional derivative of the diagram in Fig. 9(a) in momentum space will give the following expression: Due to translational invariance G 0 is diagonal in momentum space and depends only on a single momentum variable. The resulting diagram obtained from this functional differentiation is depicted in Fig. 9(b). It contributes to the Feynman diagrammatic expansion of the four-point vertex function. However, as we already assumed translational invariance before taking the derivative we can obtain information only in the subspace The contribution for arbitrary external momentum variables is given by the functional derivative of the diagram in Fig. 9(c). Translational invariance is not taken into account yet and therefore G 0 depends on two momentum variables. In this case, differentiating the diagram gives which is exactly the analytic expression for Fig. 9 This result can be generalized to more complicated symmetries such as interacting fermions with spin σ on a lattice with a translational invariant hopping matrix h and local onsite Hubbard type interaction. The free part S 0 of the action S = S 0 + S int can than be written as and the interacting part S int describes a local interaction, which does not break U(1) symmetry nor SU(2) spin symmetry nor does it violate translational invariance in space and imaginary time.
We focus in the following on the U(1) symmetry leading to particle number conservation but note that exactly the same argument can be used for SU(2) spin symmetry.
A second-order diagram contributing to the self-energy for this model is given in Fig. 10(a). We only consider the labels of the diagram due to the U(1) symmetry where incoming lines denote ψ and outgoing lines denoteψ. Due to the U(1) symmetry particle number is conserved and the number of ingoing lines equals the number of outgoing lines. Taking the derivative of this diagram with respect to Gψ ψ 0 we obtain the diagram in Fig. 10(b), which contributes to the four-point correlation function. The diagram in Fig. 10(a) is not the only diagram which upon cutting a single propagator line leads to the diagram in Fig. 10(b). If we allow for U(1) symmetry breaking we can also consider the diagram in Fig. 10(c), which contributes to the off-diagonal part of the self-energy. The particle number is not conserved for this diagram as there are two incoming lines but no outgoing lines. The diagram includes a off-diagonal G 0 line, Gψψ 0 , which is indicated by two arrows on that line pointing in different directions. Cutting this G 0 line also leads to the diagram in Fig. 10(b), which preserves particle number conservation. This example shows, in analogy with translational invariance, that in an intermediate stage of the calculation U(1) symmetry can be broken and only after the derivative has been taken all symmetries of the action have to be respected.
Concluding, we showed that in general G −1 0 has to be defined outside the subspace of physical inverse noninteracting two-point correlation functions. The full space is given by a collective index space such that the free part of the action can be written as The collective index i summarizes all possible field labels of the considered model.
The inverse noninteracting two-point correlation function G −1 0;i,j is therefore defined in this collective index space with properties discussed in the following.

Dyson-Schwinger equations
We start with the most general form of the action of a many-body system with arbitrary two-body interactions: The fields φ i are either complex numbers if they represent bosonic degrees of freedom or anticommuting Grassmann numbers for fermionic degrees of freedom. They are labeled by a single collective index i, which summarizes all possible labels for the fields as introduced in Sec. 1. The action can be formally defined on a lattice or in continuous space, where x is taken to be a discrete variable or a continuous variable, respectively. The inverse of the noninteracting Greens function G −1 0 should be thought of as a symmetric/antisymmetric matrix (for bosons/fermions) in the space of the collective index i and the matrix elements of the two-body interaction are given by a fully symmetric/antisymmetric tensor V . Summation over repeated indices is implicitly assumed. The Dyson-Schwinger equation for model (B12) can be derived by introducing the generating functional of the n-point correlation functions: For bosons the source fields J i are complex numbers, while for fermions they are Grassmann numbers anticommuting with itself and with the fields φ i . The measure of the functional integral is denoted by d(φ). All n-point correlation functions can be generated by functional derivatives of the generating functional with respect to the source fields, The Dyson-Schwinger equation is derived by considering a linear shift of a single field variable in the functional integral for the generating functional. For fermions, has to be a Grassmann number, while for bosons, is a complex number. The elements in the functional integral for the generating functional transform under this shift as Therefore the generating functional is given by Here and in the following, the upper sign in ± or ∓ is for the bosonic case where the lower sign is for the fermionic case. The differential form of the Dyson-Schwinger equation is obtained by equating the above expression in powers of , This is a functional integro-differential equation The relations between the expansion coefficients, i.e., the n-point correlation functions, can be obtained by successive differentiation of (B17) with respect to the source fields J i and setting the sources to J i = 0 afterwards. This yields an infinite tower of integral equations for the n-point correlation functions. The equation obtained by differentiating (B17) once is given by This equation relates the two-point correlation function with the four-point correlation function. The four-point correlation function can be split into a disconnected part and a connected part which on the other hand can be factorized into a contribution coming from the two-point correlation function and a genuine four-point contribution, the four-point vertex function. Therefore we rewrite (B19) such that it relates the two-point correlation function G with the four-point vertex function . In principle, this can be done by introducing further generating functionals for one-particle irreducible correlation functions but as we only need the connection between correlation functions on the four-point level we will directly introduce the relations between them. The connected four-point correlation function G (4) c is given by Here and in the following, we use the short hand notation i n = n for the collective index i n . Plugging this relations into (B19) and solving for the inverse of the two-point correlation function, we obtain The self-energy 1,2 has been introduced as a short hand notation for the contributions coming from the interaction part. The above equation is a single equation for two unknown correlation functions and therefore constitutes an underdetermined set of equations. In the following section, we derive a closed set of equations.

Functional closure of Dyson-Schwinger equations
Based on the Dyson-Schwinger equation (B22), which relates the inverse two-point correlation function with the four-point vertex function, we derive a closed set of functional integro-differential equations. The solution of this set of differential equations gives direct excess to the two-point correlation function G and to the four-point vertex function .
We start the derivation with the identity We have used that This functional derivative identity together with the relation between the correlation functions (B20) yields The four-point vertex function can also be written as a functional derivative with respect to G −1 0 by plugging the above identity (B25) into the definition of the four-point vertex function (B21). The final result is where we have used the functional chain rule . Therefore, in the following, we will use the self-energy 1,2 obtained from the Dyson-Schwinger equation (B22) and plug this expression into the identity (B26). We obtain the following expression after a long but straightforward calculation: The last term can be simplified by noting that second order functional derivatives commute. The following commutator identity can be derived using this property, This relations transforms (B30) into a definition of the four-point vertex function as a universal functional [G] in the language of a functional integro-differential equation. The complete set of closed Dyson-Schwinger equations is given by The equations are diagrammatically depicted in Fig. 2.

APPENDIX C: NUMERICAL IMPLEMENTATION
In this section, we discuss the numerical implementation for the solution of the functional integro-differential equation in (B33) for the Z 2 symmetric φ 4 model in 2D, cf. (B4). Thus, in the following, we consider model (B12) for the case where the collective index i just consists of the labels of the discrete lattice sites i = (x i ) of the 2D square lattice. We will see in the following that adding additional field labels to the collective index i in principle does not change the algorithm but complicates the already complex bookkeeping of indices in the algorithm.
For model (12) G −1 0;i,j = − i,j + m 2 δ i,j , where is the discretized Laplace operator in D dimensions and the fully symmetric tensor is given by the on-site interaction V i,j,k,l = λ δ i,j δ i,k δ i,l . We consider the construction of the universal functional [G] through the solution of the vertex equation in (B33) in momentum space. In the following, the momentum variables are denoted by p i = i and the linear combination of momentum variables by, e.g., p i + p j = i + j . In this representation, the integro-differential equation in (B33) can be written as

Homotopy analysis method
In this section, the semianalytic HAM is used to solve the functional integro-differential equation (C1). As discussed in Sec. III the starting point of the HAM is the construction of the homotopy The equation for the m-th-order deformation equation is the starting point for the tree expansion developed in Ref. [17] for a truncated functional [G]. In the next section, we will summarize the main ideas of the tree expansion and extend these ideas to account for the functional derivatives in (C4).

Tree expansion
The mth-order deformation equation (C4) gives the mth term in the series solution (C3) of the HAM. In order to calculate the mth deformation, all previous deformations and their functional derivatives have to be known. Neglecting the functional dependence on the unknown G for a moment, already the full storage of a single deformation is a formidable task as the deformations themselves depend on four external indices, which are D-dimensional vectors, i.e., u ,m (p) is a rank-4D tensor. Including the functional derivatives seems to make a numerical calculation of (C4) an impossible task. The tree expansion developed in Ref. [17] solves this apparently impossible task by a stochastic interpretation. The basic idea of the tree expansion is that in order to calculate the mth-order deformation with (C4), (C4) is used again on the right-hand side to get rid off the explicit dependence of all deformations and functional derivatives with k < m. This procedure is recursively repeated until all deformations u ,m with m > 0 are eliminated on the right-hand side. The result is a complicated expression, the tree expansion, which only depends explicitly on the deformation u ,0 . We will develop a diagrammatic language which captures all possible terms in the tree expansion and introduce a Monte Carlo algorithm to stochastically sum all terms. Each term in the tree expansion will be represented in the diagrammatic language of rooted trees which will be introduced in the following. Therefore the tree expansion for u ,m is given by the sum of all possible rooted tree diagrams, which can be constructed with respect to a certain fixed set of rules.
We start the discussion with the rooted tree diagram in Fig. 11(a). It corresponds to a single term in the tree expansion. For the moment, we are neglecting the labeling of the rooted tree and consider only the overall structure of the tree and its corresponding analytic expression. After the general structure is introduced, we discuss its labeling, i.e., how momentum variables and additional labels, e.g., external leg permutations, are taken into account. The rooted tree and the structure of the corresponding analytic expression is obtained from the following procedure.
The uppermost circle, the root of the rooted tree, corresponds to the left-hand side of (C4). In the case of Fig. 11(a), u ,3 should be expanded. The branch, indicated by a straight line, growing from the root leading to two leafs, represented by the circles, diagrammatically depicts that u ,3 is expanded in the first stage with respect to the third line in (C4), i.e., . 11. Examples of rooted tree diagrams, which are constructed by expanding the HAM in the tree expansion. The different diagrammatic elements are explained in the main text and listed in Thus, the branch corresponds to the integral kernel K (2) c and the two leafs corresponds to the product of deformations u ,k u ,m−1−k . A single term in the sum over deformations is chosen. In Fig. 11(a), k = 1 has been picked leading to u ,1 u ,1 . In the second stage, the above procedure is repeated by considering each leaf as a new individual root. In case of Fig. 11(a), the left leaf has been expanded with respect to the second line in (C4), i.e., and the right leaf corresponds to the first line, i.e., For that case, there is no contribution from an integral kernel. This is indicated by drawing the leaf as a star instead of a circle. Gathering all contributions, the structure of the analytic expression corresponding to the rooted tree in the example of Fig. 11(a) is This is a single term in the tree expansion and illustrates that the rooted trees depend explicitly only on the kernel functions and on the initial guess u ,0 , which is the starting point of the series solution (C3). The discussion also shows that in (C4) there are five different possibilities to expand u ,m , which lead to the five different branch types. All branch types with their respective leafs are shown in Fig. 12. In the diagrammatic expressions the branch type is determined by the leafs on the lower end of the branch. The branch type in Figs. 12(d) and 12(e) correspond to the expansion with respect to the functional derivative terms. The leafs with functional derivatives are drawn by boxes and the additional dangling line is introduced. It indicates the action of the functional derivative on a consecutive element in the rooted tree. On the basis of Fig. 11(c), we discuss the structure of the analytic expression corresponding to a rooted tree with functional derivatives. In the first stage of the tree expansion, the root u ,4 is expanded with respect to the branch type 12. The possible branch types of the rooted trees. Each branch type corresponds to a single line on the right-hand side of (C4). The straight line, the branch itself, represents the corresponding kernel function K (i ) c , i ∈ {1, 2} or K (i ) , i ∈ {3, 4} for the respective branch type. The star, circle, and box indicate the leaf grown from the branch. They represent the contribution from the deformation u ,m for the specific branch type. The star contributes no integration variable to the tree, cf. first line in (C4). The box corresponds to a leaf at which a functional derivative is introduced. The branch type can be read off from the combination of leafs grown on the lower end of the branch. The leaf types (c) and (e) consists of two individual single leafs because each can be individually expanded further in the tree expansion. The dangling lines on the leafs of (d) and (e) are representing the not yet completed action of the functional derivative on a consecutive branch. Fig. 12(d), which includes a single leaf with a functional derivative. The structure of the analytic expression in this expansion stage is given by In the next stage, u ,3 is expanded with respect to the branch type Fig. 12(e). This gives the analytic expression Due to the product rule of functional derivatives there are three possibilities for the action of the outer functional derivative. The first possibility is that the functional derivative is acting on the kernel function of the second expansion stage, i.e., the branch directly after the deformation u ,3 is differentiated. The second possibility is that the outer derivative is acting on the leaf with the functional derivative and therefore a secondorder functional derivative is produced. The third possibility is that it acts on the deformation u ,1 . From Fig. 11(c), we see that the dangling line is completed to an arc line, which indicates that the derivative introduced in the first stage of the tree expansion is differentiating a branch grown after the leaf without a functional derivative. Therefore the third possibility is represented by the rooted tree in Fig. 11(c) and the analytic expression for the rooted tree in this expansion stage is In the next two expansion stages, the deformations u ,1 are further expanded. In both cases with respect to the branch type Fig. 12(b). This gives the analytic expression Here and in the following, we consider only initial guesses u ,0 , which have no explicit G dependence, i.e., Taking this into account the structure of the analytic expression corresponding to the fully grown root tree in Fig. 11(c) is given by From this discussion, we see that functional derivatives can only act on the branches of a fully grown rooted tree. Nevertheless, in intermediate stages of the tree expansion the functional derivatives can act on the deformations with u ,m , m > 0. These have to be expanded further in later stages of the tree expansion and therefore the analytic expression corresponding to a fully grown rooted tree involves only u ,0 , kernel functions and functional derivatives of the kernel functions. Moreover, as we see from the above discussion the derivatives can only act on branches, which are grown from a leaf consecutive to a functional derivative leaf. Examples of valid rooted trees with functional derivatives are shown in Figs. 11(b) and 11(c). In these fully grown rooted trees, the dangling lines, which indicate the action of the functional derivative form closed arc lines connecting a functional derivative leaf to a branch. This denotes which branch is differentiated by the functional derivative. We will see later that each branch can only be differentiated a finite amount of times and therefore all possible functional derivatives on all possible branches can be tabulated, cf. Appendix E.
We have introduced the general structure of the rooted trees in the framework of the tree expansion on the basis of the examples in Fig. 11. In the following, we discuss the labeling and the specific analytic expressions of the diagrammatic elements in the rooted tree diagrams. The universal functional [G] is a functional with respect to an arbitrary G ∈ R L×L , which is not necessary defined inside the physical, translational invariant subspace. Nevertheless, we finally want to evaluate [G] in the subspace of physical G. We showed that [G] can be constructed by drawing all possible rooted tree diagrams in the tree expansion of the HAM. Therefore, in order to evaluate [G] at a physical G, the diagrammatic elements of the rooted tree diagrams, i.e., the kernel functions and their derivatives have to be build up from the physical, translational invariant two-point correlation functions. In the following, we show how this leads to the correct labeling of the rooted tree diagrams.
Consider the rooted tree diagram in the tree expansion for u ,3 depicted in Fig. 11(d). We fix the external momentum of the root, u ,3 (p), to be p = (p 1 , p 2 , p 3 , −p 1 − p 2 − p 3 ), i.e., the rooted tree diagram is evaluated inside the translational invariant subspace. The first stage in the tree expansion for the rooted tree in Fig. 11(d) is graphically illustrated in Fig. 13(a). In the example, the expansion is performed in the t-permutation channel, cf. (C4). The analytic expression corresponding to this expansion stage is diagrammatically illustrated on the right-hand side of the arrow in Fig. 13(a) and can be read off as  13. Two consecutive stages of expansions in the tree expansion leading to the full rooted tree in Fig. 11(d). The two-point correlation functions contribution to the kernel are not yet evaluated in the physical, translational invariant subspace. Therefore G still depends on two momentum variables, which are indicated by two arrows on top of the G lines. The new external momentum carried by the leafs is determined by the kernel functions in (C4) and the specific random permutation of external legs picked for this expansion stage. Diagrammatically, the contributing kernel functions together with the respective deformations are depicted on the right-hand side of the arrow. (a) The deformation u ,3 is expanded with respect to the second line in (C4) leading to a new leaf, u ,2 . The branch corresponds to the kernel function depending on two two-point correlation functions, which lead to a evaluation of the lower-order deformation at a new external momentum. In this example, the t channel is considered. (b) In the next expansion stage, the leaf u ,2 is expanded further with respect to the third line in (C4) in the u permutation. In this case, the kernel function depends on four two-point correlation functions and the leafs corresponds to two lower order deformations each evaluated at a new external momentum.
Here and in the following, summation over repeated indices is assumed implicitly. Obviously, the kernel function can be evaluated in the physical subspace of translational invariant G, G phys,1,2 = G 1 δ 1,−2 . Consequently, the new external momentum for u ,2 (p ) is p = (p 1 , p 2 , p 3 , p 4 ) = (p 1 + p 3 − p 5 , p 5 , p 2 , −p 1 − p 2 − p 3 ) and therefore momentum conservation is fulfilled and the procedure of the tree expansion can be recursively repeated without leaving the physical subspace.
The same holds for the expansion with respect to the branch type Fig. 12(c). For the example of Fig. 11(d), the second stage of the tree expansion is illustrated in Fig. 13(b). In this case, assuming that the u permutation has been chosen, × u ,0 (−6, 4 , −9, −8) The external momentum variables for the leafs u ,0 (p ) and u ,1 (p ) are after the evaluation of G in the physical subspace, p = (p 5 , p 4 , −p 6 , −p 4 − p 5 + p 6 ) and p = (p 1 + p 4 − p 6 , p 6 , p 3 , p 2 ). Therefore the evaluation of the deformations u ,0 , u ,1 is again in the physical subspace. We are choosing the convention that the leaf on the left-hand side of the grown branch of type Fig. 12(c) corresponds to the deformation depending on two external variables. In the example, the deformation u ,1 , depending on two external variables, is represented by the left leaf, whereas u ,0 , depending on one external variables, is represented by the right leaf, cf. Fig. 13(b).
In conclusion, we see from this example that there are no further complications in evaluating the rooted tree diagram at G = G phys if the diagram contains no elements with functional derivatives. In this case, the only information needed in each expansion stage is (1) the external momentum variables of the leaf which is expanded, (2) the branch type into which the leaf is expanded, and (3) the internal momentum variables over which is integrated. This provides enough information to compute in each expansion stage the new external momentum variables for the newly grown leafs and therefore the tree expansion can be carried out recursively.
Rooted tree diagrams with elements of type Figs. 12(d) and 12(e) need some further consideration. As an example consider Fig. 11(b), which contains a functional derivative.
As in the previous example, the rooted tree is evaluated in the physical subspace by fixing the external momentum of the root to be p = (p 1 , p 2 , p 3 , −p 1 − p 2 − p 3 ) and evaluate the kernel function corresponding to the branch of the corresponding branch type, Fig. 12(d), at translational invariant G.This is illustrated in Fig. 14(a). 5,...,9 K (3) (1, 5, . . . , 9) δu ,1 (−7, 2, 3, 4) δG 8,9 = G 1−5−6,7 G 5,8 G 6,9 δu ,1 (−7, 2, 3, −1−2−3) δG 8,9 The additional arrowheads on the two-point correlation function lines in Fig. 14(a) are indicating the connection to the functional derivative. Note that only the G lines contributing to the kernel function are evaluated in the physical subspace. The functional derivative still depends on two momentum variables. We will show later how this has to be understood by considering the differentiation of a kernel function. The G lines without additional arrowheads are connected to the deformation such that the momentum carried by this line is con- correct momentum conservation at the functional derivative terms can be found from the following argument.
The functional derivative of the four-point vertex function with respect to G gives a term which account for correlation functions on the six-point level. We write this as In the subspace of physical correlation functions, O (6) depends only on the relative distances between the lattice sites i 1 , . . . , i 6 . Therefore, in order to evaluate rooted trees with functional derivatives in the translational invariant subspace, the correct momentum conservation at the functional derivative is δ(p 1 + p 2 + p 3 + p 4 − p 5 − p 6 ). We will also need the more general result for the momentum conservation with multiple functional derivatives The functional derivative terms explicitly break momentum conservation on the four-point level but on a higher correlation function level it must always be satisfied with respect to the extended momentum conservation rule (C19). The additional arrowheads on the G lines, cf. Fig. 14(a), can also be viewed as indicating that these momenta are contributing with a minus sign in the extended momentum conservation rule. This consideration shows that (C17) is a valid projection into the physical subspace and the tree expansion can now be carried on if the momentum conservation rule (C19) is taken into account.
Diagrammatically the excess momenta carried by the functional derivatives are depicted as the previously introduced additional arcs which originate from leafs with functional derivatives and close on a branch, cf. Fig. 11(c). We show in the following that once an arc closes on a branch, i.e., after the functional derivative acted on a kernel function, momentum conservation will be automatically restored on the four-point level.
We consider the example in Fig. 11(c). The analytic expression obtained in the first stage of the expansion of u ,4 at fixed external momentum variables compatible with translational invariance p = (p 1 , p 2 , p 3 , −p 1 − p 2 − p 3 ) is already given in (C17) and illustrated in Fig. 14(a). In the second stage in Fig. 11(c), the leaf u ,3 is further expanded with respect to the branch of type Fig. 12(e). This expansion process is illustrated in Fig. 14(b). The external momentum variables for u ,3 (p ) are given by p = (p 1 , p 2 , p 3 , p 4 ) = (1−5−6, 2, 3, −1−2− 3). The functional derivative with respect to the momentum variable (p 5 , p 6 ) = (−p 5 , −p 6 ) does not act directly on the branch grown after the leaf with the functional derivative u ,3 , but it acts on a branch grown from a consecutive leaf. If the functional derivative would act on the branch grown from the functional derivative leaf, represented by the box, it would yield a second-order functional derivative. However, from Fig. 11(c), we see that it acts on the branch grown from the leaf without functional derivative. Therefore we obtain an expression with two first order functional derivatives, which is diagrammatically depicted in Fig. 14(b). The functional derivative with respect to G(p 10 , p 12 ) originates from u ,1 , i.e., from the current branch type, whereas the derivative with respect to G(p 5 , p 6 ) originates from u , 3 . We consider the extended momentum conservation rule (C19). The analytic expression corresponding to this expansion step is given by G=G phys where we have established the generalized momentum conservation rule for p 12 = p 1 + p 2 + p 3 + p 4 − p 5 − p 6 − p 10 and renamed in the last equality the integration variables in order to distinguish them from the integration variables of the previous expansion stage. This example shows that even though the functional derivative with respect to G(5 , 6 ) only acts on a branch consecutive to the leaf u ,1 which for itself is not a leaf with a functional derivative the information about the excess momentum running over this leaf has to be stored. Therefore an additional line through the rooted tree is drawn, which indicates that the excess momentum has to be accounted for on the intermediate leafs by Finally, in the next two expansion stages, the functional derivatives act on the branches. The functional derivative acting on a branch can be graphically depicted as just cutting out one of the two-point correlation function lines. This is shown in Fig. 15 for a particular example. It corresponds to the last expansion stage in Fig. 11(c) where the right leaf u ,1 is expanded. It is assumed that the t-permutation channel is picked. The external momentum of the leaf u ,1 (p ) is p = (7, 9, 1 + 2 + 3 + 4 − 7 − 8 − 9, 8). The analytic expression for the differentiated branch is given by δ δG (5 , 6 ) G 5,6 G 1+3−5,7 u ,0 (−6, −7, 2 , 4 ) The new external variable for the leaf grown from the differentiated branch u ,0 (p ) is p = (−6 , 1 +3 −10, 2 , 4 ). Inserting the expression for the primed variables we find that u ,0 is evaluated at p = (6,−6−8 −9, 9, 8). Therefore, after the functional derivative has acted on the branch, the excess momentum is annihilated and momentum conservation at the four-point level is restored. The discussion of the examples in Fig. 11 provides enough information to write down a set of rules to construct and label all possible rooted trees. A second set of rules assigns to each diagrammatic element in a rooted tree an analytic weight. This makes it possible to design a Monte Carlo algorithm in the space of rooted trees which directly evaluates [G] in the subspace of translational invariant G. In the following, we first write down the set of rules to construct and label a random rooted tree. The complete list of diagrammatic elements with their respective weights are tabulated in Appendixes D and E.
We would also like to note that the above examples only discuss the case where the physical subspace is defined by translational invariance. With the same line of arguments it is possible to obtain the construction and evaluation of [G] in a physical subspace for general symmetries defined by the collective index space, cf. (B12) and (B33).
The HAM gives a series solution for [G], which sums over all deformations u ,m . Thus rooted trees of arbitrary expansion order m have to be considered. The rules to generate a random rooted tree in the tree expansion of an arbitrary u ,m with external momentum variables p are. (1) Grow a random branch from the root, i.e., select the expansion of the root into one of the branch types in Fig. 12. For this choice, we assign the labels (k,b,d). The label k stands for the number of leafs grown from the branch, k=1, 2. The label b = bare, bold for whether the analytic expression for the branch involves contributions from two-point correlation functions and d = True, False whether one of the leafs is a functional derivative. Thus the respective labeling for the branch types in Fig. 12 are (a) (1,bare,False), (b) (1,bold,False), (c) (2,bold,False), (d) (1,bold,True), and (e) (2,bold,True). The probability to select the chosen (k,b,d) is stored in p a priori in order to obtain the a priori probability for the fully grown rooted tree which is needed for the Monte Carlo process. (2) According to the randomly picked branch type add the respective number of leafs with their given type to the tree. If a branch type with k=2 was picked one term in the sum m−1 i=0 has to be chosen randomly, cf. (C4), and the probability for that choice is multiplied to p a priori . If a branch type with b=bold and d = False was pick one channel of external leg permutations in the sum c has to be chosen randomly, cf. (C4), and the probability for that choice is multiplied to p a priori . Depending on the branch type new momentum variables have to be seeded randomly. These are D-dimensional random vectors in the first Brillouin zone. The number of new momentum variables for the possible types are: (a) 0, (b) 1, (c) 2, (d) 2, and (e) 3. The probability for the choice of momentum variables is multiplied to p a priori .
For the branch types (c) and (e), there are two newly grown leafs. For (c), there are two new momentum variables. We assign the left leaf to represent the deformation evaluated at two external momentum variables and the right leaf depending on only one. We choose the convention that the first of the two momentum variables is assigned to the left leaf whereas the right leaf depends on both variables.
For (e), there are three momentum variables. The leaf with the functional derivative depends only on two of the three randomly seeded momentum variables. We choose the convention that the first two of these variables are assigned to the functional derivative leaf and the right leaf depends on all three variables.
(3) For each of the new leafs, u ,m , if m > 0 consider this leaf as a new root and go to 2. If m = 0, a new branch can not be grown.
(4) Start a search through the fully grown rooted tree. For each functional derivative encountered, choose randomly a branch consecutive to the functional derivative leaf. This branch will be differentiated by the functional derivative and therefore the additional bookkeeping line is established. The probability for each choice is multiplied to p a priori .
(5) Start a search through the fully grown rooted tree. If a differentiated branch is encountered pick randomly a term in the respective list of differentiated kernel functions, cf. Appendix E. Correspondingly, reduce the number of momentum variables of the differentiated branch. The probability for each choice is multiplied to p a priori .
(6) Start a search through the fully grown rooted tree. For each expansion stage, determine the new external momentum variables of the leafs grown in this expansion stage. This is done with respect to the extended momentum conservation rule (C19), cf. also Appendixes D and E.
(7) Start a search through the fully grown rooted tree. As the momentum variables are now correctly assigned the complete weight of the rooted tree can be determined by multiplying the contributions from the kernel functions and from u ,0 at each expansion into a variable w tree . Therefore, after having applied the above rules, a randomly grown rooted tree is obtained with a corresponding weight w tree . Together with the a priori generation probability p a priori this is sufficient to perform a direct sampling of the tree expansion by the rules of detailed balance. Therefore we obtain a stochastic summation of the HAM series solution for the functional integro-differential equation defining [G].
We find it convenient to extend the direct sampling of diagram topologies with a Markov chain sampling of the integration variables. For that, we randomly pick a single integration variable of a rooted tree and suggest to update this integration variable by adding a random shift in that momentum variable. After the integration variable has been changed, we start from step 2 in the above set of rules to generate a new weight for the tree with the new integration variable. Together with detailed balance, this is used to perform a Markov chain sampling of the integration variables. In order to achieve good acceptance ratios for the direct sampling of rooted tree topologies, we restrict the randomly seeded integration variables in each stage of the tree expansion to belong to a small region around zero momentum. Therefore we rely on the Markov chain sampling of integration variables to obtain the integration over the complete first Brillouin zone. This leads to sampling problems for higher deformation orders because the sampling can get stuck in a random rooted tree topology. A Markov chain Monte Carlo algorithm has to be used, which also samples diagram topologies and not only integration variables. Therefore a better starting point to include the functional derivatives, compared to randomly suggesting a rooted tree, is the algorithm developed in Ref. [17].

APPENDIX D: KERNEL FUNCTIONS
In this section, we list the analytic expressions for the different branch types and the new external momentum variables for the corresponding leafs grown from these branches. The analytic expressions give directly the projection into the physical subspace of translational invariant G. We denote the external momentum variables at the leaf from which the branch is grown as p = (p 1 , p 2 , p 3 , p 4 ) and assume that this leaf corresponds to a deformation u ,m . We again use the shorthand notation p i = i and p i + p j = i + j and implicitly assume summation over repeated indices. Figure 12 lists all possible branch types where each branch type accounts for a single line in (C4). Furthermore, we consider functional derivatives in intermediate stages which directly act on the leafs and which have to be accounted for by the momentum conservation rule (C19). These additional functional derivatives are highlighted by the momentum labels l i . Diagrammatically the branch types are depicted in Fig. 16. The analytic expressions for the different branch types, including the intermediate functional derivatives, are given by (a) the analytic expression for the branch type in Fig. 12 .