Expansion of Onsager's Exact Partition Function for 2D Ising Model We have a question where we are given the exact expression for the 2D Ising model partition function:
$$\frac{1}{N}\ln Z ~=~ \ln(2 \cosh^2(\beta J)) $$$$+ \frac{1}{2}\int_{-\pi}^{\pi}\frac{dq_1}{2\pi}\int_{-\pi}^{\pi}\frac{dq_2}{2\pi} \ln \left( (1+x^2)^2 - 2x(1-x^2)(\cos(q_1) + \cos(q_2)\right),$$
where $x=\tanh(\beta J)$. In the first parts of the question we show that the argument in the logarithm is non-negative, and that there is a critical point at $q_1 = q_2 = 0$ and $x=\sqrt{2}-1$ where the argument is zero. We are then asked to:
"Expand the logarithm about both $q_1=q_2=0$ and this critical value of x and evaluate the resulting integral to extract the leading singular behaviour of $\ln Z$ (and hence the free energy) near the transition. What is the nature of the singularity in the heat capacity?
Now I'm not really sure how to go about expanding this. Clearly a standard multivariate taylor expansion is not going to be of any use here as $\ln(0)$ goes to $-\infty$. I thought maybe just do the expansion with $q_1$ and $q_2$ and then evaluate the integral to start but then I don't know what to do from there really. How does one generally go about expanding around a singularity in a function? Any pointers in the right direction would be much appreciated!
 A: If we expand the argument of $\ln$ in $q_1, q_2$ around $\left(q_1,q_2\right)=\left(0,0\right)$ and set $q_1,q_2$ to zero, we find $$(1+x^2)^2-2x(1-x^2)2=(1+x^2)^2-4x(1-x^2)=(x^2+2x-1)^2$$ which is obviously minimized at the root of $x^2+2x-1$, giving $x_c = \sqrt{2}-1$. 
Because the model is taken in the $N\rightarrow\infty$ limit, the integral for $F/N$ will be dominated by the neighborhood of $x_c$. 
Let's now expand the argument with respect to $x$ about $x_c$, i.e., $x\rightarrow x_c + \epsilon$. The expansion yields $$8\epsilon^2+4\sqrt{2}\epsilon^3 + \epsilon^4$$ which we will cut off at order $\epsilon^2$ since $\epsilon$ is taken to be small. 
Thus, once we put back the terms from non-zero $q_1,q_2$, the argument becomes $$8\epsilon^2 - k\left(q_1^2+q_2^2\right)$$ where $k$ is a constant resulting from the expansion of $2x(1-x^2)$ (the multiplier on the cosine terms). From here I'll rescale things so that there are no pointless constants in front of our arguments.
Now have $$F/N \sim \int dq_1dq_2 \ln\left(\epsilon^2 - q_1^2+q_2^2\right)$$
Change to polar coordinates so that $q_1^2+q_2^2 \rightarrow r^2$ and $dq_1dq_2\rightarrow rdrd\theta$ yields
$$\int \ln\left(\epsilon^2 - r^2\right)rdrd\theta = 2\pi\int \ln\left(\epsilon^2 - r^2\right)rdr$$
This integral is somewhat less daunting than the original (you can definitely evaluate it in Mathematica!) and gives
$$F/N = \epsilon^2\ln\epsilon - \ln\left(1-\epsilon^2\right)+\ln\left(\epsilon^2-1\right)$$
On taking the second derivative ($C\sim N^{-1}\frac{\partial^2F}{\partial\epsilon^2}$) we find that the singular term comes from $\epsilon^2\ln \epsilon$ which upon differentiating yields
$$\epsilon^2 \ln \epsilon \stackrel{\partial/\partial\epsilon}{\rightarrow} \epsilon\ln\epsilon + \epsilon \stackrel{\partial/\partial\epsilon}{\rightarrow} \ln\epsilon + \textrm{constants} \approx \ln \epsilon = \ln \left(x-x_c\right)$$
which shows that the specific heat diverges logarithmically near $x_c$
