The Dimer Model and Pfaffians

The dimer model is one of the most beautiful exactly solvable models in two-dimensional statistical mechanics. It begins with an elementary combinatorial question: in how many ways can one cover a finite graph by disjoint edges, so that every vertex is covered exactly once? But behind this innocent question lies a deep analytic mechanism. The partition function is related to a Pfaffians and determinants. Its local statistics are controlled by the inverse of a discrete operator. The purpose of this post is to explain the entire mechanism carefully. I will try to keep the algebra explicit, because the dimer model is one of those subjects where the conceptual statement is very short, but the real understanding comes from doing the calculation.

The fundamental dictionary is this:

\displaystyle \text{perfect matchings} \longleftrightarrow \text{pairings} \longleftrightarrow \text{Pfaffians}.

The only obstruction is sign. A Pfaffian is a signed sum over pairings, whereas a dimer partition function is a positive weighted sum over perfect matchings. Kasteleyn’s theorem says that for a planar graph one can orient the edges so that all the Pfaffian signs agree. Once this is done, the partition function is a Pfaffian, and all local statistics follow by differentiating that Pfaffian.

The finite dimer model

Let G=(V,E) be a finite graph with an even number of vertices, |V|=2n. A dimer configuration or perfect matching is a subset M\subset E such that every vertex of G is incident to exactly one edge of M . Thus each matching is a decomposition of the vertex set into n adjacent pairs. Assign each edge e\in E a positive weight w_e>0 . The weight of a matching is

\displaystyle w(M)=\prod_{e\in M} w_e.

The dimer partition function is the sum of these weights

\displaystyle Z_G = \sum_{M} \prod_{e\in M} w_e,

where the sum is over all perfect matchings.

The Gibbs probability of a matching is

\displaystyle \mathbb P(M) = \frac{1}{Z_G} \prod_{e\in M} w_e.

The basic statistical questions are the probabilities and covariances:

\displaystyle \mathbb P(e\in M), \quad \mathbb P(e_1,\dots,e_k\in M), \quad   \text{Cov}(\mathbf 1_e,\mathbf 1_f) = \mathbb E(\mathbf 1_e\mathbf 1_f) -\mathbb E(\mathbf 1_e)\mathbb E(\mathbf 1_f).

The first important observation is that the logarithmic derivative of Z_G already gives edge probabilities. Indeed,

\displaystyle w_e\frac{\partial Z_G}{\partial w_e} = \sum_{M\ni e} \prod_{f\in M} w_f.

Therefore

\displaystyle \mathbb P(e\in M) =\frac{w_e}{Z_G} \frac{\partial Z_G}{\partial w_e} = w_e\frac{\partial}{\partial w_e}\log Z_G.

Thus if one can compute \log Z_G , one can compute one-edge statistics by differentiation. Higher derivatives give higher correlations. The Pfaffian solution is powerful precisely because it turns \log Z_G into a logarithmic determinant.

Pfaffians

Let A be a 2n\times 2n skew-symmetric matrix. Thus A_{ij}=-A_{ji}, ~ A_{ii}=0. The Pfaffian of A is defined by

\displaystyle {\text{Pf}}(A) = \frac{1}{2^n n!} \sum_{\sigma\in S_{2n}} {\text{sgn}}(\sigma) \prod_{r=1}^{n} A_{\sigma(2r-1),\sigma(2r)}.

At first sight this formula looks complicated, but the idea is very simple. A permutation \sigma of the set {1,\dots,2n} arranges the indices in a row: \sigma(1),\sigma(2),\dots,\sigma(2n). Then we group this row into consecutive pairs:

\displaystyle \{\sigma(1),\sigma(2)\}, \{\sigma(3),\sigma(4)\}, \dots, \{\sigma(2n-1),\sigma(2n)\}.

So each permutation produces a pairing of the 2n indices. The product

\displaystyle \prod_{r=1}^{n} A_{\sigma(2r-1),\sigma(2r)}

assigns to that pairing the product of the corresponding matrix entries. The factor {\text{sgn}}(\sigma) records the sign of the permutation needed to put the indices into that paired order. However, the same unordered pairing is produced many times. First, the n pairs may be permuted among themselves in n! different ways. Second, inside each pair, the two elements may be interchanged, giving 2^n possibilities. This explains the normalizing factor \frac{1}{2^n n!}. Thus the Pfaffian is best thought of as a signed sum over all ways of partitioning {1,\dots,2n} into unordered pairs.

In other words,

\displaystyle {\text{Pf}}(A)= \sum_{\text{pairings } \pi} {\text{sgn}}(\pi) \prod_{ \{i,j\}\in\pi } A_{ij}.

This formula is often the most useful conceptual form. The Pfaffian is to pairings what the determinant is to permutations. The determinant sums over ways of matching rows to columns; the Pfaffian sums over ways of pairing indices with each other.

For example, take 2n=4 . Then a 4\times4 skew-symmetric matrix has the form

\displaystyle A= \begin{pmatrix} 0&a_{12}&a_{13}&a_{14}\\ -a_{12}&0&a_{23}&a_{24}\\ -a_{13}&-a_{23}&0&a_{34}\\ -a_{14}&-a_{24}&-a_{34}&0 \end{pmatrix}.

There are only three ways to pair the set {1,2,3,4} :

\displaystyle (12)(34), \quad (13)(24), \quad (14)(23).

Correspondingly,

\displaystyle {\text{Pf}}(A) = a_{12}a_{34} -a_{13}a_{24} + a_{14}a_{23}.

The signs are important. The pairings (12)(34) and (14)(23) occur with positive sign, while (13)(24) occurs with negative sign. This is the first place where one sees the essential feature of the Pfaffian: it counts pairings, but it counts them with signs.

This is exactly why Pfaffians enter the dimer model. A perfect matching of a graph is a pairing of the vertices, with the extra condition that each pair must be an edge of the graph. If the vertices of the graph are labeled 1,2, \dots,2n, then a perfect matching is a partition of these labels into n pairs

\displaystyle \{i_1,j_1\},\dots,\{i_n,j_n\},

where every \{i_r,j_r\} is an edge. Therefore a Pfaffian is naturally adapted to the enumeration of perfect matchings.

The only difficulty is that the dimer partition function is a positive sum,

\displaystyle Z_G = \sum_M \prod_{e\in M} w_e,

whereas the Pfaffian is a signed sum. Kasteleyn’s insight was that, for planar graphs, one can orient the edges so that the signs of all perfect matching terms become the same. Then the absolute value of the Pfaffian is exactly the dimer partition function.

The second basic fact about Pfaffians is the identity

\displaystyle {\text{Pf}}(A)^2=\det A.

This identity is the reason the Pfaffian method is so powerful analytically. The Pfaffian itself is combinatorial: it is a signed sum over pairings. But the square of the Pfaffian is a determinant, and determinants can be studied by linear algebra, eigenvalues, and Fourier analysis. Let us explain why the identity is natural. The determinant of a skew-symmetric matrix is a polynomial of degree 2n in the entries of A . The Pfaffian is a polynomial of degree n , since each term contains exactly n matrix entries. Thus {\text{Pf}}(A)^2 has degree 2n , the same degree as \det A . Moreover, both objects transform in the same way under a change of basis. If B is an invertible 2n\times2n matrix and we replace A by B A B^T, then

\displaystyle \det(BAB^T) = \det(B)\det(A)\det(B^T) = \det(B)^2\det(A).

On the Pfaffian side, one has the transformation law

\displaystyle {\text{Pf}}(BAB^T)= \det(B){\text{Pf}}(A).

Squaring this gives

\displaystyle {\text{Pf}}(BAB^T)^2 =\det(B)^2{\text{Pf}}(A)^2,

which is exactly the same transformation behavior as the determinant. Thus the identity {\text{Pf}}(A)^2=\det A is compatible with every linear change of coordinates. Now every skew-symmetric matrix can be put, over a suitable field such as \mathbb C , into a block-diagonal normal form consisting of 2\times2 skew blocks. That is, after a change of basis, one obtains a matrix of the form

\displaystyle A =\begin{pmatrix} 0&\lambda_1&&&&\\ -\lambda_1&0&&&&\\ &&0&\lambda_2&&\\ &&-\lambda_2&0&&\\ &&&&\ddots&\\ &&&&&\ddots \end{pmatrix}.

For this matrix, the Pfaffian is especially simple. Since the only nonzero pairings must pair the two indices inside each block, we get

\displaystyle {\text{Pf}}(A) = \lambda_1\lambda_2\cdots\lambda_n.

The determinant of each 2\times2 block

\displaystyle \begin{pmatrix} 0&\lambda_j\\ -\lambda_j&0 \end{pmatrix}

is \lambda_j^2. Therefore the determinant of the whole block-diagonal matrix is \det A = \lambda_1^2\lambda_2^2\cdots\lambda_n^2. Hence, in this normal form,

\displaystyle \det A = \big(\lambda_1\lambda_2\cdots\lambda_n\big)^2 ={\text{Pf}}(A)^2.

Since both sides transform compatibly under change of basis, the identity holds for every skew-symmetric matrix:

\displaystyle {\text{Pf}}(A)^2=\det A.

This identity is the algebraic bridge between combinatorics and analysis. In the dimer model, the Pfaffian appears because perfect matchings are pairings. But once the partition function has been written as a Pfaffian, one can square it and obtain a determinant. Then one can compute large-volume limits by diagonalizing operators, taking logarithms of eigenvalues, and replacing sums by integrals. This is precisely how the exact solution of the square-lattice dimer model passes from finite combinatorics to Fourier analysis.

In short, Pfaffian is a signed sum over pairings, while its square is equal to the determinant. The first statement explains why Pfaffians count dimers. The second explains why dimer models are analytically computable.

Matchings and Pfaffians

Label the vertices of G by 1,2,\dots,2n. We now want to build a matrix whose Pfaffian remembers the perfect matchings of G . Since a Pfaffian is naturally associated with a skew-symmetric matrix, the first step is to convert the weighted graph into a skew-symmetric weighted adjacency matrix. To do this, choose an orientation for every edge of G . Thus every undirected edge ij is assigned one of the two possible directions,

\displaystyle i\to j \quad\text{or}\qquad j\to i.

Once these choices have been made, define a 2n\times 2n matrix K by

\displaystyle K_{ij} = \begin{cases} w_{ij}, & \text{if the edge } ij \text{ is oriented } i\to j,\\ -w_{ij}, & \text{if the edge } ij \text{ is oriented } j\to i,\\ 0, & \text{if } i \text{ and } j \text{ are not adjacent}. \end{cases}

This matrix is automatically skew-symmetric. Indeed, if the edge ij is oriented i\to j , then K_{ij}=w_{ij}, K_{ji}=-w_{ij}, so K_{ji}=-K_{ij}. If i and j are not adjacent, then both entries are zero. Thus in all cases, K_{ij}=-K_{ji}. At this stage K is simply an oriented weighted adjacency matrix.

It becomes a Kasteleyn matrix only when the orientation satisfies the special Kasteleyn sign condition. For now, we are only using the fact that K is skew-symmetric, so that its Pfaffian is defined. Now expand the Pfaffian:

\displaystyle {\text {Pf}} (K) = \frac{1}{2^n n!} \sum_{\sigma\in S_{2n}} {\text {sgn} }(\sigma) \prod_{r=1}^{n} K_{\sigma(2r-1),\sigma(2r)}.

Every term in this sum chooses n pairs of vertices, \{\sigma(1),\sigma(2)\}, \{\sigma(3),\sigma(4)\}, \dots,\{\sigma(2n-1),\sigma(2n)\}. In other words, every term of the Pfaffian corresponds to a pairing of all vertices. The associated product is

\displaystyle K_{\sigma(1),\sigma(2)} K_{\sigma(3),\sigma(4)} \cdots K_{\sigma(2n-1),\sigma(2n)}.

This product is nonzero exactly when every pair is an edge of the graph. Thus the only nonzero terms in the Pfaffian expansion are precisely those pairings of the vertex set in which every pair is an edge of G . But that is exactly the definition of a perfect matching. Therefore the nonzero terms in {\text{Pf}}(K) are in one-to-one correspondence with perfect matchings of G . Thus each perfect matching contributes its usual dimer weight, but possibly with a sign. Therefore

\displaystyle {\text{Pf}}(K) = \sum_{M} \epsilon(M) \prod_{e\in M} w_e,

where the sum is over all perfect matchings of G , and where \epsilon(M)\in \{+1,-1\} is a sign depending on two things: the Pfaffian sign of the pairing and the orientation signs of the edges in M .

This is very close to the dimer partition function. The difference is only the signs. The partition function is a positive weighted sum over matchings. The Pfaffian is a signed weighted sum over the same matchings. Therefore the central problem is not to make the Pfaffian produce the correct terms; it already does that automatically. The central problem is to make the signs of those terms agree.

In other words, we want to choose the orientation of the edges so that \epsilon(M) is independent of M . If this can be done, then there is a fixed sign \epsilon_0\in \{+1,-1\} such that \epsilon(M)=\epsilon_0 for every matching M . Then

\displaystyle {\text {Pf}}(K) =\epsilon_0 \sum_M \prod_{e\in M} w_e =\epsilon_0 Z_G.

Taking absolute values gives

\displaystyle Z_G=|\text{Pf}(K)|.

This is the goal of the Kasteleyn orientation. One seeks an orientation of the edges for which the Pfaffian signs of all perfect matchings line up instead of canceling. This is why the Kasteleyn problem is a sign problem. The graph weights w_e already enter correctly. The Pfaffian already sums over exactly the right combinatorial objects. What remains is to arrange the orientation signs so that the Pfaffian becomes a positive enumeration rather than a signed enumeration.

For a completely arbitrary graph this cannot always be done in such a simple local way. The remarkable theorem of Kasteleyn is that for planar graphs it can be done. More precisely, if G is planar, one can orient the edges so that every face satisfies a certain parity condition. That local parity condition forces the relative Pfaffian sign between any two perfect matchings to be +1 . Consequently all matchings appear in \text{Pf}(K) with the same sign, and the dimer partition function becomes the absolute value of a Pfaffian.

The Kasteleyn orientation is the special choice of orientation for which this signed sum becomes, up to one global sign, the positive dimer partition function:

\displaystyle Z_G=|\text{Pf}(K)|.

This is the precise point at which the combinatorics of perfect matchings becomes linear algebra.

Kasteleyn orientation

Let M and M' be two perfect matchings of G . Their symmetric difference is

\displaystyle M\triangle M' = (M\setminus M')\cup(M'\setminus M).

This operation keeps exactly those edges that occur in one matching but not in the other. It is the natural object to study because the edges that belong to both matchings contribute in the same way to both Pfaffian terms; the only possible difference between the two signs comes from the edges on which the matchings disagree.

Now look at the graph formed by M\triangle M' . At any vertex, the matching M uses exactly one incident edge, and the matching M' also uses exactly one incident edge. If these two edges are the same, then that edge belongs to both matchings and is removed from the symmetric difference. If they are different, then both edges remain in M\triangle M' . Therefore every vertex that appears in M\triangle M' has degree exactly 2 : one incident edge comes from M and one incident edge comes from M' . A finite graph in which every vertex has degree 2 is a disjoint union of cycles. Moreover, these cycles must have even length, because as one walks around such a cycle, the edges alternate between M and M' . Thus

\displaystyle M\triangle M' =C_1\cup C_2\cup\cdots\cup C_r,

where each C_j is an even cycle, and along each cycle the edges alternate between the two matchings. These are called alternating cycles. This observation is fundamental. It tells us that to compare the Pfaffian signs of two arbitrary matchings, it is enough to compare two matchings that differ on one alternating cycle. Once the one-cycle comparison is understood, the general case follows by multiplying the contributions over all cycles in M\triangle M' .

So let C=(v_1,v_2,\dots,v_{2m},v_1) be one alternating cycle of length 2m . One matching uses the alternating edges (v_1v_2),(v_3v_4),\dots,(v_{2m-1}v_{2m}), while the other uses (v_2v_3),(v_4v_5),\dots,(v_{2m}v_1). We want to compare the two corresponding Pfaffian terms. There are two kinds of signs involved. First, there is the Pfaffian permutation sign. This is the sign coming from the order in which the vertices are paired in the Pfaffian expansion. Second, there is the orientation sign. This is the sign coming from the skew-symmetric entries of the matrix K . If an edge is oriented in the same direction as the ordered pair in the Pfaffian term, its matrix entry is +w_e ; if it is oriented in the opposite direction, its matrix entry is -w_e .

The Kasteleyn orientation is precisely designed so that these two sign effects cancel in the right way. To see this carefully, fix the cyclic order v_1,v_2,\dots,v_{2m}. For the first matching on the cycle, the local Pfaffian contribution may be written as K_{v_1v_2}K_{v_3v_4} \cdots K_{v_{2m-1}v_{2m}}. For the second matching, take the ordered pairs (v_2,v_3),(v_4,v_5),\dots,(v_{2m},v_1). The corresponding local product is K_{v_2v_3} K_{v_4v_5} \cdots K_{v_{2m}v_1}. The permutation (v_1,v_2,\dots,v_{2m}) \mapsto (v_2,v_3,\dots,v_{2m},v_1) is a cyclic shift of 2m objects. Its sign is (-1)^{2m-1}=-1. Thus, before considering the orientations of the edges, the second local pairing has an extra factor -1 relative to the first.

Now define the orientation sign along the cycle. Let

\displaystyle \eta_i = \begin{cases} +1, & \text{if the edge } v_iv_{i+1} \text{ is oriented } v_i\to v_{i+1},\\ -1, & \text{if the edge } v_iv_{i+1} \text{ is oriented } v_{i+1}\to v_i. \end{cases}

Here indices are understood modulo 2m , so v_{2m+1}=v_1 . Define

\displaystyle \eta(C) = \prod_{i=1}^{2m}\eta_i.

This number records whether the cycle has an even or odd number of edges oriented against the chosen cyclic direction. Since the cycle has even length, it is also equivalent to recording whether the number of edges oriented with the cyclic direction is even or odd. The ratio of the orientation signs of the two alternating matchings is exactly \eta(C) . Therefore the total relative sign between the two Pfaffian terms is -\eta(C). For the two matchings to contribute with the same sign, we want this relative sign to be +1 . Hence we need -\eta(C)=1. Equivalently,

\displaystyle \eta(C)=-1.

This is the key sign condition on an alternating cycle: the product of the orientation signs around the cycle must be -1 . In words, the cycle must have an odd number of edges oriented in the chosen cyclic direction. Such a cycle is often called oddly oriented. Thus the real goal is to have every alternating cycle should be oddly oriented. If this holds, then whenever two perfect matchings differ along an alternating cycle, their Pfaffian signs agree. Since any two perfect matchings differ by a disjoint union of alternating cycles, it follows that all perfect matchings contribute with the same Pfaffian sign. This is the heart of Kasteleyn’s method.

Now suppose G is planar and embedded in the plane. A Kasteleyn orientation is an orientation of the edges such that each bounded face is oddly oriented. More explicitly, if f is a bounded face and we traverse its boundary clockwise, then the number of boundary edges oriented clockwise is odd.

Equivalently, if \eta(f)=\prod_{e\in\partial f}\eta_e, where \eta_e=+1 if e is oriented in the clockwise boundary direction and \eta_e=-1 otherwise, then the Kasteleyn condition is \eta(f)=-1 for every bounded face f . For example, for a square face, this means that an odd number of the four boundary edges must be oriented clockwise. Thus the number of clockwise-oriented edges may be 1 or 3. For a hexagonal face, with six boundary edges, the same convention says again that an odd number of edges must be oriented clockwise. Thus the number may be 1, 3, 5.

Let us now explain why this local face condition implies the desired condition on alternating cycles. Let C be an alternating cycle that appears in the symmetric difference of two perfect matchings. Since C comes from two perfect matchings, the vertices strictly inside C are matched among themselves by the edges common to the two matchings. Hence the number of interior vertices is even. Let V_{\mathrm{int}} be the number of vertices strictly inside C , E_{\mathrm{int}} the number of edges strictly inside C , and F the number of bounded faces inside C . Euler’s formula for the planar region enclosed by C gives V_{\mathrm{int}}-E_{\mathrm{int}}+F=1. Since V_{\mathrm{int}} is even, this implies, modulo 2 , E_{\mathrm{int}}+F\equiv 1 \pmod 2. Now multiply the Kasteleyn face condition \eta(f)=-1 over all faces inside C . The left side becomes the product of all orientation signs along the boundaries of these faces. Every strictly interior edge is counted twice, once from each of its adjacent faces. The two clockwise boundary directions along that shared edge are opposite. Therefore an interior edge contributes a factor -1 to the total product. The boundary edges of C are counted once and together contribute precisely \eta(C) .

Therefore \prod_{f\text{ inside }C}\eta(f) = \eta(C)(-1)^{E_{\mathrm{int}}}. But each face satisfies \eta(f)=-1 , so the left side is (-1)^F. Thus \eta(C)(-1)^{E_{\mathrm{int}}} =(-1)^F. Hence \eta(C) =(-1)^{F+E_{\mathrm{int}}}. Using E_{\mathrm{int}}+F\equiv1\pmod2, we obtain \eta(C)=-1. So every alternating cycle is oddly oriented. This proves exactly what we needed. If two perfect matchings differ on one alternating cycle, their relative Pfaffian sign is -\eta(C). But for a Kasteleyn orientation we showed, \eta(C)=-1. Therefore -\eta(C)=1 and the two matchings have the same Pfaffian sign. Since any two perfect matchings differ by a disjoint union of alternating cycles, it follows that all perfect matchings have the same sign in the Pfaffian expansion. Thus there exists a fixed sign \epsilon_0\in \{+1,-1\} such that

\displaystyle \text{Pf}(K) = \epsilon_0 \sum_M \prod_{e\in M} w_e.

Therefore \text{Pf}(K)=\epsilon_0 Z_G. Taking absolute values gives the Kasteleyn formula:

\displaystyle Z_G=|\text{Pf}(K)|=\sqrt{|\det K|}..

This is the fundamental bridge from planar dimer combinatorics to linear algebra. The partition function, originally defined as a sum over exponentially many perfect matchings, becomes the square root of a determinant.

Statistics

Let e={i,j} be an edge of G . Suppose that, in our chosen orientation, this edge is oriented from i to j . Then the corresponding entries of the Kasteleyn matrix are K_{ij}=w_e, K_{ji}=-w_e. All dependence of K on the weight w_e occurs in these two entries. This simple observation is what makes edge probabilities computable by differentiating the Pfaffian.

If we differentiate Z_G with respect to w_e , only those matchings containing e contribute. Indeed,

\displaystyle \frac{\partial}{\partial w_e} = \prod_{f\in M}w_f =0 \quad \text{if } e\notin M,

whereas if e\in M , then

\displaystyle \frac{\partial}{\partial w_e} = \prod_{f\in M}w_f = \frac{1}{w_e} \prod_{f\in M}w_f.

Therefore

\displaystyle w_e\frac{\partial Z_G}{\partial w_e} = \sum_{M\ni e} \prod_{f\in M}w_f.

Dividing by Z_G , we obtain the elementary but very important identity

\displaystyle \mathbb P(e\in M) = \frac{w_e}{Z_G} \frac{\partial Z_G}{\partial w_e} = w_e\frac{\partial}{\partial w_e}\log Z_G.

This formula says that edge occupation probabilities are logarithmic derivatives of the partition function. In statistical mechanics language, inserting the observable \mathbf 1_{{e\in M}} is the same as differentiating the free energy with respect to the coupling \log w_e .

Now use the Kasteleyn formula Z_G=|\text{Pf}(K)|. Since a Kasteleyn orientation makes all perfect matchings appear with one common sign, the absolute value only removes a global sign. If the weights vary in a small positive neighborhood, that global sign remains fixed. Thus, for differentiation, we may compute with \log \text{Pf}(K) instead of \log |\text{Pf}(K)|.

Using the identity \text{Pf}(K)^2=\det K, we have

\displaystyle 2\log\text{Pf}(K)=\log\det K.

Differentiating with respect to w_e gives

\displaystyle 2\frac{\partial}{\partial w_e}\log\text{Pf}(K) = \frac{\partial}{\partial w_e}\log\det K.

The usual determinant identity is

\displaystyle \frac{\partial}{\partial t}\log\det K(t) = \text{Tr}\big(K(t)^{-1}\frac{\partial K(t)}{\partial t}\big).

Applying this with t=w_e , we obtain

\displaystyle \frac{\partial}{\partial w_e}\log\det K = \text{Tr}\big(K^{-1}\frac{\partial K}{\partial w_e}\big).

Hence

\displaystyle \frac{\partial}{\partial w_e}\log{\text{Pf}}(K) = \frac12 {\text{Tr}}\big(K^{-1}\frac{\partial K}{\partial w_e} \big).

We now compute this trace explicitly. Since only the entries K_{ij} and K_{ji} depend on w_e , we have

\displaystyle \frac{\partial K_{ij}}{\partial w_e}=1, \quad \frac{\partial K_{ji}}{\partial w_e}=-1,

and all other entries of \partial K/\partial w_e are zero. Let D_e=\frac{\partial K}{\partial w_e}. Then (D_e)_{ij}=1, (D_e)_{ji}=-1, and all other entries vanish. The trace is

\displaystyle \text{Tr}(K^{-1}D_e)=\sum_{a,b}(K^{-1})_{ab}(D_e)_{ba}.

Since D_e is nonzero only at (i,j) and (j,i) , only two terms survive:

\displaystyle \text{Tr}(K^{-1}D_e) =(K^{-1})_{ji}(D_e)_{ij} + (K^{-1})_{ij}(D_e)_{ji}.

Using (D_e)_{ij}=1, (D_e)_{ji}=-1, we get

\displaystyle \text{Tr}(K^{-1}D_e) =(K^{-1})_{ji} + (K^{-1})_{ij}.

Because K is skew-symmetric, its inverse is also skew-symmetric: (K^{-1}){ij}=-(K^{-1}){ji}. Therefore

\displaystyle \text{Tr}(K^{-1}D_e) =(K^{-1})_{ji} + (K^{-1})_{ji} =2(K^{-1})_{ji}.

Substituting this into the derivative of the Pfaffian gives

\displaystyle \frac{\partial}{\partial w_e} \log\text{Pf}(K)=(K^{-1})_{ji}.

Finally we obtain

\displaystyle \mathbb P(e\in M) = K_{ij}(K^{-1})_{ji}.

This is one of the main formulas in the theory. It says that a one-edge probability is obtained by multiplying the Kasteleyn matrix entry on that edge by the corresponding reversed entry of the inverse Kasteleyn matrix. One should not think of K^{-1} as a probability matrix by itself. Its entries may be negative or complex, depending on the orientation convention. The probability appears only after combining the inverse entry with the original Kasteleyn weight. The product K_{ij}(K^{-1})_{ji} is the invariant probabilistic quantity. This formula also explains why the inverse Kasteleyn matrix is central. Once K^{-1} is known, all one-edge densities are known.

We now pass to several edges.

Let e_r=\{i_r,j_r\},\quad r=1,\dots,k, be pairwise disjoint edges. Suppose each e_r is oriented as i_r\to j_r. Let U=\{i_1,j_1,\dots,i_k,j_k\} be the set of all vertices incident to these edges. We want the probability that all the edges e_1,\dots,e_k occur in the random matching:

\displaystyle \mathbb P(e_1,\dots,e_k\in M).

If these edges are forced to appear, then their endpoints are already matched. The remaining problem is to match all vertices not in U . Thus the remaining contribution is the partition function of the graph with those vertices deleted. Therefore

\displaystyle \mathbb P(e_1,\dots,e_k\in M) = \frac{ \big(\prod_{r=1}^k w_{e_r}\big) Z_{G\setminus U}}{Z_G}.

This is the direct probabilistic formula. The numerator says: first pay the weights of the forced edges, then sum over all possible matchings of the remaining graph. By Kasteleyn’s theorem, Z_G=|\text{Pf}(K)|. Similarly, Z_{G\setminus U} =|\text{Pf}(K_{\widehat U})|, where K_{\widehat U} is the skew-symmetric matrix obtained by deleting the rows and columns indexed by U . Thus the probability becomes a ratio of Pfaffians:

\displaystyle \mathbb P(e_1,\dots,e_k\in M) = \big(\prod_{r=1}^k w_{e_r}\big) \frac{|\text{Pf}(K_{\widehat U})|}{|\text{Pf}(K)| }.

The ratio can be expressed using the inverse matrix. This is the Pfaffian analogue of the Jacobi complementary minor identity for determinants. The Pfaffian Jacobi identity says that if K is invertible and U is an even subset of indices, then

\displaystyle \text{Pf}(K_{\widehat U}) = \pm \text{Pf}(K) \text{Pf}\big((K^{-1})_U\big).

Here (K^{-1})_U denotes the submatrix of K^{-1} formed by keeping only the rows and columns indexed by U . The sign depends only on the order in which the indices of U are listed. Once an ordering convention is fixed, this sign is fixed and can be absorbed into the same orientation convention used for the forced edges. Substituting the Pfaffian Jacobi identity into the probability formula gives

\displaystyle \mathbb P(e_1,\dots,e_k\in M) = \big(\prod_{r=1}^k w_{e_r}\big) \big | \text{Pf}\big((K^{-1})_U\big)\big |.

If we order the indices in U as i_1,j_1,i_2,j_2,\dots,i_k,j_k, and keep track of the signs consistently, then the formula may be written without absolute values as

\displaystyle \mathbb P(e_1,\dots,e_k\in M) = \big (\prod_{r=1}^{k} K_{i_r j_r} \big ) {\text{Pf}} ( (K^{-1})_U).

This is the finite-dimensional Pfaffian formula for dimer correlations. Let us check that this agrees with the one-edge formula. If k=1 , then U=\{i,j\}. the formula reduces to

\displaystyle \mathbb P(e\in M)=K_{ij}(K^{-1})_{ji},

which is the formula obtained by logarithmic differentiation. For two edges, the formula becomes more concrete. Suppose e_1=\{i_1,j_1\},  e_2=\{i_2,j_2\}. Then

\displaystyle \mathbb P(e_1,e_2\in M) = K_{i_1j_1}K_{i_2j_2} \text{Pf} \big( (K^{-1})_{{i_1,j_1,i_2,j_2}}\big).

The Pfaffian of a 4\times4 skew matrix has three terms, so this probability is explicitly a quadratic expression in entries of K^{-1} . This is the beginning of the general Pfaffian point process structure. The terminology Pfaffian point process means precisely that every joint inclusion probability is a Pfaffian built from a fixed kernel. In the dimer model, the kernel is essentially the inverse Kasteleyn matrix. Therefore all local statistics reduce to the analytic study of K^{-1} .

This is the main transition from combinatorics to analysis:

\displaystyle \text{partition function} \quad\longrightarrow\quad \text{Pf}(K),

\displaystyle \text{edge probabilities} \quad\longrightarrow\quad K^{-1},

\displaystyle \text{correlations} \quad\longrightarrow\quad \text{Pfaffians of submatrices of }K^{-1}.

Thus, after Kasteleyn’s theorem, the statistical mechanics of dimers is governed by the inverse of a discrete skew symmetric operator.

Bipartite Graphs

Now suppose that G is bipartite. Thus the vertex set is divided into two classes, V=B\sqcup W, and every edge joins a black vertex to a white vertex. There are no edges between two black vertices and no edges between two white vertices. A perfect matching of a bipartite graph pairs every black vertex with exactly one white vertex, and every white vertex with exactly one black vertex. Therefore, if a perfect matching exists, the two vertex classes must have the same size. We write |B|=|W|=n. The bipartite case is especially important because the Pfaffian formula simplifies to an ordinary determinant. This is not just a notational convenience. It is the reason that bipartite dimer models become determinental point processes, and it is the reason Fourier analysis becomes particularly clean on periodic bipartite graphs.

Order the vertices so that all black vertices come first and all white vertices come second: B=\{b_1,\dots,b_n\},~W=\{w_1,\dots,w_n\}. Since there are no black-black or white-white edges, the skew-symmetric Kasteleyn matrix has the block form

\displaystyle K= \begin{pmatrix} 0&A\\ -A^T&0 \end{pmatrix}.

Here A is an n\times n matrix whose rows are indexed by black vertices and whose columns are indexed by white vertices. Its entry A_{bw} is the signed weighted adjacency entry between b and w . More explicitly,

\displaystyle A_{bw} = \begin{cases} \pm w_{bw}, & \text{if } b \text{ and } w \text{ are adjacent},\\ 0, & \text{otherwise}. \end{cases}

The sign depends on the chosen Kasteleyn orientation. If the oriented edge agrees with the convention from black to white, one gets +w_{bw} ; otherwise one gets -w_{bw} . In many periodic examples one also allows complex signs, such as i , because this can make the Fourier symbol simpler. The underlying principle is the same: A is the signed weighted bipartite adjacency matrix.

The Pfaffian of the block matrix is related to the determinant of A by the standard identity

\displaystyle \text{Pf} \begin{pmatrix} 0&A\\ -A^T&0 \end{pmatrix} =(-1)^{n(n-1)/2}\det A.

Let us briefly explain where the sign comes from. The Pfaffian pairs the 2n indices. Since the upper-left and lower-right blocks are zero, a nonzero Pfaffian term can only pair black indices with white indices. Thus a nonzero term is exactly a bijection B\to W. That is precisely the kind of object summed over by the determinant of A . The Pfaffian expansion gives the same products, the only difference is the universal sign caused by the convention that all black vertices were listed first and all white vertices second. This universal sign is (-1)^{n(n-1)/2}. Since the dimer partition function uses an absolute value, this global sign is irrelevant. Therefore Kasteleyn’s formula becomes

\displaystyle Z_G = |\text{Pf}(K)| =|\det A|.

Thus, in the bipartite planar case, the Pfaffian solution reduces to the determinant formula

\displaystyle  Z_G=|\det A|.

This is the first major simplification. The second major simplification concerns probabilities. In the general planar case, edge correlations are Pfaffians of submatrices of K^{-1} . In the bipartite case, they become determinants of submatrices of A^{-1} . Let e=(b,w) be an edge, with b\in B and w\in W . The one-edge probability is

\displaystyle \mathbb P(e\in M) = A_{bw}(A^{-1})_{wb}.

Here one must pay attention to the order of indices. The matrix A has black rows and white columns, so A_{bw} is natural. But A^{-1} maps in the opposite direction: its rows are indexed by white vertices and its columns by black vertices. Thus the inverse entry is written (A^{-1})_{wb}.

This is the bipartite analogue of the general formula. Again, A_{bw} and (A^{-1})_{wb} separately may have signs or complex phases, depending on the Kasteleyn convention. But their product is the real positive probability that the edge appears in the matching.

Now take several pairwise disjoint edges. The probability that all these edges occur is

\displaystyle \mathbb P(e_1,\dots,e_k\in M)=\frac{ \big(\prod_{r=1}^{k} A_{b_rw_r}\big)\det A_{\widehat B,\widehat W}}{\det A},

up to the same harmless global sign convention. Here A_{\widehat B,\widehat W} is the matrix obtained by deleting the black rows b_1,\dots,b_k and the white columns w_1,\dots,w_k . The determinant Jacobi identity says that complementary minors of A are expressed by minors of A^{-1} . In the present setting, this gives

\displaystyle \frac{\det A_{\widehat B,\widehat W}}{\det A} = \pm \det\big [(A^{-1})_{w_rb_s}\big ]_{r,s=1}^{k}.

With the consistent ordering convention, the sign is absorbed into the product of the signed edge weights. Therefore

\displaystyle \mathbb P(e_1,\dots,e_k\in M) = \big (\prod_{r=1}^{k} A_{b_rw_r}\big )\det\big [(A^{-1})_{w_rb_s}\big ]_{r,s=1}^{k}.

This is the determinantal form of the bipartite dimer correlations.

It is worth emphasizing what this says. To know the probability that the k edges (b_1,w_1),\dots,(b_k,w_k) all appear, one forms the k\times k matrix \big [(A^{-1}){w_rb_s}\big ]_{r,s=1}^{k}. The row index remembers the white endpoint of the forced edge e_r , and the column index remembers the black endpoint of the forced edge e_s . Taking the determinant of this matrix and multiplying by the edge weights gives the joint probability. This is why the bipartite dimer model is called determinantal. Every finite joint edge probability is a determinant built from the same kernel A^{-1} .

Let us now write out the two-edge case explicitly, because this is the formula that controls correlations. Let e=(b,w), f=(b',w') be two disjoint edges. Then the two-edge probability is

\displaystyle \mathbb P(e,f\in M) = A_{bw}A_{b'w'} \det \begin{pmatrix} (A^{-1})_{wb} & (A^{-1})_{wb'}\\(A^{-1})_{w'b} & (A^{-1})_{w'b'} \end{pmatrix}.

On the other hand,

\displaystyle \mathbb P(e\in M)\mathbb P(f\in M) = A_{bw}A_{b'w'} (A^{-1})_{wb}(A^{-1})_{w'b'}.

Subtracting this from the two-edge probability, we obtain

\displaystyle \text{Cov}(\mathbf 1_e,\mathbf 1_f)= A_{bw}A_{b'w'} (A^{-1})_{wb'}(A^{-1})_{w'b}.

This is one of the central formulas of dimer statistics. It says that the connected two-point function factors into a product of two inverse Kasteleyn entries. Consequently, the decay of dimer-dimer correlations is governed by the decay of A^{-1} . If (A^{-1})(x,y)=O(|x-y|^{-\alpha}), then the covariance of two distant edge occupation variables typically decays like O(|x-y|^{-2\alpha}).

For the critical square-lattice dimer model, one has A^{-1}(x,y)=O(|x-y|^{-1}), and therefore

\displaystyle \text{Cov}(\mathbf 1_e,\mathbf 1_f)=O(|x-y|^{-2}).

Thus the determinant reduction does much more than simplify notation. It turns the entire local statistical theory of bipartite dimers into the analysis of the inverse matrix A^{-1} . In periodic models, A^{-1} is computed by Fourier transform, so the study of correlations becomes the study of a Fourier integral.

Rectangular Lattice

We now specialize the general Pfaffian method to the first completely explicit example: the ordinary rectangular square grid. This is the best place to see the calculation because the graph is planar, the Kasteleyn sign issue is real but still completely local, and the diagonalization is just the elementary sine diagonalization of the one-dimensional path. Only after this case is understood should one pass to the torus, where the same local operator appears but the topology introduces four global sign sectors.

Let R_{m,n} be the m\times n rectangle of unit squares. We regard the squares themselves as the vertices of a graph. Two vertices are adjacent if the corresponding squares share an edge. A domino tiling is then exactly a perfect matching of this graph: each domino chooses one edge, and the condition that every square is covered exactly once says that every vertex is incident to exactly one chosen edge. If horizontal dominoes have weight a>0 and vertical dominoes have weight b>0 , then the partition function is

\displaystyle Z_{m,n}(a,b)=\sum_D a^{h(D)}b^{v(D)},

where D runs over domino tilings, h(D) is the number of horizontal dominoes, and v(D) is the number of vertical dominoes. In the unweighted model one sets a=b=1 , and then Z_{m,n}(1,1) is simply the number of domino tilings. If mn is odd, then no tiling exists, since every domino covers two squares. Thus Z_{m,n}=0 in that case. From now on assume mn is even.

Color the squares like a chessboard: the square (x,y) is black if x+y is even and white if x+y is odd. Adjacent squares always have opposite colors, because moving one step horizontally or vertically changes x+y by 1 . Hence every domino covers one black square and one white square. Since mn is even, the rectangle has equally many black and white squares. We write these two sets as B and W , with |B|=|W|=mn/2. A domino tiling is now the same thing as a bijective matching from black squares to adjacent white squares.

Because the graph is bipartite, the Pfaffian reduces to a determinant. We build a matrix A with rows indexed by black squares and columns indexed by white squares. If b\in B and w\in W are not adjacent, put A_{bw}=0. If they are horizontally adjacent, put A_{bw}=a. If they are vertically adjacent, put A_{bw}=ib. Thus the matrix is

\displaystyle A_{bw}=\begin{cases} a, & b,w\text{ horizontally adjacent},\\ ib, & b,w\text{ vertically adjacent},\\ 0, & b,w\text{ not adjacent}. \end{cases}

The factor i is the whole sign device. It is not a probabilistic weight. The actual vertical domino weight is b . The complex phase i is inserted so that the determinant signs agree with the positive dimer weights. This is the concrete square-grid version of the Kasteleyn orientation. One can replace these complex signs by a real Kasteleyn orientation after a transformation, but for computation the complex convention is cleaner.

Now expand the determinant. Since the rows are black squares and the columns are white squares,

\displaystyle \det A=\sum_{\pi} {\text{sgn}}(\pi)\prod_{b\in B}A_{b,\pi(b)},

where \pi runs over all bijections from B to W. A term is nonzero exactly when every chosen pair b,\pi(b) is an adjacent black-white pair. Since \pi uses every white square exactly once, a nonzero term is exactly a domino tiling. Thus the determinant is already summing over the correct objects. The only issue is that the determinant gives signs, whereas the dimer partition function is a positive sum. More explicitly, if D is the tiling corresponding to \pi_D , then its term in the determinant is {\text{sgn}}(\pi_D)i^{v(D)}a^{h(D)}b^{v(D)}. Therefore the determinant will equal the partition function up to one global phase if the quantity {\text{sgn}}(\pi_D)i^{v(D)} is independent of D.

The smallest calculation already shows exactly why the i is needed. Consider a 2\times2 rectangle. There are two tilings: two horizontal dominoes or two vertical dominoes. With a suitable ordering of the two black and two white squares, the matrix is

\displaystyle A=\begin{pmatrix} a&ib\\ ib&a \end{pmatrix}.

The determinant is a\cdot a-(ib)(ib). Since (ib)(ib)=i^2b^2=-b^2, we get \det A=a^2+b^2. This is exactly the weighted partition function of the 2\times2 square. Without the factor i , the determinant would be a^2-b^2, which is wrong. The determinant naturally inserts a minus sign between the two possible pairings. The vertical pair contributes an additional minus sign because i^2=-1. The two signs cancel.

The same cancellation holds on the whole rectangle. To see why, compare two tilings D and D'. Their symmetric difference D\triangle D' is a disjoint union of even alternating cycles. Along such a cycle, the two tilings use alternating edges. Therefore it is enough to compare two tilings that differ only on one alternating cycle. Suppose that cycle contains r black vertices and r white vertices. Passing from one alternating choice of edges to the other changes the associated bijection B\to W by an r -cycle, and hence changes the determinant sign by (-1)^{r-1}. On the other hand, because each elementary square face has the local sign property just computed, switching across the region enclosed by the cycle changes the product of the i -phases by the same factor. Equivalently, the complex square-grid convention satisfies the Kasteleyn face condition: every elementary square is oddly signed. Since every alternating cycle in a rectangle bounds a planar region made from elementary faces, the face condition forces the cycle sign to be exactly the one needed to cancel the determinant sign. Hence the total phase {\text{sgn}}(\pi_D)i^{v(D)} is the same for every tiling D.

This is the crucial planar point. The determinant already knows the matchings. The only question is whether the signs cancel or reinforce. For the rectangle, the i on vertical edges makes every local square correct, and because the rectangle is simply connected, every alternating cycle is built from these local faces. Thus there is one global phase c , with |c|=1, such that \det A=cZ_{m,n}(a,b). Taking absolute values gives the exact identity

\displaystyle Z_{m,n}(a,b)=|\det A|.

It remains to compute this determinant.

Rather than diagonalize A directly, it is cleaner to double it into an operator on all squares. Let T be the mn\times mn matrix indexed by all squares of the rectangle, defined by putting weight a between horizontally adjacent squares and weight ib between vertically adjacent squares. Thus T is the signed adjacency matrix of the rectangular grid, with the same horizontal and vertical weights used above. If the vertices are ordered with all black squares first and all white squares second, then T has block form

\displaystyle T=\begin{pmatrix}0&A\\ A^T&0\end{pmatrix}.

If N=mn/2 , then a standard block determinant calculation gives \det T=(-1)^N(\det A)^2. Therefore |\det T|=|\det A|^2, and since Z_{m,n}=|\det A|, we have

\displaystyle Z_{m,n}(a,b)=|\det T|^{1/2}.

Now the advantage is that T separates into horizontal and vertical one-dimensional parts. Let H_m be the adjacency matrix of the path with m vertices, and let H_n be the adjacency matrix of the path with n vertices. Thus H_m has 1 immediately above and below the diagonal and zeros elsewhere. Then

\displaystyle T=aH_m\otimes I_n+ibI_m\otimes H_n.

This formula just says that horizontal moves are governed by the path matrix in the x -direction, while vertical moves are governed by the path matrix in the y -direction.

We now diagonalize the path matrix. For j=1,\dots,m, define u_j(x)=\sin\frac{j\pi x}{m+1}. These functions vanish at the artificial boundary points x=0 and x=m+1, which is exactly the Dirichlet boundary condition appropriate to an open rectangle. If \theta_j=j\pi/(m+1), then

\displaystyle u_j(x+1)+u_j(x-1)=2\cos\theta_j ~u_j(x).

This follows from the elementary identity \sin(\alpha+\beta)+\sin(\alpha-\beta)=2\sin\alpha\cos\beta. Therefore u_j is an eigenvector of H_m with eigenvalue 2\cos\frac{j\pi}{m+1}. Similarly, the eigenvectors of H_n are v_k(y)=\sin\frac{k\pi y}{n+1} with eigenvalues 2\cos\frac{k\pi}{n+1}.

The product functions u_j(x)v_k(y) form a basis of eigenvectors for T. Indeed, applying T=aH_m\otimes I_n+ibI_m\otimes H_n gives one contribution from the x -direction and one contribution from the y -direction. Thus the eigenvalue corresponding to (j,k) is

\displaystyle \mu_{jk}=2a\cos\frac{j\pi}{m+1}+2ib\cos\frac{k\pi}{n+1}.

Since these product eigenvectors form a basis, the determinant of T is the product of all these eigenvalues. Taking absolute values gives

\displaystyle |\mu_{jk}|=\big(4a^2\cos^2\frac{j\pi}{m+1}+4b^2\cos^2\frac{k\pi}{n+1}\big)^{1/2}.

Finally Z_{m,n}=|\det T|^{1/2}. Hence the exact formula is

\displaystyle  Z_{m,n}(a,b)=\prod_{j=1}^{m}\prod_{k=1}^{n}\big(4a^2\cos^2\frac{j\pi}{m+1}+4b^2\cos^2\frac{k\pi}{n+1}\big)^{1/4}.

In the unweighted case this becomes

\displaystyle  Z_{m,n}=\prod_{j=1}^{m}\prod_{k=1}^{n}\big(4\cos^2\frac{j\pi}{m+1}+4\cos^2\frac{k\pi}{n+1}\big)^{1/4}.

The exponent 1/4 has a simple origin. First, the absolute value of each complex eigenvalue produces a square root, because |X+iY|=(X^2+Y^2)^{1/2}. Second, the partition function is |\det T|^{1/2} rather than |\det T|. These two square roots combine to produce the fourth root in the final product.

The exact product immediately gives the thermodynamic limit. Take logarithms first:

\displaystyle \log Z_{m,n}(a,b)=\frac14\sum_{j=1}^{m}\sum_{k=1}^{n}\log\big(4a^2\cos^2\frac{j\pi}{m+1}+4b^2\cos^2\frac{k\pi}{n+1}\big).

After division by mn, this is a double Riemann sum. The points j\pi/(m+1) fill (0,\pi) , and the points k\pi/(n+1) fill (0,\pi). Therefore

\displaystyle \lim_{m,n\to\infty}\frac{1}{mn}\log Z_{m,n}(a,b)=\frac{1}{4\pi^2}\int_0^\pi\int_0^\pi \log\big(4a^2\cos^2\theta+4b^2\cos^2\phi\big) d\theta d\phi.

The logarithm has an integrable singularity at (\theta,\phi)=(\pi/2,\pi/2) when a,b>0. Near that point the expression inside the logarithm is comparable to a positive constant times (\theta-\pi/2)^2+(\phi-\pi/2)^2. Thus the singularity is of the form \log r in two dimensions, and \int_0^\epsilon r|\log r| dr<\infty. Hence the Riemann-sum limit is legitimate.

In the unweighted case, the free energy per square is

\displaystyle f=\frac{1}{4\pi^2}\int_0^\pi\int_0^\pi\log\big(4\cos^2\theta+4\cos^2\phi\big) d\theta d\phi.

This integral evaluates to G/\pi, where G is Catalan’s constant,

\displaystyle G=\sum_{r=0}^{\infty}\frac{(-1)^r}{(2r+1)^2}.

Thus the number of domino tilings of a large rectangle satisfies

\displaystyle Z_{m,n}(1,1)=\exp\big(\frac{G}{\pi}mn+o(mn)\big ).

Numerically, G/\pi\approx0.291560904. Thus the number of tilings grows exponentially like \exp(0.291560904,mn) up to subexponential factors.

Periodic Lattice

We now move from the rectangle to periodic boundary conditions. The local operator is almost the same, but the topology changes the finite formula. On a rectangle, every alternating cycle bounds a planar region, and the local Kasteleyn sign condition is enough to make all matching signs agree. On a torus, some cycles wind around the two fundamental directions and do not bound disks. Those cycles are the reason one Pfaffian is no longer sufficient.

Use a fundamental cell containing one black vertex and one white vertex. Write the black vertex in cell (m,n) as b_{m,n} and the white vertex as w_{m,n}. This is just a convenient bookkeeping convention for the bipartite square lattice. The Kasteleyn operator acts from functions on white vertices to functions on black vertices. With horizontal weight a and vertical weight b , define

\displaystyle (Af)(m,n)=a f(m+1,n)+a f(m-1,n)+ib f(m,n+1)+ib f(m,n-1).

In the uniform case a=b=1, this becomes

(Af)(m,n)=f(m+1,n)+f(m-1,n)+i f(m,n+1)+i f(m,n-1).

Again, the vertical factor i is a sign device, not a probability. It is chosen so that every elementary square has the correct Kasteleyn sign and so that the Fourier symbol is especially simple.

We now explain carefully why the actual torus partition function is not just one determinant. This is the first genuinely new phenomenon that appears when one passes from the rectangle to periodic boundary conditions. On a rectangle, every closed alternating cycle bounds a planar region. Therefore the local Kasteleyn sign condition around elementary faces forces all matching signs to agree. On a torus this is no longer true. There are closed cycles which do not bound disks. They go once around the torus in the horizontal direction, or once around the torus in the vertical direction, or around both. These non-contractible cycles are exactly what forces the four-Pfaffian formula.

Let the finite periodic graph be drawn on a torus. Concretely, start with an M\times N periodic square grid and identify the left side with the right side and the bottom side with the top side. A path on this graph may return to its starting point in the torus even though, if we lift the torus to the infinite plane, the lifted path has moved by some multiple of the horizontal period and some multiple of the vertical period. This is the simplest way to understand winding. A closed curve on the torus may lift to a curve in the plane whose endpoint is shifted by rM units horizontally and sN units vertically. The integers r and s are the winding numbers of the curve. For the Pfaffian sign problem we only need their parities, so we only remember them modulo 2. Equivalently, one can cut the torus open into a rectangle. The two cuts are called seams: one vertical seam where the left and right sides are identified, and one horizontal seam where the bottom and top sides are identified. A closed cycle has odd horizontal winding if it crosses the vertical seam an odd number of times. It has odd vertical winding if it crosses the horizontal seam an odd number of times. This gives the same mod-two information as the lifting picture. Thus a closed cycle has two winding parities, one horizontal and one vertical.

Now fix once and for all a reference matching M_0. This reference matching is not meant to be special probabilistically. It is just a bookkeeping device. For any other matching M, consider the symmetric difference

\displaystyle M\triangle M_0=(M\setminus M_0)\cup(M_0\setminus M).

As before, every vertex that appears in this symmetric difference has degree 2, because one incident edge comes from M and one incident edge comes from M_0. Therefore M\triangle M_0 is a disjoint union of even alternating cycles. Along each such cycle, the edges alternate between the matching M and the reference matching M_0.

On the plane, each of these cycles bounds a region, so it has no global winding. On the torus, some of the cycles in M\triangle M_0 may wind around the torus. Add up their winding numbers modulo 2. This gives a pair (p,q)\in \{0,1\} \times \{0,1\}. Here p=1 means that the total horizontal winding of M\triangle M_0 is odd, and p=0 means it is even. Similarly, q=1 means that the total vertical winding is odd, and q=0 means it is even. In the seam picture, p is the parity of the number of times the alternating cycles cross the vertical seam, and q is the parity of the number of times they cross the horizontal seam. This is what winding parity means. It is a topological label attached to the matching M relative to the reference matching M_0. There are four possible labels:

\displaystyle (0,0),\quad (1,0),\quad (0,1),\quad (1,1).

Let Z_{pq} be the positive weighted sum of all matchings whose winding parity relative to M_0 is (p,q). Thus Z_{00} is the total weight of matchings with even horizontal and even vertical winding relative to M_0, Z_{10} is the total weight of matchings with odd horizontal and even vertical winding, and so on. Since every matching belongs to exactly one of these four classes, the true partition function is

\displaystyle Z=Z_{00}+Z_{10}+Z_{01}+Z_{11}.

The point is that a single torus Pfaffian usually does not give this positive sum. It gives a signed sum in which the sign depends on the winding class. The local Kasteleyn condition still fixes all contractible cycles, just as on the rectangle. Therefore all matchings in the same winding class have the same relative sign. But the local condition cannot force agreement between different winding classes, because a non-contractible cycle does not bound a union of faces. So the Pfaffian sign can still distinguish the four classes (0,0),(1,0),(0,1),(1,1).

For a standard genus-one Kasteleyn orientation, the remaining topological sign is (-1)^{pq}. This formula says the following. The class (0,0) has sign +1, the class (1,0) has sign +1, the class (0,1) has sign +1, and the class (1,1) has sign -1. In other words, for this convention the only class that appears with the opposite sign is the class that winds oddly in both fundamental directions. This is a standard normalization of the torus Kasteleyn orientation. A different initial orientation may move the signs around, but it gives an equivalent four-Pfaffian formula after relabeling signs.

Now we introduce the four Pfaffian sectors. These should not be confused with the four winding classes. The winding classes label matchings. The Pfaffian sectors label four different Kasteleyn matrices. To obtain them, cut the torus open along the two seams. Along each seam, we are allowed to flip the signs of all Kasteleyn entries crossing that seam. We may either flip or not flip in the horizontal direction, and we may either flip or not flip in the vertical direction. Therefore there are four choices. Write them as

\displaystyle (\epsilon,\eta)\in \{0,1\}\times \{0,1\}.

Here \epsilon=1 means that we have inserted an extra minus sign across the horizontal period, and \epsilon=0 means we have not. Similarly, \eta=1 means that we have inserted an extra minus sign across the vertical period, and \eta=0 means we have not. In Fourier language, these are the periodic and antiperiodic boundary conditions. In geometric language, they are sign twists around the two fundamental cycles of the torus.

Now ask what this twist does to a matching of winding class (p,q). If \epsilon=1, every edge crossing the horizontal seam has its sign changed. A matching whose symmetric difference from M_0 crosses that seam an odd number of times picks up an extra minus sign relative to M_0. Therefore the horizontal twist contributes (-1)^{\epsilon p}. Similarly, the vertical twist contributes (-1)^{\eta q}. Thus the total twist factor is

\displaystyle (-1)^{\epsilon p+\eta q}.

Combining this twist factor with the original torus Kasteleyn sign (-1)^{pq}, the Pfaffian in sector (\epsilon,\eta) has the form

\displaystyle P_{\epsilon,\eta}=\sum_{p,q\in{0,1}}(-1)^{pq+\epsilon p+\eta q}Z_{pq}.

This formula is the whole four-Pfaffian mechanism. The variables p,q classify matchings by winding parity. The variables \epsilon,\eta classify the four boundary-sign choices of the Kasteleyn matrix. The exponent pq+\epsilon p+\eta q records the Pfaffian sign of that winding class in that sector.

Let us write out the four equations. When (\epsilon,\eta)=(0,0), there is no seam twist, so the only topological sign is (-1)^{pq}. Hence

\displaystyle P_{00}=Z_{00}+Z_{10}+Z_{01}-Z_{11}.

When (\epsilon,\eta)=(1,0), the horizontal twist changes the sign of the classes with p=1. Thus

\displaystyle P_{10}=Z_{00}-Z_{10}+Z_{01}+Z_{11}.

When (\epsilon,\eta)=(0,1), the vertical twist changes the sign of the classes with q=1. Thus

\displaystyle P_{01}=Z_{00}+Z_{10}-Z_{01}+Z_{11}.

Finally, when (\epsilon,\eta)=(1,1), both twists are present, and the signs are

\displaystyle P_{11}=Z_{00}-Z_{10}-Z_{01}-Z_{11}.

These four equations should be read as a linear system. The four Pfaffians P_{00},P_{10},P_{01},P_{11} are four signed measurements of the four positive quantities Z_{00},Z_{10},Z_{01},Z_{11}. The actual partition function is the sum of those four positive quantities. We now solve for that sum. Add the first three Pfaffians and subtract the fourth:

\displaystyle P_{00}+P_{10}+P_{01}-P_{11}.

The coefficient of Z_{00} is 1+1+1-1=2. The coefficient of Z_{10} is 1-1+1-(-1)=2. The coefficient of Z_{01} is 1+1-1-(-1)=2. The coefficient of Z_{11} is -1+1+1-(-1)=2. Therefore

\displaystyle P_{00}+P_{10}+P_{01}-P_{11}=2(Z_{00}+Z_{10}+Z_{01}+Z_{11}).

Dividing by 2 gives the exact torus partition function:

\displaystyle  Z=\frac12\big(P_{00}+P_{10}+P_{01}-P_{11}\big).

Equivalently, since the coefficient in front of P_{\epsilon,\eta} is (-1)^{\epsilon\eta}, one can write the same formula more compactly as

\displaystyle  Z=\frac12\sum_{\epsilon,\eta\in{0,1}}(-1)^{\epsilon\eta}P_{\epsilon,\eta}.

This is the four-Pfaffian formula. It is important to understand what it is saying. On the plane, one Kasteleyn Pfaffian is enough because all alternating cycles are contractible and the face sign condition controls them. On the torus, alternating cycles can have four possible winding parities. One Pfaffian gives one signed combination of these four classes. By changing the boundary signs in the two fundamental directions, we get four different signed combinations. The above linear combination cancels the unwanted signs and gives every winding class coefficient +1. Thus the four Pfaffians are exactly the correction for the topology of the torus.

In the bipartite square-lattice case, each P_{\epsilon,\eta} is, up to a harmless global sign depending on the ordering convention, the determinant of the corresponding bipartite Kasteleyn matrix A^{\epsilon,\eta}. Thus the finite torus computation is still done by Fourier products. What changes is that the correct finite partition function is not one product; it is the above signed combination of four Pfaffians, or equivalently four determinant sectors with the appropriate Pfaffian signs.

Now let us connect this topological discussion with the Fourier calculation. In sector (\epsilon,\eta), the functions satisfy

\displaystyle f(m+M,n)=(-1)^\epsilon f(m,n),\qquad f(m,n+N)=(-1)^\eta f(m,n).

Thus \epsilon=0 means periodic and \epsilon=1 means antiperiodic in the first direction; similarly \eta=0 and \eta=1 mean periodic and antiperiodic in the second direction. The allowed Fourier modes are

\displaystyle f_{j,k}(m,n)=\exp(i\theta_j^{(\epsilon)}m+i\phi_k^{(\eta)}n),

where

\displaystyle \theta_j^{(\epsilon)}=\frac{(2j+\epsilon)\pi}{M},\qquad \phi_k^{(\eta)}=\frac{(2k+\eta)\pi}{N}.

These formulas are just the equations e^{i\theta M}=(-1)^\epsilon and e^{i\phi N}=(-1)^\eta. Therefore the antiperiodic sector simply shifts the allowed momenta by half a lattice spacing. Applying the square-lattice Kasteleyn operator to such a mode gives the eigenvalue

\displaystyle \lambda_{j,k}^{\epsilon,\eta}=2a\cos\frac{(2j+\epsilon)\pi}{M}+2ib\cos\frac{(2k+\eta)\pi}{N}.

Hence in a fixed sector,

\displaystyle \det A^{\epsilon,\eta}=\prod_{j=0}^{M-1}\prod_{k=0}^{N-1}\big(2a\cos\frac{(2j+\epsilon)\pi}{M}+2ib\cos\frac{(2k+\eta)\pi}{N}\big).

This determinant is the Fourier product for one Pfaffian sector. The exact finite torus partition function remembers the Pfaffian signs and combines the four sectors as above. For asymptotic purposes it is often useful to introduce the positive quantity

\displaystyle \mathcal Z_{\epsilon,\eta}=|\det A^{\epsilon,\eta}|=\prod_{j=0}^{M-1}\prod_{k=0}^{N-1}\big(4a^2\cos^2\frac{(2j+\epsilon)\pi}{M}+4b^2\cos^2\frac{(2k+\eta)\pi}{N}\big)^{1/2}.

One should be careful here. The exact finite formula is a signed Pfaffian formula, not simply the same expression with absolute values inserted everywhere. The absolute values are harmless only when one is studying the leading exponential growth, because all four sectors have the same thermodynamic limit. Their momentum grids differ only by shifts of size one half of a mesh, and after division by the area those shifts disappear.

Consider one sector. Taking logarithms of the absolute determinant and dividing by MN gives

\displaystyle \frac{1}{MN}\log|\det A^{\epsilon,\eta}|=\frac{1}{MN}\sum_{j=0}^{M-1}\sum_{k=0}^{N-1}\log\big |2a\cos\frac{(2j+\epsilon)\pi}{M}+2ib\cos\frac{(2k+\eta)\pi}{N}\big |.

As M,N\to\infty, the allowed momenta become equidistributed on the Fourier torus. Thus the limit is

\displaystyle f_{\mathrm{cell}}(a,b)=\frac{1}{(2\pi)^2}\int_0^{2\pi}\int_0^{2\pi}\log|2a\cos\theta+2ib\cos\phi| d\theta d\phi.

This is the free energy per fundamental cell. In the convention used here, one fundamental cell contains one black vertex and one white vertex. Therefore the free energy per graph vertex is one half of f_{\mathrm{cell}}. In the unweighted case,

\displaystyle f_{\mathrm{cell}}=\frac{1}{(2\pi)^2}\int_0^{2\pi}\int_0^{2\pi}\log|2\cos\theta+2i\cos\phi | d\theta d\phi.

Since | 2\cos\theta+2i\cos\phi |=2(\cos^2\theta+\cos^2\phi)^{1/2}, this integral evaluates to

\displaystyle f_{\mathrm{cell}}=\frac{2G}{\pi}.

Thus the entropy per graph vertex is

\displaystyle f_{\mathrm{vertex}}=\frac{G}{\pi}.

This agrees with the rectangle calculation, because the rectangle free energy above was measured per square of the board, and those squares are exactly the graph vertices in the domino-tiling representation.

The inverse Kasteleyn kernel

The partition function is governed by the determinant or Pfaffian of the Kasteleyn matrix, but local statistics are governed by the inverse Kasteleyn matrix. On the infinite square lattice, the operator is translation-invariant, so Fourier transform turns A into multiplication by its symbol P(e^{i\theta},e^{i\phi}). Therefore the inverse operator is obtained by multiplying by 1/P in Fourier space and transforming back. The inverse entry at displacement (m,n) is

\displaystyle A^{-1}(m,n)=\frac{1}{(2\pi)^2} \int_0^{2\pi}\int_0^{2\pi}\frac{e^{-i(m\theta+n\phi)}}{2\cos\theta+2i\cos\phi }  d\theta d\phi.

This is just the Fourier inversion formula. The large-distance behavior is controlled by the singularities of 1/P. In the square-lattice case, P vanishes at the four points (\pi/2,\pi/2),(\pi/2,3\pi/2),(3\pi/2,\pi/2),(3\pi/2,3\pi/2). Near the zero (\pi/2,\pi/2), write \theta=\pi/2+p and \phi=\pi/2+q. Since \cos(\pi/2+p)=-p+O(p^3) and similarly for \phi, we have

\displaystyle P(e^{i\theta},e^{i\phi})=-2p-2iq+O(|p|^3+|q|^3).

Thus the singular part of the inverse multiplier has the local form -\frac12(p+iq)^{-1}. A two-dimensional Fourier integral with such a Cauchy-type singularity decays like reciprocal distance. Consequently

\displaystyle A^{-1}(m,n)=O\left((m^2+n^2)^{-1/2}\right).

By the determinantal formula for bipartite dimers, the covariance of two distant edge occupation variables is a product of two inverse Kasteleyn entries. Hence, for two edges separated by distance r,

\displaystyle {\text{Cov}}(\mathbf 1_e,\mathbf 1_f)=O(r^{-2}).

This is the basic critical exponent of the square-lattice dimer model. The determinant gives the partition function, the inverse determinant kernel gives probabilities, the zeros of the Fourier symbol give the singularities of the inverse kernel, and those singularities give the power-law decay.

Analytic phases: liquid, gaseous, frozen

For a general periodic bipartite dimer model, the central analytic object is the characteristic polynomial. If the fundamental domain has several black and white vertices, the Fourier-transformed Kasteleyn operator is a finite matrix A(z,u), and one defines

\displaystyle P(z,u)=\det A(z,u).

This Laurent polynomial is the characteristic polynomial of the periodic dimer model. The square lattice is the simplest case, where the fundamental domain contains one black vertex and one white vertex, and the characteristic polynomial is simply P(z,u)=z+z^{-1}+i(u+u^{-1}).

The corresponding spectral curve is

\displaystyle P(z,u)=0.

The large-distance behavior of the dimer model is governed by how this curve intersects the unit torus |z|=|u|=1. This is the analytic origin of the usual distinction between liquid, gaseous, and frozen regimes.

In the liquid phase, the characteristic polynomial has zeros on the unit torus. Generically these zeros are simple. Since the inverse Kasteleyn matrix is obtained by Fourier inversion of 1/P, zeros of P on the unit torus become singularities of the Fourier integrand. Those singularities prevent exponential decay. Instead, the inverse Kasteleyn matrix has power-law decay. In the square-lattice case, the zero is first order and the local form of 1/P is Cauchy-like, so

\displaystyle A^{-1}(x,y)=O(|x-y|^{-1}).

The determinantal formula for bipartite dimers then gives edge-edge covariance as a product of two inverse Kasteleyn entries. Hence, for two distant edges e and f,

\displaystyle {\text{Cov}}(\mathbf 1_e,\mathbf 1_f)=O(|x-y|^{-2}).

This is the critical regime. Correlations decay polynomially, not exponentially. The uniform square-lattice dimer model is the standard example.

In the gaseous phase, the characteristic polynomial has no zeros on the unit torus. Then 1/P is analytic in a complex neighborhood of the real Fourier torus. Analytically, this is the massive situation. Since the Fourier integrand is analytic, one can shift the contour slightly into the complex domain. The oscillatory factor then gains exponential decay. Consequently,

\displaystyle A^{-1}(x,y)=O(e^{-c|x-y|})

for some c>0. Edge correlations also decay exponentially. Thus the gaseous phase is noncritical: it has a finite correlation length, and distant local events become independent exponentially fast.

In the frozen phase, the Gibbs measure becomes nearly deterministic at large scale. One or more edge types dominate, and local randomness disappears. In slope language, frozen phases occur at the boundary of the Newton polygon of P. In such a region, changing the microscopic configuration costs too much entropy or weight, and the tiling locks into a rigid pattern.

The whole story can be read as one chain of ideas:

\displaystyle \text{planarity}\Rightarrow\text{Kasteleyn signs}\Rightarrow\text{Pfaffian partition function}\Rightarrow\text{inverse Kasteleyn kernel}\Rightarrow\text{Fourier symbol}\Rightarrow\text{phases and limit shapes}.

This is the analytic heart of the dimer model. The same algebraic object that counts matchings gives the free energy. The same inverse matrix that gives edge probabilities gives correlation decay. The same Laurent polynomial that appears in Fourier diagonalization determines the liquid, gaseous, and frozen regimes. Exact solvability is not a single trick; it is the coherence of this entire structure.

Leave a comment