How are the coefficients determined in the high temperature expansion of the 2D Ising model? I have been studying the 2D Ising model lately and have been looking at high and low temperatures. But I'm having problems when trying to understand the high temperature one. The final expansion looks like this:
$$Z =(\cosh K)^{2N}2^{N}\sum (\tanh K)^{l}$$
with $$K = \beta J$$
I understand the part inside the sigma sums for all the possible closed loops, with $l$ being the length of the loop. When computing the expansion to the 8th order of $l$ the answer is (I'll use $\tanh K = \epsilon$ ):
 $$Z =(\cosh K)^{2N}2^{N}(1+ N\epsilon^4+2N\epsilon^6+N\frac{N+9}{2}\epsilon^8 + ...)$$
What I don't understand is how are the closed loops in the lattice counted. 
 A: I am certainly not the person on Physics SE with most expertise concerning lattice models, but since nobody has offered an answer yet, here is mine.
As you have indicated, the partition function can be expressed as a high-temperature expansion involving closed loops or polygons (of nearest-neighbour interaction terms) on the square lattice. The essential point is that all the diagrams that do not involve an even number of lines at each vertex will cancel out to zero, once the summation over spins is carried out. The loops can be categorized according to number of lines: larger loops contribute at a higher order in the expansion.
The counting of loops is explained clearly here and probably in several textbooks as well. I'm just going to reproduce some of the material in the table on p9 of that link (in case it disappears in future). For the smaller loops, you can do the counting "by hand". For larger ones, it's a numerical exercise best tackled on a computer.
The smallest loops are of order 4 and 6. Here they are:

For an $N$-spin system, assuming periodic boundaries so there are no "edge effects" which might restrict the placement of the loop, there are $N$ possibilities for the location of the 4-line loop. Just consider the number of options for the bottom left corner, for instance. This is the coefficient of the $\epsilon^4$ term. For the 6-line loop, there are again $N$ possible locations, but also it may take either of two orientations (horizontal or vertical). So the coefficient of the $\epsilon^6$ term is $2N$.
For the 8-line loops, there are several arrangements. Here are two of them.

The one on the left is two squares. Having placed the first one in any of the allowed $N$ positions, the second one can be placed in any of $N-5$ remaining positions. The $5$ excluded ones are directly on top of the first one, and in $4$ adjacent positions (north, south, east and west). Then there is a factor $2$ to account for the fact that the two loops are identical. So this gives a contribution $N(N-5)/2$. For the one on the right, there are $2N$ possibilities, similar to the 6-line loop seen before. There are still two more shapes to consider:

The one on the left has $N$ possible locations, but for each one there are $4$ possible orientations, so the contribution is $4N$. The one on the right just has $N$ possible locations. Adding all these up gives the coefficient of $\epsilon^8$ as
$$
\frac{N(N-5)}{2} + 2N + 4N + N = \frac{N(N+9)}{2}
$$
The table in the linked chapter also gives the calculation for $\epsilon^{10}$, but I'll stop here.
