Grassmannization of classical 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 $\{{\alpha }_{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 ${\alpha }_{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}\equiv f({\alpha }_{b})$. The site factor, g_j, is a function that depends on all variables residing on links attached to the site j, denoted as ${\{{\alpha }_{b}\}}_{j}$. Then, ${g}_{j}\equiv g({\{{\alpha }_{b}\}}_{j})$. Solely for the purpose of avoiding heavy notations, we consider translational invariance when $f({\alpha }_{b})\equiv {f}_{b}$ is the same function on all links and g_j site independent, ${g}_{j}\equiv 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 ${\{{\alpha }_{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 ${\{{\alpha }_{b}\}}_{j}$ (e.g., the zero-divergency constraint in J-current models [28], or the even-number constraint in the high-temperature expansion of Z₂ models), in which case ${g}_{j}({\{{\alpha }_{b}\}}_{j})$ is identically zero for each non-physical state of ${\{{\alpha }_{b}\}}_{j}$.

For each label $\alpha \ne 0$ of the link b, introduce four Grassmann variables: ${\xi }_{\alpha ,b}$, ${\xi }_{\alpha ,b}^{\prime }$, ${\bar{\xi }}_{\alpha ,b}$, and ${\bar{\xi }}_{\alpha ,b}^{\prime }$. For a textbook introduction to Grassmann variables, we refer to [29]. For $\alpha =0$ we assume that ${\xi }_{0,b}={\xi }_{0,b}^{\prime }={\bar{\xi }}_{0,b}={\bar{\xi }}_{0,b}^{\prime }=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=\mathrm{Tr}\prod {A}_{b}\prod {B}_{j}$—by the following rules

$$\begin{eqnarray}{A}_{b} & = & \exp \left\{\displaystyle \sum _{\alpha \ne 0}\left[\displaystyle \frac{{\bar{\xi }}_{\alpha ,b}^{\prime }{\xi }_{\alpha ,b}^{\prime }}{\sqrt{f(\alpha )}}+\displaystyle \frac{{\bar{\xi }}_{\alpha ,b}{\xi }_{\alpha ,b}}{\sqrt{f(\alpha )}}\,\right]\right\}\\ & = & \displaystyle \prod _{\alpha \ne 0}\exp \left\{\displaystyle \frac{{\bar{\xi }}_{\alpha ,b}^{\prime }{\xi }_{\alpha ,b}^{\prime }}{\sqrt{f(\alpha )}}+\displaystyle \frac{{\bar{\xi }}_{\alpha ,b}{\xi }_{\alpha ,b}}{\sqrt{f(\alpha )}}\right\},\end{eqnarray} \tag{ 2.1 }$$

$$\begin{eqnarray}{B}_{j} & = & \displaystyle \sum _{{\{{\alpha }_{b}\}}_{j}}g({\{{\alpha }_{b}\}}_{j})\displaystyle \prod _{b\in \{b\}{}_{j}}{\breve{\xi }}_{{\alpha }_{b},b}\,{\breve{\xi }}_{{\alpha }_{b},b}^{* }\\ & = & 1\,+\displaystyle \sum _{{\{{\alpha }_{b}\}}_{j}\ne {0}_{j}}g({\{{\alpha }_{b}\}}_{j})\displaystyle \prod _{b\in \{b\}{}_{j}}{\breve{\xi }}_{{\alpha }_{b},b}\,{\breve{\xi }}_{{\alpha }_{b},b}^{* }.\end{eqnarray} \tag{ 2.2 }$$

Here $\{b\}{}_{j}$ stands for the set of all links incident to the site j, and variables ${\breve{\xi }}_{{\alpha }_{b},b}$ and ${\breve{\xi }}_{{\alpha }_{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 figure 1 for an illustration)

$$\begin{eqnarray}{\breve{\xi }}_{{\alpha }_{b},b} & = & {\xi }_{{\alpha }_{b},b},\,\,{\breve{\xi }}_{{\alpha }_{b},b}^{* }={\xi }_{{\alpha }_{b},b}^{\prime }\,\,{\rm{(for}}\,{\rm{incoming}}\,{\rm{end)}},\\ {\breve{\xi }}_{{\alpha }_{b},b} & = & {\bar{\xi }}_{{\alpha }_{b},b}^{\prime },\,\,{\breve{\xi }}_{{\alpha }_{b},b}^{* }={\bar{\xi }}_{{\alpha }_{b},b}\,\,{\rm{(for}}\,{\rm{outgoing}}\,{\rm{end)}}.\end{eqnarray} \tag{ 2.3 }$$

Zoom In Zoom Out Reset image size

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.2) for the sites of the incoming (j = 1) and outgoing (j = 2) ends of the link should match each other: ${\alpha }_{1}={\alpha }_{2}$. Indeed, at ${\alpha }_{1}\ne {\alpha }_{2}$, it is not possible to find an appropriate term in the expansion of the link exponential (2.1) such that—upon multiplying by the site factors ${\breve{\xi }}_{{\alpha }_{1},b}\,{\breve{\xi }}_{{\alpha }_{1},b}^{* }$ and ${\breve{\xi }}_{{\alpha }_{2},b}\,{\breve{\xi }}_{{\alpha }_{2},b}^{* }$—all powers of the Grassmann variables ${\xi }_{{\alpha }_{1},b}$, ${\xi }_{{\alpha }_{1},b}^{\prime }$, ${\bar{\xi }}_{{\alpha }_{1},b}$, ${\bar{\xi }}_{{\alpha }_{1},b}^{\prime }$, ${\xi }_{{\alpha }_{2},b}$, ${\xi }_{{\alpha }_{2},b}^{\prime }$, ${\bar{\xi }}_{{\alpha }_{2},b}$, ${\bar{\xi }}_{{\alpha }_{2},b}^{\prime }$ are exactly equal to 1 to ensure that the Grassmann integral is non-zero. For ${\alpha }_{1}={\alpha }_{2}\equiv \alpha$, we need to consider two cases: $\alpha =0$ and $\alpha \ne 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 (2.1):

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \prod _{\gamma \ne 0}\displaystyle \int { \mathcal D }{[\bar{\xi }^{\prime} \xi ^{\prime} \bar{\xi }\xi ]}_{\gamma }\exp \left\{\displaystyle \frac{{\bar{\xi }}_{\gamma }^{\prime }{\xi }_{\gamma }^{\prime }}{\sqrt{f(\gamma )}}+\displaystyle \frac{{\bar{\xi }}_{\gamma }{\xi }_{\gamma }}{\sqrt{f(\gamma )}}\right\}\\ & & \quad =\displaystyle \displaystyle \prod _{\gamma \ne 0}\int { \mathcal D }{[\bar{\xi }^{\prime} \xi ^{\prime} \bar{\xi }\xi ]}_{\gamma }\left[1+\displaystyle \frac{{\bar{\xi }}_{\gamma }^{\prime }{\xi }_{\gamma }^{\prime }}{\sqrt{f(\gamma )}}\right]\left[1+\displaystyle \frac{{\bar{\xi }}_{\gamma }{\xi }_{\gamma }}{\sqrt{f(\gamma )}}\right]\\ & & \quad =\displaystyle \prod _{\gamma \ne 0}\displaystyle \frac{1}{f(\gamma )}\equiv \displaystyle \frac{1}{{f}_{* }},\end{array}\end{eqnarray} \tag{ 2.4 }$$

where we defined f_* in the last step. In the second case, the two end sites contribute the factor ${\xi }_{\alpha ,b}\,{\xi }_{\alpha ,b}^{\prime }\,{\bar{\xi }}_{\alpha ,b}^{\prime }\,{\bar{\xi }}_{\alpha ,b}={\bar{\xi }}_{\alpha ,b}^{\prime }\,{\xi }_{\alpha ,b}^{\prime }\,{\bar{\xi }}_{\alpha ,b}\,{\xi }_{\alpha ,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 (2.4)

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \prod _{\gamma \ne 0}\displaystyle \int { \mathcal D }{[\bar{\xi }^{\prime} \xi ^{\prime} \bar{\xi }\xi ]}_{\gamma }{\bar{\xi }}_{\alpha }^{\prime }{\bar{\xi }}_{\alpha }\left[1+\displaystyle \frac{{\bar{\xi }}_{\gamma }^{\prime }{\xi }_{\gamma }^{\prime }}{\sqrt{f(\gamma )}}\right]\left[1+\displaystyle \frac{{\bar{\xi }}_{\gamma }{\xi }_{\gamma }}{\sqrt{f(\gamma )}}\right]{\xi }_{\alpha }{\xi }_{\alpha }^{\prime }\\ & & \mathop{\displaystyle \quad =\prod _{\gamma \ne 0,\alpha }}\limits_{}\displaystyle \frac{1}{f(\gamma )}=\displaystyle \frac{f(\alpha )}{{f}_{* }}.\end{array}\end{eqnarray} \tag{ 2.5 }$$

We see that, apart from the irrelevant global factor ${\prod }_{b}1/{f}_{* }$, we reproduce the configuration space and weight factors of the original link model.

To generate the Feynman diagrammatic expansion, we need to represent the Grassmann weight factor in the exponential form. The link factors (2.1) have the form of Gaussian exponentials already. Hence, it is only the site factors that need to be rewritten identically as

$$\begin{eqnarray}&&{B}_{j}=\exp \left[\displaystyle \sum _{{\{{\alpha }_{b}\}}_{j}}\lambda ({\{{\alpha }_{b}\}}_{j})\displaystyle \prod _{b\in \{b\}{}_{j}}{\breve{\xi }}_{{\alpha }_{b},b}\,{\breve{\xi }}_{{\alpha }_{b},b}^{* }\right].\end{eqnarray} \tag{ 2.6 }$$

The constants $\lambda ({\{{\alpha }_{b}\}}_{j})$ are readily related to the site factors $g({\{{\alpha }_{b}\}}_{j})$ by simple algebraic equations obtained by expanding the exponential and equating sums of similar terms to their counterparts in the rhs of equation (2.2).

By expanding the non-Gaussian part of the exponential (2.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 (2.6) takes us back to equation (2.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.

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 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\equiv f(\alpha )$. 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 ${\{{\alpha }_{p}\}}_{j}$ for the edge j); it encodes, if necessary, constraints on the allowed sets of ${\{{\alpha }_{p}\}}_{j}$.

Without loss of generality (up to a global normalization factor), we identify the 'ground state' as ${\alpha }_{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 $\nu \equiv {\nu }_{p}=0,1,2,3$ enumerates four vertices in such a way that $\nu \pm 1$ modulo 4 stands for the next/previous vertex with respect to the vertex ν in the clockwise direction.

For each state $\alpha \ne 0$ of the plaquette p, we introduce eight Grassmann variables: ${\xi }_{\alpha ,p,{\nu }_{p}}$, ${\bar{\xi }}_{\alpha ,p,{\nu }_{p}},{\nu }_{p}=0,1,2,3$. As before, for $\alpha =0$ the variables ξ and $\bar{\xi }$ are not Grassmannian, ${\xi }_{0,p,\nu }=0$, ${\bar{\xi }}_{0,p,\nu }=1$. The corresponding plaquette weight in the Grassmann partition function reads

$$\begin{eqnarray}&&{A}_{p}=\exp \left\{\displaystyle \sum _{\alpha \ne 0}\,{[-f(\alpha )]}^{-1/4}\displaystyle \sum _{{\nu }_{p}=0}^{3}\,{\bar{\xi }}_{\alpha ,p,{\nu }_{p}}\,{\xi }_{\alpha ,p,{\nu }_{p}}\right\}.\end{eqnarray} \tag{ 3.1 }$$

Note a close analogy with equation (2.1). Site weights equation (2.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 ${\alpha }_{p}=0$, we write

$$\begin{eqnarray}&&{B}_{j}=1+\displaystyle \sum _{{\{{\alpha }_{p}\}}_{j}\ne {0}_{j}}g({\{{\alpha }_{p}\}}_{j})\displaystyle \prod _{p\in \{p\}{}_{j}}{\xi }_{\alpha ,p,({\nu }_{p}^{(j)}+1)}\,{\bar{\xi }}_{\alpha ,p,{\nu }_{p}^{(j)}},\end{eqnarray} \tag{ 3.2 }$$

where ${\nu }_{p}^{(j)}$ 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 $({\nu }_{p}^{(j)}+1)$.) In what follows, we will associate ${\nu }_{p}^{(j)}$ 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 ${\alpha }_{p}$ is used in all edge weights. The $\alpha =0$ contribution comes from the term

$$\begin{eqnarray}&&-\displaystyle \prod _{\gamma }\displaystyle \frac{1}{f(\gamma )}\,\displaystyle \prod _{{\nu }_{p}=0}^{3}\,{\bar{\xi }}_{\gamma ,p,{\nu }_{p}}\,{\xi }_{\gamma ,p,{\nu }_{p}}\qquad ({\rm{at}}\,\,\alpha =0)\end{eqnarray} \tag{ 3.3 }$$

in the expansion of the exponential (3.1). It contributes a factor $1/{q}_{* }$, where

$$\begin{eqnarray}&&{q}_{* }=\displaystyle \prod _{\gamma }(-1)f(\gamma ).\end{eqnarray} \tag{ 3.4 }$$

The $\alpha \ne 0$ contribution comes from the plaquette term

$$\begin{eqnarray}&&\displaystyle \prod _{\gamma \ne \alpha }\displaystyle \frac{1}{f(\gamma )}\,\displaystyle \prod _{{\nu }_{p}=0}^{3}\,{\bar{\xi }}_{\gamma ,p,{\nu }_{p}}\,{\xi }_{\gamma ,p,{\nu }_{p}}\qquad ({\rm{at}}\,\,\alpha \ne 0)\end{eqnarray} \tag{ 3.5 }$$

multiplied by the product ${\prod }_{{\nu }_{p}=0}^{3}\,{\bar{\xi }}_{\alpha ,p,{\nu }_{p}}\,{\xi }_{\alpha ,p,{\nu }_{p}}$ originating from the boundary edge terms ${\xi }_{\alpha ,p,({\nu }_{p}^{(j)}+1)}\,{\bar{\xi }}_{\alpha ,p,{\nu }_{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(\alpha )$ in equation (3.1). Upon Grassmann integration, the contribution to the partition function of the resulting term equals to $f(\alpha )/{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 equation (3.1) with four pairs of Grassmann fields for every non-zero plaquette state. The interaction part of the Grassmann action is contained in edge weights (3.2) after they are written in an exponential form

$$\begin{eqnarray}&&{B}_{j}=\exp \left[\displaystyle \sum _{{\{{\alpha }_{p}\}}_{j}}\lambda ({\{{\alpha }_{p}\}}_{j})\displaystyle \prod _{p\in \{p\}{}_{j}}{\xi }_{\alpha ,p,({\nu }_{p}^{(j)}+1)}\,{\bar{\xi }}_{\alpha ,p,{\nu }_{p}^{(j)}}\right],\end{eqnarray} \tag{ 3.6 }$$

with the constants $\lambda ({\{{\alpha }_{b}\}}_{j})$ unambiguously related to the edge factors $g({\{{\alpha }_{b}\}}_{j})$.

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

$$\begin{eqnarray}&&W(\{{\alpha }_{b}\})=\displaystyle \prod _{{b}_{1},{b}_{2}}F({\alpha }_{{b}_{1}},{b}_{1};{\alpha }_{{b}_{2}},{b}_{2}),\end{eqnarray} \tag{ 3.7 }$$

based on products of two-link factors. Without loss of generality, these factors can be cast into the exponential form

$$\begin{eqnarray}&&W(\{{\alpha }_{b}\})=\displaystyle \prod _{{b}_{1},{b}_{2}}\,{{\rm{e}}}^{-(1/2){\eta }_{{\alpha }_{{b}_{1}},{b}_{1};{\alpha }_{{b}_{2}},{b}_{2}}},\end{eqnarray} \tag{ 3.8 }$$

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 (3.8) identically (up to normalization) as

$$\begin{eqnarray}&&W(\{{\alpha }_{b}\})=\int { \mathcal D }X\displaystyle \prod _{b}\,{{\rm{e}}}^{{{\rm{i}}{X}}_{{\alpha }_{b},b}}{W}_{{\rm{G}}}(\{{X}_{{\alpha }_{b},b}\}).\end{eqnarray} \tag{ 3.9 }$$

Here $\{{X}_{{\alpha }_{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

$$\begin{eqnarray}&&{\eta }_{{\alpha }_{{b}_{1}},{b}_{1};{\alpha }_{{b}_{2}},{b}_{2}}=\langle {X}_{{\alpha }_{{b}_{1}},{b}_{1}}{X}_{{\alpha }_{{b}_{2}},{b}_{2}}\rangle .\end{eqnarray} \tag{ 3.10 }$$

What we achieve for a fixed set of X variables is a link model that contains only single-link factors

$$\begin{eqnarray}&&\forall b\,:\qquad {f}_{b}({\alpha }_{b})={{\rm{e}}}^{{{\rm{i}}{X}}_{{\alpha }_{b},b}}.\end{eqnarray} \tag{ 3.11 }$$

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

$$\begin{eqnarray}&&\displaystyle \sum _{{\alpha }_{b}}{f}_{b}({\alpha }_{b})\,\to \,{{ \mathcal W }}_{b}^{(G)}=\exp \left[{\bar{\xi }}_{b}{\xi }_{b}\left(\displaystyle \sum _{{\alpha }_{b}}\,{{\rm{e}}}^{{{\rm{i}}{X}}_{{\alpha }_{b},b}}\right)\right].\end{eqnarray} \tag{ 3.12 }$$

Note that here Grassmann variables have nothing to do with the discrete index ${\alpha }_{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 ${\alpha }_{s}$ residing on lattice sites when the configuration weight is given by

$$\begin{eqnarray}&&W(\{{\alpha }_{s}\})=\displaystyle \prod _{{s}_{1},{s}_{2}}\,{{\rm{e}}}^{-(1/2){\eta }_{{\alpha }_{{s}_{1}},{s}_{1};{\alpha }_{{s}_{2}},{s}_{2}}}.\end{eqnarray} \tag{ 3.13 }$$

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 non-negative 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 $\zeta =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 $\zeta =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 ζ.

As an illustration, consider the theory (3.7) and (3.8) and its Grassmann counterpart (3.12). Introduce the ζ-dependence by the replacement

$$\begin{eqnarray}&&{{\rm{e}}}^{{{\rm{i}}{X}}_{{\alpha }_{b},b}}\,\to \,{{\rm{e}}}^{{\rm{i}}\zeta {X}_{{\alpha }_{b},b}}.\end{eqnarray} \tag{ 3.14 }$$

In terms of the original theory, the replacement (3.14) means $\eta \to {\zeta }^{2}\eta$, for all η's in equation (3.8). If amplitudes of all η values in (3.8) 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 (3.12), the expansion of the exponential ${{\rm{e}}}^{{\rm{i}}\zeta {X}_{s,{\alpha }_{b}}}$ in powers of ζ generates an infinite series of interaction vertexes (the zeroth-order term defines the harmonic action):

$$\begin{eqnarray}&&{\bar{\xi }}_{b}{\xi }_{b}\displaystyle \sum _{{\alpha }_{b}}\left({\rm{i}}\zeta {X}_{{\alpha }_{b},b}-\displaystyle \frac{1}{2}{\zeta }^{2}{X}_{{\alpha }_{b},b}^{2}-\displaystyle \frac{{\rm{i}}}{3!}{\zeta }^{3}{X}_{\,{\alpha }_{b},b}^{3}+\cdots \right).\end{eqnarray} \tag{ 3.15 }$$

Higher-order vertexes in X come with a higher power of ζ and this sets unambiguously the rules for defining the diagram order.

Consider the 2D Ising model on the square lattice with the Hamiltonian

$$\begin{eqnarray}&&-H/T=\beta \displaystyle \sum _{\langle i,j\rangle }{\sigma }_{i}{\sigma }_{j}+\displaystyle \sum _{i}{h}_{i}{\sigma }_{i}.\end{eqnarray} \tag{ 4.1 }$$

The Ising variables $\sigma =\pm 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\geqslant 0$ without loss of generality. The partition function of the Ising model reads

$$\begin{eqnarray}&&Z=\displaystyle \sum _{\{{\sigma }_{{i}}\}}\underset{\langle {i},j\rangle }{\overset{}{{\rm{\Pi }}}}{{\rm{e}}}^{\beta {\sigma }_{{i}}{\sigma }_{j}}\,\displaystyle \prod _{{i}}{{\rm{e}}}^{{h}_{{i}}{\sigma }_{{i}}}.\end{eqnarray} \tag{ 4.2 }$$

The most typical observable of the Ising model is the spin–spin correlation function ${\rho }_{{ij}}$

$$\begin{eqnarray}&&{\rho }_{{ij}}=\langle {\sigma }_{i}{\sigma }_{j}\rangle =\displaystyle \frac{1}{Z}{\displaystyle \frac{{\partial }^{2}Z}{\partial {h}_{i}\partial {h}_{j}}| }_{{h}_{i}={h}_{j}=h.}\end{eqnarray} \tag{ 4.3 }$$

Using the well-known identities

$$\begin{eqnarray}{{\rm{e}}}^{\beta {\sigma }_{i}{\sigma }_{j}} & = & \cosh \beta (1+{\sigma }_{i}{\sigma }_{j}\tanh \beta )\\ {{\rm{e}}}^{h{\sigma }_{i}} & = & \cosh h(1+{\sigma }_{i}\tanh h),\end{eqnarray} \tag{ 4.4 }$$

the partition function can be written as $Z={Z}_{0}Z^{\prime}$ with ${Z}_{0}={(\cosh \beta )}^{2N}{(\cosh h)}^{N}$ for a lattice of N sites and $2N$ links. With the notation $\zeta =\tanh \beta$ and $\eta =\tanh h$ the remaining factor is given by

$$\begin{eqnarray}&&Z^{\prime} =\displaystyle \sum _{\{{\sigma }_{i}\}}\,\displaystyle \prod _{\langle i,j\rangle }(1+{\sigma }_{i}{\sigma }_{j}\zeta )\,\displaystyle \prod _{i}(1+{\sigma }_{i}\eta ).\end{eqnarray} \tag{ 4.5 }$$

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+{\sigma }_{i}{\sigma }_{j}\zeta )$. In the partition function, terms with an odd power of ${\sigma }_{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

$$\begin{eqnarray}&&f(0)=1,\,\,f(1)={f}_{* }\,=\,\zeta ,\end{eqnarray} \tag{ 4.6 }$$

$$\begin{eqnarray}&&g(0)=g(2)=g(4)=1,\,\,g(1)=g(3)=\eta .\end{eqnarray} \tag{ 4.7 }$$

Here we label site factors using the total sum of incident link variables, ${\sum }_{b\in \{b\}{}_{j}}{\alpha }_{b}$, to avoid unnecessary rank-4 tensor notations. If we further redefine ${Z}_{0}\to {Z}_{0}{2}^{N}{f}_{* }^{2N}$, then the Grassmann representation of the partition function $Z^{\prime}$ is given by

$$\begin{eqnarray}Z^{\prime} & = & \displaystyle \int { \mathcal D }{[\bar{\xi }^{\prime} \xi ^{\prime} \bar{\xi }\xi ]}_{\{{\alpha }_{b}\}}\displaystyle \prod _{\{{\alpha }_{b}\}}\exp \left(\displaystyle \frac{1}{\sqrt{\zeta }}{\bar{\xi }}_{{\alpha }_{b}}^{\prime }{\xi }_{{\alpha }_{b}}+\displaystyle \frac{1}{\sqrt{\zeta }}{\bar{\xi }}_{{\alpha }_{b}}{\xi }_{{\alpha }_{b}}\right)\\ & & \times \exp \left(\displaystyle \sum _{j}{\lambda }_{{\alpha }_{\{b\}}}\displaystyle \prod _{{b}_{j}}{\breve{\xi }}_{{\alpha }_{b}}{\breve{\xi }}_{{\alpha }_{b}}^{* }\right).\end{eqnarray} \tag{ 4.8 }$$

We now compute the factors λ. To this end, we first introduce notations (for a fixed site j and suppressing the site index for clarity)

$$\begin{eqnarray}\begin{array}{rcl}{V}_{1} & = & {\breve{\xi }}_{{R}}{\breve{\xi }}_{{R}}^{* }+{\breve{\xi }}_{{U}}{\breve{\xi }}_{{U}}^{* }+{\breve{\xi }}_{{L}}{\breve{\xi }}_{{L}}^{* }+{\breve{\xi }}_{{D}}{\breve{\xi }}_{{D}}^{* }={n}_{{R}}+{n}_{{U}}+{n}_{{L}}+{n}_{{D}},\\ {V}_{2} & = & {n}_{{R}}{n}_{{U}}+{n}_{{R}}{n}_{{L}}+{n}_{{R}}{n}_{{D}}+{n}_{{U}}{n}_{{L}}+{n}_{{U}}{n}_{{D}}+{n}_{{L}}{n}_{{D}},\\ {V}_{3} & = & {n}_{{R}}{n}_{{U}}{n}_{{L}}+{n}_{{R}}{n}_{{U}}{n}_{{D}}+{n}_{{R}}{n}_{{L}}{n}_{{D}}+{n}_{{U}}{n}_{{L}}{n}_{{D}},\\ {V}_{4} & = & {n}_{{R}}{n}_{{U}}{n}_{{L}}{n}_{{D}}\end{array}\end{eqnarray} \tag{ 4.9 }$$

and then Taylor expand

$$\begin{eqnarray}&&\exp [{\lambda }_{1}{V}_{1}+{\lambda }_{2}{V}_{2}+{\lambda }_{3}{V}_{3}+{\lambda }_{4}{V}_{4}].\end{eqnarray} \tag{ 4.10 }$$

The only non-zero terms generated by this expansion are ${V}_{1}^{2}=2{V}_{2},{V}_{1}^{3}=6{V}_{3},{V}_{1}^{4}=24{V}_{4},{V}_{1}{V}_{2}=3{V}_{3},{V}_{1}{V}_{3}=4{V}_{4}$ and ${V}_{2}^{2}=6{V}_{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

$$\begin{eqnarray}&&\begin{array}{l}\exp [{\lambda }_{1}{V}_{1}+{\lambda }_{2}{V}_{2}+{\lambda }_{3}{V}_{3}+{\lambda }_{4}{V}_{4}]\\ \quad =1+{\lambda }_{1}{V}_{1}+{\lambda }_{2}{V}_{2}+{\lambda }_{3}{V}_{3}+{\lambda }_{4}{V}_{4}\\ \quad \quad +\displaystyle \frac{1}{2}({\lambda }_{1}^{2}2{V}_{2}+{\lambda }_{2}^{2}6{V}_{4}+2{\lambda }_{1}{\lambda }_{2}3{V}_{3}+2{\lambda }_{1}{\lambda }_{3}4{V}_{4})\\ \quad \quad +\displaystyle \frac{1}{6}({\lambda }_{1}^{3}6{V}_{3}+3{\lambda }_{1}^{2}{\lambda }_{2}12{V}_{4})+\displaystyle \frac{1}{24}{\lambda }_{1}^{4}24{V}_{4}\end{array}\end{eqnarray} \tag{ 4.11 }$$

Term-by-term matching with equation (4.7) then leads to

$$\begin{eqnarray}&&{g}_{1}=\eta ={\lambda }_{1},\end{eqnarray} \tag{ 4.12 }$$

$$\begin{eqnarray}&&{g}_{2}=1={\lambda }_{2}+{\lambda }_{1}^{2},\end{eqnarray} \tag{ 4.13 }$$

$$\begin{eqnarray}&&{g}_{3}=\eta ={\lambda }_{3}+3{\lambda }_{1}{\lambda }_{2}+{\lambda }_{1}^{3},\end{eqnarray} \tag{ 4.14 }$$

$$\begin{eqnarray}&&{g}_{4}=1={\lambda }_{4}+3{\lambda }_{2}^{2}+4{\lambda }_{1}{\lambda }_{3}+6{\lambda }_{2}{\lambda }_{1}^{4}+{\lambda }_{1}^{4}.\end{eqnarray} \tag{ 4.15 }$$

The solution is immediate

$$\begin{eqnarray}&&{\lambda }_{1}=\eta ,\end{eqnarray} \tag{ 4.16 }$$

$$\begin{eqnarray}&&{\lambda }_{2}=1-{\eta }^{2},\end{eqnarray} \tag{ 4.17 }$$

$$\begin{eqnarray}&&{\lambda }_{3}=-2\eta +2{\eta }^{3},\end{eqnarray} \tag{ 4.18 }$$

$$\begin{eqnarray}&&{\lambda }_{4}=-2+8{\eta }^{2}-6{\eta }^{4}.\end{eqnarray} \tag{ 4.19 }$$

In zero external field the only vertexes with non-zero coupling in the partition function are V₂ and V₄

$$\begin{eqnarray}&&\displaystyle \prod _{j}\exp ({V}_{2}^{(j)}+{V}_{4}^{(j)})=\exp \left\{\displaystyle \sum _{j}({V}_{2}^{(j)}-2{V}_{4}^{(j)})\right\}\end{eqnarray} \tag{ 4.20 }$$

and the spin–spin correlation function is given by

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{ij}} & = & \displaystyle \frac{1}{Z^{\prime} }\displaystyle \prod _{b}\displaystyle \int { \mathcal D }{[...]}_{b}\exp \left(\displaystyle \frac{1}{\sqrt{\zeta }}{\bar{\xi }}_{b}^{\prime }{\xi }_{b}^{\prime }+\displaystyle \frac{1}{\sqrt{\zeta }}{\bar{\xi }}_{b}{\xi }_{b}\right)\\ & & \times \exp \left\{\displaystyle \sum _{j}({V}_{2}^{(j)}-2{V}_{4}^{(j)})\right\}\\ & & \times ({V}_{1}^{(i)}-2{V}_{3}^{(i)})({V}_{1}^{(j)}-2{V}_{3}^{(j)}).\end{array}\end{eqnarray} \tag{ 4.21 }$$

In order to derive the Feynman rules generating the diagrammatic series, we write the partition function in the form

$$\begin{eqnarray}&&Z^{\prime} ={Z}_{\updownarrow }\left(\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{{x}_{1},\ldots ,{x}_{n}}\displaystyle \frac{{(+1)}^{n}}{n!}\langle V({x}_{1})\ldots V({x}_{n}){\rangle }_{0}\right),\end{eqnarray} \tag{ 4.22 }$$

where ${Z}_{\updownarrow }$ is the partition function of the Gaussian part (it is the product of local link contributions), ${Z}_{\updownarrow }={\prod }_{b}\int { \mathcal D }[...]\exp \left({\bar{\xi }}_{b}^{\prime }{\xi }_{b}^{\prime }+\tfrac{1}{\zeta }{\bar{\xi }}_{b}{\xi }_{b}\right)={(1+\zeta )}^{(2N)}$. The Feynman rules for the correlation function of the 2D Ising model now follow from the textbook considerations:

• (i)  
  The bare propagators ${G}^{(0)}=\sqrt{\zeta }$ for primed and non-primes variables are local and reside on the links of the original lattice. In the correlation function they always occur in pairs of conjugate Grassmann variables and each pair contributes a factor ζ. The propagation lines do not have arrows. The bare interaction vertexes (or pre-diagrams, see figure 2) are also local and live on the sites of the lattice. There are different types belonging to the V₂ and V₄ classes with weight 1 and −2, respectively (see equation (4.20)). On the first (and last) site of the correlator we have a vertex belonging to the class V₁ or V₃ (see figures 3–7) with weight 1 and −2, respectively (see (4.21)).
• (ii)  
  Draw in order n all topologically distinct connected diagrams with n pairs of bi-grassmann variables living on the links of the lattice. The number of interaction vertexes, excluding the end points, is at most $n-1$.
• (iii)  
  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.
• (iv)  
  The total weight of the diagram in order n is hence ${(-1)}^{P}{(-2)}^{q}{\zeta }^{n}$ with P the signature of the exchange permutation and q the sum of all type-3 and type-4 vertexes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 figure 8). We check the connectivity of a diagram by the breadth-first algorithm. We have implemented the Feynman rules in two different appraoches: an algorithm that enumerates and evaluates all possible diagrams, as well as a DiagMC approach. Details about the implementation can be found in the appendix A

Zoom In Zoom Out Reset image size

For illustrative purposes, let us focus on the first element of the correlation function connecting the sites $(0,0)$ and $(1,0)$ (using translational invariance, any two neighboring sites $\langle {{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2}\rangle$ can be taken). To first order, we put a V₁ 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₂ (RD) vertex on the site $(0,1)$ and a V₂ 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{\zeta }^{3}$. These diagrams contributing in first and third order are shown in figure 3.

In fifth order, there are four diagrams with a V₃ vertex on one of the endpoints, yielding a contribution $-8{\zeta }^{5}$. There are 14 diagrams consisting of only V₁ and V₂ vertexes and single pair-lines, yielding a contribution $14{\zeta }^{5}$. The contributions to fifth order are shown in figures 4, 5, 7, and 8. There are however additional diagrams with two pairs of Grassmann variables living on the same link, as is shown in figure 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₁ and a V₂ (RU) vertex, and on the target site $(1,0)$ a V₁ and a V₂ (UL) vertex. On the site $(1,1)$ there is a V₂ (LD) and on site $(0,1)$ a V₂ (RD) vertex. Let us look more carefully at the link between the origin and target site:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\zeta }^{2}}\int { \mathcal D }[\bar{\xi }^{\prime} \bar{\xi }\xi \xi ^{\prime} ]\bar{\xi }^{\prime} \bar{\xi }\xi \xi ^{\prime} \bar{\xi }^{\prime} \bar{\xi }\xi \xi ^{\prime} .\end{eqnarray} \tag{ 4.23 }$$

The origin is associated with $\bar{\xi }^{\prime} \bar{\xi }$ and the target with $\xi \xi ^{\prime}$ by our convention. Applying Wick's theorem, there are four possible ways to pair the Grassmann variables:

• (i)  
  The pairing combination comes with the sign +1 and leads to a connected diagram (this is the lower right panel in figure 8).
• (ii)  
  The pairing combination comes with the sign −1 and leads to a connected diagram (this is the upper right panel in figure 8).
• (iii)  
  The pairing combination comes with the sign −1 and leads to a connected diagram (this is the lower left panel in figure 8).
• (iv)  
  The pairing combination leads to a disconnected diagram and does not contribute to the correlation function (this is the upper left panel in figure 8).

The net contribution of these four 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){\zeta }^{5}=4{\zeta }^{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 (4.23); and, in particular, from arbitrarily high powers of one and the same interaction vertex.

Our results for the spin–spin correlation function are shown in table 1. The correlation function is known recursively from [30–32]. It is also known as a Painlevé-VI nonlinear differential equation [33] 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 [26]), which are convenient for series expansions

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{(\mathrm{1,0})} & = & \coth (2\beta )\left[\displaystyle \frac{1}{2}+\displaystyle \frac{{\cosh }^{2}2\beta }{\pi }(2\,{\tanh }^{2}\,2\beta -1)K({k}_{\gt })\right]\\ & \to & \zeta +2{\zeta }^{3}+4{\zeta }^{5}+12{\zeta }^{7}+42{\zeta }^{9}+\cdots ,\end{array}\end{eqnarray} \tag{ 5.1 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{(\mathrm{1,1})} & = & \displaystyle \frac{2}{\pi {k}_{\gt }}[K({k}_{\gt })+({k}_{\gt }^{2}-1)K({k}_{\gt })]\\ & \to & 2{\zeta }^{2}+4{\zeta }^{4}+10{\zeta }^{6}+32{\zeta }^{8}+118{\zeta }^{10}+\cdots \mathrm{.}\end{array}\end{eqnarray} \tag{ 5.2 }$$

with ${k}_{\gt }={\sinh }^{2}(2\beta )$, $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 two rows in table. 1.

Table 1.  Expansion coefficients for the correlation function up to order 11.

Site/order	ζ	${\zeta }^{2}$	${\zeta }^{3}$	${\zeta }^{4}$	${\zeta }^{5}$	${\zeta }^{6}$	${\zeta }^{7}$	${\zeta }^{8}$	${\zeta }^{9}$	${\zeta }^{10}$	${\zeta }^{11}$
(1, 0)	1	0	2	0	4	0	12	0	42	0	164
(1, 1)	0	2	0	4	0	10	0	32	0	118	0
(2, 0)	0	1	0	6	0	16	0	46	0	158	0
(2, 1)	0	0	3	0	11	0	31	0	97	0	351
(2, 2)	0	0	0	0	6	0	24	76	0	248	0
(3, 0)	0	0	1	0	12	0	48	0	152	0	506
(3, 1)	0	0	0	4	0	26	0	92	0	298	0
(3, 2)	0	0	0	0	10	0	55	0	201	0	684
(3, 3)	0	0	0	0	0	20	0	120	0	480	0
(4, 0)	0	0	0	1	0	20	0	118	0	7 452	0
(4, 1)	0	0	0	0	5	0	52	0	244	0	885
(4, 2)	0	0	0	0	0	15	0	118	0	521	0
(4, 3)	0	0	0	0	0	0	25	0	259	0	1176
(4, 4)	0	0	0	0	0	0	0	70	0	560	0
(5, 0)	0	0	0	0	1	0	30	0	250	0	1200
(5, 1)	0	0	0	0	0	6	0	92	0	574	0
(5, 2)	0	0	0	0	0	0	21	0	231	0	1266
(5, 3)	0	0	0	0	0	0	0	56	0	532	0
(5, 4)	0	0	0	0	0	0	0	0	126	0	1176
(5, 5)	0	0	0	0	0	0	0	0	0	252	0
(6, 0)	0	0	0	0	0	1	0	42	0	474	0
(6, 1)	0	0	0	0	0	0	7	0	149	0	1215
(6, 2)	0	0	0	0	0	0	0	28	0	416	0
(6, 3)	0	0	0	0	0	0	0	0	84	0	1026
(6, 4)	0	0	0	0	0	0	0	0	0	210	0
(6, 5)	0	0	0	0	0	0	0	0	0	0	462
(7, 0)	0	0	0	0	0	0	1	0	56	0	826
(7, 1)	0	0	0	0	0	0	0	8	0	226	0
(7, 2)	0	0	0	0	0	0	0	0	36	0	699
(7, 3)	0	0	0	0	0	0	0	0	0	120	0
(7, 4)	0	0	0	0	0	0	0	0	0	0	330
(8, 0)	0	0	0	0	0	0	0	1	0	72	0
(8, 1)	0	0	0	0	0	0	0	0	9	0	326
(8, 2)	0	0	0	0	0	0	0	0	0	45	0
(8, 3)	0	0	0	0	0	0	0	0	0	0	165
(9, 0)	0	0	0	0	0	0	0	0	1	0	90
(9, 1)	0	0	0	0	0	0	0	0	0	10	0
(9, 2)	0	0	0	0	0	0	0	0	0	0	55
(10, 0)	0	0	0	0	0	0	0	0	0	1	0
(10, 1)	0	0	0	0	0	0	0	0	0	0	11
(11, 0)	0	0	0	0	0	0	0	0	0	0	1

The spin susceptibility is related to the zero momentum value of the Green function by ${\beta }^{-1}\chi =1+\rho ({\boldsymbol{p}}=0)$. We can hence sum over the entire lattice to obtain

$$\begin{eqnarray}\begin{array}{rcl}{\beta }^{-1}\chi & = & 1+4\zeta +12{\zeta }^{2}+36{\zeta }^{3}+100{\zeta }^{4}+276{\zeta }^{5}\\ & & +740{\zeta }^{6}+1972{\zeta }^{7}+5172{\zeta }^{8}+13492{\zeta }^{9}\\ & & +34876{\zeta }^{10}+89764{\zeta }^{11}+229628{\zeta }^{12}\\ & & +585508{\zeta }^{13}+1486308{\zeta }^{14}+\cdots \mathrm{.}\end{array}\end{eqnarray} \tag{ 5.3 }$$

To this order the series expansion agrees with the ones from [34, 35]. For a library of high-temperature series expansions, see [36]. Currently, the series is known (at least) up to order 2000 and still topic of active research [32, 34]. 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 $\gamma =7/4$. We plot in figure 9 the susceptibility versus β for different expansion orders, and also plot the asymptotic behavior for comparison.

Zoom In Zoom Out Reset image size

The critical temperature and the exponent γ can be found from a study of the convergence radius of the series. Since

$$\begin{eqnarray}{\beta }^{-1}\chi & = & \displaystyle \sum _{n}{\chi }_{n}{\zeta }^{n}\propto {(1-\zeta /{\zeta }_{c})}^{-\gamma }\\ & = & 1+\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{\gamma (\gamma -1)\cdots (\gamma +n-1)}{n!}{\left(\displaystyle \frac{\zeta }{{\zeta }_{c}}\right)}^{n}\end{eqnarray} \tag{ 5.4 }$$

the ratio of coefficients asympotically behaves as

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{n}}{{\chi }_{n-1}}=\displaystyle \frac{1}{{\zeta }_{c}}+\displaystyle \frac{\gamma -1}{{\zeta }_{c}}\displaystyle \frac{1}{n}.\end{eqnarray} \tag{ 5.5 }$$

In figure 10 we extract the critical point ${\zeta }_{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 [37], γ can be determined independently from ${\zeta }_{c}$ as $\gamma \approx 1.751949$ on the square lattice when the series is known up to 14th order, i.e., an accuracy of $0.5 \%$.

Zoom In Zoom Out Reset image size

In 1D, the only allowed vertex is RL (the last one in the second line of figure 2). It has weight +1. The only allowed endpoints are L and R (the second and fourth vertexes shown in the first line of figure 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 ${\zeta }^{n}$ (up to the lattice symmetry). The susceptibility is hence

$$\begin{eqnarray}&&T\chi =1+2(\zeta +{\zeta }^{2}+\cdots )=1+2\displaystyle \frac{\zeta }{1-\zeta },\end{eqnarray} \tag{ 6.1 }$$

reproducing the exact solution with asymptotic behavior $\chi \propto \beta \exp (2\beta )$ as $T\to 0$.

The ${G}^{2}W$ skeleton expansion becomes exact already in 0th order

$$\begin{eqnarray}&&{\rm{\Pi }}={{\rm{\Pi }}}^{0}=\zeta \end{eqnarray} \tag{ 6.2 }$$

$$\begin{eqnarray}&&{\rm{\Sigma }}=0\end{eqnarray} \tag{ 6.3 }$$

which yields $G={G}_{0}=\sqrt{\zeta }$, $W=V/(1-V{\rm{\Pi }})=1/(1-\zeta )$, and also ${\rm{\Pi }}=\zeta /(1-\zeta )$. This immediately leads to the same result as in equation (6.1) when adding the end-point vertexes L and R to Π.

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:

• (i)  
  MOVE–RETRACT. We choose one of the four 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₁ type of the tail turns into a V₂, otherwise the head goes back and changes the previous V₂ into a V₁ type (unless the diagram order is 0 or 1, when only V₁ 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₁.
• (ii)  
  SWAP VERTEX. Swaps between the vertexes ${V}_{1}+{V}_{2}\leftrightarrow {V}_{3}$ (for head and/or tail) and ${V}_{2}+{V}_{2}\leftrightarrow {V}_{4}$ (anywhere in the diagram). This update is its own reciprocal.
• (iii)  
  RELINK. On a given link, relink primed and non-primed 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.

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₃ and V₄ vertexes as well as possible fermionic exchanges. All diagrams in order n can be found as follows:

• (i)  
  Write down all possible words of the form ${X}_{1}{X}_{2}\ldots {X}_{n}$ with the alphabet ${X}_{j}\in \{0,3\}$ corresponding to the four directions on the square lattice. Make sure that subsequent directions are not backtracking. For example, if X₄ is in the positive $+\hat{x}$ direction, then X₅ cannot be in the negative $-\hat{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.
• (ii)  
  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).
• (iii)  
  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₃ ends), then for the head alone, and finally for the tail alone.
• (iv)  
  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.
• (v)  
  We compute the product of all the vertex weights, according to the Feynman rules.
• (vi)  
  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.
• (vii)  
  We check the connectivity of the diagram using the breadth-first algorithm. Disconnected diagrams contribute 0.
• (viii)  
  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}\approx {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⁻⁸ 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₂ 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.

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 figure B1 . In this setup, the expansion order is defined by the number of W-lines (obviously, the discussion of section 3.3 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₄ vertexes to account for all contributions. We start with neglecting V₄ vertexes, and discuss their role later.

Zoom In Zoom Out Reset image size

In more detail, the formalism of the ${G}^{2}W$ expansion in the absence of V₄ vertexes is as follows:

• (i)  
  There are six bare two-body interaction vertexes V₂ $({RU},{RL},{RD},{LU},{UD},{LD})$, see the second line in figure 2. They reside on the sites of the original square lattice and all have weight 1. Symbolically, we encode the tensor structure of V₂ using a convenient short hand notation ${V}_{2}={\sum }_{\alpha ,\gamma =1}^{4}V(\alpha ,\gamma ){n}_{\alpha }{n}_{\gamma }$, where

  $$\begin{eqnarray}V(\alpha ,\gamma )=\left[\begin{array}{cccc}0 & 1 & 1 & 1\\ 1 & 0 & 1 & 1\\ 1 & 1 & 0 & 1\\ 1 & 1 & 1 & 0\end{array}\right].\end{eqnarray} \tag{ B.1 }$$

  The row index represents the first leg enumerated according to the convention $(R,U,L,D)\to (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.
• (ii)  
  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

  $$\begin{eqnarray}&&G(\alpha ,\gamma )={G}^{(0)}(\alpha ,\gamma )+\displaystyle \sum _{\mu ,\nu }{G}^{(0)}(\alpha ,\mu ){\rm{\Sigma }}(\mu ,\nu )G(\nu ,\gamma ).\end{eqnarray} \tag{ B.2 }$$

  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 equation (B.2): 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₀₂; it can be alternatively denoted as G_x and, by 90^o rotation symmetry of the square lattice, is the same for all links.
• (iii)  
  The matrix structure of polarization Π is similar to that of V. The 0th order expression based on bare Green functions is given by

  $$\begin{eqnarray}{{\rm{\Pi }}}_{({\rm{x}},{\rm{y}})}^{(0)}(\alpha ,\gamma )=\zeta \left[\begin{array}{cccc}0 & 0 & {\delta }_{x,1}{\delta }_{y,0} & 0\\ 0 & 0 & 0 & {\delta }_{x,0}{\delta }_{y,1}\\ {\delta }_{x,-1}{\delta }_{y,0} & 0 & 0 & 0\\ 0 & {\delta }_{x,0}{\delta }_{y,-1} & 0 & 0\end{array}\right].\end{eqnarray} \tag{ B.3 }$$
• (iv)  
  The effective potential W is defined through the Dyson equation in momentum representation

  $$\begin{eqnarray}&&{W}_{{\boldsymbol{q}}}(\alpha ,\gamma )=V(\alpha ,\gamma )+\displaystyle \sum _{\mu ,\nu }V(\alpha ,\mu ){{\rm{\Pi }}}_{{\boldsymbol{q}}}(\mu ,\nu ){W}_{{\boldsymbol{q}}}(\nu ,\gamma ).\end{eqnarray} \tag{ B.4 }$$

  We expect to see signatures of the ferromagnetic transition in matrix elements of ${W}_{{\boldsymbol{q}}=0}$ because they directly relate to the divergent uniform susceptibility χ.

To obtain the 0th order result, we replace Π with ${{\rm{\Pi }}}^{(0)}$ in equation (B.4). For any ζ we compute ${W}_{{\boldsymbol{q}}=0}$ from equation (B.4) by matrix inversion. We find a divergence at ${\zeta }_{c}=1/3$ (shown in figure B2 ) that can be also established analytically. We see that W diverges as ${({\zeta }_{c}-\zeta )}^{-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 ${\zeta }_{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.

Zoom In Zoom Out Reset image size

We now include the diagrams with one W line for the selfenergy and the polarization. In real space we find

$$\begin{eqnarray}{{\rm{\Sigma }}}_{x}^{(1)} & = & {{\rm{\Sigma }}}^{(1)}=-{G}_{x}{W}_{(\mathrm{1,0})}(2,0)=-{{GW}}_{(\mathrm{1,0})}(2,0)\\ {{\rm{\Pi }}}^{(1)} & = & {G}^{4}{W}_{(\mathrm{0,0})}(0,0)+\mathrm{cycl}.\end{eqnarray} \tag{ B.5 }$$

The matrix structure of ${{\rm{\Pi }}}^{(1)}$ is identical to that of ${{\rm{\Pi }}}^{(0)}$ and is not shown here explicitly. Coupled self-consistent equations (B.2), (B.4) and (B.5), are solved by iterations.

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 third-order graphs for Ψ. The corresponding expressions in real space are

$$\begin{eqnarray}{{\rm{\Sigma }}}^{(2)} & = & -{W}_{(\mathrm{0,0})}(0,0){W}_{(\mathrm{0,0})}(2,2){G}^{3}\\ {{\rm{\Pi }}}_{(0,0)}^{(2)}(0,0) & = & {G}^{6}{W}_{(\mathrm{0,0})}(2,2){W}_{(\mathrm{1,0})}(0,2)\\ {{\rm{\Pi }}}_{(1,0)}^{(2)}(0,2) & = & {G}^{6}{W}_{(1,0)}^{2}(0,2)\\ & & +{G}^{6}{W}_{(\mathrm{0,0})}(0,0){W}_{(\mathrm{0,0})}(2,2).\end{eqnarray} \tag{ B.6 }$$

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 figure B3 . The transition point has slightly shifted to larger values of ζ compared to the zeroth-order result, and the exponent has also slightly increased.

Zoom In Zoom Out Reset image size

The ${G}^{2}W$-expansion scheme treats different bare vertexes (see figure 2) on unequal footing: the V₂ vertexes are fully dressed, but the V₄ vertexes are included perturbatively (we neglected them so far). These higher-rank vertexes have a weight of comparable magnitude to the V₂ vertexes (–2 for V₄ versus +1 for V₂). 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 ${\zeta }_{c}$ increases towards its exact value).

Formally, there is no valid reason for neglecting the V₄ vertexes altogether. Let us show how they can be taken care of in the spirit of the shifted action approach [40]. 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 $\zeta \ll 1$. By assuming that the skeleton sequence (without V₄) is solved, we introduce the full polarization function $\bar{{\rm{\Pi }}}(\hat{\alpha },\hat{\beta })$ through the Dyson equation

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{{\boldsymbol{q}}}(\alpha ,\gamma )={{\rm{\Pi }}}_{{\boldsymbol{q}}}(\alpha ,\gamma )+\displaystyle \sum _{\mu ,\nu }{{\rm{\Pi }}}_{{\boldsymbol{q}}}(\alpha ,\mu )V(\mu ,\nu ){\bar{{\rm{\Pi }}}}_{{\boldsymbol{q}}}(\nu ,\gamma ).\end{eqnarray} \tag{ B.7 }$$

To be specific, we focus on the n.n. element ${\rho }_{(\mathrm{1,0})};$ similar manipulations hold for any other distance. Now consider all diagrams for this correlator without the V₄ vertexes within the ${G}^{2}W$ formulation (see [40]):

• Put one V₁ vertex on the origin site $(0,0)$ and the other V₁ vertex on the target site ${\bf{r}}=(1,0)$, see equation (4.21). There are $4\times 4=16$ different ways of doing that depending on the directions of legs. Connect the legs with ${\bar{{\rm{\Pi }}}}_{{\boldsymbol{r}}}(\alpha ,\gamma )$. For example, in the limit of $\zeta \ll 1$, choosing the ($\alpha =0$)-leg on site $(0,0)$ and the $\gamma =2$-leg on site $(1,0)$ results in the contributions $\zeta -4{\zeta }^{5}+\cdots$. Similarly, the choice of $\alpha =1$ and $\gamma =1$ leads to the contribution ${\zeta }^{3}$.
• Put V₃ on $(0,0)$ and V₁ on $(1,0)$, and connect all legs with $\bar{{\rm{\Pi }}}$ lines. There are four ways to orient the V₃ vertex and for each one there are two choices for connecting legs with $\bar{{\rm{\Pi }}}$ propagators. The leading contribution to ${\rho }_{(\mathrm{1,0})}$ goes hence as $-8{\zeta }^{5}$.
• Putting V₁ on $(0,0)$ and V₃ on $(1,0)$ gives the same contribution by symmetry.
• Put one V₃ vertex on $(0,0)$ and the other V₃ vertex on $(1,0)$. Now there are 16 ways of orienting both V₃ vertexes, and for each orientation there are 15 choices for connecting the legs. These contributions start at order $\propto {\zeta }^{9}$.

Next, we repeat the above procedure of connecting legs by adding one V₄ vertex, which can be put on any site, after that we can add two V₄ vertexes etc to generate a perturbative expansion in the number of V₄ terms Compared to the original bare series in powers of ζ, we have reordered the series: the effective potential is summing up all V₂ vertexes, whereas we expand (and sample in a Monte Carlo framework) in powers of ${\lambda }_{4}$.

To illustrate this framework, let us take $\zeta =0.01$ and recall that in the bare series ${\rho }_{(\mathrm{1,0})}=\zeta +2{\zeta }^{3}+4{\zeta }^{5}+12{\zeta }^{7}+\cdots$. The first three terms can be reproduced without V₄ vertexes and with only 1 V₃ on either the origin or the target site, see figures 3–8. The fifth order coefficient originates from 16 'simple' diagrams containing just V₁ and V₂ vertexes without any exchange. The diagrams containing a V₃ vertex yield a coefficient −8, and the exchange diagrams yield a coefficient −4. On a 16 × 16 lattice, the propagators obtained in section B.2 (i.e., to zeroth order) are

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}(0,2)=1.00000002\times {10}^{-02},\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}(1,1)=1.00010011\times {10}^{-06},\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}(1,2)=1.00080057\times {10}^{-10},\end{eqnarray} \tag{ B.10 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(0,0)}^{(0)}(1,2)=1.00020021\times {10}^{-08}.\end{eqnarray} \tag{ B.11 }$$

We do not mention explicitly other symmetry-related elements. The sum of all matrix elements for ${\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}$ is 0.01 000 200 160. One clearly recognizes the coefficients 1, 2 and 16 for the first-, third- and fifth-order contributions to the bare series. Contributions from the V₃ vertexes can be estimated from multiplying ${\bar{{\rm{\Pi }}}}_{(x,y)=(0,0)}^{(0)}(1,2)\times {\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}(0,2)$ which yields $\approx {10}^{-10}$. There are four different diagrams, each with weight −2, resulting in the above-mentioned coefficient −8.

On a 16 × 16 lattice, the propagators obtained in section B.4 (i.e., to second order) are

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(\mathrm{1,0})}(0,2)=9.99999980\times {10}^{-03},\end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(\mathrm{1,0})}(1,1)=1.00009999\times {10}^{-06},\end{eqnarray} \tag{ B.13 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(\mathrm{1,0})}(1,2)=1.00120089\times {10}^{-10},\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{\bar{{\rm{\Pi }}}}_{(x,y)=(\mathrm{0,0})}(1,2)=1.00020005\times {10}^{-08}.\end{eqnarray} \tag{ B.15 }$$

The sum of all matrix elements for ${\bar{{\rm{\Pi }}}}_{(x,y)=(1,0)}^{(0)}$ is 0.01 000 200 120. 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₃ diagrams in the way described above we recover the correct result to this order in ζ (which is +4).

The first instance of a V₄ vertex occurs in order ${\zeta }^{6}$ in the bare series. The relevant bare diagrams are the ones for ${\rho }_{(\mathrm{1,1})}$ with a V₄ 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₄ vertex on this site and connect all unpaired legs with $\bar{{\rm{\Pi }}}$ propagators. However, with the propagators obtained in section B.4 we are not supposed to account for all possible diagrams in the bare series to order 6 because our bold expansion in section B.4 is only accurate up to order ${\zeta }^{3}$: consider again ${\rho }_{(\mathrm{1,1})}$ and the bare diagrams where exchanges are possible on the links between the sites $(0,0)\mbox{--}(1,0)$ and $(1,0)\mbox{--}(1,1)$. Then there are irreducible non-local contributions that are not accounted for in section B.4 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 ${\zeta }^{6}$ contribution to ${\rho }_{(\mathrm{1,1})}$, which is below the correct value of 10.