Grassmannization of classical models

Applying Feynman diagrammatics to non-fermionic strongly correlated models with local constraints might seem generically impossible for two separate reasons: (i) the necessity to have a Gaussian (non-interacting) limit on top of which the perturbative diagrammatic expansion is generated by Wick's theorem, and (ii) the Dyson's collapse argument implying that the expansion in powers of coupling constant is divergent. We show that for arbitrary classical lattice models both problems can be solved/circumvented by reformulating the high-temperature expansion (more generally, any discrete representation of the model) in terms of Grassmann integrals. Discrete variables residing on either links, plaquettes, or sites of the lattice are associated with the Grassmann variables in such a way that the partition function (and correlations) of the original system and its Grassmann-field counterpart are identical. The expansion of the latter around its Gaussian point generates Feynman diagrams. A proof-of-principle implementation is presented for the classical 2D Ising model. Our work paves the way for studying lattice gauge theories by treating bosonic and fermionic degrees of freedom on equal footing.


I. INTRODUCTION
Feynman diagrammatic technique is a powerful tool of statistical mechanics. Among the hallmarks of the method are the ability to deal-both analytically and numerically-with the thermodynamic limit rather than a finite-size cluster, the possibility of partial summations up to infinite order, and the fully self-consistent formulation in terms of renormalized (dressed) quantities. The latter properties allow one to go beyond the Taylor expansion in terms of the coupling constant or any other parameter.
Advantages of the diagrammatic technique come at a price. The most serious issue is the divergence of expansions in powers of the coupling constant for systems prone to Dyson's collapse 1 (i.e., pathological system behavior when the coupling constant is rotated in the complex plane). For partial summation techniques to work, the non-perturbed part of the theory has to be Gaussian (in terms of either real, or complex, or Grassmann variables) to ensure the validity of Wick's theorem. These issues are often related: for example, Ising and XY mod-els formulated in terms of original spin variables do not suffer from Dyson's collapse but lack the Gaussian (noninteracting) limit, while their classical (lattice) field counterparts with the well-defined Gaussian limit are subject to Dyson's collapse. It would be a mistake, however, to think that meaningful diagrammatic series are only possible for a very limited class of Hamiltonians, namely, when the original system is that of interacting lattice fermions. As already clearly explained by Samuel in a series of papers, 2-4 a broad class of classical spin and dimer models can be reformulated in terms of familiar interacting fermions and studied with field-theoretical techniques. Similarly, rather arbitrary quantum spin/boson lattice models can be rigorously mapped onto fermionic field theories. [5][6][7] As expected, grassmannian formulations of spin/link/boson models with local constraints are generically strongly-coupled theories at low temperature, and even the most advanced self-consistent treatments based on the lowest-order graphs are not supposed to provide quantitatively (and often qualitatively) accurate answers. Moreover, these theories may contain arbitrary multi-particle interaction vertexes, which further complicate the structure of the diagrammatic expansion. One of the promising numerical techniques currently under development for strongly correlated systems is diagrammatic Monte Carlo (DiagMC). It is based on the stochastic evaluation of irreducible Feynman graphs up to some high order and can be implemented in a number of ways, from perturbative expansions in powers of the coupling constant to various self-consistent skeleton schemes based on fully renormalized one-or two-body propagators. In such contexts as resonant fermions, 8 frustrated magnetism, 9,10 and out-of-equilibrium impurity-like models 11,12 the method was recently shown to be able to go significantly beyond the state of the art. Also, significant progress has been made in understanding superfluid properties of the Hubbard-type models. [13][14][15] Notably, the infamous signproblem preventing conventional Monte Carlo methods from simulating fermionic system with sizes large enough for reliable extrapolation to the thermodynamic limit, is absent as such in DiagMC. Instead, the computational complexity is now linked to the number of diagrams growing factorially with their order.
Nevertheless, millions of diagrams can be accounted for and the approach is flexible enough to deal with an arbitrary interaction Hamiltonian/action.
The current paradigm for generic lattice gauge models, as they occur in lattice-QCD as well as in solid state and ultra-cold atomic physics, is to work with finite-size systems and to treat link variables separately from the fermionic sector. More precisely, link variables are simulated using classical Monte Carlo techniques (with local updates), and fermions (quarks) are described by determinants. This approach suffers from a severe signproblem for finite density of fermions (non-zero chemical potential). 16,17 If link variables are straightforwardly represented by bosonic fields, then the thermodynamic limit can be addressed within the diagrammatic approach that treats bosonic and fermionic degrees of freedom on equal footing. However, in this formulation the bosonic fields pose a fundamental problem, which manifests itself in a zero convergence radius. It is thus desirable to have a generic scheme for replacing link variables with Grassmann fields to ensure that the diagrammatic expansion has proper analytic properties around the Gaussian point.
In this paper, we introduce a general procedure of grassmannization for classical lattice models. It is by no means a unique one, and in certain specific cases more compact/simpler representations can be found. There is a strong connection to the anti-commuting variables approach introduced by S. Samuel, 2-4 which can solve the 2D Ising model exactly (free fermion operators to solve the Ising model exactly were first found by Kaufman 18 and refined by Schultz, Mattis and Lieb 19 ) and provides a good starting point for field-theoretic studies of the 3D Ising model. For the latter system our approach amounts to an alternative but equally complicated field theory.
Our prime goal is to build on these ideas and develop a scheme that is flexible enough to apply to a broader class of link models with arbitrary multi-bond interactions and local constraints.
The idea of grassmannization is to represent the partition function of the model as a Grassmann integral from the exponential of a Grassmann functional. The Feynman rules then emerge by Taylor-expanding the non-Gaussian part of the exponential and applying Wick's theorem to the Gaussian averages. Paradigmatic lattice systems are link and plaquette models featuring discrete degrees of freedom-integer numbers-residing on links (plaquettes) of square lattices and subject to certain local constraints in terms of the allowed values of the sum of all link (plaquette) variables adjacent to a given site (edge). It turns out that it is these constraints that require special tricks involving multiple Grassmann variables for each value of each discrete variable. Link models often emerge as high-temperature expansions of lattice systems 20 in Ising, XY, O(3), etc. universality classes no matter whether the original degrees of freedom are discrete or continuous (e.g., classical vector-field variables). Link models may also emerge as dual (low-temperature) expansions, and specific examples are provided by the 2D Ising model 21 and the 3D |ψ| 4 model (the latter case leads to the so-called J-current model with long-range interactions). Similarly, plaquette models emerge as a high-temperature expansion of lattice gauge theories, but sometimes they represent the dual (low-temperature) expansion, as in the case of the 3D Ising model. Finally, it is worth mentioning how the models with the same general structure are generated by strong-coupling expansions in lattice-QCD. 22 The paper is structured as follows. In Sec. II we explain how a partition function of a discrete link model can be written as a Grassmann integral. The equivalence between the two formulations is readily proved through term-by-term comparison. Standard properties of Grassmann variables then immediately allow one to express the Grassmann weight in the exponential form in order to define the field-theory. In Sec. III we discuss generalizations of the proposed grassmannization scheme. We start by describing the procedure for a broad class of plaquette models. Next we show a simple way to introduce Grassmann variables for non-local link models with pairwise interactions between the link variables. The construction is further simplified when constraints are replaced with statistical penalties for certain configurations of link (plaquette) variables. We conclude this section with defining the meaning of the term "order of expansion" for the resulting field theory. In Sec. IV we deliberately choose the most general grassmannization scheme for the 2D Ising model to illustrate and test how our construction works in practice. We stress that our goal is not to solve the 2D Ising model exactly 21,23 or determine a series expansion for it 20 but to develop a general framework-including numeric component-for applying Grassmann variables to link and plaquette models and show that its evaluation can be done realistically. After determining all field-theoretic parameters, characterizing various interaction terms and source operators for calculating correlation functions (with and without magnetic field), and explaining Feynman rules for constructing the perturbative expansion, we proceed with the description of algorithms to compute them (Monte Carlo and deterministic) in Sec. V. Results are presented and discussed in Sec. VI. By comparing with the exact solution we show that the critical exponent γ for magnetic susceptibility could be determined with an accuracy of about 5%, while the critical point could be located with sub-percent accuracy. In Sec. VII we discuss the implementation of the self-consistent skeleton technique within the so-called G 2 W -expansion 24 which computes irreducible (skeleton) diagrams for the self-energy and "polarization" function and uses them in the Dyson equations in order to find the renormalized propagators and screened interactions. We also present results that emerge when this technique is based on several low-order diagrams. We briefly comment in Sec. VIII that both bare-series and the G 2 Wexpansion methods readily solve the 1D Ising model exactly. We conclude with prospects for future work in Sec. IX.

A. Local link models
For the purposes of this article, we mean by a link model a classical statistical model with states labeled by a set of discrete variables {α b } residing on links (bonds) of a certain lattice. In addition, we require that the ground state is unique. Without loss of generality, it can be chosen to be the state with α b = 0 on each link b.
We further narrow the class of link models-to which we will refer to as local link models-by the requirement that the statistical weight of a state factors into a product of link and site weights (to be referred to as link and site factors, respectively). A link factor, f b , is a function of the corresponding link variable, f b ≡ f (α b ). The site factor, g j , is a function that depends on all variables residing on links attached to the site j, denoted as {α b } j . Then, g j ≡ g({α b } j ). Solely for the purpose of avoiding heavy notations, we consider translational invariance when f (α b ) ≡ f b is the same function on all links and g j site independent, g j ≡ g. Given that only the relative weights of the states matter, we set f (0) = 1 and g(0 j ) = 1, where 0 j stands for the {α b = 0} j set.
The site factors play the key role in link models. They describe interactions between (otherwise independent) link degrees of freedom. In particular, this interaction can take the extreme form of a constraint on the allowed physical configurations of {α b } j (e.g., the zerodivergency constraint in J-current models, 25 or the evennumber constraint in the high-temperature expansion of

B. Grassmannization
For each label α = 0 of the link b, introduce four Grassmann variables: ξ α,b , ξ α,b ,ξ α,b , andξ α,b . For a textbook introduction to Grassmann variables, we refer to Ref. 26. For α = 0 we assume that ξ 0,b = ξ 0,b =ξ 0,b =ξ 0,b = 1. In terms of these variables, define the Grassmann weight-a product of link, A b , and site, B j , factors such that tracing over all degrees of freedom yields the partition function Z = Tr A b B j -by the following rules, Here {b} j stands for the set of all links incident to the site j, and variablesξ α b ,b andξ * α b ,b are defined differently for different links. We first introduce the notion of direction (on each link) so that one of the two link ends becomes "incoming" and its counterpart "outgoing" (with respect to the site adjacent to the end). Next, we assign (see Fig. 1 for an illustration) The claim is that the Grassmann integral of the weight over all variables reproduces the partition function of the original link model. For a link b to yield a non-zero contribution to the integral the link labels in (2) for the sites of the incoming (j = 1) and outgoing (j = 2) ends of the link should match each other: α 1 = α 2 . Indeed, at α 1 = α 2 , it is not possible to find an appropriate term in the expansion of the link exponential (1) such that-upon multiplying by the site factorsξ α1,bξ ξ α2,bξ * α2,b -all powers of the Grassmann variables ξ α1,b , ξ α1,b ,ξ α1,b ,ξ α1,b , ξ α2,b , ξ α2,b ,ξ α2,b ,ξ α2,b are exactly equal to 1 to ensure that the Grassmann integral is non-zero. For α 1 = α 2 ≡ α, we need to consider two cases: α = 0 and α = 0. In the first case, the non-zero contribution to the integral comes from the product of second terms in the expansion of the link exponentials (1): where we defined f * in the last step. In the second case, the two end sites contribute the factor ξ α,b ξ α,bξ α,bξ α,b = ξ α,b ξ α,bξ α,b ξ α,b . Now we have to consider the first term in the expansion of the link exponential for state α, while for other variables the calculation is repeated as in (4) γ =0 We see that, apart from the irrelevant global factor b 1/f * , we reproduce the configuration space and weight factors of the original link model.

C. Field-theoretical formulation
To generate the Feynman diagrammatic expansion, we need to represent the Grassmann weight factor in the exponential form. The link factors (1) have the form of Gaussian exponentials already. Hence, it is only the site factors that need to be rewritten identically as The constants λ({α b } j ) are readily related to the site factors g({α b } j ) by simple algebraic equations obtained by expanding the exponential and equating sums of similar terms to their counterparts in the r.h.s. of Eq. (2). By expanding the non-Gaussian part of the exponential (6) and applying Wick's theorem, we arrive at Feynman rules for the diagrammatic series. The reader should avoid confusion by thinking that an expansion of the exponential (6) takes us back to Eq. (2). Recall that connected Feynman diagrams are formulated for the free energy density, not the partition function, and summation over all lattice sites is done for a given set of interaction vertexes in the graph, as opposite to the summation over all vertex types for a given set of lattice points. Therefore, the "coupling constants" in Feynman diagrams are λ's, not g's.

D. Absorbing link factors into site factors
The separation of the weight factors into link and site ones is merely a convention. Indeed, each link factor can be ascribed to one of the two site factors at its ends. This leads to a slightly different Grassmannization protocol. This trick may prove convenient for generalization to non-local models considered below.

A. Plaquette models
A plaquette model can be viewed as a certain generalization of the local link model. States (configurations) of a plaquette model are indexed by a set of discrete labels residing on (oriented) plaquettes of a hyper-cubic lattice. The plaquette label α takes on either a finite or countably infinite number of values. The statistical weight of each state factors into a product of plaquette and edge weights (to be referred to as plaquette and edge factors, respectively). A plaquette factor, f , is a function of the corresponding plaquette variable, f ≡ f (α). An edge factor, g, is a function which depends on the labels of all plaquettes sharing this edge (this set of labels will be denoted as {α p } j for the edge j); it encodes, if necessary, constraints on the allowed sets of {α p } j .
Without loss of generality (up to a global normalization factor), we identify the "ground state" as α p = 0 for all plaquettes, and set f (0) = 1. The orientation of the plaquette (for some models it is merely a matter of convenience) is enforced by an ordered enumeration of sites at its boundary. For a plaquette p, the vertex label ν ≡ ν p = 0, 1, 2, 3 enumerates four vertices in such a way that ν ± 1 modulo 4 stands for the next/previous vertex with respect to the vertex ν in the clockwise direction.
For each state α = 0 of the plaquette p, we introduce eight Grassmann variables: ξ α,p,νp ,ξ α,p,νp , ν p = 0, 1, 2, 3. As before, for α = 0 the variables ξ andξ are not Grassmannian, ξ 0,p,ν = 0,ξ 0,p,ν = 1. The corresponding plaquette weight in the Grassmann partition function reads Note a close analogy with Eq. (1). Site weights Eq. (2) are now replaced with edge weights B j . Using the notation {p} j for the set of all plaquettes sharing the edge j, and 0 j for the state when all plaquettes in this set have α p = 0, we write where ν (j) p is the site enumeration index within the plaquette p, with respect to which the edge j is outgoing.
[Accordingly, the edge j is incoming with respect to site (ν In what follows, we will associate ν (j) p not only with the site, but also with the corresponding edge.
The proof that the classical and Grassmannian partition functions are identical (up to a global factor) is similar to the one for the link model after we notice that a non-zero contribution from plaquette p is possible only if the same plaquette label α p is used in all edge weights. The α = 0 contribution comes from the term in the expansion of the exponential (7). It contributes a factor 1/q * , where The α = 0 contribution comes from the plaquette term multiplied by the product 3 νp=0ξ α,p,νp ξ α,p,νp originating from the boundary edge terms ξ α,p,(ν (j) Because of the Grassmann anticommutation rules, this four-edge factor yields an additional minus sign, explaining the use of the negative sign in front of f (α) in Eq. (7). Upon Grassmann integration, the contribution to the partition function of the resulting term equals to f (α)/q * .
Feynman diagrammatics for the plaquette model is obtained by following the same basic steps as for the link models. The Gaussian part is given by Eq. (7) with four pairs of Grassmann fields for every non-zero plaquette state. The interaction part of the Grassmann action is contained in edge weights (8) after they are written in an exponential form with the constants λ({α b } j ) unambiguously related to the edge factors g({α b } j ).

B. Unconstrained discrete models with pair-wise interaction
The hallmark of the considered link (plaquette) models is the non-trivial interaction introduced via site (edge) factors. It is due to this type of interaction-and, in particular, its extreme form of a constraint on allowed combinations of discrete variables-that we had to introduce multiple Grassmann variables for each state of the link (plaquette). The situation simplifies dramatically if we are dealing with unconstrained discrete degrees of freedom with pair interactions between them. Consider a link model defined by the statistical weight based on products of two-link factors. Without loss of generality, these factors can be cast into the exponential form We assume that all factors in the product are bounded and properties of the η-matrix are well-conditioned. Grassmannization of this model can be done by taking advantage of properties of Gaussian integrals that allow one to express (14) identically (up to normalization) as Here {X α b ,b } is a collection of auxiliary real continuous variables. For briefness, we do not show explicitly the Gaussian weight W G that is uniquely defined by the values of all pairwise averages performed with this weight What we achieve for a fixed set of X variables is a link model that contains only single-link factors ∀b : For models with site constraints, link factors can be attributed to site factors at the incoming (or outgoing) ends with subsequent Grassmannization of the latter as discussed above. For unconstrained models, Grassmannization is accomplished by replacing sums over link variables with Note that here Grassmann variables have nothing to do with the discrete index α b , in contrast with previous considerations. The resulting formulation contains both Grassmann and real-number integrations. Clearly, all considerations can be repeated identically (up to a trivial change in notations) for a model based on discrete variables α s residing on lattice sites when the configuration weight is given by C. Order of expansion The notion of the order of expansion is absolutely central for practical applications when diagrammatic series are truncated. Normally, it is defined as an integer nonnegative power of a certain dimensionless parameter ζ playing the role of a generalized coupling constant, such that the diagrammatic expansion corresponds to a Taylor expansion in ζ about the point ζ = 0. Without loss of generality, we can always select ζ (by an appropriate rescaling) in such a way that the physical value of ζ is 1. This is especially convenient in cases when there is more than one interaction vertex, and ascribing different powers of ζ to them results in (re-)grouping of different terms in the series. A reasonable guiding principle behind such a (re-)grouping is the requirement to end up with Taylor series having finite convergence radius around ζ = 0. The latter is guaranteed if the theory is analytic in ζ at the origin; the necessary condition for this to be true is the absence of Dyson's collapse when changing the sign (more generally, the phase) of ζ.
Introduce the ζdependence by the replacement In terms of the original theory, the replacement (20) means η → ζ 2 η, for all η's in Eq. (14). If amplitudes of all η values in (14) are bounded, we expect that such a dependence on ζ is analytic not only for a finite system, but also in the thermodynamic limit at finite temperature. In the Grassmann action (18), the expansion of the exponential e iζXs,α b in powers of ζ generates an infinite series of interaction vertexes (the zeroth-order term defines the harmonic action): Higher-order vertexes in X come with a higher power of ζ and this sets unambiguously the rules for defining the diagram order.  (24). The contribution of these diagrams is ζ + 2ζ 3 .

A. Model and observables
Consider the 2D Ising model on the square lattice with the Hamiltonian The Ising variables σ = ±1 live on the sites of the 2D square lattice and interact ferromagnetically with their nearest neighbors, as is represented by the first term in the Hamiltonian. We write the dimensionless coupling as β in units of the temperature T . Additionally, every spin feels a dimensionless magnetic field h i = h, which can be taken h ≥ 0 without loss of generality. The partition function of the Ising model reads The most typical observable of the Ising model is the spin-spin correlation function ρ ij , B. Grassmannization of the high-temperature expansion Using the well-known identities the partition function can be written as Z = Z 0 Z with Z 0 = (cosh β) 2N (cosh h) N for a lattice of N sites and 2N links. With the notation ζ = tanh β and η = tanh h the remaining factor is given by Upon summation over spin variables we are left with a link model, where link variables take only two values, 0 or 1, to specify whether we are dealing with the first or the second term in the sum (1 + σ i σ j ζ). In the partition function, terms with an odd power of σ i on any of the sites yield zero upon spin summation. The remaining terms depend on link variables in a unique way. The formalism of the previous section can be straightforwardly applied, and we obtain Here we label site factors using the total sum of incident link variables, b∈{b}j α b , to avoid unnecessary rank-4 tensor notations. If we further redefine Z 0 → Z 0 2 N f 2N * , then the Grassmann representation of the partition function Z is given by

C. Vertex coefficients
We now compute the factors λ. To this end, we first introduce notations (for a fixed site j and suppressing the site index for clarity) and then Taylor expand The only non-zero terms generated by this expansion are and V 2 2 = 6V 4 . All other powers and multiplications of operators yield zero. Note that operators from different sites commute and may be excluded from consideration here. The final result is Term-by-term matching with Eq. (28) then leads to In what follows we will discuss the η = 0 case (zero external field) when the only vertexes with non-zero coupling in the partition function are V 2 and V 4 , The expansion of Z in powers of ζ then goes as and the spin-spin correlation function is given by

D. Feynman rules
In order to arrive at the Feynman perturbative expansion we need to write the partition function in the form where Z is the partition function of the Gaussian part (it is the product of local link contributions), 3. For links with multiple occupancy, a minus sign occurs when swapping 2 Grassmann variables. The minus sign can also be found by counting all closed fermionic loops.
4. The total weight of the diagram in order n is hence (−1) P (−2) q ζ n with P the signature of the exchange permutation and q the sum of all type-3 and type-4 vertexes.
Disconnected diagrams are defined with respect to both the primed and non-primed Grassmann variables simultaneously. Thus, a link can lead to a disconnected diagram only if the primed and non-primed variables simultaneously lead to disconnected pieces (such as the upper left panel in Fig. 8). We check the connectivity of a diagram by the breadth-first algorithm.

E. Example: the first element of the spin correlation function
Let us focus on the first element of the correlation function connecting the sites (0, 0) and (1, 0) (using translational invariance, any 2 neighboring sites r 1 , r 2 can be taken). To first order, we put a V 1 vertex on the origin and target site. There is one way to combine them, thus the total contribution is ζ. By the symmetry of the lattice, even expansion orders do not contribute. In third order, we can construct a diagram by putting a V 2 (RD) vertex on the site (0, 1) and a V 2 vertex (LD) on the site (1, 1). The mirror symmetry of this diagram about the x-axis is also a valid diagram. Hence, the contribution is 2ζ 3 . These diagrams contributing in first and third order are shown in Fig. 3.
In fifth order, there are 4 diagrams with a V 3 vertex on one of the endpoints, yielding a contribution −8ζ 5 . There are 14 diagrams consisting of only V 1 and V 2 vertexes and single pair-lines, yielding a contribution 14ζ 5 . The contributions to fifth order are shown in Figs. 4, 5, 7, and 8. There are however additional diagrams with 2 pairs of Grassmann variables living on the same link, as is shown in Fig. 8 (there are equivalent diagrams obtained by mirror symmetry around the x-axis which are not shown). They all have on the origin a V 1 and a V 2 (RU ) vertex, and on the target site (1, 0) a V 1 and a V 2 (U L) vertex. On the site (1, 1) there is a V 2 (LD) and on site (0, 1) a V 2 (RD) vertex. Let us look more carefully at the link between the origin and target site: The origin is associated withξ ξ and the target with ξξ by our convention. Applying Wick's theorem, there are 4 possible ways to pair the Grassmann variables: 1. The pairing combinationξ ξ ξ ξ ξ ξ ξ ξ comes with the sign +1 and leads to a connected diagram (this is the lower right panel in Fig. 8).
2. The pairing combinationξ ξ ξξ ξ ξ ξ ξ comes with the sign -1 and leads to a connected diagram (this is the upper right panel in Fig. 8) 3. The pairing combinationξ ξ ξ ξ ξ ξ ξ ξ comes with the sign -1 and leads to a connected diagram (this is the lower left panel in Fig. 8) 4. The pairing combinationξ ξ ξ ξ ξ ξ ξ ξ leads to a disconnected diagram and does not contribute to the correlation function (this is the upper left panel in Fig. 8).
The net contribution of these 4 distinct diagrams is hence −1 (also the diagrams obtained by mirror symmetry around the x-axis yield −1, so the total contribution to fifth order is (−8 + 14 − 2)ζ 5 = 4ζ 5 . It is instructive to notice that the sum of all diagrams in which multiple Grassmann pairs live on the same link always produces zero in case all diagrams are connected, in line with the nilpotency of Grassmann variables. Wick's theorem splits however these contributions in connected and disconnected diagrams, where the disconnected diagrams cancel against the denominator of the Feynman expansion. It is this non-trivial regrouping imposed by Wick's theorem that can yield non-zero contributions from terms like (45); and, in particular, from arbitrarily high powers of one and the same interaction vertex.

V. IMPLEMENTATION
We explored two ways of evaluating the (bare) series for the spin-correlator: a stochastic Monte Carlo approach and a deterministic full evaluation of all diagrams.

A. Monte Carlo sampling
In order to perform a Monte Carlo sampling over all Feynman diagrams, we introduce a head and a tail that represent the endpoints of the correlation function. By moving them around the lattice and changing the diagrammatic elements in between the head and tail, we are able to reach an ergodic sampling. The algorithm can be formulated as follows: The tail remains stationary at the origin whereas the head can move around the lattice. When the head and tail are on the same site and the expansion order is 0, the value of the correlation function is 1 which can be used for normalization of the Monte Carlo process. A Monte Carlo measurement contributes +1 or −1 depending on the sign of the diagram weight. The simplest Monte Carlo procedure samples according to the absolute weights of the diagrams and consists of the following pairs of reciprocal updates: 1. MOVE-RETRACT. We choose one of the 4 directions randomly, and attempt to place the head on the site adjacent to the current head site according to this direction. In case this direction does not correspond to backtracking, the current V 1 type of the tail turns into a V 2 , otherwise the head goes back and changes the previous V 2 into a V 1 type (unless the diagram order is 0 or 1, when only V 1 types are possible). When moving forward, the way of pairing primed and non-primed variables is always unique which in turns implies that we can only retract when the head is connected via a "straight pair connection" to the previous vertex (both primed and non-primed Grassmann variables of the head are connected to the same vertex on the previous site). We only allow the MOVE-RETRACT updates if the end vertex types are V 1 .
2. SWAP VERTEX. Swaps between the vertexes V 1 + V 2 ↔ V 3 (for head and/or tail) and V 2 + V 2 ↔ V 4 (anywhere in the diagram). This update is its own reciprocal.
3. RELINK. On a given link, relink primed and nonprimed Grassmann variables. This can change the sign of the weight only. This update is its own reciprocal.
The second and third type of updates may lead to disconnected diagrams. In such cases, the configuration is unphysical. We opt to allow such configurations, but a Monte Carlo measurement is forbidden and type-1 updates remain impossible until the diagram is connected again. For small values of ζ the sign problem is nearly absent, but only low expansion orders can be reached. For higher values of ζ (close to and above the critical one) an increasing number of orders contributes significantly, consequently more time is spent in higher orders and the sign problem significantly worsens.

B. Deterministic full evaluation
For the case of the 2D Ising model, a Monte Carlo approach offers no advantages over a full series expansion approach. With this we mean the explicit listing and evaluation of all possible diagrams as opposed to the stochastic sampling over all topologies. This is because all diagrams in a given expansion order contribute a number of order unity (times the same power of ζ), often with alternating sign, leading to huge cancellations. Only the exact cancellation has physical information, and this requires that every diagram is evaluated multiple times before the correct convergence can be seen. A Monte Carlo approach makes much more sense if the dominant contributions to the total weight are coming from a narrow parameter region, which is usually the case if there are additional integrals over internal momenta.
We therefore wrote a code that evaluates all diagrams for the correlation function up to a maximum order. The construction is based on the fact that there is an easy way to construct all the "easy" diagrams (the ones that formally look like originating from a high-temperature series expansion). These can serve as parent diagrams, from which further offspring diagrams can be constructed which have one or multiple V 3 and V 4 vertexes as well as possible fermionic exchanges. All diagrams in order n can be found as follows: 1. Write down all possible words of the form X 1 X 2 . . . X n with the alphabet X j ∈ {0, 3} corresponding to the 4 directions on the square lattice. Make sure that subsequent directions are not backtracking. For example, if X 4 is in the positive +x direction, then X 5 cannot be in the negative −x direction. From this word we also know all sites and links that are visited, as well as all type-1 and type-2 vertexes that are used to make this diagram.
2. Such a parent diagram is added to a list of different topologies only if it has a unique topology. To store the topological information of a bare vertex, we need to store a pair consisting of a site index and a vertex type. The diagram is then stored as an ordered map where the "key" values are given first by the lattice site index and second by the vertex type (in binary format). The ordered map may have multiple entries with the same key if multiple vertexes reside on the same site and if they are of the same type (e.g., two RL vertexes on the same site).
3. We iterate over this configuration list and check if the tail and head sites can be merged into a type-3 vertex by combining them with type-2 vertexes that reside on the same lattice site. If so, and if the resulting topology is unique, the diagram is added to the list. This step is performed in three parts: first for the head and tail together (in order to find all diagrams with 2 V 3 ends), then for the head alone, and finally for the tail alone.
4. We iterate again over the full configuration list and check if 2 type-2 vertexes that live on the same site can be merged into a type-4 vertex. This last step has to be repeated until no further merges are possible (since it may happen that a diagram has multiple type-4 vertexes or even multiple type-4 vertexes on the same site). Diagrams thus created are also added the configuration list if their topology is unique. After completion of this step, all possible topologies have been generated.

5.
We compute the product of all the vertex weights, according to the Feynman rules.
6. From this list of parent diagrams we need to generate all offspring diagrams which feature all possible fermionic permutations for multiply occupied links. This first requires that we know how the vertexes are connected in the parent diagram, which is stored in the configuration list. The parent diagram always has permutation sign +1 (because the connections of the primed and non-primed Grassmann variables are always the same). Next we generate all possible permutations by relinking the primed and/or the non-primed Grassmann variables using Heap's algorithm. If a link has occupation number m, then there are (m!) 2 combinations to be generated (and there may be more than one multiply occupied link). The permutation signature is also stored.
7. We check the connectivity of the diagram using the breadth-first algorithm. Disconnected diagrams contribute 0.
8. Finally, we compute the isomorphism factor: if m identical vertexes on the same site are found, a factor 1/m! must be taken into account. This is a consequence of how we construct the diagrams: topology checks were only performed on the parent diagrams (and based on vertexes only), not on offsprings obtained by fermionic exchange. (It would be prohibitively expensive to add the offspring diagrams to the list of all possible diagrams.) Hence, just as we generate illegal disconnected diagrams, we also have a double counting problem when identical vertexes occur in the list.
In order 14, there were about 140,000 parent diagrams contributing to the first entry on the diagonal of the correlator. The hugest number of permutations was (4!) 4 (3!) 4 ≈ 10 8 . Since the sum of these permutations has a net contribution of order 1, Monte Carlo has roughly a sign problem of the order of 10 −8 for these diagrams. The first time a nontrivial isomorphism factor is seen is in order 6 for the first element on the diagonal of the spin correlator: There are diagrams in which two links are doubly occupied, and those links are connected by an identical V 2 vertex, hence the isomorphism factor 1/2. More efficient ways of evaluating and storing the diagrams can probably be devised and implemented, but the above scheme is sufficient to check the validity of the technique and study the transition.

A. Spin-spin correlation function
Our results for the spin-spin correlation function are shown in Table. I. The correlation function is known recursively from Refs. 27-29. It is also known as a Painlevé-VI nonlinear differential equation 30 but this is not so well suited to obtain the series coefficients. Along the principal axes and the diagonal it can also be expressed as a Toeplitz determinant. The first element along and the axis and the diagonal can be recast in terms of complete elliptic integrals (see pp. 200-201 in Ref. 21), which are convenient for series expansions, → ζ + 2ζ 3 + 4ζ 5 + 12ζ 7 + 42ζ 9 + . . .
with k > = sinh 2 (2β), K(.) and E(.) the complete elliptic K and E functions, respectively. The above-cited recursion relations could be initialized with these expansions and shown to yield the same results as the top 2 rows in Table. I.
site/order ζ ζ 2 ζ 3 ζ 4 ζ 5 ζ 6 ζ 7 ζ 8 ζ 9 ζ 10 ζ 11  The series is convergent for any finite expansion order, i.e., in the thermodynamic limit the infinite series will diverge first at the phase transition point. It is hence possible to study the critical behavior of the susceptibility, which is governed by the critical exponent γ = 7/4. We plot in Fig. 9 the susceptibility versus β for different expansion orders, and also plot the asymptotic behavior for comparison. The critical temperature and the exponent γ can be found from a study of the convergence radius of the series. the ratio of coefficients asympotically behaves as In Fig. 10 we extract the critical point ζ c from the intercept and the critical exponent γ from a linear fit through the ratio of the coefficients. The critical point could be determined with an accuracy of 0.5%, whereas the error on γ is of the order of 5%. However, according to more advanced extrapolation techniques discussed in Ref. 34, γ can be determined independently from ζ c as γ ≈ 1.751949 on the square lattice when the series is known up to 14th order, i.e., an accuracy of 0.5%.

VII. THE G 2 W SKELETON SCHEME
The expansion of susceptibility in terms of ζ is, of course, identical to the one found by the hightemperature series expansion method. To make the distinction between the high-temperature series formalism and Grassmannization approach clear, we discuss the skeleton formulation of the interacting fermionic fieldtheory based on dressed (or "bold") one-body propagators (G) and bold interaction lines (W ). This leads to the so-called G 2 W skeleton scheme (see for instance Refs. 24 and 35 for the terminology): all lines in all diagrams are assumed to be fully renormalized propagators and effective potentials, but vertex functions remain bare. In Sec. VIII we show that the G 2 W -expansion scheme offers a very simple way to solve the 1D Ising model exactly.

A. Objects and notation
The key objects in the standard skeleton scheme are the selfenergy (Σ) and the polarization function (Π). They are related to the Green function (G) and the effective potential (W ) by their respective Dyson equations. The diagrams for Π and Σ are obtained by removing one W -or G-line, respectively, from connected graphs for the generalized Luttinger-Ward functional Ψ, shown to second order in Fig. 11. In this setup, the expansion order is defined by the number of W -lines (obviously, the discussion of Sec. III C does not apply to the self-consistent skeleton sequence). All objects of interest are tensors; they have a coordinate (or momentum) dependence, as well as the legs orientation dependence for the incoming and outgoing parts. This conventional scheme has to be supplemented with Ψ-graphs involving V 4 vertexes to account for all contributions. We start with neglecting V 4 vertexes, and discuss their role later.
In more detail, the formalism of the G 2 W expansion in the absence of V 4 vertexes is as follows: 1. There are six bare two-body interaction vertexes V 2 (RU, RL, RD, LU, U D, LD), see the second line in Fig. 2. They reside on the sites of the original square lattice and all have weight 1. Symbolically, we encode the tensor structure of V 2 using a convenient short hand notation The row index represents the first leg enumerated according to the convention (R, U, L, D) → (0, 1, 2, 3), and the column index represents the second leg. By doing so, we artificially double the number of vertexes from 6 to 12. For example, the element (0, 2) corresponds to n L n R whereas (2, 0) corresponds to n R n L , which is exactly the same term.
2. The selfenergies Σ for the primed and non-primed Grassmann variables take the same value. Thus, we have to compute only one of them and we can suppress the index that distinguishes between the two Grassmann fields. The selfenergy defines the Green function through the Dyson equation For a link going from site i to site j, the first index α refers to site i (in the above-defined sense), and the second index γ refers to site j. Note the absence of the momentum dependence in Eq. (52): The bold Green function remains local on the links in any order of renormalization. It means, in particular, that the only non-zero element for a link between sites (0, 0) and (1, 0) is G 02 ; it can be alternatively denoted as G x and, by 90 o rotation symmetry of the square lattice, is the same for all links.
3. The matrix structure of polarization Π is similar to that of V . The 0th order expression based on bare Green functions is given by  (53)

The effective potential W is defined through the Dyson equation in momentum representation
We expect to see signatures of the ferromagnetic transition in matrix elements of W q=0 because they directly relate to the divergent uniform susceptibility χ.

B. Zeroth order result
To obtain the 0th order result, we replace Π with Π (0) in Eq. (54). For any ζ we compute W q=0 from Eq. (54) by matrix inversion. We find a divergence at ζ c = 1/3 (shown in Fig. 12) that can be also established analytically. We see that W diverges as (ζ c − ζ) −1 . We get the same power law behavior for the (0, 1) matrix element as well as for the Frobenius norm-they just differ by a constant factor. It is not surprising that our ζ c is below the exact value for this model; the skeleton approach at 0th order is based exclusively on simple "bubble" diagrams in terms of bare Green functions that are all positive, leaving to an overestimate of the critical temperature. Fermionic exchange cycles and vertexes with negative weights do not contribute at this level of approximation.

C. First order result
We now include the diagrams with one W line for the selfenergy and the polarization. In real space we find The matrix structure of Π (1) is identical to that of Π (0) and is not shown here explicitly. Coupled self-consistent Eqs. (55), (52), and (54) are solved by iterations.

D. Second order result
As mentioned previously, to account for second-order terms in Σ, one goes to the second order graphs for Ψ and removes a G line, whereas the second-order terms for Π are obtained by removing one W -line from the thirdorder graphs for Ψ. The corresponding expressions in real space are The remaining non-zero contributions are obtained by invoking discrete lattice symmetries. Note that to this order the polarization function is extremely local and contains only same site and n.n. terms. Again, coupled self-consistent GW-equations are solved by fixed-point iterations. The resulting behavior for W is analyzed in Fig. 13. The transition point has slightly shifted to larger values of ζ compared to the zeroth-order result, and the exponent has also slightly increased.

E. Relating Π to the spin correlation function
The G 2 W -expansion scheme treats different bare vertexes (see Fig. 2) on unequal footing: the V 2 vertexes are fully dressed, but the V 4 vertexes are included perturbatively (we neglected them so far). These higher-rank vertexes have a weight of comparable magnitude to the V 2 vertexes (-2 for V 4 vs +1 for V 2 ). In addition, the difference in sign between the weights is expected to result in important cancellations between the diagrams and better convergent series for the spin correlation function (this is how ζ c increases towards its exact value).
Formally, there is no valid reason for neglecting the V 4 vertexes altogether. Let us show how they can be taken care of in the spirit of the shifted action approach. 36 This discussion also gives us the opportunity to explain how the spin correlator is related to the G 2 W skeleton expansion, which is most easily understood in the limit ζ 1. By assuming that the skeleton sequence (without V 4 ) is solved, we introduce the full polarization function Π(α,β) through the Dyson equation To be specific, we focus on the n.n. element ρ (1,0) ; similar manipulations hold for any other distance. Now consider all diagrams for this correlator without the V 4 vertexes within the G 2 W formulation (see Ref. 36): • Put one V 1 vertex on the origin site (0, 0) and the other V 1 vertex on the target site r = (1, 0), see Eq. (43). There are 4 × 4 = 16 different ways of doing that depending on the directions of legs. Connect the legs withΠ r (α, γ). For example, in the limit of ζ 1, choosing the (α = 0)-leg on site (0, 0) and the γ = 2-leg on site (1, 0) results in the contributions ζ − 4ζ 5 + . . .. Similarly, the choice of α = 1 and γ = 1 leads to the contribution ζ 3 .
• Put V 3 on (0, 0) and V 1 on (1, 0), and connect all legs withΠ lines. There are four ways to orient the V 3 vertex and for each one there are two choices for connecting legs withΠ propagators. The leading contribution to ρ (1,0) goes hence as −8ζ 5 .
• Put one V 3 vertex on (0, 0) and the other V 3 vertex on (1, 0). Now there are 16 ways of orienting both V 3 vertexes, and for each orientation there are 15 choices for connecting the legs. These contributions start at order ∝ ζ 9 .
Next, we repeat the above procedure of connecting legs by adding one V 4 vertex, which can be put on any site, after that we can add two V 4 vertexes etc. to generate a perturbative expansion in the number of V 4 terms. Compared to the original bare series in powers of ζ, we have reordered the series: the effective potential is summing up all V 2 vertexes, whereas we expand (and sample in a Monte Carlo framework) in powers of λ 4 .
To illustrate this framework, let us take ζ = 0.01 and recall that in the bare series ρ (1,0) = ζ + 2ζ 3 + 4ζ 5 + 12ζ 7 + . . .. The first 3 terms can be reproduced without V 4 vertexes and with only 1 V 3 on either the origin or the target site, see Figs. 3-8. The fifth order coefficient originates from 16 "simple" diagrams containing just V 1 and V 2 vertexes without any exchange. The diagrams containing a V 3 vertex yield a coefficient −8, and the exchange diagrams yield a coefficient −4. On a 16 × 16 lattice, the propagators obtained in Sec. VII B (i.e., to zeroth order) arē We do not mention explicitly other symmetry-related elements. The sum of all matrix elements forΠ The sum of all matrix elements forΠ (0) (x,y)=(1,0) is 0.01000200120. One clearly recognizes the coefficients 1, 2 and 12 for first, third and fifth order contributions to the bare series. For the fifth order contribution, we now obtain 12 instead of 16 thanks to the Grassmann exchange contribution that is accounted for properly at this level of approximation. By adding the V 3 diagrams in the way described above we recover the correct result to this order in ζ (which is +4).
The first instance of a V 4 vertex occurs in order ζ 6 in the bare series. The relevant bare diagrams are the ones for ρ (1,1) with a V 4 vertex on site (1, 0) (and all cases related by the lattice symmetry). Our bold expansion can correctly account for this contribution if we put a V 4 vertex on this site and connect all unpaired legs withΠ propagators. However, with the propagators obtained in Sec. VII D we are not supposed to account for all possible diagrams in the bare series to order 6 because our bold expansion in Sec. VII D is only accurate up to order ζ 3 : Consider again ρ (1,1) and the bare diagrams where exchanges are possible on the links between the sites (0, 0) − (1, 0) and (1, 0) − (1, 1). Then there are irreducible non-local contributions that are not accounted for in Sec. VII D with a positive weight that involves exchanges on both links in a correlated fashion. These contributions would obviously be accounted for in higher order corrections to Ψ, when Π becomes non-local. This is also seen in the numerics: the G 2 W approach to second order yields a coefficient of 6 for ζ 6 contribution to ρ (1,1) , which is below the correct value of 10.

VIII. THE ISING MODEL IN ONE DIMENSION
Let us show that the proposed approach solves the 1D Ising model exactly, both in the bare formulation as well as in the G 2 W skeleton formulation.

A. Bare series
In 1D, the only allowed vertex is RL (the last one in the second line of Fig. 2). It has weight +1. The only allowed endpoints are L and R (the second and fourth vertexes shown in the first line of Fig. 2). As expected, this means that there are no loops, no fermionic exchanges, and no minus signs in 1D. At order n of the expansion for the spin correlator there is only one contributing diagram with weight ζ n (up to the lattice symmetry). The susceptibility is hence T χ = 1 + 2(ζ + ζ 2 + . . .) = 1 + 2 reproducing the exact solution with asymptotic behavior χ ∝ β exp(2β) as T → 0.

IX. CONCLUSION
We have developed a general scheme for mapping a broad class of classical statistical link (plaquette) models onto interacting Grassmann-field theories that can be studied by taking full advantage of the diagrammatic technique. This mapping, in particular, would allow to formulate an all-diagrammatic approach to (d + 1)dimensional lattice gauge theories with finite density of fermions. The resulting field-theory looks very complex because it contains a large number of Grassmann variables with numerous multi-point interaction vertexes. Moreover, it is generically strongly-coupled at low temperature meaning that an accurate solution using diagrammatic methods is only possible when calculations are performed to high order and extrapolated to the infinite-order limit.
The complexity of the problem should not be taken as an indication that the entire idea is hopeless. Monte Carlo methods were designed to deal with configuration spaces of overwhelming size and complexity and arbitrary weights. In this sense, diagrammatic Monte Carlo methods simulating the configuration space of irredicuble connected Feynman graphs are based on the same general principles and one should not be surprised that they can evaluate the sum of millions of bare (or skeleton) graphs, enough to attempt an extrapolation to the infinite-order limit. What makes diagrammatic Monte Calro distinctly unique (apart from working with everchanging number of continuous variables without systematic errors) is the radical transformation of the sign problem. It is completely eliminated in conventional sense because the thermodynamic limit is taken first. Given that the number of diagrams increases factorially with their order, finite convergence radius in ζ is only possible if same-order diagrams cancel each other to such a degree that at high order their combined contribution is not increasing factorially. In other words, non-positive weights are required for the entire approach to work and we call it the "sign-blessing" phenomenon. Diagram weights for Grassmann/fermion fields alternate in sign depending on the diagram topology; this leads to the sign-blessing phenomenon for lattice models.
We illustrated the proposed approach by considering the 2D Ising model as a prototypical example. We have deliberately chosen to work with the generic formulation to avoid model specific simplifications because our goal was not to solve the model but to demonstrate how one would proceed in the general case. The ultimate goal is to explore how this field-theoretical approach can help with understanding properties of lattice gauge models.

X. ACKNOWLEDGEMENT
We are grateful to A. J. Guttmann for providing us with references to the high-temperature series expansions and U. Wolff for drawing our attention to the work