Expansion of the effective action around non-Gaussian theories

This paper derives the Feynman rules for the diagrammatic perturbation expansion of the effective action around an arbitrary solvable problem. The perturbation expansion around a Gaussian theory is well known and composed of one-line irreducible diagrams only. For the expansions around an arbitrary, non-Gaussian problem, we show that a more general class of irreducible diagrams remains in addition to a second set of diagrams that has no analogue in the Gaussian case. The effective action is central to field theory, in particular to the study of phase transitions, symmetry breaking, effective equations of motion, and renormalization. We exemplify the method on the Ising model, where the effective action amounts to the Gibbs free energy, recovering the Thouless-Anderson-Palmer mean-field theory in a fully diagrammatic derivation. Higher order corrections follow with only minimal effort compared to existing techniques. Our results show further that the Plefka expansion and the high-temperature expansion are special cases of the general formalism presented here.


Introduction
Many-particle systems are of interest in various fields of physics. Field theory offers a versatile unique language to describe such systems and powerful methods to treat diverse problems arising in classical statistical mechanics, quantum mechanics, quantum statistics, quantum field theory, and stochastic dynamical systems [1][2][3][4][5]. Diagrammatic techniques in particular are convenient and efficient to organize practical calculations that arise in the context of systematic perturbation expansions and fluctuation expansions around a solvable problem. But most techniques rely on the solvable part being Gaussian: The basic connecting elements in Feynman diagrams, or Mayer graphs [6], are lines. The purpose of this paper is to extend the diagrammatic computation of a central quantity, the effective action Γ, beyond this Gaussian case.
The effective action, vertex generating functional, or Gibbs free energy, is the first Legendre transform of the generating functional of connected Green's functions or cumulants. Originally developed in statistical mechanics to study non-ideal gases in the 1930s, its modern formulation as a Legendre transform was introduced in the context of quantum statistics [7], but has soon been applied in quantum field theory [8] (similar concepts had been developed in parallel in quantum field theory [9]).
The effective action is central to the study of phase transitions, because its stationary points determine the values of the order parameters that characterize the states in which matter may exist; it therefore reduces the problem of a phase transition to finding bifurcation points in a variational problem [2,3,5]. This formulation yields self-consistency equations for the mean value of a field that incorporate fluctuation corrections. At lowest order, at the tree level of diagrams, one obtains the mean-field theory of the Curie-Weiss type. Higher Legendre transforms generalize this concept to finding self-consistent solutions for renormalized higher order Green's functions [7] [5, chapter 6] [10,11]; the second transform leading, at lowest order, to a self-consistency equation for the propagator, the Hartree-Fock approximation [12][5, i.p. chapter 6 for a review].
The second reason for the importance of the effective action is that it compactly encodes all dynamical or statistical properties of a system: Connected Green's functions decompose into tree diagrams, whose vertices are formed by derivatives of Γ. Conversely, relatively fewer diagrams contribute to Γ, making its computation favorable: Its perturbative expansion around a Gaussian solvable theory, or around Gaussian fluctuations with regard to a background field, requires only one-line (particle) irreducible graphs (1PI) [2,3,7,13].
Third, treating critical points by renormalization group methods in classical [14,15] and quantum systems [16] leads to the study of flow equations for the effective action; most explicitly shown in its functional formulation [10,17,18].
Perturbative calculations are essential for most applied problems. But practically all techniques in field theory are based on expansions around the Gaussian case. The topic of this article is a systematic, diagrammatic, and perturbative expansion of the effective action in the general case, where the solvable problem is not of Gaussian type; in statistics: the unperturbed problem has cumulants of orders three or higher; in field theory: the bare theory has non-vanishing connected Green's functions of orders three or higher; in diagrams: we have bare propagators that tie together three or more field points instead of or besides the usual propagator lines.
It may be obvious that the linked cluster theorem in this case remains intact without change, shown below for completeness. Also the converse problem, the decomposition of Green's functions into connected Green's functions by Wick's theorem is immediately replaced by its general counterpart, the factorization of moments into all product sets of cumulants [19]. But diagrammatic rules for the direct expansion of the effective action do not follow as simply. Only for particular non-Gaussian cases, diagrammatic rules have been derived earlier: for the classical non-ideal gas [6,20] and for the Ising model [21,22] (compare [5, end of chap. 6.3.1.]). For the classical gas, one finds that only so called star-graphs contribute [5, 6, 21, chap. 6.3.2.], which is known as Mayer's second theorem. Star-graphs are diagrams that do not fall apart by removing one cumulant. For the Ising model one has to introduce a special notation for the sums, the so called N -form, which only includes contributions from star graphs [22]. At the end of Section 2.6, we will give more details on this result and show its relation to our work.
We first derive a recursive diagrammatic algorithm to compute the effective action perturbatively around a non-Gaussian solvable theory. Then we show that the correction terms produced by the iteration satisfy a set of Feynman rules: Only two topologically defined classes of diagrams remain. The first being irreducible diagrams in a more general sense than in the known Gaussian case. The second being diagrams that contain vertex functions of the solvable theory. The Feynman rules for the former directly follow from the linked cluster theorem. The rules for the latter graphs follow from a set of problem-independent skeleton diagrams that we present subsequently. This set of rules then allows the computation of corrections at arbitrary order in perturbation theory in a non-iterative manner. We find, however, that the iterative procedure itself proves highly efficient. Applied to the special case of the Ising model, it affords only a few more calculations than the method that is specific to the problem [21,22]. In this particular case, the iterative procedure is closely related to the known high temperature expansion [23].
The remainder of the paper is organized as follows: In Section 2.1 we introduce the notation and the minimal extension of the diagrammatic language required. Section 2.2 sets up the general perturbative problem and shows why the expansion of the partition function and the generating functional of connected Green's functions straight forwardly extends to the non-Gaussian setting. These preparations are required to derive in Section 2.3 our main result, a recursive algebraic equation that determines Γ. Section 2.4 interprets the algebraic expressions in diagrammatic language, leading to a set of Feynman rules to calculate all graphs contributing to Γ at arbitrary given order in a non-iterative manner. In Section 2.5 we exemplify the algorithm for the simplest non-Gaussian extension of a ϕ 3 theory. In Section 2.6 we demonstrate its application to the Ising model (which turns into the Sherrington-Kirkpatrick spin glass model for Gaussian random couplings [24]) to derive the Thouless-Anderson-Palmer mean-field theory [25] and its higher order corrections [21-23, 26, 27] in a purely diagrammatic manner.

Results
We here chose a notation that should be transparent with regard to the nature of the problem (following to some extent [13]): our results easily transfer to classical statistical mechanics, quantum mechanics, quantum statistics, or quantum field theory. For clarity, we here stick to the language of classical statistical field theory, in particular to bosonic fields. In the following Section 2.1 we set up the language and define elementary quantities.

Definition of a field theory
We assume that the physical system has degrees of freedom x ∈ E. We may think of x as a scalar, a vector (as in the example in Section 2.6) or a (possibly multi-component) field. The domain E must be chosen accordingly. The system is described by the action S S : E → R.
with the property that a particular configuration of x appears with probability p(x) ∝ exp(S(x)). The partition function Z is then given as where we denote as j T x the inner product on the space E. The symbol x stands for the sum over all configurations of x ∈ E, which technically may be a sum, an integral, or a path-integral, depending on the space E. We call j the source field. The properly normalized moment-generating function is denoted as Z(j). Its logarithm is the cumulant generating function also called generating functional of the connected Green's functions in field theory; in statistical mechanics it is related to the Helmholtz free energy F (j) = ln Z(j) by W (j) = F (j) − F (0). We denote the n-th derivative of a function f by f (n) and the n-th cumulant as x n .
To define the effective action Γ(x * ), we eliminate the dependence on the source field j in favor of the mean value x by a Legendre-Fenchel transform We provide the derivation of (4) in Appendix 4.1, because we need an intermediate step as the starting point of the perturbative expansion. To this end we employ the reciprocity property of Legendre transforms which, as derived in the appendix (with (53)), leads to the integral-differential equation [13,Sec. 3.23.6] The latter expression is the starting point for the perturbative expansion of Γ. It is also the natural starting point of a loopwise expansion. The derivative operators that lead to W (1) and Γ (1) are given according to the space E as either partial derivatives or, in the case of fields, functional derivatives.

Perturbative problems
The equation of state (5) provides us with a self-consistency equation for the mean value x * . The general strategy is therefore to obtain an approximation of Γ that includes fluctuation corrections and then use the equation of state to get an approximation for the mean value including these very corrections. We will here derive a perturbative procedure to calculate approximations of Γ and will find the graphical rules for doing so. To solve a problem perturbatively we decompose the action as with a part S 0 that can be solved exactly, that is we know its corresponding cumulant generating function W 0 (j), and the remaining terms collected in ǫV (x). For the special case of S 0 being quadratic in the fields x, this leads to the well known result that corrections to Γ are composed of one-line irreducible graphs only; our algorithm includes but extends this case to non-quadratic S 0 . An example of this situation is depicted in Figure 1, where the solvable theory W 0 already contains cumulants of order one to three. We use the common notation for the interaction or bare vertices, which we will call just "vertices" in the following: the n-th power of the field x is denoted as a vertex with n emerging lines. The connecting elements in the graphs are the cumulants contained in W 0 . These are symbolized by circles, where the number of legs corresponds to the order of the cumulant (following e.g. [1,3]). In ordinary Feynman diagrams, W 0 would only have second order cumulants, which are typically denoted by straight lines connecting the vertices.
The proof for the diagrammatic expansion closely follows the proof of the linked cluster theorem (sometimes also called first Mayer theorem [6]), the connectedness of contributions to the cumulant generating function W . We therefore briefly recapitulate the most important steps of the latter proof (adapted from [2, Sec. 6.1.1]) here and provide all details in Appendix 4.3 for completeness and consistency of notation.
The perturbative corrections W V (j) = W (j) − W 0 (j) to the cumulant generating function obey the operator equation (following from (55)) where the latter constant factor is immaterial for the determination of cumulants. Expanding in ǫ, the function W can be constructed iteratively, where each step produces additional diagrams. The proof of connectedness then proceeds by induction under the assumption that all graphs contained in W l are connected. In the following we will call a "component" a term that appears in W V that is composed of a product of vertices, the Taylor coefficients V , and of cumulants W (n) 0 , the derivatives of the cumulant generating function of the unperturbed theory. Representing V in its Taylor series, every vertex of V ties together components already contained in W l . The Expansion around a theory with first three non-vanishing cumulants. a) Bare vertex, here corresponding to a term ∝ x 3 . b) Cumulants of the solvable part W 0 of the theory: the first three cumulants are non-vanishing, corresponding to circles with one, two, and three legs, respectively.
cumulant generating function is then given as Combinatorics shows that only those graphs survive in the limit that pick up each of their k vertices in a different step in (9), from which the factor ǫ So we see that W l at any intermediate step l may contain connected components with arbitrary numbers of external legs. By the inductive nature of the proof it is therefore clear that the linked cluster theorem holds independently of the number of connected components in W 0 ; that is to say it is generally true, even if already W 0 contains connected components with arbitrary numbers of legs; if it is the cumulant generating function of a non-Gaussian theory.

Recursion for the effective action
We now proceed to derive our main result, the diagrammatic expansion of the effective action Γ. We here explain the main steps and provide all details in Appendix 4.4 and Appendix 4.5.
To lowest order in perturbation theory, setting ǫ = 0 in (7), we get W (j) = W 0 (j); the leading order term in Γ is the corresponding Legendre transform We first need to derive a recursive equation to obtain approximations of the form where we defined Γ V (x * ) to contain all correction terms that have at least one interaction vertex due to the interaction potential V to some order ǫ k in perturbation theory. We will use the iteration in a second step to proof Feynman rules for the diagrams.
It is a priori not clear that the decomposition (11) is useful. We show in Appendix 4.4 that this indeed so and, moreover, that it leads to the recursive operator = exp(−W 0 (j)) exp ǫV (∂ j ) + Γ which has a similar form as (8), defining the linked cluster expansion of W . The term Γ can hence be seen as a monopole vertex. We want to solve the latter equation iteratively order by order in the number of vertices k, defining Γ V,k . Analogous to the proof of the linked cluster theorem, we arrive at a recursion by writing the exponential of the differential operator in (12) as a limit and by expanding the logarithm to obtain the iteration with initial condition We obtain the perturbation correction to Γ (11) as the limit To obtain the final result (17), we need to express j = Γ (1) 0 (x * ) in g l (j). The latter step is crucial to obtain Γ V (x * ) as a function of x * . Because the mean value is given by (3) and (5)), this step effectively expresses all cumulants x n 0 ≡ W (n) 0 of the unperturbed theory in terms of the first cumulant x 0 (x * ) = x * , where x * is the value of the first cumulant of the full theory. The relation j 0 (x * ) can practically be obtained either by inverting 0 (j 0 ) or directly as j 0 = Γ (1) 0 (x * ), if Γ 0 (x * ) is known explicitly. From the last point follows that the differentiation in (15) 0 (x * )) produces an inner derivative Γ This property will become important in the proof that the iteration (13) generates only one-line irreducible graphs. To see this, we have to generalize the notion of irreducibility known from the Gaussian case, which we will do in the following section.

Proof of a generalized irreducibility
The term "one-line irredicibility" in the Gaussian case refers to the absence of diagrams that can be disconnected by cutting a single second order bare propagator (a line in the original language of Feynman diagrams). In the slightly generalized graphical notation introduced in Figure 1, these graphs would have the form where two sub-graphs of k and k ′ vertices, respectively, are joined by a bare second order cumulant . Before proceeding to the proof, we need to define irreducibility in a more general sense than used in the Gaussian case. If a graph can be decomposed into a pair of sub-graphs both containing vertices by disconnecting the end point of a single vertex, we call this graph reducible. In the Gaussian case, this definition is identical to one-line reducibility, because all end points of vertices necessarily connect to a second order propagator, a line. This is not necessarily the case if the bare theory has higher order cumulants. We may have components of graphs, such as , where the three-point cumulant (could also be of higher order) connects to two third (or higher) order interactions on either side. Disconnecting a single leg, either to the left or to the right, decomposes the diagram into two parts, each of which contains at least one interaction vertex. We call such a diagram reducible and diagrams without this property irreducible here. Note that a single leg of an interaction vertex that ends on a first order cumulant does not make a diagram reducible; both components need to contain at least one interaction.

2.4.1.
Cancellation of reducible diagrams We employ the following graphical notation: Since g l (j) = : g l depends on j only indirectly by the j-dependence of the contained bare cumulants, we denote the derivative by attaching one leg, which is effectively attached to one of the cumulants of W 0 contained in g l j g l := ∂ j g l := ∂ j g l (j).
We first note that (13) generates two kinds of contributions to g l+1 , corresponding to the lines (14) and (15), respectively. The first line causes contributions that come from the vertices of ǫV (∂ j ) alone. These are similar as in the linked cluster theorem (9). Determining the first order correction yields with g 0 = 0 which contains all graphs with a single vertex from V and connections formed by cumulants of W 0 . These graphs are trivially irreducible, because they only contain a single vertex. The proof of the linked cluster theorem (see Appendix 4.3) shows how the construction proceeds recursively: correspondingly the l + 1-st step (14) generates all connected graphs from components already contained in W 0 + g l . These are tied together with a single additional vertex from ǫV (x). In each step, we only need to keep those graphs where the new vertex in (14) joins at most one component from g l to an arbitrary number of components of W 0 , hence we maximally increase the number of vertices in each component by one. This is so, because comparing the combinatorial factors (60) and (61), contributions formed by adding more than a single vertex (joining two or more components from g l by the new vertex) in a single step are suppressed with at least L −1 , so they vanish in the limit (17).
The second term (15) is similar to (14) with three important differences: • The differential operator appears in the form ∂ j − x * . As a consequence, when setting j 0 = Γ (1) 0 (x * ) in the end in (17), all terms cancel where ∂ j acts directly on W 0 (j), because W (1) 0 (j 0 ) = x * ; non-vanishing contributions only arise if the ∂ j acts on a component contained in g l . Since vertices and cumulants can be composed to a final graph in arbitrary order, the diagrams produced by ∂ j − x * acting on g l are the same as those in which ∂ j − x * first acts on W 0 and in a subsequent step of the iteration another ∂ j acts on the produced W (1) 0 . So to construct the set of all diagrams it is sufficient to think of ∂ j as acting on g l alone; the reversed order of construction, where ∂ j first acts on W 0 and in subsequent steps of the iteration the remainder of the diagram is attached to the resulting W (1) 0 , is contained in the combinatorics.
• The single appearance of the differential operator ∂ j acts like a monopole vertex: the term therefore attaches an entire sub-diagram contained in Γ  V (x * ) do not depend on j; the j-dependence of all contained cumulants is fixed to the value j = Γ (1) 0 (x * ), as seen from (12). As a consequence, these sub-graphs cannot form connections to vertices in subsequent steps of the iteration.
Considering that Γ (1) V (x * ) is represented by (19), in step l + 1 the line (15) contributes graphs of the form Since by their definition as a pair of Legendre transforms we have we notice that the subtraction of the graphs (22) may cancel certain connected graphs produced by the line (14). In the case of a Gaussian solvable theory W 0 this cancellation is the reason why only one-line irreducible contributions remain.
We here obtain the more general result, that these contributions cancel all reducible components, according to the definition above. We will prove this central result in the following.
To see the cancellation, we note that a reducible graph by our definition has at least two components joined by a single leg of a vertex. Let us first consider the case of a diagram consisting of exactly two irreducible sub-diagrams joined by a single leg, as it is generated by (15). This leg may either belong to the part g (1) or to g (1) l in (22), so either to the left or to the right sub-diagram. In both cases, there is a second cumulant W (2) 0 either left or right of Γ (2) 0 . This is because if the two components are joined by a single leg, this particular leg must have terminated on a W The second point to check is the combinatorial factor of graphs of the form (22). To construct a graph of order k, where the left component has k ′ bare vertices and the right has k − k ′ , we can choose one of L steps within the iteration in which we may pick up the left term by (15). The remaining k − k ′ vertices are picked up by (14), which are ones. Every addition of a component to the graph comes with The symmetry factors s 1 , s 2 of the two sub-graphs generated by (22) enter the symmetry factor s = s 1 · s 2 · c of the composed graph as a product, where c is the number of ways in which the two sub-graphs may be joined. But the factor s, by construction, excludes those symmetries that interchange vertices between the two sub-graphs. Assuming, without loss of generality, a single sort of interaction vertex, ways of choosing k ′ of the k vertices to belong to the left part of the diagram. Therefore the symmetry factor s is smaller by the factor s ′ than the symmetry factor of the corresponding reducible diagram constructed by (14) alone, because the latter exploits all symmetries, including those that mix vertices among the sub-graphs. Combining the defect s ′ with the combinatorial factor (23) yields k! , which equals the combinatorial factor of the reducible graph under consideration.
The generalization to the case of an arbitrary number of sub-diagrams of which M are irreducible and connected to the remainder of the diagram by exactly one link, called "leaves", is straightforward: we can always pick one of these sub-diagrams and replace it by its corresponding contribution to Γ (1) (x * ). Then, by the same arguments as before, we produce the same diagram with the same prefactor, only with opposite sign. In summary we conclude that all reducible graphs are canceled by (22). We therefore obtain the first part of our Feynman rules: All connected, irreducible diagrams that are also contained in the perturbation expansion of W V also contribute to Γ V with a negative sign and the same combinatorial factor.
But there is a second sort of graphs produced by (22) that does not exist in the Gaussian case: If the connection between the two sub-components by Ô ends on a third or higher order cumulant. These graphs cannot be produced by (14), so they remain with a minus sign. We show an example of such graphs in Figure 4c). One might wonder why this contribution does not cancel while making a "reducible impression", if one interprets Γ as a kind of two-point interaction. The solution is that for diagrams of this kind, indeed contributions of opposite signs and same absolute values are generated, but the contribution with the one sign is generated once and the contribution with the other sign twice, therefore the total contribution does not vanish. We address this issue more thoroughly in 4.7.
Moreover, subsequent application of ∂ x * on such components produces graphs that contain higher order derivatives of Γ (n) 0 . By the arguments given in the proof for the generalized irreducibility, we deduce that all diagrams cancel that have at least one leaf from W (1) V (that includes only unperturbed cumulants and interactions, but no vertices Γ (n) (x * )) and that is connected to the remainder of the diagram by a single leg of an interaction (and not via a leg of Γ (n) (x * ), see 4.7). Reducible non-standard diagrams that violate these prerequisites can occur. In the end of the following subsection, we will also give an example for this case. To enumerate all all non-standard diagrams and to find their Feynman rules, we also derive a set of skeleton diagrams in the following subsection that allow the computation of all contributions at any given order including their combinatorial factor.

Perturbative diagrammatics derived from skeleton diagrams
To obtain a better understanding of the types of diagrams that contribute to the perturbation expansion, in particular the non-standard diagrams not contained in W V , it is useful to derive a non-iterative approach based on a set of skeleton diagrams. We use the term "skeleton diagram" for a diagram containing dressed as opposed to bare cumulants and vertices. We start with a Taylor expansion around a adroitly chosen point x 1 where we wrote the first two terms explicitly. We now choose a particular point x 1 by making use of the involutive property of the Legendre transform, which is given for all smooth and convex cumulant generating functions, irrespective of the form of the underlying theory. Concretely, we choose As an immediate consequence, we see that Using these relations and the definition of the Legendre transform, the two terms in the first line of (24) can be rewritten as We now see that the first two terms in (24), besides Γ 0 (x * ), also contain all diagrams from W V (j 0 ), but with the opposite sign. We note that by their dependence on x * and j 0 , respectively, these contributions can be calculated: The relation (25) allows us to express all bare cumulants W The second line in (24) produces additional diagrams that are of second order in the interaction or higher; this is because the "leaves" (26) of these tree diagrams are −W V , which, by definition, contain at least one interaction vertex. These terms include those diagrams that cancel the reducible diagrams included in W V (j 0 ), as proven above. In addition, these terms contain the non-standard diagrams. We will describe in the following, how the exact form and combinatorial factors of these diagrams can be obtained.
To express the dressed vertex functions Γ (n) (x 1 ) that appear in the second line of (24), we first derive a perturbative expansion for Γ (2) (x 1 ), which is basically Dyson's equation. We here include the derivation to see that it holds beyond the Gaussian case. Decomposing W into the solvable part and its perturbative corrections, the reciprocity relation Multiplying from left with Γ 0 (j 0 ), we obtain by rearranging The latter term contains at least one interaction vertex in W V , the former is of order zero. So we may iterate this equation to obtain an inverted Dyson's equation for Γ (2) that reads

It is easily checked by insertion that (30) solves (29).
So we have expressed the second order vertex at x 1 by means of quantities that are directly accessible or can be calculated perturbatively in a straight-forward way, by their dependence on x * and j 0 , respectively. Higher order vertices are determined in the standard way by differentiating the reciprocity relation (28) multiple times with respect to j 0 , and removing the legs W (2) (j 0 ) that arise from inner derivatives by multiplication with the corresponding inverse Γ (2) (x 1 ). The result is the well-known decomposition of vertex in terms of tree skeleton diagrams [1,2] Diagrammatically, the first of these relations reads We can hence express vertices of arbitrary order at the point x 1 by cumulants at the point j 0 (which is given by x * by means of (25)) and Γ (2) (x 1 ), which is again given by the same cumulants and Γ (2) 0 (x * ) by Dyson's equation (30). The important point to note here is that the set of rules to translate skeleton diagrams into their perturbative expansion is problem independent; it just results from the properties of the Legendre transform and the decomposition into solvable part and perturbation.
Because both factors in every term of the sum in the second line of (24) depend on the interaction V , the order of each such term is the sum of the orders of its factors. We derive two additional rules from this view. First, since also W 0 (j 0 ) and Γ 0 (x * ) form a pair of Legendre transforms, their reciprocity relation W 0 (x * ) = 1 gives rise to a relation corresponding to (31) . This means that the lowest order terms for the n-th derivative in (31) resum to Γ (n) 0 (x * ). The final set of skeleton diagrams derived from (24), the main result of this section therefore reads , which contains all diagrams from (30) and (31) that have at least one interaction vertex, so which are of order O(ǫ) or higher.
The second rule follows from W V (j 0 ) being of first order or higher, so that at the n-th order in the interaction, we have to consider at most the n-th term in the sum in (32). In this last term we must in addition drop theΓ-term. The only term to consider at second order, for example, is This term cancels the contributions of the diagrams b) and c) in Figure 3. However, that all reducible diagrams cancel can be seen more systematically by the constructive arguments given in Section 2.4.1.
But eq. (27) is useful to get a handle on the non-standard diagrams. If we want to know all such diagrams of a certain order k in the interaction, we first assign the number of interactions to be contained in each factor in any term 0 (x * ), as said above. For non-zero numbers of vertices inΓ (n) , we decompose the latter via (30) and (31) into an expression only containing full cumulants and Γ (2) 0 (x * ) and distribute again the number of interactions assigned toΓ (n) among these cumulants. Each full cumulant is then broken down into all its perturbative expansion, for which the standard Feynman rules hold according to the linked cluster theorem (see Appendix 4.3). Throughout this expansion, we skip all diagrams that are also contained in −W V , because these are just the reducible ones that we know to cancel each other.
Another insight we gain by (32) in combination with (30) is that non-standard diagrams are not necessarily irreducible -in fact, they can even be reducible in the sense that they fall into two parts both containing interactions if we cut a second-order cumulant. We show an example for this phenomenon in a toy model at the end of Section 2.5.

Summary of Feynman rules for the non-Gaussian case
We now summarize the algorithmic rules derived from the above observations to obtain Γ: explicitly by finding j 0 that extremizes the right hand side. At this order g 0 = 0. (ii) At order k in the perturbation expansion: (a) add all irreducible graphs in the sense of the definition above that have k vertices; (b) add all graphs containing derivatives Γ (n) 0 as connecting elements that cannot be reduced to the form of a graph contained in the expansion of W V (j 0 ); the graphs left out are the counterparts of the reducible ones in W V (j 0 ). The topology and combinatorial factors of these non-standard contributions are generated iteratively by (13) from the previous order in perturbation theory; this iteration, by construction, only produces diagrams, where at least two legs of each Γ (n) 0 connect to a third or higher order cumulant. We can also directly leave out diagrams, in which a subdiagram contained in W V is connected to the remainder of the diagram by a single leg of an interaction vertex. Alternatively, the additional diagrams can be constructed directly for any order k by instantiating the skeleton diagrams generated by (32) to the desired order, as explained in Section 2.4.2. Note that here diagrams may appear in which Γ (n>2) couples to less than two third or higher order cumulants. After decomposing all Γ (n>2) into Γ (2) 0 and diagrams from W (n) by (31) and (30), we can discard diagrams, in which at least one subdiagram from W V forms a leave that is connected by a single link of an interaction vertex; these cancel by the same argument as given in the iterative approach.
(iii) assign the factor ǫ k r1!···r l+1 ! to each diagram with r i -fold repeated occurrence of vertex i; assign the combinatorial factor that arises from the possibilities of joining the connecting elements as usual in Feynman diagrams (see examples below and e.g. [3]); (iv) express the j-dependence of the n-th cumulant x n (x * ) in all terms by the first cumulant 0 (j 0 ); this can be done, for example, by inverting the last equation or directly by using j 0 = Γ (1) 0 (x * ); express the occurrence of Γ This set of Feynman rules constitutes the central result of our work. Note that the rules include the well-known Gaussian case, because irreducibility in the here-defined sense is identical to one-line-irreduciblity in the expansion around a Gaussian theory: every leg of an interaction vertex necessarily connects to a line there. Non-standard diagrams hence cannot appear in this case.
The Feynman rules hold for any order k. So one may compute corrections at order k directly without having computed any lower order. The combinatorial factor is fixed unambiguously, too. To get an intuitive understanding of the cancellation, it is still useful to illustrate the recursive construction of graphs in Section 2.5 in the application to a minimal, non-Gaussian setting. In Section 2.6 we demonstrate the application to systems of Ising spins, recovering the TAP approximation [25], the high temperature expansion [23], and the Plefka expansion [28]. It turns out that the recursive algorithmic procedure leads to a dramatic decrease of required computations.

Illustrative example for the graphical rules
As a first example let us consider a zero-dimensional field theory, that is to say, a probability distribution of a scalar variable. By assumption, S 0 constitutes the solvable theory, so that W 0 (j) can be determined exactly by (2). We here illustrate the method for a solvable theory that has non-vanishing cumulants of orders one, two, and three, hence W 0 is a polynomial of order three where the cumulants x n 0 with n ∈ {1, 2, 3} appear as Taylor coefficients. Its graphical representation is shown in Figure 1b. As perturbation we assume a three point vertex ǫV (x) = ǫ x 3 , shown in Figure 1a.
We now calculate the diagrams recursively according to step (ii). At order l = 0 we have g 0 = 0. At first order l = 1, we therefore get from (21) the graphs shown in Figure 2.
All diagrams contributing to g 1 at first order l = 1 generated by (21) for the theory of Figure 1.
Each graph has a single vertex, thus we get a factor ǫ 1! . The combinatorial factors for each graph are stated explicitly above. So the diagrams in Γ V at first order are identical to all connected diagrams with a single vertex; in this regard, the expansion is identical to the well known Gaussian case, except that the diagram Figure 2 c) would not appear in the absence of third order cumulants of the solvable theory.
At order l = 2, the graphs in W 0 ∪ g 1 contain those that have been produced in the previous step. The step (14) therefore produces additional diagrams out of these components, some of which are shown in Figure 3. Figure 3. Some diagrams generated by (14) contributing to g 2 − g 1 at order l = 2. Diagrams b) and c) are reducible, so they will drop out.
The diagrams a) and b) in Figure 3 are composed of the diagram a) in Figure 2 contained in g 1 and one additional vertex. The diagram c) in Figure 3 is composed of diagram b) in Figure 2; its combinatorial factor 3 · 3 is the combined factor of the first order diagram and the factor 3 due to three possibilities to select a leg of the vertex. We here skipped further diagrams; in particular all diagrams, that are generated by the first order diagram c) in Figure 2.
The line (15) produces additional diagrams with negative sign. We here indicate the simplification Ô = 1 by an overbrace. Some diagrams are shown in Figure 4.
additional non−cancelling diagram We observe that the diagram a) in Figure 4 cancels the diagram b) in Figure 3. This cancellation is of the ordinary type; the diagram b) in Figure 3 is one-line irreducible in the original sense: its two components are connected by a second order cumulant, which would be denoted by a line in the original graphical language of Feynman diagrams; here denoted as .
An example of the more general cancellation appears between diagram b) in Figure 4 and c) in Figure 3: the two components in Figure 3c) are not connected by a second order cumulant, hence the diagram is not one-line irreducible in the original sense, but it is reducible in the sense we defined above. We see that it is canceled by Figure 4b), consistent with the general proof. If we want to obtain the diagrams of the next order, note that we have to keep reducible diagrams, because the cancellation occurs only after setting j = Γ (1) 0 (x * ). The diagram Figure 4c) produced by (15) cannot be generated by (14), because the connecting proper (effective) vertex Γ (2) 0 = Ô connects to two third order cumulants. It is generated in two ways; with two Γ (1) (x * )-components or only one of them. These two contributions have opposite signs, but different prefactors, therefore they do not cancel (see also 4.7 for details). We now derive the appearance of this latter diagram by the application of the method using skeleton diagrams. The relevant skeleton diagram is the one that appears at second order on the Taylor expansion in the second line of (24) 1 2! Ô .
We now instantiate this diagram with perturbative expansions of the desired order, so that the resulting diagram is of order two: The two point vertex function, by Dyson's equation (30) to lowest order is Γ (2) (x 1 ) = Γ V attached to either side, to first order, contain contributions of the form where all elements on the right are cumulants of the unperturbed system and bare interaction vertices. The combinatorial factor of this latter sub-diagram is computed with the standard rules for connected diagrams, which yields a factor 3 from the possible ways of attaching the first order cumulant to the interaction vertex. There is only a single possibility to attach the additional external leg to the third order cumulant and only a single way to join the elements in the skeleton diagram in (33), besides the explicit factor 1/2!. Finally, each remaining diagram must be interpreted algebraically by steps (iii) and (iv). These steps are straight forward here. We note that this step is possible at all in the technique using the skeleton diagrams, because all cumulants appearing are W 0 (x * )) and in Dyson's equation (30) in addition only Γ (2) 0 (x * ) appears; so all elements, by assumption of S 0 being the solvable theory, are known.
There is yet another observation that we can make in the framework of of this toy model, namely that there are reducible non-standard diagrams that contribute to Γ V . One example is the following fourth-order cumulant emerging from the n = 2 term in the sum in (27) and using the first-order correction from (30) (with second order correction in the interaction for W

Application to the Ising model
We here illustrate the method on a simple system, the classical Ising model with the action where x ∈ {−1, 1} N is a vector of Ising spins and J a symmetric matrix with J ii = 0.
Note that there is no additional constraint on J; it could, for example, be drawn from a random distribution, as done in spin-glass models. In particular there is no restriction to nearest-neighbour interactions. We are here interested in the weak coupling limit with small ǫ.
There are essentially two different ways to arrive at an approximation of the effective action Γ(x * ). The first approach represents the pairwise coupling term as the result of a Gaussian average over auxiliary variables (see e.g. [1, Chapter 4.3 eq. 4.50a], [5, chapter (5.2), eq. (5.49)], [29, eq. (4), (5)], or [30]). This Hubbard-Stratonovich transform reduces the calculation of Z to the summation of n-th moments of the Gaussian employing Wick's theorem, weighted by the Taylor coefficients of V (x) = ln xi e jixi = i ln 2 cosh j i ; the latter play the role of vertices here. Consequently, standard Feynman diagrammatic rules apply. This approach is only applicable to this model because the interaction is pairwise; interactions of higher order cannot be written by help of Gaussian auxiliary fields.
The second way, which we will follow here, considers as the solvable part the single-spin term of the action which is directly summable, yielding a cumulant-generating function W 0 whose decomposition into a sum over individual sites shows their statistical independence. For each i, the solvable theory therefore has the infinitely many non-vanishing cumulants of a binary variable. The interaction is in turn treated as the perturbation This approach has been followed in quite a number of variations over decades (see e.g. [ (37) is treated as a perturbation and (36) as the solvable part. In particular this approach is identical to Plefka's method [28]. Both routes of course lead to identical results. To second order in ǫ one obtains the mean-field theory presented by Thouless, Anderson, and Palmer (TAP) [25] as a "fait acompli", without proof, but mentioning a previous diagrammatic derivation. Indeed, Vasiliev and Radzhabov ( [21,22], summarized in [5, section 6.3.1, 6.3.2, 6.3.4]) had derived this result diagrammatically with the help of the analog to the Dyson Schwinger equation for the effective action [35] before (see [5,

section 5.2] for a review of this method).
If one performed a perturbation expansion directly on the level of Z (55) with f (x) 0 = x f (x) exp(j T x), one would quickly obtain unwieldy expressions; but most of these terms cancel by the subsequent transformations Z ln → W L → Γ in the final result, making a direct expansion of Γ desirable. The calculation of Γ at higher orders in ǫ by this direct approach has prompted for computer algebra systems [26]. Georges and Yedidia [23,33] developed a sequence of shortcuts and tricks for this problem by which they managed the approximation up to forth order (reviewed in [34, chapter 3]). As mentioned, Γ has been derived diagrammatically by Vasiliev and Radzhabov [21, eq. (21)]. Using Dyson-Schwinger equations, they explicitly presented all orders up to and including the third.
We here chose this model to exemplify the application of the general framework exposed in previous sections. For one, because it allows the comparison to a wealth of methods developed over decades, as mentioned above. And also because the Ising model has proven useful in many other applications than micromagnetism, for example artificial neural networks [36]. Especially the TAP-approximation and its higher order corrections are employed to derive analytical approximations [37] to contrastive divergence [38], a learning rule commonly employed to train restricted Boltzmann machines [39,40]. Another major field, in which the Ising model is applied, are inference problems, where the inverse Ising problem has to be solved, that is, the sources h i and couplings J ij have to be computed for given means and pairwise covariances. Depending on the problem, the spins represent activities of neural units [41][42][43], active or inactive genes in the case of gene regulatory networks or participants in financial markets (see [44] for an excellent review of the inverse Ising problem).
The most natural quantity for this type of problem would be the second Legendre transform, the usual entropy of the multi-variate binary distribution for given first two moments, because the couplings and sources turn out to be the extrema of this quantity. We present an iterative method to compute the second Legendre transform in 4.6. Still, the first Legendre seems to be commonly used also for the inverse problem [27,45]. The reason for this might be that calculating the second Legendre transform is technically challenging; nevertheless, specific approximations exist that are valid for small correlations, first computed by a technique specialized to the Ising model [46], later obtained by techniques borrowed from the functional renormalization group [47], using the Wetterich equation [17].. An extension to more than pairwise coupling (37) is desirable in these fields [48][49][50]. In particular the Hubbard-Stratonovich method mentioned above would not work in these cases, while the here proposed approach does.
Performing the here presented diagrammatic expansion, we observe a drastic decrease of computations and their transparent organization by diagrams. The calculation up to fourth order indeed takes only minor effort.
We start with the unperturbed theory, given by (36), in which all spins decouple and the cumulants up to fourth order read To zeroth order the Legendre transform of W 0 is the entropy of a binary variable where (1 ± m i )/2 are the probabilities for x i = ±1. The result (39) follows from step (i) of the algorithm by a short calculation from (36); the form is moreover clear, because the distribution (36) maximizes the entropy and the Legendre transform to Γ just fixes j so that the mean is x i = m i . To obtain corrections to (39) we represent the n-th cumulant by an empty circle with n legs and the n-th derivative of Γ 0 by a hatched circle with n legs. An interaction J ij is denoted by an edge: To evaluate the diagrams, we will only need the second derivative Γ 0 , either directly given by differentiating (39) twice or by using the relation W . Within this language, the perturbative corrections up to the third order to be added to (39) are readily constructed from steps ii-iv at the end of Section 2.4: For the third order diagrams, the symmetry factors 2 3 and 2 2 are noted separately in front of the respective terms. For all diagrams, they are determined as usual for Feynman diagrams (see Appendix 4.8). Until this order only the first sort of ordinary diagrams composed of vertices and cumulants contribute. All diagrams are irreducible in the general sense defined above.
Going to fourth order, we first consider the four diagrams that are formed out of vertices and cumulants alone. These would also contribute to an expansion of W .
In the appendix, Appendix 4.8, we show that the symmetry factors are given by S R = 48, S TM = 96 and S dC = 8. The fourth standard diagram is shown below. It is similar to the first and second non-standard contribution. To indicate the origin of the different contributions, we denote sub-graphs originating from Γ (1) by filled circles -they are, however, translated in the same way as the empty ones: With the symmetry factors given by S a = 48, S b = 24 and S c = 4 (see appendix, sec. Appendix 4.8 for the derivation), we can add up the contributions of the three diagrams, which gives, leaving out factors that are equal in all diagrams for simplicity: The last term looks like a second cumulant and, interestingly, this contribution is indeed exactly canceled by the contribution of the ring-diagram of fourth order (45) in the case that exactly one pair of indices, which belong to cumulants represented at opposite sites of the ring, are equal. This has to be so, as shown by Vasiliev and Radzhabov [22], because, according to their terminology, diagram a is a "nonstar graph", that is, not a star graph -it decomposes into two parts by removing a cumulant (a "vertex" in their words). Note that the star graph property implies irreducibility, but not vice versa. As Vasiliev and Radzhabov have shown, only star diagrams contribute to the effective action of the Ising model, if we introduce so called compensating graphs. These graphs are constructed as follows: Take all star graphs of a certain order and draw contractions (in [22] depicted by dashed lines) in all possible ways between circles that are not connected by an interaction and therefore represent indices that could take the same value (keep in mind that J ii = 0). If two or more circles are contracted, they are considered as a single circle, are associated with a common index and are translated as the product of all cumulants associated to this point. Then substract all compensating graphs that are non-star-graphs and sum up original and compensating graphs [22, eq. (8)]. We may check that our result is consistent with this rule. In the ring-diagram (45), the only contraction that leads to a non-star graph is the one that identifies the indices of two opposite circles. This contribution therefore cancels from the final expression for the effective action, analogous to the summed up contribution of the diagrams a, b, and c. In our framework not only the double Citroën diagram (47) contributes to the sum over only two unequal indices, but also the ring diagram in the case that the indices of pairwise opposite circles and the diagrams a, b and c in the case that the indices of the outmost circles are equal (see Appendix 4.8 for details). In summary, the general method that we present here is consistent with the specific result known for the Ising model.

Discussion
We presented a systematic diagrammatic scheme to calculate the effective action for any problem that can be decomposed into a solvable part with known correlators and a perturbing part. We have proved that corrections are composed of two types of diagrams: • irreducible diagrams in a more general sense than in the Gaussian case: those that cannot be decomposed by detaching a single leg of a vertex; • diagrams of special form that are neither contained in the perturbation expansions of Z nor W ; they are composed of sub-graphs that are connected by either a vertex Γ (n) 0 , n ≥ 2 that, by at least two legs, connects to a third or higher order bare cumulant. This set of diagrams is found by instantiating a known sequence of skeleton diagrams to the desired order in perturbation theory.
The appearance of the latter diagrams can be regarded a as generalization of the amputation in the Gaussian case, because the inverse propagator Γ (2) 0 may attach to cumulants of the solvable problem at any order, not only to second order cumulants; only in the latter case it "amputates" the lines.
The presented inductive proof in addition yields an iterative equation which allows the algebraic construction of all graphs and their combinatorial factors from elementary rules of calculus.
One may wonder why we derived the recursion for Γ based on the ideas of the proof of the linked cluster theorem (see Appendix 4.3, extending the proof in [2, Sec. 6.1.1] to the non-Gaussian case). It might seem more direct to modify the corresponding proof of one-line irreducibility to the non-Gaussian setting considered here. The latter, however, appeared impossible to us: The elegant proof by Zinn-Justin [2, section 6.5] rests on the assumption that the underlying theory is Gaussian; in their eq. 6.59 each line is disconnected in all possible ways and the result is shown to remain connected. The proof hence requires that the only connecting elements of the bare theory be lines; this is precisely the restriction we lift here.
The proof by Weinberg [51, section 16.1] shows elegantly that W decomposes into tree graphs, whose vertices -which one could call "effective vertices" in this case -are generated by Γ. This statement remains of course true also in the non-Gaussian case, because W L ↔ Γ form a pair of Legendre transforms. We use this fact to derive the decomposition into skeleton diagrams. This decomposition also follows directly from ‡ In comparing the expressions to these works, note the different sum conventions: Georges et al. and Nakanishi et al. are only summing over those tuples of distinct indices that lead to different terms, while we allow multiple occurrence of terms already in the definition of the action (35) and correct this by using J ij 2 instead of J ij as the interaction. Note that the former sum convention in the second-to-last line of eq. (11) in Georges et al. amounts to a summation over all tuples of three different indices and all even permutations thereof, whereas in eq. (15) of Nakanishi et al, these permutations are written out explicitly, therefore the summation is only over the tuples with three different indices. Therefore, even if these two expressions seem to disagree by a factor 3 at first sight, they do not, because, depending on the form of the sum, the interpretation of their sum convention changes. We thank Adam Rançon for clarifying this summation convention. the reciprocity relation 1 = Γ (2) W (2) of the Hessians [2, section 6.2]. The connecting elements in these trees are the full propagators W (2) , which Dyson's equation expresses in terms of bare propagators and the self energy. The reciprocity relation implies the recursion for the self-energy Γ 0 . Since W V is composed of all connected graphs, including those containing higher order cumulants of W 0 , we see that the terms · · · Γ (2) 0 · · · appearing in the iteration produce those unusual terms of the form W An alternative approach is the analog of the Dyson Schwinger equation for the effective action [35,52,53], reviewed in [5, sec. 1.8, 6.2]. It leads to equations of motion for Γ that enable an iterative expansion, for example in the interaction strength. An iterative solution leads to a proof of the 1PI property in the Gaussian case [5, sec. 1.8.3]. We believe that the latter approach could be generalized to obtain the same results as presented here. The equation of motion can be derived from the invariance of the integral measure with respect to translations [5, sec. 1.8.3], but also extends to problems on discrete state spaces, such as the Ising model [5,21,22,32]. From this equation of motion, eq. (6.136) in [5] is derived, similar to our (32) with the difference that Vasiliev expands in terms of unperturbed cumulants and not full vertices. Probably, both equations are suitable starting points to rederive our results in Section 2.4. The Feynman rules would then be deduced from this step. Possibly these approaches could render the taxonomy of non-standard diagrams clearer. However, it will unlikely be as simple as in the Gaussian case because, as pointed out, reducible non-standard diagrams do not necessarily cancel. We leave this question for future work. In any case, more refined diagrammatic rules previously had to be obtained in a case by case manner so far, depending on the problem at hand (compare also the last paragraph of [5, chap. 6.3.1.]). The algorithm presented in Section 2.4, in contrast, is the same for any model.
One notes that eq. (21) in the work by Vasiliev and Radzhabov [21] is identical to the third order derivation recovered much later by other means [23,[26][27][28]; in particular the result already includes the TAP approximation [25], if only the first three terms are considered. In principle, these early works [21,22] also derive the Feynman rules for the Ising model; however, without giving any concrete expression for orders higher than three. In contrast to these results, the set of Feynman rules that we present here are applicable to general non-Gaussian theories. In case of the Ising model these had been sought for some time (see [34, p. 28]).
We hope that the presented technique may prove useful in finding new approximations around known limiting cases. Examples may include expansions of the Hubbard model around the atomic limit [23]. An extension of our theory to higher order Legendre transforms, as broadly discussed in [5], could lead to a diagrammatic formulation of the results derived in the context of the inverse Ising problem [46,47] and to an alternative field-theoretic formulation of the extended Plefka-expansion for stochastic systems, that has recently been developed [54]. The appendix 4.6 presents an iterative algorithm as a first step towards this goal, an iterative equation to compute corrections to the second Legendre transform. Further work is required to derive a set of Feynman rules from this equation.
The application of the here presented method is possible whenever a model admits a closed-form solution. An interesting regime of application may therefore be spherical models [55]. In the thermodynamic limit, the free energy, the cumulant-generating function of the model, is known. Extensions of the spherical model that include four point coupling terms for example appear in the field of random lasers [56]. If the quartic term is small compared to the quadratic term, the here proposed method could be applied to obtain approximate self-consistency equations. Such quartic spherical models, moreover, appear in inverse problems of diverse systems [57]. Future work is needed to see if the perturbative results offered by the current work may help at obtaining approximate solutions to such inverse problems.
We also note that the procedure as presented in Section 2.4.2 yields a non-iterative way to generate all non-standard diagrams. We believe that this algorithm should be amenable to numerical implementation. The diagrammatic Monte Carlo technique, for example, presents an effective method to calculate the kernels of Schwinger-Dyson equations -for the one-particle Greens function this is the self-energy [58][59][60]. One way to derive Schwinger-Dyson equations relies on multiple Legendre transforms, expressing the generating functionals in terms of correlation functions instead of potentials, as pioneered by de Dominicis and Martin [7]. The equations of state then constitute self-consistency equations for these correlation functions, such as eq. 6.11 in [5]. For example one needs to consider the second Legendre transform to determine the self-energy and the fourth to determine the two-particle Green's function selfconsistently. We do a first step towards developing diagrammatic rules for the second Legendre transform in the appendix 4.6, whereas higher order Legendre transforms are beyond the scope of the current work. Their diagrammatics is already quite involved if the underlying theory is Gaussian [5, chap. 6.28.f]; the parquet equations form an approximation resulting from the fourth order Legendre transform. A different way to obtain a set of self-consistency equations, however, is by a direct resummation of diagrams, as in the derivation of the Bethe-Salpeter equation [61]. The diagrammatic rules derived here could be useful in such an approach.
In general, the methods also seems particularly promising for hierarchical problems: assuming that a problem can be decomposed into small, but strongly interacting clusters that can be solved exactly, the method may be used to systematically expand in the interaction strength between such clusters.

Definition of the effective action
To define the effective action Γ(x * ) we eliminate the dependence on the source field j in favor of the mean value x * (j) := x = ∂ j W (j) by using (1) and (2) and by following the usual background field method [13, Chapter 3.23.6], briefly summarized here: We express W as the integral and then separate the fluctuations δx = x − x * from the background value x * to get For given x * , we now choose j in a way that x * = x (j) becomes the mean value of the field, so that δx has vanishing mean where we used (53) in the last step. Since the exponential function is monotonic, exp(x) ′ > 0 ∀x, the latter expression vanishes at the point where the exponent is stationary The condition (54) has the form of a Legendre transform (4) from the function W (j) to the effective action Γ(x * ). The supremum follows from stationarity (54) and because W is convex down, its Hessian W (2) , the covariance matrix, is positive definite (cf. Section 4.2 or [2, p. 166]). Therefore, −W (j) + j T x * , as a function of j, is convex up (concave), so we may define the Legendre transform (4) by the supremum over j.

Convexity of W
W is convex, because W (2) is the covariance matrix: it is symmetric and therefore has real eigenvalues. For covariance matrices these are in addition always positive [2, p. 166]. This can be seen from the following argument. Let us define the bilinear form A positive definite bilinear form has the property f (η) > 0 ∀η. With δx := x − x we can express the covariance as W (2) kl = δx k δx l , so we may explicitly write f (η) as which is the expectation value of a positive quantity.

Linked cluster theorem
The following proof of connectedness of all diagrammatic contributions to W does not rely on the solvable part S 0 being Gaussian. We here start from the general expression (1) to derive an expansion of W (j), using the definition (2) to write Taking the logarithm, the latter factor turns into an additive constant ln Z0(0) which ensures W (0) = 0. Since we are ultimately interested in the derivatives of W , namely the cumulants, we may drop the constant and ensure W (0) = 0 by finally dropping the zeroth order Taylor coefficient.
The idea to prove connectedness follows to some extent [2]. The proof is by induction, dissecting the operator exp (ǫV (∂ j )) appearing in (55) into infinitesimal operators using the definition of the exponential function as the limit For large L given and fixed and 0 ≤ l ≤ L, we define It fulfills the trivial recursion exp(W l+1 (j)) = (1 + ǫ L V (∂ j )) W l (j) from which follows an iteration by multiplication with exp(−W l (j)) and taking the logarithm The desired result W (j) then follows as the limit W (j) = lim L→∞ W L (j). Expanding We start the induction by noting that for l = 0 we have W l=0 = W 0 , the cumulant generating function of the solvable system. At this order, W hence does not contain any diagrammatic corrections; so in particular no disconnected ones. We assume that the assumption is true until some 0 ≤ l ≤ L. Stated precisely, we assume that all perturbative corrections to W l (j) with k vertices are connected and are ∝ ǫ L k ; the sub-leading terms O ǫ L 2 in (58) vanish in the limit L → ∞, as shown below. Representing the potential V as a Taylor series, we see that each step adds terms of the form shown in the second line where ∂ i is used in short for ∂ ji . The Taylor coefficient V (n 1 ,...,n N ) n1!···nN ! is graphically represented by a vertex (see Figure 1). Noting that the two exponential factors cancel each other after the differential operator has been applied to the latter factor, what remains is a set of connected components of W l (j) tied together by the vertex V (n 1 ,...,n N ) n1!···nN ! . Disconnected components cannot appear, because there is only a single vertex; each of its legs belongs to one ∂ i , which, by acting on W l (j), attaches to one leg of the components in W l .
The iteration (59) shows that the second term in each step adds to W l (j) a set of diagrams to obtain W l+1 (j). It is clear from the single appearance of V (n1,...,nN ) that in each iteration only one such additional vertex is added. We will show now that we, moreover, only need to consider such additional diagrams, where the new vertex V (n1,...,nN ) connects to a perturbative correction contained in W l (with k ≥ 1 vertices), while all of its remaining legs connect to a cumulant of the unperturbed theory W 0 (with k = 0 vertices). Stated differently, a perturbative correction with k vertices picks up each of its vertices in a different iteration step l in (59); contributions where a single iteration step increases the number of vertices in a component by more than one vanish in the L → ∞ limit.
To understand why this is so, we consider the overall factor in front of a resulting diagram with k vertices after L iterations of (59). In each step of (59) the first term copies all diagrams from W l . The second term adds those formed by help of the additional vertex V (n1,...,nN ) . Following how one particular graph is generated by the iteration, in each step we have the binary choice to either leave it as it is or to combine it with other components by help of an additional vertex.
We first consider the case that each of the k vertices is picked up in a different step (at different l) in the iteration. Each such step comes with a factor ǫ L and in there are L k ways to select k out of the L iteration steps to pick up a vertex to construct this particular diagram. So in total we get a factor which is independent of L. Now consider the case that we pick up the k vertices along the iteration (59) such that in one step we combined two sub-components with each one or more vertices already. Consequently, to arrive at k vertices in the end, we only need k ′ < k iteration steps in which the second rather than the first term of (59) acted on the component. The overall factor therefore is ǫ L We can hence neglect the latter option and conclude that W (j) is composed of all connected graphs, where a perturbative correction with k vertices comes with the factor ǫ k k! given by (60). By the same reasoning we may neglect the O ǫ L 2 term in (58), because also here in a single step we would increase the order of the diagram by more than one factor L −1 .

Operator equation for Γ V
Let us first see why the decomposition into a sum in (11) holds. To this end, we consider (6) with δx = x − x * and use the decomposition (7) of the action as well as the decomposition (11) of Γ to obtain Collecting the terms depending on x * on the left hand side we get with (5) where we moved the perturbing potential in front of the integral, making the replacement x → ∂ j and we identified the unperturbed cumulant generating function exp(W 0 (j)) = Z −1 (0) x exp S 0 (x) + j T x from the second to the third line. With the term Γ (1) 0 (x * ) , which follows from the definition (10). Multiplying with exp(−W 0 (j)) j=Γ (1) 0 (x * ) from left then leads to a recursive equation (12) for Γ V , which shows that our ansatz (11) was indeed justified and that we may determine Γ V recursively, since Γ V appears again on the right hand side.

Recursion for g l
To construct the diagrams iteratively we write the perturbing term in (12) as Inserted into (12) we assume L fixed but large and choose some 0 ≤ l ≤ L. We define G l (j) as the result after application of l factors of the right hand side of (62) Obviously we have G 0 ≡ W 0 .
(63) For l = L → ∞, we obtain the desired result as .
By its definition, G l (j) obeys the trivial iteration exp(G l+1 (j)) Multiplying from left with exp(−G l (j)) and taking the logarithm on both sides while using ln(1 + 1 L x) = 1 L x + O(L −2 ) we get the recursion for the additional diagrams produced in step l + 1 . By the initial condition (63) and the form of the additive iteration (64) it is clear that all graphs of W 0 are also contained in G l for any step l. We may therefore define only the perturbative corrections as g l := G l − W 0 .
We first note that indeed (64) yields a closed iteration: Constructing the graphs of up to order l + 1 in G l+1 by (64) we only need the graphs in G l on the right hand side of (64), which, by construction, are of order ≤ l . This is because V contains exactly one bare vertex and Γ

Second Legendre transform
We here extend the iterative procedure to compute perturbative corrections to the second Legendre transform. To this end, we define an effective action that is a function of the first moment x * and the second moment c * = x 2 . We express W as the integral Here k T x 2 must be understood as a bilinear form in x, hence il k il x i x l . We have We would like to define an effective action that is a function of these latter coordinates. So we define Γ(x * , c * ) as with the additional constraints that ∂/∂j and ∂/∂k of the right hand side vanishes. Consequently, the so-defined function Γ fulfills the two equations of state ∂Γ ∂x * = j, ∂Γ ∂c * = k, which can be obtained by application of the chain rule. The second of these equations, for example, follows as Hence we may write the definition of Γ also as . We next decompose S(x) = S 0 (x) + ǫV (x), where we assume that we may compute the cumulant generating function W 0 (j, k) = ln x exp(S 0 (x) + j T x + k T x 2 ) exactly. Correspondingly, we assume a decomposition of the effective action into the solvable part Γ 0 and the perturbative corrections Γ V as With the notation ∂Γ ∂x * =: Γ (1,0) and ∂Γ ∂c * =: Γ (0,1) we may express the integral equation (66) as or, bringing all terms that are independent of the integration variable x to the left, we get with j 0 = Γ (1,0) 0 We hence obtain the operator form of the equation as . Using the representation of the exponential function as a series and the series expansion of the logarithm to first order, we get the iteration with the initial condition g 0 = 0. The perturbative corrections to the effective action then result as the limit .
One could use this iterative procedure to compute perturbative corrections. It seems that a proof of a generalized irreducibility should follow along similar lines as in the case of a first order Legendre transform. In particular, the term Γ (0,1) V (∂ k − c * ) establishes a double link between a component contained in g l and one in Γ V , thus suggesting a more general 2PI property. We will leave a more careful consideration open for subsequent works.

Taxonomy of reducible diagrams
One might wonder why the diagram in Figure 4c) remains in the perturbative expansion of Γ despite being reducible in the sense that it is divided into two parts by cutting one of the legs of Γ (2) 0 , which otherwise acts similarly as a leg of a bare interaction vertex. We here provide an explanation in the framework of the iterative construction of Γ V . There are four different types of diagrams that are reducible in the sense that they fall apart if one detaches one leg from a cumulant. We classify them by the type of subdiagrams (leaves) they are composed of and how these leaves are connected. We can distinguish five cases: I Two irreducible diagrams from W V connected by a single leg that belongs to a bare interaction vertex. Ia Multiple irreducible diagrams from W V that are connected to the remainder by the same junction as in case I.
II Two irreducible diagrams from W where empty and hatched circles in this section always denote the unperturbed quantities W (n) 0 and Γ (2) 0 , respectively. We therefore see that we can replace both leaves by a part of Γ (1) V (x * ). So, together with the original diagram, we have four contributions that all contribute with the same magnitude, two of them positive, two of them negative: The sum is 0. Schematically, the situation is as follows: Here the circles filled with squares denote some arbitrary continuation of the respective diagram. We conclude that in second order in the interactions, no diagrams with irreducible W (1) V -leaves contribute. We will proof by induction in the number of vertices n that no diagram described in case Ia (with at least one reducible link) contributes to Γ V . Case I is our induction assumption (n = 2). Assume that the assumption holds ∀ k ≤ n. We observe that if the whole diagram contains n + 1 interactions we can pick an arbitrary leaf that is connected to the rest by a single link, so that the (reducible) rest contains at most n interactions. By the induction assumption, this rest is not contained in the diagrammatic expansion of Γ V . Therefore, we have just one contribution, which can be replaced by a part of Γ This concludes the induction. We conclude from these observations that no diagrams with a single irreducible W V -leaf contribute.
Considering case II, it is obvious that we can choose one of the leaves to be provided by Γ (1) V and the other one to be composed of bare vertices. We can also insert a unity (68) to the left of the given Γ (2) 0 , so that we can interpret both leaves to come from Γ (1) V (x * ). Because the diagram composed solely of bare vertices is missing, we only have three contributions. They all contribute with the same absolute value, but two of them with one sign, one with the other sign and therefore, the contributions do not sum to zero and the whole diagram contributes to the expansion of Γ V , as claimed in the main text. Diagrammatically, this is illustrated by In point II, we have proved that there are diagrams contributing to Γ V , but not to W V (for the example Figure 4c)), therefore the situations of points III and IV can appear.
The diagrams described under point III cancel, because the irreducible components within the W V -part or can be left as it is. Both contribute with same absolute value, but opposite sign and cancel. Because the other part cannot be composed just by bare vertices and unperturbed cumulants, it occurs only once, namely as part of Γ (1) V . Therefore, the contribution from this diagram remains zero.
Point IV is even simpler, because by construction we only have a unique way to construct the composed diagram, namely by joining two Γ (1) V -parts and therefore, the whole diagram is not canceled.
Using these rules to determine if a given diagram contributes to Γ V or not, we conclude from case Ia that, as mentioned in the main text, diagrams with at least one component from W V do not contribute. However, one has to be careful applying case III to rule out the occurrence of a certain diagram, that can be decomposed according to III, because this decomposition is not necessarily unique and so the contribution of the whole diagram is actually not zero.

Calculation of symmetry factors in the diagrammatic solution of the Ising model
We determine here the symmetry factors of the diagrams given in the main text using the following scheme: First, label all legs of interactions by indices. Then, count the possible combinations of swaps of the two legs at an interaction and permutations of interactions that lead to a new labeled diagram.
For the second-order-diagram (42), this is just 2 coming from the fact that flipping a single vertex produces a new labeled diagram, but flipping both vertices leaves the diagram invariant.
For the ring diagram in third order (43), permuting interactions produces the same diagram ("at most" it mirrors the diagram). The same holds for the other third order diagram (44).
Swapping legs in the last diagram (44) yields a factor 2 2 , because switching all vertex legs does not yield a new labeled diagram, whereas for the ring diagram, it does, leading to the symmetry factor 2 3 .
Proceeding to fourth order, we first compute the symmetry factor S R of the ring diagram in (45). For a given interaction, we have 3 possibilities to choose another interaction not to pair it with and every vertex flip produces a new labeled diagram. This gives S R = 3 · 2 4 = 48.
For the second diagram (46), labeled by "TM" (because it is depicted in a hexagonal shape reminiscent of the elements of the game board of "Terra Mystica"), we may choose 4 2 vertices for the upper position and flip every of the four vertices: S TM = 6 · 2 4 = 96.
The diagram (47) that looks like a doubled Citroën logo (dC), whose symmetry factor is given -analogous to its third-order-counterpart -by S dC = 2 3 = 8. Including the numerical factors contained in the third-oder-cumulants into the prefactor, we end up with the following contribution to the effective action: Finally, let us determine the symmetry factors for the "non-standard" fourthorder-contributions in (48), S i for i ∈ {a, b, c}: For diagram a, which we term "glasses"diagram for obvious reasons, there are three ways to pair, say, the first vertex with any of the three others. The other pair is then fixed, too. Then, every flip of a vertex produces a new labeled diagram because it matters which vertex is attached to the four-point-cumulant. Together, this gives S a = 3 · 2 4 .
For diagram S b , flipping each of the two single vertices in the left part produces a new labeled diagram, because it matters which vertex is attached to the three-pointcumulant (2 2 ). The leaf brings in a factor 2, because also here, flipping both vertices brings in a factor 2 each and a factor 1/2 because the leaf is of second order in the interaction. We count the leaf as a different interaction type, therefore the prefactor is given by 2 3 2!·1! . Finally, there is only one way to construct diagram c and we are left with the intrinsic symmetry factors of the clusters: S c = 2 2 . Here again, Γ V acts as a special kind of interaction, which occurs twice, which leads to the prefactor 2 2 2 . Adding up the contributions of the three diagrams leads us to the following expression for the term inside the square bracket in (48): In conclusion, the diagrams in (48) add up to −4 2 4 i =j =k The term in curly braces, despite emerging from a sum of three different diagrams, looks like a second cumulant. Interestingly, the contribution (72) for the case i = k is exactly canceled by that of the ring-diagram in the case that the indices of exactly two of the opposite cumulants in the ring are equal. We see this as follows: There are two ways how exactly two indices of opposite circles in the ring can be identified (i = k i − 3m 2 j + 15m 2 i m 2 j . Taken together this leads to (49).

Acknowledgements
We are indebted to Adam Rançon who commented on an earlier version of this manuscript and especially for his hints on Georges' and Yedidia's summation convention, which allowed us to see that our result is consistent with earlier works. Furthermore, we would like to acknowledge the two anonymous reviewers for many thoughtful comments that helped us to considerably improve the manuscript. In particular, we thank them for pointing out the importance of a non-iterative algorithm, that lead to the idea to use skeleton diagrams, and for pointing us towards literature on spherical models and laser physics as a potential field of application.