How can I derive the analog of the susceptibility sum rule for the specific heat? How can I derive the analog of the susceptibility sum rule for the specific heat? Does an infinite correlation length imply an infinite specific heat?
$$
\chi = \frac{\partial M}{\partial H}
= \frac{1}{N}\sum_{i,j} \Gamma(i,j)
$$
 A: I'll restrict my answer to the nearest neighbour Ising model with $N$ spins $s_i=\pm1$,
for which the energy and magnetization are given by 
$$
E=-J\sum_{\langle j,k\rangle} s_j s_k, \qquad M=\sum_j s_j
$$
and the notation $\langle j,k\rangle$ denotes nearest-neighbour pairs of spins.
Recapping, the susceptibility formula comes from standard linear response theory
$$
\chi = \frac{\partial M}{\partial H} 
\propto \langle M^2\rangle - \langle M\rangle^2
$$
where $\chi$ is the susceptibility of the whole system.
Inserting the definition of $M$
$$
\chi 
\propto \sum_j\sum_k \left(\langle s_j s_k\rangle - \langle s_j\rangle\langle s_k\rangle \right)
\propto
N \sum_k \left(\langle s_0 s_k\rangle - \langle s_0\rangle\langle s_k\rangle \right)
$$
where we use translational invariance to replace one of the sums by a factor $N$;
now $j$ has become an arbitrarily chosen spin $0$ which we can take as
the origin of coordinates.
The factor of $N$ gives the expected extensive behaviour of $\chi$,
and we can define the quantity inside the sum as the 
spin-spin correlation function $c$,
which is expected to depend only on the 
vector between the spins: $c(\mathbf{r}_k-\mathbf{r}_0)=c(\mathbf{r}_k)$.
So the susceptibility per spin is
$$
\frac{\chi}{N} 
\propto
\sum_{\mathbf{r}} c(\mathbf{r})
$$
where now we sum over all vectors from a lattice site at the origin
to all other lattice sites.
We can approximate the sum as an integral,
and can often assume that $c$ has a finite range
with a correlation length $\xi$, so $c\sim \exp(-r/\xi)$,
and the integral will give a finite result even if we let $N\rightarrow\infty$.
However, if the correlation length diverges,
for instance near the critical point,
the integrand will not decay (fast enough) with distance $r$,
the integral will diverge,
and we expect a diverging susceptibility.
Formally we can go through a similar procedure for the heat capacity.
Now
$$
C_V = \frac{\partial E}{\partial T} 
\propto \langle E^2\rangle - \langle E\rangle^2
$$
where $C_V$ is the heat capacity of the whole system.
This will be related to a correlation function of a different variable.
Define, for each bond $b$ between a nearest neighbour pair $\langle j,k\rangle$,
the quantity
$$
\varepsilon_b = -J s_j s_k
$$
so that
$$
E = \sum_b \varepsilon_b
$$
The derivation follows exactly the same pattern,
and we end up with the heat capacity per spin
$$
\frac{C_V}{N}
\propto
\sum_b \left(\langle \varepsilon_0 \varepsilon_b\rangle - \langle \varepsilon_0\rangle\langle \varepsilon_b\rangle \right)
\propto
\sum_{\mathbf{r}} c'(\mathbf{r})
$$
We used translational invariance to sum over all the bond vectors
that we could use as an origin: that's equal to the number of spins $N$
multiplied by $q/2$ where $q$ is the coordination number of the lattice,
and this $q/2$ factor has been absorbed into the proportionality constant.
The remaining sum is over all vectors connecting the centres of bonds,
to the centre of the arbitrary bond that we have chosen as origin.
This again can be treated as an integral.
The correlation function $c'(\mathbf{r})$ is different from $c(\mathbf{r})$,
but again we expect it to have a characteristic correlation length $\xi'$,
and if this diverges, we expect to see a divergent $C_V$.
Naturally,
near a critical point,
the same phenomenon is giving rise to all these divergences.
So we expect to see critical exponents for the correlation length(s),
for $\chi$, and for $C_V$,
which are all related to each other.
That's the area covered by the scaling hypothesis,
which is another story.
For a more general physical model, the analysis may not be so direct,
but the idea will be the same.
Provided the interactions are short range, it should be possible
to define a "local" energy, or energy density,
which will have a correlation function in space.
The magnitude of the total energy fluctuations will tend to diverge
if the corresponding correlation length diverges,
and hence the heat capacity will diverge.
I'm not saying that all second-order phase transitions
are accompanied by a diverging heat capacity,
but this is what is expected for the Ising universality class
on approach to a
critical point.
