# Coarse graining procedure for the Ising model

To make things more clear, let us see how the coarse graining procedure works for the Ising model. If we call ${\displaystyle m_{i}=\left\langle S_{i}\right\rangle }$ the local magnetization at the ${\displaystyle i}$-th site and ${\displaystyle d}$ the dimensionality of the system, every "block" will have volume ${\displaystyle \ell ^{d}}$; we define for every block of the system centered in ${\displaystyle {\vec {r}}}$ the coarse grained magnetization as:

${\displaystyle m_{\ell }({\vec {r}})={\frac {1}{N_{\ell }}}\sum _{i\in {\vec {r}}}m_{i}}$
where ${\displaystyle N_{\ell }=(\ell /a)^{d}}$ is the number of spins (degrees of freedom in general) which belong to the block centered in ${\displaystyle {\vec {r}}}$; this definition is reasonable as long as ${\displaystyle N_{\ell }}$ is large. Since it has been built as an average, ${\displaystyle m_{\ell }}$ does not fluctuate much on microscopic scales but varies smoothly in space. Of course, in general we need to specify ${\displaystyle \ell }$ in order to determine ${\displaystyle m_{\ell }}$, but the coarse graining procedure we are applying will be useful only if the final results are independent of ${\displaystyle \ell }$ (at least in the spatial scales considered).

We now must express the partition function in terms of ${\displaystyle m_{\ell }({\vec {r}})}$, and as we have stated before:

${\displaystyle Z=\int {\mathcal {D}}[m_{\ell }({\vec {r}})]e^{-\beta {\mathcal {H}}_{\text{eff.}}[m_{\ell }({\vec {r}})]}}$
so we must now compute ${\displaystyle {\mathcal {H}}_{\text{eff.}}}$. Since we now have a system made up of "blocks" this effective Hamiltonian will be composed of two parts: a bulk component relative to the single blocks and an interface component relative to the interaction between the blocks; let us consider them individually.

bulk component
Suppose that every block of volume ${\displaystyle \ell ^{d}}$ is separate from the rest of the system; inside every one of them the magnetization is uniform (since the linear dimension of the blocks is much smaller than the correlation length), so we can use Landau theory for uniform systems. In the case of the Ising model, it led to:

${\displaystyle {\mathcal {L}}={\frac {\overline {a}}{2}}tm^{2}+{\frac {\overline {b}}{4}}m^{4}}$
The total bulk energy is thus obtained summing over all the blocks:
${\displaystyle \beta {\mathcal {H}}_{\text{eff.}}^{\text{bulk}}=\sum _{\vec {r}}\left({\frac {\overline {a}}{2}}tm^{2}({\vec {r}})+{\frac {\overline {b}}{4}}m^{4}({\vec {r}})\right)}$

interaction component
We now must take into account the fact that adjacent blocks do interact. In particular since as we have stated ${\displaystyle m}$ does not vary much on microscopic scales, the interaction between the blocks must be such that strong variations of magnetization between neighbouring blocks is energetically unfavourable. If we call ${\displaystyle {\vec {\delta }}}$ a vector of magnitude ${\displaystyle \ell }$ that points from one block to a neighbouring one, the most simple analytic expression that we can guess for such a term can be a harmonic one[1]:

${\displaystyle \beta {\mathcal {H}}_{\text{eff.}}^{\text{int.}}=\sum _{\vec {r}}\sum _{\vec {\delta }}{\frac {\overline {k}}{2}}\left(m({\vec {r}})-m({\vec {r}}+{\vec {\delta }})\right)^{2}}$
(the factor ${\displaystyle 1/2}$ multiplying ${\displaystyle {\overline {k}}}$, just like the numeric factors multiplying ${\displaystyle {\overline {a}}}$ and ${\displaystyle {\overline {b}}}$, have been inserted for future convenience). We can also think of this as a first approximation of a general interaction between the blocks, namely as the first terms of a Taylor expansion of the real interaction energy.

Now, since the linear dimension of the blocks ${\displaystyle \ell }$ is much smaller than the characteristic length ${\displaystyle L}$ of the system we can treat ${\displaystyle {\vec {r}}}$ as a continuous variable and thus substitute the sum over ${\displaystyle {\vec {r}}}$ with an integral:

${\displaystyle \sum _{\vec {r}}\longrightarrow {\frac {1}{\ell ^{d}}}\int d^{d}{\vec {r}}}$
(while the sum over ${\displaystyle {\vec {\delta }}}$ remains a sum, since for every ${\displaystyle {\vec {r}}}$ there is only a finite number of nearest neighbours). Therefore:
${\displaystyle \beta {\mathcal {H}}_{\text{eff.}}^{\text{bulk}}={\frac {1}{\ell ^{d}}}\int \left({\frac {\overline {a}}{2}}tm^{2}({\vec {r}})+{\frac {\overline {b}}{4}}m^{4}({\vec {r}})\right)d^{d}{\vec {r}}}$
${\displaystyle \beta {\mathcal {H}}_{\text{eff.}}^{\text{int.}}={\frac {1}{\ell ^{d}}}\int \sum _{\vec {\delta }}{\frac {\overline {k}}{2}}\left(m({\vec {r}})-m({\vec {r}}+{\vec {\delta }})\right)^{2}d^{d}{\vec {r}}}$

Keeping in mind that ${\displaystyle |{\vec {\delta }}|=\ell }$, the interaction term can be rewritten in terms of ${\displaystyle {\vec {\nabla }}m}$:

{\displaystyle {\begin{aligned}{\frac {1}{\ell ^{d}}}\int \sum _{\vec {\delta }}{\frac {\overline {k}}{2}}\left(m({\vec {r}})-m({\vec {r}}+{\vec {\delta }})\right)^{2}d^{d}{\vec {r}}={\frac {1}{\ell ^{d-2}}}\int {\frac {\overline {k}}{2}}\sum _{\vec {\delta }}\left({\frac {m({\vec {r}})-m({\vec {r}}+{\vec {\delta }})}{\ell }}\right)^{2}d^{d}{\vec {r}}=\\=\int {\frac {\overline {k}}{2\ell ^{2-d}}}\sum _{\vec {\delta }}\left({\frac {\partial m}{\partial r_{\delta }}}\right)^{2}d^{d}{\vec {r}}=\int {\frac {\overline {k}}{2\ell ^{2-d}}}\left({\vec {\nabla }}m\right)^{2}d^{d}{\vec {r}}\end{aligned}}}
where we have called ${\displaystyle r_{\delta }}$ the components of ${\displaystyle {\vec {\delta }}}$. Thus, if we now define for the sake of simplicity:
${\displaystyle a:={\frac {\overline {a}}{\ell ^{d}}}\qquad \qquad b:={\frac {\overline {b}}{\ell ^{d}}}\qquad \qquad k:={\frac {\overline {k}}{\ell ^{2-d}}}}$
we will have:
${\displaystyle \beta {\mathcal {H}}_{\text{eff.}}=\int \left[{\frac {a}{2}}tm^{2}({\vec {r}})+{\frac {b}{4}}m^{4}({\vec {r}})+{\frac {k}{2}}\left({\vec {\nabla }}m({\vec {r}})\right)^{2}\right]d^{d}{\vec {r}}}$
Therefore, the (functional) partition function of the system will be:
${\displaystyle Z=\int e^{-\int \left({\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}+{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}\right)d^{d}{\vec {r}}}{\mathcal {D}}[m_{\ell }({\vec {r}})]}$
Let us now make a couple of considerations:

• If ${\displaystyle m({\vec {r}})=m={\text{const.}}}$ the energy of the system has the same structure of the one used in Landau theory
• The term proportional to ${\displaystyle ({\vec {\nabla }}m)^{2}}$ is completely new but we could have introduced it intuitively to a Landau-like mean field functional, since the introduction of spatial variations in the order parameter has an energetic cost which must depend on how it varies in space, i.e. it depends on the gradient of ${\displaystyle m}$. In particular, it must involve ${\displaystyle ({\vec {\nabla }}m)^{2}}$ because of the symmetries of the model: since the system is isotropic and ${\displaystyle \mathbb {Z} _{2}}$-invariant, we must use combinations of derivatives that are invariant under rotations and parity, and ${\displaystyle ({\vec {\nabla }}m)^{2}}$ is the simplest of them[2].

If there is also an external magnetic field ${\displaystyle {\vec {h}}({\vec {r}})=\beta {\vec {H}}({\vec {r}})}$, we must add to the Hamiltonian the term:

${\displaystyle -\int {\vec {h}}({\vec {r}})\cdot {\vec {m}}({\vec {r}})d^{d}{\vec {r}}}$
so that the partition function becomes:
${\displaystyle Z=\int e^{-\int \left({\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}+{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}-hm\right)d^{d}{\vec {r}}}{\mathcal {D}}[m_{\ell }({\vec {r}})]}$
which is a functional of ${\displaystyle m({\vec {r}})}$ and ${\displaystyle h({\vec {r}})}$. As usual, all the thermodynamics of the system can be obtained from ${\displaystyle Z}$, provided that now we take functional derivatives instead of usual derivatives.

## Saddle point approximation: Landau theory

We can now compute ${\displaystyle Z}$, as a first approach, using The saddle point approximation; as we will see this will reproduce a Landau-like mean field theory which will also take into account the presence of inhomogeneities. In particular thanks to the new term involving ${\displaystyle {\vec {\nabla }}m}$ we will be able to compute the fluctuation correlation function[3] and so also to determine the critical exponents ${\displaystyle \eta }$ and ${\displaystyle \nu }$.

Therefore we approximate ${\displaystyle Z}$ with the leading term of the integral, i.e. we must determine the function ${\displaystyle m_{0}}$ that maximizes the exponent, namely minimizes:

${\displaystyle {\mathcal {L}}[m({\vec {r}})]:=\beta {\mathcal {H}}_{\text{eff.}}-\int hmd^{d}{\vec {r}}=\int \left({\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}+{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}-hm\right)d^{d}{\vec {r}}}$
and then compute ${\displaystyle Z}$ as:
${\displaystyle Z\approx e^{-{\mathcal {L}}[m_{0}({\vec {r}})]}}$
where ${\displaystyle m_{0}}$ is determined imposing the stationarity of the functional ${\displaystyle {\mathcal {L}}}$ with respect to ${\displaystyle m}$:
${\displaystyle {\frac {\delta {\mathcal {L}}}{\delta m}}_{|m_{0}}=0}$
This leads to the state equation of the system:
${\displaystyle h({\vec {r}})={\frac {\delta {\tilde {\mathcal {H}}}}{\delta m}}}$
where we have defined ${\displaystyle {\tilde {\mathcal {H}}}:=\beta {\mathcal {H}}_{\text{eff.}}}$ for brevity. If we now call ${\displaystyle {\mathfrak {h}}}$ the integrand of ${\displaystyle {\tilde {\mathcal {H}}}}$, since[4]:
${\displaystyle {\frac {\delta {\tilde {\mathcal {H}}}}{\delta m}}={\frac {\partial {\mathfrak {h}}}{\partial m}}-{\vec {\nabla }}{\frac {\partial {\mathfrak {h}}}{\partial ({\vec {\nabla }}m)}}}$
we have:
${\displaystyle h({\vec {r}})=-k\nabla ^{2}m({\vec {r}})+atm({\vec {r}})+bm^{3}({\vec {r}})}$
Note that if ${\displaystyle h={\text{const.}}}$ and ${\displaystyle m_{0}={\text{const.}}}$ (i.e. the system is uniform) we get the same state equation that we have obtained with Landau theory:
${\displaystyle h=atm_{0}+bm_{0}^{3}}$

### Correlation function in the saddle point approximation

We can now proceed to compute the correlation function within our approximations. In order to do that, we take the (functional) derivative of the state equation with respect to ${\displaystyle h({\vec {r}}{\text{ }}')}$, so that ${\displaystyle \chi _{T}=\delta m/\delta h}$ appears:

${\displaystyle \delta ({\vec {r}}-{\vec {r}}{\text{ }}')=\left[-k\nabla ^{2}+at+3bm^{2}\right]\chi _{T}({\vec {r}}-{\vec {r}}{\text{ }}')}$
Now, from fluctuation-dissipation theorem we know that:
${\displaystyle G({\vec {r}}-{\vec {r}}{\text{ }}')=k_{B}T\chi _{T}({\vec {r}}-{\vec {r}}{\text{ }}')}$
so that:
${\displaystyle \beta \left[-k\nabla ^{2}+at+3bm^{2}\right]G({\vec {r}}-{\vec {r}}{\text{ }}')=\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$
Note that this means that the correlation function ${\displaystyle G({\vec {r}}-{\vec {r}}{\text{ }}')}$ can be interpreted as the Green's function of the operator written between the square brackets.

In case of translationally invariant (i.e. uniform) systems, ${\displaystyle m}$ is constant and equal to the equilibrium values given by the Landau theory for the Ising model; in particular, depending on the sign of ${\displaystyle t}$ there are two possible situations:

${\displaystyle {\boldsymbol {t>0}}}$
In this case ${\displaystyle m({\vec {r}})=m_{0}=0}$, so the last equation becomes:

${\displaystyle \left(-k\nabla ^{2}+at\right)G({\vec {r}}-{\vec {r}}{\text{ }}')=k_{B}T\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$
Defining:
${\displaystyle \xi _{>}:={\sqrt {\frac {k}{at}}}}$
this can be rewritten as:
${\displaystyle \left(-\nabla ^{2}+\xi _{>}^{-2}\right)G({\vec {r}}-{\vec {r}}{\text{ }}')={\frac {k_{B}T}{k}}\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$

${\displaystyle {\boldsymbol {t<0}}}$
In this case the magnetization is:

${\displaystyle m_{0}=\pm {\sqrt {-{\frac {at}{b}}}}}$
so the differential equation for ${\displaystyle G}$ becomes:
${\displaystyle \left(-k\nabla ^{2}-2at\right)G({\vec {r}}-{\vec {r}}{\text{ }}')=k_{B}T\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$
This can be rewritten in a form similar to the previous case; in fact, if we define:
${\displaystyle \xi _{<}={\sqrt {-{\frac {k}{2at}}}}}$
we get:
${\displaystyle \left(-\nabla ^{2}+\xi _{<}^{-2}\right)G({\vec {r}}-{\vec {r}}{\text{ }}')={\frac {k_{B}T}{k}}\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$

We will shortly see that ${\displaystyle \xi _{>}}$ and ${\displaystyle \xi _{<}}$ are the expressions of the correlation length for ${\displaystyle T>T_{c}}$ and ${\displaystyle T, respectively. We can therefore see that in both cases we get:

${\displaystyle \nu ={\frac {1}{2}}}$

Thus, for both the cases ${\displaystyle t>0}$ and ${\displaystyle t<0}$ the correlation function ${\displaystyle G({\vec {r}}-{\vec {r}}{\text{ }}')}$ can be obtained by solving the differential equation:

${\displaystyle \left(-\nabla ^{2}+\xi ^{-2}\right)G({\vec {r}}-{\vec {r}}{\text{ }}')={\frac {k_{B}T}{k}}\delta ({\vec {r}}-{\vec {r}}{\text{ }}')}$
which can be done with Fourier transforms. If we use the following conventions for the Fourier transform:
${\displaystyle {\hat {f}}({\vec {q}})=\int f({\vec {r}}-{\vec {r}}{\text{ }}')e^{-i{\vec {q}}\cdot ({\vec {r}}-{\vec {r}}{\text{ }}')}d^{d}({\vec {r}}-{\vec {r}}{\text{ }}')}$
${\displaystyle f({\vec {r}}-{\vec {r}}{\text{ }}')={\frac {1}{(2\pi )^{d}}}\int {\hat {f}}({\vec {q}})e^{i{\vec {q}}\cdot ({\vec {r}}-{\vec {r}}{\text{ }}')}d^{d}{\vec {q}}}$

then transforming both sides of the equation we get:

${\displaystyle (q^{2}+\xi ^{-2}){\hat {G}}({\vec {q}})={\frac {k_{B}T}{k}}\quad \Rightarrow \quad {\hat {G}}({\vec {q}})={\frac {k_{B}T}{k}}{\frac {1}{q^{2}+\xi ^{-2}}}}$
where ${\displaystyle q=|{\vec {q}}|}$[5]. From this last equation we can also foresee that when ${\displaystyle T=T_{c}}$, since ${\displaystyle \xi =\infty }$ we have ${\displaystyle {\hat {G}}({\vec {q}})\sim q^{-2}}$ and so ${\displaystyle G({\vec {r}})\sim 1/|{\vec {r}}|^{2-d}}$, from which we have that the critical exponent ${\displaystyle \eta }$ is null (we will see that explicitly once we have computed ${\displaystyle G}$). Therefore, renaming ${\displaystyle {\vec {x}}={\vec {r}}-{\vec {r}}{\text{ }}'}$ we can now determine ${\displaystyle G({\vec {x}})}$ with the Fourier antitransform:
${\displaystyle G({\vec {x}})=\int {\frac {d^{d}{\vec {q}}}{(2\pi )^{d}}}{\frac {e^{i{\vec {q}}\cdot {\vec {x}}}}{q^{2}+\xi ^{-2}}}}$
This integral is a bit tedious to compute, and in general its result depends strongly on the dimensionality ${\displaystyle d}$ of the system; the general approach used to solve it is to shift to spherical coordinates in ${\displaystyle \mathbb {R} ^{d}}$ and then complex integration for the remaining part, which involves ${\displaystyle |{\vec {q}}|}$. In order to do some explicit computations, let us consider the case ${\displaystyle d=3}$; we will then have:
${\displaystyle G({\vec {x}})=\int {\frac {d^{3}{\vec {q}}}{(2\pi )^{3}}}{\frac {e^{i{\vec {q}}\cdot {\vec {x}}}}{q^{2}+\xi ^{-2}}}={\frac {1}{(2\pi )^{3}}}\int _{0}^{\infty }{\frac {q^{2}}{q^{2}+\xi ^{-2}}}dq\int _{-1}^{1}e^{iqx\cos \theta }d(\cos \theta )\int _{0}^{2\pi }d\varphi }$
Therefore:
${\displaystyle G({\vec {x}})={\frac {1}{(2\pi )^{2}}}\int _{0}^{\infty }{\frac {q^{2}}{q^{2}+\xi ^{-2}}}dq\left({\frac {e^{iqx}}{iqx}}\right)_{-1}^{1}={\frac {1}{(2\pi )^{2}x}}\int _{0}^{\infty }{\frac {q\sin(qx)}{q^{2}+\xi ^{-2}}}dq}$
This last integral can be computed, using the residue theorem, extending it to the complex plane:
{\displaystyle {\begin{aligned}\int _{0}^{\infty }{\frac {q\sin(qx)}{q^{2}+\xi ^{-2}}}dq={\frac {1}{2}}\int _{-\infty }^{+\infty }{\frac {q\sin(qx)}{q^{2}+\xi ^{-2}}}dq={\frac {1}{2}}\operatorname {Im} \int {\frac {ze^{izx}}{z^{2}+\xi ^{-2}}}dz=\\={\frac {1}{2}}\operatorname {Im} \int {\frac {ze^{izx}}{(z+i\xi ^{-1})(z-i\xi ^{-1})}}dz\end{aligned}}}
Now, the integrand exhibits two poles at ${\displaystyle \pm i\xi ^{-1}}$; we choose as the contour of integration ${\displaystyle \Gamma }$ the upper semicircle in the complex plane, which contains only the pole at ${\displaystyle i\xi ^{-1}}$ and so using the residue theorem we will have:
${\displaystyle {\frac {1}{2}}\operatorname {Im} \int _{\Gamma }{\frac {ze^{izx}}{(z+i\xi ^{-1})(z-i\xi ^{-1})}}dz={\frac {1}{2}}\operatorname {Im} \left[2\pi i{\frac {i\xi ^{-1}e^{-x/\xi }}{2i\xi ^{-1}}}\right]={\frac {1}{2}}\operatorname {Im} \left[2\pi i{\frac {e^{-x/\xi }}{2}}\right]={\frac {\pi }{2}}e^{-x/\xi }}$
Therefore, in the end we have:
${\displaystyle G({\vec {x}})={\frac {1}{8\pi }}{\frac {e^{-|{\vec {x}}|/\xi }}{|{\vec {x}}|}}}$

We see now clearly that the correlation function has indeed an exponential behaviour (as we have stated also in Long range correlations) and that ${\displaystyle \xi }$ is really the correlation length; furthermore, ${\displaystyle G({\vec {x}})\sim 1/|{\vec {x}}|}$ and from the definition of the exponent ${\displaystyle \eta }$ we have ${\displaystyle G({\vec {x}})\sim 1/|{\vec {x}}|^{d-2+\eta }}$, so since ${\displaystyle d=3}$ we indeed have ${\displaystyle \eta =0}$.

Therefore, we have seen that for the Ising model ${\displaystyle \nu =1/2}$. If we also consider the values of the other critical exponents we see that the upper critical dimension for this model is ${\displaystyle d_{c}=4}$. In other words, mean field theories are actually good approximations for the Ising model if ${\displaystyle d\geq 4}$. We will later see some other confirmations of this fact.

## Gaussian approximation

Until now even if we have introduced Ginzburg-Landau theory we are still neglecting the effects of the fluctuations since we are regarding the mean field theory approximation for non-homogeneous systems as a saddle point approximation of a more general theory; in other words, since we are approximating ${\displaystyle Z}$ as ${\displaystyle e^{-{\mathcal {L}}[m_{0}({\vec {r}})]}}$ we are still regarding the magnetization ${\displaystyle m}$ as non fluctuating over the system. In order to include the fluctuations we must do more and go further the simple saddle point approximation. The simplest way we can include fluctuations in our description is expanding ${\displaystyle Z}$expressed as a functional integral around the stationary solution and keeping only quadratic terms; this means that we are considering fluctuations that follow a normal distribution around the stationary value. The important thing to note, however, is that in this approximation these fluctuations are independent, i.e. they do not interact with each other[6]. As we will see, with this assumption the values of some critical exponents will differ from the "usual" ones predicted by mean field theories.

Let us apply this approximation from easier cases to more complex ones (and finally to the one we are interested in).

### Gaussian approximation for one degree of freedom

Let us consider a system with a single degree of freedom ${\displaystyle q}$, and call ${\displaystyle {\mathcal {H}}(q)}$ its Hamiltonian. Supposing that ${\displaystyle q_{0}}$ is a minimum for ${\displaystyle {\mathcal {H}}}$, i.e. ${\displaystyle \partial {\mathcal {H}}/\partial q_{|q_{0}}=0}$, expanding ${\displaystyle {\mathcal {H}}}$ around ${\displaystyle q_{0}}$ we get:

${\displaystyle {\mathcal {H}}(q)={\mathcal {H}}(q_{0})+{\frac {1}{2}}{\frac {\partial ^{2}{\mathcal {H}}}{\partial q^{2}}}_{|q_{0}}\delta q^{2}+O(\delta q^{3})}$
where ${\displaystyle \delta q=q-q_{0}}$ is the fluctuation of ${\displaystyle q}$ around its stationary value ${\displaystyle q_{0}}$. The Boltzmann factor needed to compute the partition function will therefore be:
${\displaystyle e^{-\beta {\mathcal {H}}(q)}\sim e^{-\beta {\mathcal {H}}(q_{0})-{\frac {1}{2}}{\frac {\delta q^{2}}{\lambda ^{2}}}}}$
where for the sake of simplicity we have defined:
${\displaystyle {\frac {1}{\lambda ^{2}}}:=\beta {\frac {\partial ^{2}{\mathcal {H}}}{\partial q^{2}}}_{|q_{0}}}$
With this approximation, the partition function of the system results:
${\displaystyle Z=\int _{-\infty }^{+\infty }e^{-\beta {\mathcal {H}}(q)}\sim e^{-\beta {\mathcal {H}}(q_{0})}\int _{-\infty }^{+\infty }e^{-{\frac {\delta q^{2}}{2\lambda ^{2}}}}dq}$
The last one is a gaussian integral, which readily gives:
${\displaystyle \int _{-\infty }^{+\infty }e^{-{\frac {(q-q_{0})^{2}}{2\lambda ^{2}}}}dq={\sqrt {2\pi \lambda ^{2}}}}$
Therefore:
${\displaystyle Z=e^{-\beta {\mathcal {H}}(q_{0})}{\sqrt {2\pi \lambda ^{2}}}}$
and since ${\displaystyle Z=e^{-\beta F}}$ the free energy of the system is:
${\displaystyle F={\mathcal {H}}(q_{0})-{\frac {k_{B}T}{2}}\ln \left(2\pi \lambda ^{2}\right)}$
We thus see that the introduction of the fluctuations in the degree of freedom ${\displaystyle q}$ has led to the appearance of an entropic term in the free energy.

### Gaussian approximation for N degrees of freedom

This is a simple generalization of the previous case; the Hamiltonian will now be a function of the ${\displaystyle N}$-component vector ${\displaystyle {\vec {q}}=(q_{1},\dots ,q_{N})}$, and calling ${\displaystyle {\vec {q}}_{0}}$ the minimum of ${\displaystyle {\mathcal {H}}}$, expanding around ${\displaystyle {\vec {q}}_{0}}$ we get:

${\displaystyle {\mathcal {H}}({\vec {q}})={\mathcal {H}}({\vec {q}}_{0})+\sum _{\alpha ,\beta }{\frac {1}{2}}\left(q_{\alpha }-{q_{0}}_{\alpha }\right){\frac {\partial ^{2}{\mathcal {H}}}{\partial q_{\alpha }\partial q_{\beta }}}_{|{\vec {q}}_{0}}\left(q_{\beta }-{q_{0}}_{\beta }\right)+\cdots }$
Now, the Hessian matrix:
${\displaystyle {\mathcal {H}}_{\alpha \beta }={\frac {\partial ^{2}{\mathcal {H}}}{\partial q_{\alpha }\partial q_{\beta }}}_{|{\vec {q}}_{0}}}$
is of course symmetric (thanks to Schwarz's theorem), so it can be diagonalized. Calling ${\displaystyle {\hat {\lambda }}_{i}}$ its eigenvalues and defining ${\displaystyle \beta {\hat {\lambda }}_{i}:=1/\lambda _{i}^{2}}$, we can write:
${\displaystyle \beta {\mathcal {H}}({\vec {q}})=\beta {\mathcal {H}}({\vec {q}}_{0})+{\frac {1}{2}}\sum _{i=1}^{N}{\frac {{\tilde {q}}_{i}^{2}}{\lambda _{i}^{2}}}}$
where ${\displaystyle {\tilde {q}}_{i}}$ is the fluctuation from the eigenvector of the Hessian relative to the eigenvalue ${\displaystyle {\hat {\lambda }}_{i}}$. We therefore have:
${\displaystyle Z=\int _{-\infty }^{+\infty }e^{-\beta {\mathcal {H}}({\vec {q}})}d{\vec {q}}=e^{-\beta {\mathcal {H}}({\vec {q}}_{0})}\prod _{i=1}^{N}\int _{-\infty }^{+\infty }e^{-{\frac {{\tilde {q}}_{i}^{2}}{2\lambda _{i}^{2}}}}dq_{i}}$
The last integrals are all independent and of the same form as the previous case. Thus, in the end we have:
${\displaystyle F={\mathcal {H}}({\vec {q}}_{0})-{\frac {k_{B}T}{2}}\sum _{i=1}^{N}\ln \left(2\pi \lambda _{i}^{2}\right)}$

### Gaussian approximation for infinite degrees of freedom

Let us now move to the really interesting case, i.e. the case of infinite degrees of freedom. In general terms (we shall shortly see this explicitly for the Ising model) we want to compute a partition function of the form:

${\displaystyle Z=\int e^{-\beta {\mathcal {H}}(m({\vec {r}}))}{\mathcal {D}}[m({\vec {r}})]}$
and the Gaussian approximation will be obtained, in analogy with the previous cases, determining the extremal (uniform) solution ${\displaystyle m({\vec {r}})=m_{0}}$ and then expanding ${\displaystyle {\mathcal {H}}}$ around ${\displaystyle m_{0}}$ to the second order in the fluctuation ${\displaystyle \delta m({\vec {r}})=m({\vec {r}})-m_{0}}$. Thus, we will in general obtain a Hamiltonian of the form:
${\displaystyle {\mathcal {H}}(m({\vec {r}}))={\mathcal {H}}(m_{0})+{\frac {1}{2}}\int \delta m({\vec {r}})G_{0}^{-1}({\vec {r}},{\vec {r}}{\text{ }}')\delta m({\vec {r}}{\text{ }}')d^{d}{\vec {r}}d^{d}{\vec {r}}{\text{ }}'}$
where:
${\displaystyle G_{0}^{-1}({\vec {r}},{\vec {r}}{\text{ }}'):={\frac {\delta ^{2}{\mathcal {H}}}{\delta m({\vec {r}})\delta m({\vec {r}}{\text{ }}')}}_{\left|{\begin{smallmatrix}m({\vec {r}})=m_{0}\\h=0\end{smallmatrix}}\right.}}$
Therefore:
${\displaystyle F={\mathcal {H}}(m_{0})+{\frac {k_{B}T}{2}}\int \ln \left({\frac {G_{0}^{-1}({\vec {r}},{\vec {r}}{\text{ }}')}{2\pi }}\right)d^{d}{\vec {r}}}$

### Gaussian approximation for the Ising model in Ginzburg-Landau theory

Let us now apply what we have just only stated to a concrete case, i.e. the Ginzburg-Landau theory for the Ising model we were considering. In this case:

${\displaystyle \beta {\mathcal {H}}(m({\vec {r}}))=\int \left[{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}+{\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}-hm\right]d^{d}{\vec {r}}}$
and we have seen that the stationary solution ${\displaystyle m_{0}}$ is such that:
${\displaystyle h=atm_{0}+bm_{0}^{3}}$
Let us now include also the fluctuations of ${\displaystyle m}$, substituting ${\displaystyle m=m_{0}+\delta m}$. To the second order in ${\displaystyle \delta m}$ we have:

and so setting ${\displaystyle h=0}$:

{\displaystyle {\begin{aligned}\beta {\mathcal {H}}(m({\vec {r}}))=\int \left[{\frac {a}{2}}tm_{0}^{2}+{\frac {b}{4}}m_{0}^{4}+(atm_{0}+bm_{0}^{3})\delta m+\left({\frac {a}{2}}t+{\frac {3}{2}}bm_{0}^{3}\right)\delta m^{2}+\right.\\\left.+{\frac {k}{2}}\left({\vec {\nabla }}\delta m\right)^{2}\right]d^{d}{\vec {r}}\end{aligned}}}
Since ${\displaystyle atm_{0}+bm_{0}^{3}=h=0}$, defining for simplicity ${\displaystyle A_{0}:=atm_{0}^{2}/2+bm_{0}^{4}/4}$ and calling ${\displaystyle V}$ the volume of the system we get:
${\displaystyle \beta {\mathcal {H}}(m({\vec {r}}))=A_{0}V+\int \left[\left({\frac {a}{2}}t+{\frac {3}{2}}bm_{0}^{3}\right)\delta m^{2}+{\frac {k}{2}}\left({\vec {\nabla }}\delta m\right)^{2}\right]d^{d}{\vec {r}}}$
In order to compute this integral it is more convenient to shift to Fourier space.

Let us make some remarks on what happens when we apply Fourier transformations in this case. If our system is enclosed in a cubic box of volume ${\displaystyle V=L^{d}}$, we can define the Fourier components of the magnetization as:

${\displaystyle m_{\vec {k}}=\int _{V}e^{-i{\vec {k}}\cdot {\vec {r}}}m({\vec {r}})d^{d}{\vec {r}}}$
where ${\displaystyle {\vec {k}}=2\pi {\vec {n}}/L}$ and ${\displaystyle {\vec {n}}}$ is a vector whose components are integer numbers. We can therefore expand the magnetization in a Fourier series:
${\displaystyle m({\vec {r}})={\frac {1}{V}}\sum _{\vec {k}}e^{i{\vec {k}}\cdot {\vec {r}}}m_{\vec {k}}}$
Substituting this expression of ${\displaystyle m}$ in ${\displaystyle m_{\vec {k}}}$ we obtain an integral representation for the Kronecker delta; in fact:
${\displaystyle m_{\vec {k}}=\sum _{{\vec {k}}{\text{ }}'}m_{{\vec {k}}{\text{ }}'}\left({\frac {1}{V}}\int _{V}e^{i({\vec {k}}-{\vec {k}}{\text{ }}')\cdot {\vec {r}}}d^{d}{\vec {r}}\right)}$
and this is true only if:
${\displaystyle {\frac {1}{V}}\int _{V}e^{i({\vec {k}}-{\vec {k}}{\text{ }}')\cdot {\vec {r}}}d^{d}{\vec {r}}=\delta _{{\vec {k}},{\vec {k}}{\text{ }}'}\quad \Rightarrow \quad \int _{V}e^{i({\vec {k}}-{\vec {k}}{\text{ }}')\cdot {\vec {r}}}d^{d}{\vec {r}}=V\delta _{{\vec {k}},{\vec {k}}{\text{ }}'}}$
Let us now make two observations. First: since ${\displaystyle m({\vec {r}})}$ is real we have that ${\displaystyle m_{\vec {k}}^{*}=m_{-{\vec {k}}}}$. Second: our coarse graining procedure is based on the construction of blocks which have a linear dimension that cannot be smaller than ${\displaystyle a}$, the characteristic microscopic length of the system; this means that not all the ${\displaystyle {\vec {k}}}$ are allowed, and in particular we must have ${\displaystyle |{\vec {k}}|\leq \pi /a:=\Lambda }$. Now, thinking about the functional integral form of the partition function, what does the trace ${\displaystyle \int {\mathcal {D}}[m({\vec {r}})]}$ become in Fourier space? Since ${\displaystyle m({\vec {r}})}$ is expressed in terms of the Fourier modes ${\displaystyle m_{\vec {k}}}$, which are in general complex, the measure of the integral becomes:
${\displaystyle \int _{-\infty }^{+\infty }\prod _{|{\vec {k}}|<\Lambda }d\left(\operatorname {Re} m_{\vec {k}}\right)d\left(\operatorname {Im} m_{\vec {k}}\right)}$
However, since ${\displaystyle m({\vec {r}})}$ is real (i.e. ${\displaystyle m_{\vec {k}}^{*}=m_{-{\vec {k}}}}$) the real and imaginary parts of the Fourier modes are not independent, because we have:
${\displaystyle \operatorname {Re} m_{\vec {k}}=\operatorname {Re} m_{-{\vec {k}}}\quad \qquad \operatorname {Im} m_{\vec {k}}=-\operatorname {Im} m_{-{\vec {k}}}}$
This means that if we use the trace we have written above we would integrate twice on the complex plane; we must therefore change the measure so as to avoid this double counting. We can for example simply divide everything by ${\displaystyle 2}$, or restrict the integration on the region where for example the last coordinate of ${\displaystyle {\vec {k}}}$, let us call it ${\displaystyle k_{z}}$, is positive. Therefore:
{\displaystyle {\begin{aligned}\operatorname {Tr} =\int {\mathcal {D}}[m({\vec {r}})]={\frac {1}{2}}\int _{-\infty }^{+\infty }\prod _{|{\vec {k}}|<\Lambda }d\left(\operatorname {Re} m_{\vec {k}}\right)d\left(\operatorname {Im} m_{\vec {k}}\right)=\\=\int _{-\infty }^{+\infty }\prod _{|{\vec {k}}|<\Lambda ,k_{z}>0}d\left(\operatorname {Re} m_{\vec {k}}\right)d\left(\operatorname {Im} m_{\vec {k}}\right):=\int _{-\infty }^{+\infty }\prod _{\vec {k}}'dm_{\vec {k}}\end{aligned}}}
where the last step defines the symbolic notation ${\textstyle \prod _{\vec {k}}'}$. In the end, we have:
${\displaystyle Z=\int _{-\infty }^{+\infty }\prod _{\vec {k}}'e^{-\beta {\mathcal {H}}\left(m_{\vec {k}}\right)}dm_{\vec {k}}}$

Let us now compute the partition function of the system in the simpler case ${\displaystyle T>T_{c}}$[7], so that in the end we can determine the free energy of the system. In this case ${\displaystyle m_{0}=0}$ so ${\displaystyle \delta m=m}$, and therefore substituting:

${\displaystyle m_{\vec {k}}=\int _{V}e^{-i{\vec {k}}\cdot {\vec {r}}}m({\vec {r}})d^{d}{\vec {r}}}$
into:

${\displaystyle \beta {\mathcal {H}}(m({\vec {r}}))=A_{0}V+\int \left[\left({\frac {a}{2}}t+{\frac {3}{2}}bm_{0}^{3}\right)\delta m^{2}+{\frac {k}{2}}\left({\vec {\nabla }}\delta m\right)^{2}\right]d^{d}{\vec {r}}}$
we get (renaming ${\displaystyle {\vec {q}}}$ the coordinate in Fourier space, so as not to confuse it with the constant ${\displaystyle k}$):
{\displaystyle {\begin{aligned}\beta {\mathcal {H}}-A_{0}V=\\=\int d^{d}{\vec {r}}\left({\frac {k}{2V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}{\vec {q}}\cdot {\vec {q}}{\text{ }}'e^{i({\vec {q}}-{\vec {q}}{\text{ }}')\cdot {\vec {r}}}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}^{*}+{\frac {at}{2V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}e^{i({\vec {q}}-{\vec {q}}{\text{ }}')\cdot {\vec {r}}}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}^{*}\right)=\\{}\\={\frac {1}{2V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}^{*}\left(k{\vec {q}}\cdot {\vec {q}}{\text{ }}'+at\right)\underbrace {\int e^{i({\vec {q}}-{\vec {q}}{\text{ }}')\cdot {\vec {r}}}d^{d}{\vec {r}}} _{V\delta _{{\vec {q}},{\vec {q}}{\text{ }}'}}=\\={\frac {1}{2V}}\sum _{\vec {q}}m_{\vec {q}}m_{-{\vec {q}}}\left(k|{\vec {q}}|^{2}+at\right)=\\={\frac {1}{2V}}\sum _{\vec {q}}\left(k|{\vec {q}}|^{2}+at\right)|m_{\vec {q}}|^{2}\end{aligned}}}
Therefore substituting in the expression of the partition function the exponentials factorize, and in the end:
${\displaystyle Z=\prod _{\vec {q}}'\int _{-\infty }^{+\infty }e^{-{\frac {\beta }{2V}}(kq^{2}+at)|m_{\vec {q}}|^{2}+\beta A_{0}V}dm_{\vec {q}}}$
Since ${\displaystyle |m_{\vec {q}}|^{2}=\operatorname {Re} ^{2}m_{\vec {q}}+\operatorname {Im} ^{2}m_{\vec {q}}}$, changing variables to:
${\displaystyle x=\operatorname {Re} m_{\vec {q}}\quad \qquad y=\operatorname {Im} m_{\vec {q}}}$
the integration in ${\displaystyle dm_{\vec {q}}=d\left(\operatorname {Re} m_{\vec {q}}\right)d\left(\operatorname {Im} m_{\vec {q}}\right)}$ gives:
${\displaystyle \int _{-\infty }^{+\infty }dxdye^{-{\frac {kq^{2}+at}{2V}}(x^{2}+y^{2})}=\pi {\frac {2V}{kq^{2}+at}}}$
Thus:
${\displaystyle Z=\prod _{\vec {q}}'{\frac {2\pi V}{kq^{2}+at}}e^{-\beta A_{0}V}}$
We therefore have that the free energy of the system is:
${\displaystyle F=-k_{B}T\ln Z=A_{0}V-{\frac {k_{B}T}{2}}\sum _{|{\vec {q}}|<\Lambda }\ln \left({\frac {2\pi V}{kq^{2}+at}}\right)}$

We can now compute the specific heat of the system, and so determine its critical exponent ${\displaystyle \alpha }$. We therefore want to compute:

${\displaystyle C_{V}=-T{\frac {\partial ^{2}}{\partial T^{2}}}{\frac {F}{V}}}$
The derivatives are straightforward, and in the end we get:
${\displaystyle {\frac {C_{V}}{T}}={\frac {k_{B}Ta}{VT_{c}^{2}}}\sum _{|{\vec {q}}|<\Lambda }{\frac {1}{(kq^{2}+at)^{2}}}-{\frac {k_{B}a^{2}}{2VT_{c}}}\sum _{|{\vec {q}}|<\Lambda }{\frac {1}{kq^{2}+at}}}$
Let us now consider the two terms separately and study their behaviour for ${\displaystyle t\to 0^{+}}$ (we are in fact considering ${\displaystyle T>T_{c}}$). Neglecting the proportionality constants that we don't need, we can rewrite the first contribution as:
${\displaystyle I_{1}={\frac {1}{V}}\sum _{|{\vec {q}}|<\Lambda }{\frac {1}{(kq^{2}+at)^{2}}}=\int _{|{\vec {q}}|<\Lambda }{\frac {d^{d}{\vec {q}}}{(2\pi )^{d}}}{\frac {1}{(kq^{2}+at)^{2}}}}$
where we have substituted the sum with an integral[8], since the density of states in Fourier space is high (it is proportional to ${\displaystyle V}$, see the footnote). Now, using the definition of ${\displaystyle \xi :=\xi _{>}={\sqrt {k/at}}}$ that we have previously seen we have:
${\displaystyle I_{1}={\frac {1}{k^{2}}}\int _{|{\vec {q}}|<\Lambda }{\frac {d^{d}{\vec {q}}}{(2\pi )^{d}}}{\frac {1}{(q^{2}+\xi ^{-2})^{2}}}}$
In order to understand the divergent behaviour of the integral for ${\displaystyle t\to 0^{+}}$, we use a "scaling trick"; we change variable defining:
${\displaystyle {\vec {x}}=\xi {\vec {q}}}$
so that the integral becomes:
${\displaystyle I_{1}={\frac {1}{k^{2}}}\int _{|{\vec {x}}|<\xi \Lambda }{\frac {d^{d}{\vec {x}}}{(2\pi )^{d}}}\xi ^{-d}\left({\frac {x^{2}}{\xi ^{2}}}+{\frac {1}{\xi ^{2}}}\right)^{-2}={\frac {\xi ^{4-d}}{k^{2}}}\int _{|{\vec {x}}|<\xi \Lambda }{\frac {d^{d}{\vec {x}}}{(2\pi )^{d}}}{\frac {1}{(1+x^{2})^{2}}}}$
and for ${\displaystyle t\to 0^{+}}$ we know that ${\displaystyle \xi \to \infty }$, so the integral is computed for all values of ${\displaystyle x=|{\vec {x}}|}$. Now, this integral must be computed shifting to spherical coordinates, so we will have ${\displaystyle d^{d}{\vec {x}}\propto x^{d-1}dx}$ (a part from numerical factors that involve all the angles, which are integrated trivially since the integrand only depends on ${\displaystyle x}$). Therefore the integrand of ${\displaystyle I_{1}}$ is ${\displaystyle x^{d-1}/(1+x^{2})^{2}}$, and its behaviour for large or small ${\displaystyle x}$ is:
${\displaystyle {\frac {x^{d-1}}{(1+x^{2})^{2}}}\sim {\begin{cases}0&x\to 0\\x^{d-5}&x\to \infty \end{cases}}}$
where of course in the first case we have ${\displaystyle d>1}$; for large ${\displaystyle x}$, since in general:
${\displaystyle \int _{\overline {x}}^{\infty }{\frac {dx}{x^{a}}}<\infty \quad \quad {\text{if }}a>1}$
then the integral will converge if ${\displaystyle 5-d>1}$, i.e. ${\displaystyle d<4}$. To sum up, the integral that appears in ${\displaystyle I_{1}}$ converges for ${\displaystyle 1 (since in this case the integrand does not diverge in the domain of integration); of course, when this integral converges its result is simply a number. We therefore have that:
${\displaystyle I_{1}\propto \xi ^{4-d}\quad \quad d<4}$
and since ${\displaystyle \xi \to \infty }$ for ${\displaystyle t\to 0^{+}}$, we see that ${\displaystyle I_{1}}$ brings to ${\displaystyle C_{V}}$ a diverging contribution[9] for ${\displaystyle T\to T_{c}^{+}}$. We can also wonder what happens for ${\displaystyle d>4}$. From what we have stated about the rescaled form of ${\displaystyle I_{1}}$ we could think that the integral diverges, but we must also take into account the prefactor ${\displaystyle \xi ^{4-d}}$, which tends to zero for ${\displaystyle d>4}$ as the transition is approached (since ${\displaystyle \xi \to \infty }$). The net result is finite, as could also be argued from the original (unscaled) form of ${\displaystyle I_{1}}$; in fact if ${\displaystyle \xi \to \infty }$ in ${\displaystyle I_{1}}$ then (in spherical coordinates) the integrand is proportional to ${\displaystyle q^{d-1-4}=q^{d-5}}$, and since:
${\displaystyle \int _{0}^{\overline {x}}{\frac {dx}{x^{a}}}<\infty \quad \quad {\text{for }}a<1}$
then ${\displaystyle I_{1}}$ converges if ${\displaystyle 5-d<1}$, i.e. ${\displaystyle d>4}$. To sum up, we can say that the first contribution to the specific heat behaves as:
${\displaystyle {\frac {k_{B}Ta}{VT_{c}^{2}}}\sum _{|{\vec {q}}|<\Lambda }{\frac {1}{(kq^{2}+at)^{2}}}\sim {\begin{cases}t^{-\nu (4-d)}&d<4\\{\text{finite}}&d>4\end{cases}}}$
where we have also used the definition of the exponent ${\displaystyle \nu }$, i.e. ${\displaystyle \xi \sim t^{-\nu }}$. Note that this also means that in the Gaussian approximation this term brings no corrections to the exponent ${\displaystyle \alpha }$ above four dimensions.

Let us now consider the second contribution to ${\displaystyle C_{V}/T}$. In particular, as we have done before we rewrite it substituting the sum with an integral and also using the definition of ${\displaystyle \xi }$, so that:

${\displaystyle I_{2}={\frac {1}{k}}\int _{|{\vec {q}}|<\Lambda }{\frac {d^{d}{\vec {q}}}{(2\pi )^{d}}}{\frac {1}{q^{2}+\xi ^{-2}}}}$
Again, we are interested in its behaviour for ${\displaystyle \xi \to \infty }$, and as we have already noted it is only the behaviour of the integrand for ${\displaystyle q\to 0}$ that can cause any divergence. As previously done, changing variable to ${\displaystyle {\vec {x}}=\xi {\vec {q}}}$:
${\displaystyle I_{2}={\frac {\xi ^{2-d}}{k}}\int _{|{\vec {x}}|<\xi \Lambda }{\frac {d^{d}{\vec {x}}}{(2\pi )^{d}}}{\frac {1}{x^{2}+1}}}$
and using spherical coordinates:
${\displaystyle I_{2}\propto \xi ^{2-d}\int _{0}^{\xi \Lambda }{\frac {x^{d-1}}{x^{2}+1}}dx}$
The integrand behaves as:
${\displaystyle {\frac {x^{d-1}}{x^{2}+1}}\sim {\begin{cases}0&x\to 0\\x^{d-3}&x\to \infty \end{cases}}}$
where again in the first case ${\displaystyle d>1}$. Therefore, the integral in ${\displaystyle I_{2}}$ (not${\displaystyle I_{2}}$ itself) in the limit ${\displaystyle \xi \to \infty }$ converges if ${\displaystyle 3-d>1}$, i.e. ${\displaystyle d<2}$; this means that for ${\displaystyle d<2}$ and ${\displaystyle t\to 0^{+}}$ we have that ${\displaystyle I_{2}}$ behaves as:
${\displaystyle I_{2}\sim \xi ^{2-d}\sim t^{-\nu (2-d)}}$
For the case ${\displaystyle d>2}$ it is more convenient to consider the unscaled form of ${\displaystyle I_{2}}$. In this case, in the limit ${\displaystyle \xi \to \infty }$ we have (again in spherical coordinates):
${\displaystyle I_{2}\propto \int _{0}^{\Lambda }{\frac {q^{d-1}}{q^{2}}}dq}$
and the integrand converges for ${\displaystyle q\to 0}$ if ${\displaystyle d>2}$. Therefore for ${\displaystyle d<2}$ the second contribution to ${\displaystyle C_{V}/T}$ diverges, but in the same range of ${\displaystyle d}$ the divergence of the first contribution is more relevant; on the other hand, for ${\displaystyle 2\leq d<4}$ only the first contribution diverges. It is therefore the term containing ${\displaystyle I_{1}}$ that determines the divergence of the specific heat, and in particular for ${\displaystyle d<4}$ we have ${\displaystyle C_{V}\sim t^{-\nu (4-d)}}$ and so we see that in the Gaussian approximation the inclusion of the fluctuations has changed the value of the critical exponent ${\displaystyle \alpha }$ to[10]:
${\displaystyle \alpha =\nu (4-d)}$
In order to compute it, however, we still must determine ${\displaystyle \nu }$ so we now proceed to compute the two-point correlation function in order to determine both ${\displaystyle \eta }$ and ${\displaystyle \nu }$.

### Two-point correlation function in the Gaussian approximation

We know that the (simple) correlation function is defined as:

${\displaystyle G({\vec {r}},{\vec {r}}{\text{ }}')=\left\langle m({\vec {r}})m({\vec {r}}{\text{ }}')\right\rangle }$
so we first have to determine:
${\displaystyle m({\vec {r}})m({\vec {r}}{\text{ }}')={\frac {1}{V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}e^{i{\vec {q}}\cdot ({\vec {r}}+{\vec {r}}{\text{ }}')}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}}$
Shifting to Fourier space, we have:
${\displaystyle \left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle ={\frac {\int dm_{{\vec {q}}_{1}}dm_{{\vec {q}}_{2}}\cdots dm_{\vec {q}}dm_{{\vec {q}}{\text{ }}'}e^{-\beta {\mathcal {H}}}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}}{\int dm_{{\vec {q}}_{1}}dm_{{\vec {q}}_{2}}\cdots dm_{\vec {q}}dm_{{\vec {q}}{\text{ }}'}e^{-\beta {\mathcal {H}}}}}}$
where:
${\displaystyle \beta {\mathcal {H}}=A_{0}V+{\frac {1}{2V}}\sum _{{\vec {k}},{\vec {k}}{\text{ }}'}\left(k|{\vec {q}}|^{2}+at\right)|m_{\vec {q}}|^{2}}$
It is clear that in ${\displaystyle \left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle }$ all the integrals factorize since the Fourier modes are all independent (they are decoupled); therefore, all the integrals in the numerator that don't involve ${\displaystyle {\vec {q}}}$ or ${\displaystyle {\vec {q}}{\text{ }}'}$ simplify with the same integrals in the denominator, so that in the end we are left with:
${\displaystyle \left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle ={\frac {\int e^{-{\frac {k|{\vec {q}}|^{2}+at}{2V}}|m_{\vec {q}}|^{2}}e^{-{\frac {k|{\vec {q}}{\text{ }}'|^{2}+at}{V}}|m_{{\vec {q}}{\text{ }}'}|^{2}}m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}dm_{\vec {q}}dm_{{\vec {q}}{\text{ }}'}}{\int e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}|m_{\vec {q}}|^{2}}e^{-{\frac {k|{\vec {q}}{\text{ }}'|^{2}+at}{V}}|m_{{\vec {q}}{\text{ }}'}|^{2}}dm_{\vec {q}}dm_{{\vec {q}}{\text{ }}'}}}}$
(where the factor ${\displaystyle 2}$ in the denominator of the exponent has disappeared because we must remember that ${\displaystyle m_{\vec {q}}=m_{-{\vec {q}}}^{*}}$). There are now two possible cases:

${\displaystyle {\boldsymbol {{\vec {q}}\neq \pm {\vec {q}}{\text{ }}'}}}$
In this case (which can be re-expressed as ${\displaystyle |{\vec {q}}|\neq |{\vec {q}}{\text{ }}'|}$) the two coefficients ${\displaystyle m_{\vec {q}}}$ and ${\displaystyle m_{{\vec {q}}{\text{ }}'}}$ are distinct, and in the numerator the double integral factorizes into two integrals of the form:

${\displaystyle \int e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}|m_{\vec {q}}|^{2}}m_{\vec {q}}dm_{\vec {q}}=0}$
since the integrand is odd. We therefore have:
${\displaystyle G_{{\vec {q}},{\vec {q}}{\text{ }}'}=\left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle =0\qquad \qquad |{\vec {q}}|\neq |{\vec {q}}{\text{ }}'|}$

${\displaystyle {\boldsymbol {{\vec {q}}=\pm {\vec {q}}{\text{ }}'}}}$
In this case (equivalent to ${\displaystyle |{\vec {q}}|=|{\vec {q}}{\text{ }}'|}$) we can either have ${\displaystyle {\vec {q}}={\vec {q}}{\text{ }}'}$ so that ${\displaystyle \left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle =\langle m_{\vec {q}}^{2}\rangle }$, or ${\displaystyle {\vec {q}}=-{\vec {q}}{\text{ }}'}$ so that ${\displaystyle \left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle =\langle m_{\vec {q}}m_{\vec {q}}^{*}\rangle =\langle |m_{\vec {q}}^{2}|\rangle }$.

Let us first consider the case ${\displaystyle {\vec {q}}={\vec {q}}{\text{ }}'}$. Using polar coordinates, we define ${\displaystyle m_{\vec {q}}=|m_{\vec {q}}|e^{i\theta _{\vec {q}}}}$ so that the measure ${\displaystyle dm_{\vec {q}}}$ in the complex plane becomes:

${\displaystyle dm_{\vec {q}}=|m_{\vec {q}}|d|m_{\vec {q}}|d\theta _{\vec {q}}}$
Thus:
${\displaystyle \left\langle m_{\vec {q}}^{2}\right\rangle ={\frac {\int e^{-\beta {\mathcal {H}}}m_{\vec {q}}^{2}dm_{\vec {q}}}{\int e^{-\beta {\mathcal {H}}}dm_{\vec {q}}}}\propto \int _{0}^{2\pi }e^{2i\theta _{\vec {q}}}d\theta _{\vec {q}}=0}$
We are therefore left with the last case ${\displaystyle {\vec {q}}=-{\vec {q}}{\text{ }}'}$: now, shifting to polar coordinates the integrals in both the numerator and denominator involving ${\displaystyle \theta _{\vec {q}}}$ factorize and simplify, so in the end renaming ${\displaystyle x=|m_{\vec {q}}|}$ in the integrals for simplicity:
${\displaystyle \left\langle |m_{\vec {q}}|^{2}\right\rangle ={\frac {\int _{0}^{\infty }e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}x^{2}}x^{2}\cdot xdx}{\int _{0}^{\infty }e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}x^{2}}xdx}}}$
Noting that ${\displaystyle xdx=dx^{2}/2}$, if we change variable to ${\displaystyle y=x^{2}}$ and integrate we get[11]:
${\displaystyle \left\langle |m_{\vec {q}}|^{2}\right\rangle ={\frac {\int _{0}^{\infty }e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}y}ydy}{\int _{0}^{\infty }e^{-{\frac {k|{\vec {q}}|^{2}+at}{V}}}dy}}={\frac {V}{k|{\vec {q}}|^{2}+at}}}$

Therefore since the correlation function in Fourier space is non null only when ${\displaystyle {\vec {q}}=-{\vec {q}}{\text{ }}'}$, in general we can write:

${\displaystyle G_{{\vec {q}},{\vec {q}}{\text{ }}'}=\left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle ={\frac {V}{k|{\vec {q}}|^{2}+at}}\delta _{{\vec {q}},-{\vec {q}}{\text{ }}'}}$
Going back to real space we have:
{\displaystyle {\begin{aligned}G({\vec {r}},{\vec {r}}{\text{ }}')=\left\langle m_{\vec {r}}m_{{\vec {r}}{\text{ }}'}\right\rangle ={\frac {1}{V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}e^{i({\vec {q}}\cdot {\vec {r}}+{\vec {q}}{\text{ }}'\cdot {\vec {r}}{\text{ }}')}\left\langle m_{\vec {q}}m_{{\vec {q}}{\text{ }}'}\right\rangle =\\={\frac {1}{V^{2}}}\sum _{{\vec {q}},{\vec {q}}{\text{ }}'}e^{i({\vec {q}}\cdot {\vec {r}}+{\vec {q}}{\text{ }}'\cdot {\vec {r}}{\text{ }}')}{\frac {V\delta _{{\vec {q}},{\vec {q}}{\text{ }}'}}{k|{\vec {q}}|^{2}+at}}={\frac {1}{V}}\sum _{\vec {q}}{\frac {e^{i{\vec {q}}\cdot ({\vec {r}}-{\vec {r}}{\text{ }}')}}{k|{\vec {q}}|^{2}+at}}\end{aligned}}}
We see that appropriately substituting the sum with an integral (see the footnote on page ) and defining:
${\displaystyle \xi ={\sqrt {\frac {k}{at}}}}$
this correlation function acquires the same form of the one computed in mean field theory. This means that the critical exponents ${\displaystyle \nu }$ and ${\displaystyle \eta }$ now have the same values predicted by mean field theory, namely:
${\displaystyle \nu ={\frac {1}{2}}\qquad \qquad \eta =0}$

## Interaction between fluctuations: expansion to the fourth order

We have therefore seen that mean field theories can be improved including the fluctuations of the order parameter around its extremal values; in particular with the Gaussian approximation we have stopped the expansion at the second order and this led to a change in the critical exponent of the specific heat, which now really diverges instead of exhibiting a jump discontinuity as simple mean field theories predict. However, the quartic term ${\displaystyle \delta m^{4}}$ that we ignore within the Gaussian approximation (and which basically represent the interactions between the fluctuations) becomes crucial when we approach a critical point. We could thus wonder if the Gaussian approximation can be improved. In particular, reconsidering the expression of ${\displaystyle \beta {\mathcal {H}}}$ with ${\displaystyle h=0}$:

${\displaystyle \beta {\mathcal {H}}(m({\vec {r}}))=\int \left[{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}+{\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}\right]d^{d}{\vec {r}}}$
then keeping all the terms in ${\displaystyle \delta m}$ when expanding ${\displaystyle m}$ around ${\displaystyle m_{0}}$, remembering that ${\displaystyle atm_{0}+bm_{0}^{3}=0}$ and considering that the odd terms in ${\displaystyle \delta m}$ give no contributions (since we integrate an odd function over an even domain), we have:
${\displaystyle \beta {\mathcal {H}}=A_{0}V+\int \left[\left({\frac {a}{2}}t+{\frac {3}{2}}bm_{0}^{3}\right)\delta m^{2}+{\frac {k}{2}}\left({\vec {\nabla }}\delta m\right)^{2}\right]d^{d}{\vec {r}}+{\frac {b}{4}}\int \delta m^{4}d^{d}{\vec {r}}}$
A natural approach to compute the partition function ${\displaystyle Z}$ would now consist in using a perturbative method in order to expand the term ${\textstyle \exp \left(-{\frac {b}{4}}\int \delta m^{4}d^{d}{\vec {r}}\right)}$ in powers of the parameter ${\displaystyle b}$. This is of course a reasonable approach if ${\displaystyle b}$ is small; however, with a simple dimensional analysis we can show (and this is what we are going to do in the following) that for ${\displaystyle t\to 0}$ and ${\displaystyle d<4}$ this parameter diverges. The approach of the Gaussian approximation is therefore inconsistent, at least from this point of view.

### Dimensional analysis of Landau theory

We know that the partition function of the system~is:

${\displaystyle Z=\int e^{-\beta {\mathcal {H}}(m({\vec {r}}))}{\mathcal {D}}[m({\vec {r}})]}$
where, in the case ${\displaystyle h=0}$:
${\displaystyle \beta {\mathcal {H}}(m({\vec {r}}))=\int \left[{\frac {k}{2}}\left({\vec {\nabla }}m\right)^{2}+{\frac {a}{2}}tm^{2}+{\frac {b}{4}}m^{4}\right]d^{d}{\vec {r}}}$
It is now convenient to rescale the order parameter ${\displaystyle m}$ so that the term proportional to ${\displaystyle {\vec {\nabla }}m}$ has only a numerical coefficient. This can be done defining:
${\displaystyle \varphi =m{\sqrt {k}}\quad \qquad r_{0}={\frac {at}{k}}\quad \qquad u_{0}={\frac {b}{4k^{2}}}}$
so that:
${\displaystyle {\mathcal {H}}_{\text{eff.}}:=\beta {\mathcal {H}}=\int \left[{\frac {1}{2}}\left({\vec {\nabla }}\varphi \right)^{2}+{\frac {r_{0}}{2}}\varphi ^{2}+u_{0}\varphi ^{4}\right]d^{d}{\vec {r}}}$
where we have defined the dimensionless[12] effective Hamiltonian ${\displaystyle {\mathcal {H}}_{\text{eff.}}}$. Now, since ${\displaystyle {\mathcal {H}}_{\text{eff.}}}$ is dimensionless all the three integrals that appear in ${\displaystyle {\mathcal {H}}_{\text{eff.}}}$ must be so; this means that ${\displaystyle \varphi }$, ${\displaystyle r_{0}}$ and ${\displaystyle u_{0}}$ must have precise dimensions. In fact, from the first contribution we have that:
${\displaystyle \left[\int \left({\vec {\nabla }}\varphi \right)^{2}d^{d}{\vec {r}}\right]=1\quad \Rightarrow \quad L^{d}{\frac {[\varphi ]^{2}}{L^{2}}}=1\quad \Rightarrow \quad [\varphi ]=L^{1-d/2}}$
and similarly:
${\displaystyle \left[\int r_{0}\varphi ^{2}d^{d}{\vec {r}}\right]=\left[\int u_{0}\varphi ^{4}d^{d}{\vec {r}}\right]=1}$
from which we have that:
${\displaystyle [r_{0}]=L^{-2}\quad \qquad [u_{0}]=L^{d-4}}$
We can therefore use ${\displaystyle r_{0}}$ to define a length scale independent of the dimension of the system, i.e. ${\displaystyle r_{0}{}^{-1/2}}$. Since ${\displaystyle r_{0}\propto t}$ by definition and since ${\displaystyle \xi \sim t^{-1/2}}$ within mean field theories and Gaussian approximation, we see that this choice is equivalent to measuring lengths in units of the correlation length ${\displaystyle \xi }$ (which we know is independent of the dimension of the system). We now rescale all the variables, by setting:
${\displaystyle \phi ={\frac {\varphi }{\ell ^{1-d/2}}}\quad \qquad {\vec {x}}={\frac {\vec {r}}{\ell }}\quad \qquad u={\frac {u_{0}}{\ell ^{d-4}}}}$
where ${\displaystyle \ell =r_{0}{}^{-1/2}}$. This way, the partition function can be written in the form:
${\displaystyle Z=\int e^{-{\mathcal {H}}_{0}(\phi )-{\mathcal {U}}(\phi )}{\mathcal {D}}[\phi ]}$
where:
${\displaystyle {\mathcal {H}}_{0}=\int \left[{\frac {1}{2}}\left({\vec {\nabla }}\phi \right)^{2}+{\frac {1}{2}}\phi ^{2}\right]d^{d}{\vec {x}}\quad \qquad {\mathcal {U}}=\int u\phi ^{4}d^{d}{\vec {x}}}$
and ${\displaystyle {\mathcal {U}}}$, depending on the quartic term in ${\displaystyle \phi }$, is the contribution due to the interaction between the fluctuations. It is the presence of this terms that prevents us from computing ${\displaystyle Z}$ exactly and that forces us to resort to approximations.

The most standard procedure to apply in this case would be a perturbative method, namely to consider the dimensionless parameter ${\displaystyle u}$ as small, i.e. ${\displaystyle u\ll 1}$, and expanding the exponential term containing ${\displaystyle {\mathcal {U}}}$:

${\displaystyle Z=\int e^{-{\mathcal {H}}_{0}}\left(1-{\mathcal {U}}+{\frac {1}{2}}{\mathcal {U}}^{2}+\cdots \right){\mathcal {D}}[\phi ]}$
Written out explicitly, we have:
${\displaystyle u={\frac {u_{0}}{r^{(d-4)/2}}}=u_{0}\left({\frac {at}{k}}\right)^{(d-4)/2}}$
We thus immediately see that if ${\displaystyle d<4}$ then ${\displaystyle u}$ diverges when ${\displaystyle t\to 0}$, making the perturbative approach infeasible. On the other hand, if ${\displaystyle d>4}$ we indeed have ${\displaystyle u\to 0}$ when ${\displaystyle t\to 0}$, so the Gaussian approximation is actually a good one. Now, we can say that the perturbative expansion is reasonable if ${\displaystyle u\leq 1}$, which gives:
${\displaystyle t^{(d-4)/2}\geq u_{0}\left({\frac {a}{k}}\right)^{(d-4)/2}}$
This, using also the definition of ${\displaystyle u_{0}}$, can be rewritten as:
${\displaystyle t^{(4-d)/2}\geq {\frac {b}{4a^{2}\xi _{1}^{d}}}\qquad \qquad {\text{where}}\quad \xi _{1}:=\xi (t=1)={\sqrt {\frac {a}{k}}}}$
which is therefore a criterion that tells us if the perturbative approach is valid. Note that from this analysis we have determined the upper critical dimension of the system (${\displaystyle d=4}$), so we can say that this last criterion is equivalent to the Ginzburg criterion.

A final remark. The argument we have shown is not really convincing. In fact, we have only shown that every term of the perturbation theory diverges as ${\displaystyle t\to 0}$; however, this does not necessarily mean that the whole perturbation series is divergent. For example, consider the exponential series:

${\displaystyle e^{-x}=\sum _{n=0}^{\infty }{\frac {(-x)^{n}}{n!}}}$
This is convergent for ${\displaystyle x\to \infty }$ (${\textstyle \lim _{x\to \infty }e^{-x}=0}$), but each term of the sum diverges in the same limit. We thus understand that we must be careful when handling perturbation expansions, since the appropriate resummation of divergent terms can bring to a convergent series.

1. Sometimes this approximation is called elastic free energy.
2. At this point we could wonder why the interaction part of the Hamiltonian does not contain other terms, like ${\displaystyle m\nabla ^{2}m}$: in fact this is in principle perfectly acceptable since is of second order in ${\displaystyle m}$ and is invariant under rotations and parity ${\displaystyle m\to -m}$. However, we have that:
${\displaystyle m\nabla ^{2}m={\vec {\nabla }}\cdot (m{\vec {\nabla }}m)-({\vec {\nabla }}m)^{2}}$
and so when we integrate over ${\displaystyle {\vec {r}}}$, the first term vanishes in the thermodynamic limit supposing that the magnetization or its gradient goes to zero sufficiently rapidly as ${\displaystyle |{\vec {r}}|\to \infty }$. Therefore, we are left only with ${\displaystyle ({\vec {\nabla }}m)^{2}}$: the two terms are perfectly equivalent.
3. We will do our computations on the Ising model, as usual.
4. This follows from the integral form of ${\displaystyle {\tilde {\mathcal {H}}}}$. In fact, if in general a functional is of the form:
${\displaystyle F[\eta ]=\int f(\eta ,{\vec {\nabla }}\eta )d^{d}{\vec {r}}}$
then from the definition of functional derivative we have:
{\displaystyle {\begin{aligned}\int {\frac {\delta F}{\delta \eta }}\varphi ({\vec {r}})d^{d}{\vec {r}}={\frac {d}{d\varepsilon }}\left[\int f\left(\eta +\varepsilon \varphi ,{\vec {\nabla }}(\eta +\varepsilon \varphi )\right)d^{d}{\vec {r}}\right]_{|\varepsilon =0}=\\=\int \left({\frac {\partial f}{\partial \eta }}\varphi +{\frac {\partial f}{\partial {\vec {\nabla }}\eta }}\cdot {\vec {\nabla }}\varphi \right)d^{d}{\vec {r}}\end{aligned}}}
where ${\displaystyle \varphi ({\vec {r}})}$ is an arbitrary function that vanishes on the boundary of integration. Integrating by parts we get:
${\displaystyle \int {\frac {\delta F}{\delta \eta }}\varphi ({\vec {r}})d^{d}{\vec {r}}=\int \left({\frac {\partial f}{\partial \eta }}-{\vec {\nabla }}\cdot {\frac {\partial f}{\partial {\vec {\nabla }}\eta }}\right)\varphi ({\vec {r}})d^{d}{\vec {r}}}$
and so finally:
${\displaystyle {\frac {\delta F}{\delta \eta }}={\frac {\partial f}{\partial \eta }}-{\vec {\nabla }}\cdot {\frac {\partial f}{\partial {\vec {\nabla }}\eta }}}$
5. In the following, for the sake of simplicity we will indicate the magnitude of a vector simply removing the arrow sign.
6. In solid state physics this assumption is often called random phase approximation, while in field theory free field approximation.
7. We could have equivalently considered the case ${\displaystyle T, but it is a bit more complicated since ${\displaystyle m_{0}\neq 0}$ and so there is another term that contributes to the free energy. In other words, if ${\displaystyle T then in ${\displaystyle m}$ we are not considering the term with ${\displaystyle {\vec {k}}=0}$ (which is exactly equal to ${\displaystyle m}$).
8. The substitution:
${\displaystyle \sum _{|{\vec {q}}|<\Lambda }\longrightarrow {\frac {V}{(2\pi )^{d}}}\int _{0}^{\Lambda }d^{d}{\vec {q}}}$
can be justified as follows:
${\displaystyle \sum _{|{\vec {q}}|<\Lambda }=\sum _{|{\vec {q}}|<\Lambda }{\frac {\Delta {\vec {q}}}{\Delta {\vec {q}}}}=\left(\Delta {\vec {q}}\right)^{-1}\sum _{|{\vec {q}}|<\Lambda }\Delta {\vec {q}}}$
Now, ${\displaystyle {\vec {q}}}$ is quantized since ${\displaystyle V}$ is finite and we have periodic boundary conditions, and in particular ${\displaystyle \Delta {\vec {q}}=\Delta q_{1}\cdot \Delta q_{2}\cdots \Delta q_{d}=\left(2\pi /L\right)^{d}=(2\pi )^{d}/V}$. Therefore:
${\displaystyle \sum _{|{\vec {q}}|<\Lambda }={\frac {V}{(2\pi )^{d}}}\int _{|{\vec {q}}|<\Lambda }d{\vec {q}}}$
9. A little remark: the origin of this divergence does not come from the behaviour of the integral for large wavelengths. In fact, in the original definition of ${\displaystyle I_{1}}$ the integral has an upper limit, ${\displaystyle \Lambda }$, so it cannot diverge because of the large ${\displaystyle q}$ behaviour; if it diverges, it must be because of the behaviour for ${\displaystyle q\to 0}$, which corresponds to large wavelengths (this is also why this divergence is sometimes called infrared divergence).
10. We stress again that the same calculations could have been done in the case ${\displaystyle T, but we have not done so only for the sake of simplicity.
11. In order to compute the integral in the numerator we have used a standard "trick":
${\displaystyle \int _{0}^{\infty }e^{-kx}xdx=-{\frac {\partial }{\partial k}}\int _{0}^{\infty }e^{-kx}dx=-{\frac {\partial }{\partial k}}\left({\frac {e^{-kx}}{-k}}\right)_{0}^{+\infty }=-{\frac {\partial }{\partial k}}{\frac {1}{k}}={\frac {1}{k^{2}}}}$
12. Remember that ${\displaystyle \beta }$ has the dimension of the inverse of an energy.