My question is: what is the fundamental reason why power laws should occur at phase transition specifically? What assumptions on the system should we make to obtain this result?
IMHO this is answered by reversing another sentence from the Q.:
It is often stated that this is a (the?) reason why the system is scale-free in critical regime.
When the system is in scale-free regime, power laws come naturally, just as exponents come naturally when we assume that correlations are local.
Indeed, the exponent arises, e.g., if we expect that a state at distance $1$ from the origin is influence to degree $1/l$ by the state at the origin. Then the state at distance $2$ is influenced to degree $1/l$ by the state at $1$, and to degree $1/l^2$ by the state at the origin, and so on.
This reasoning is inapplicable when the system does not have a specific scale of "influence" $\alpha$. Power laws is the simplest way of modeling
such scale-invariant state - they are not exact, just like correlation functions are only approximately exponential (note that Boltzmann exponent is also really a factorial in Stirling approximation, neglecting lots of minor details.)
Why there is no scale? Because we are talking about transition between two different phases - e.g., two states characterized by drastically different scales. We could model it, e.g., by a logistic function:
$$
l(T) = l_0 + \frac{\Delta l}{1+e^{N(T-T_0)/T_0}}
$$
Here for $|T-T_0|/T_0\gg 1$ we have scale $l_0$, if $T>T_0$ and scale $l_0+\Delta l$, if $T<T_0$. For finite $N$ the change may happen very abruptly, but for $T$ close to $T_0$ the scale is still finite. However, when we take thermodynamic limit $N\rightarrow\infty$ (assuming that for our temperature resolution we always have $N|T - T_0|/T_0\gg 1$), then the scale is completely undefined at $T=T_0$ and a different (scale-free) description invites itself.
Note also that in practice correlation lengths in different phases also have completely different nature - in a paramagnetic phase of ferromagnet we can meaningfully talk about the correlation between neighboring spins, whereas in a ferromagnetic phase we are more likely to talk of spin waves, characterized by wave number and group velocity, but extending throughout the whole crystal. Likewise, we can talk about distance correlations between molecules in a liquid, but phonons and Bloch waves in a solid - in other words, we are not simply talking about sudden change in the correlation length, as the above equation may suggest, but about one parameter vanishing and another one emerging.
Disclaimer:
Strongly correlated systems are not my field, but I find this question very interesting. So please do not treat the above as an expert/definitive opinion.
Update: Mechanical analogy
I think a very instructive example is a particle in a double-well potential,
$$V(x)=\lambda (x^2-d^2)^2$$
To avoid any confusion: I am not talking here about Landau theory of phase transitions, Kramers escape, or instantons - where one often sees such pictures. I discuss a non-linear oscillator, whose motion is described by Newton's equation:
$$
m\ddot{x}(t)=-\frac{d}{dx}V(x)=-4\lambda x(x^2 - d^2).
$$
In general case the period of motion can be calculated from the energy conservation:
$$
T=\sqrt{2m}\int_{x_L}^{x_R}\frac{dx}{\sqrt{E-V(x)}},
$$
where $x_L, x_R$ are the turning points, i.e., the roots of equation $V(x)=E$.
We can readily see the existence of two limiting regimes:
- $E<\lambda d^4$ For energies below the barrier energy the motion of the particle is confined to one of the wells, approaching at low energies the motion of a Harmonic oscillator with frequency:
$$
V(x)\longrightarrow=4\lambda d^2(x-d)^2=\frac{m\Omega^2x^2}{2} \text{ as } {x\rightarrow d},\\
\Omega = \sqrt{\frac{8\lambda d^2}{m}}.
$$
- $E>\lambda d^4$ For energies above the barrier the motion is close to that of a quartic oscillator with potential $V_4(x)=\lambda x^4$:
$$
T=\sqrt{\frac{8m}{\sqrt{\lambda E}}}\int_0^1\frac{dy}{\sqrt{1-y^4}}, \text{ where } \int_0^1\frac{dy}{\sqrt{1-y^4}}=\frac{\left[\Gamma\left(\frac{1}{4}\right)\right]^2}{4\sqrt{2\pi}}
$$
the value of the last integral can be found in Gradshtein&Ryzhik or evaluated by reducing it to a beta-function with replacement $y^4=t$ - the main point is that it is a constant, and the oscillations occur at finite frequency, even if energy dependent.
However, we see that there is no smooth transition between the two regimes described above: when the energy approaches the barrier, the motion of the particle becomes very slow, and its period diverges. Indeed, for $E=\lambda d^4$ there is an unstable equilibrium state at the top of the barrier, which means that the particle could stay there infinitely long (provided no noise or other external influences.)
One could show this by solving the Newton's equation, approximating the potential by an inverted parabola - since for $|E-\lambda d^4|\ll\lambda d^4$ the particle spends most of its time near the barrier, we can consider only a small interval around it, $[-x_0,x_0]$, where cutoff distance $x_0$ is an arbitrary point in the available region - that our result will be insensitive to its choice is the manifestation of scaling (for simplicity I treat only the case $E$ is slightly above the barrier):
$$
V(x)\approx\lambda d^4-2\lambda d^2x^2=V(0)-\frac{m\Omega^2x^2}{2},\\
T=\frac{4}{\Omega}\log\left[\frac{x_0}{a}+\sqrt{\left(\frac{x_0}{a}\right)^2+1}\right],
\text{ where } a=\frac{2[E-V(0)]}{m\Omega^2}=\frac{[E-V(0)]}{2\lambda d^2}
$$
since close to the barrier $a\rightarrow 0$, we have:
$$
T\approx \frac{4}{\Omega}\log\left[\frac{2\sqrt{2}x_0}{d}\sqrt{\frac{V(0)}{E-V(0)}}\right]=
\text{const}-\frac{2}{\Omega}\log\left(\frac{E-V(0)}{V(0)}\right)
$$
We see here scaling as $\log\left(\frac{E-V(0)}{V(0)}\right)$. Note in passing that logarithm can be viewed as a power law with extremely small exponent:
$$
\frac{x^\alpha - 1}{\alpha}\rightarrow\log x \text{ for }\alpha\rightarrow 0.
$$
This property is known to physicists as a replica trick, and to statisticians as Box-Cox transformation.
Since the logarithm diverges for energies close to the to of the barrier, the value of the constant, dependent on the chosen cutoff is irrelevant, i.e., we are free from this scale.
In fact, using Gradshtein&Ryzhik, one can obtain a closed expression for the period in terms of the complete elliptic integrals of the first kind:
$$
T=\begin{cases}
2\sqrt{\frac{m}{\lambda d^2}\sqrt{\frac{\lambda d^4}{E}}}
K\left(\sqrt{\frac{1+\sqrt{\frac{\lambda d^4}{E}}}{2}}\right), \text{ if }
E>\lambda d^4,\\
\sqrt{\frac{2m}{\lambda d^2\left(1+\sqrt{\frac{E}{\lambda d^4}}\right)}}K\left(\sqrt{\frac{2}{1+\sqrt{\frac{E}{\lambda d^4}}}}\right)
\end{cases}
$$
This allows us to plot the dependence of the period of energy, using any software familiar with special functions:
The spike at $E=\lambda d^4$ is the region where the period of the oscillations diverges, as we move from a quasi-harmonic to quartic oscillations.
