Differences between critical temperatures in Ising Model The critical temperature $T_c$ for the 2D square lattice ferromagnetic Ising Model is known exactly to be $$T_c=\frac{2}{\ln(1+\sqrt{2})}$$.
For other geometries, say, a 3D cubic lattice or a 2D triangular lattice $T_c$ is significantly raised. I was thinking about the causes for this to happen. Considering the geometry is what has changed, I assume a big factor is the contribution of the number of neighbours ($4$ for square lattice, $6$ for cubic and triangular lattice). The reasoning:
"Increasing the number of neighbours means that the spin interaction contribution is higher with respect to the temperature fluctuation contribution. Therefore, it will be more difficult to go from an ordered state to a disordered state."
Would you say my reasoning above is correct? And if so, are there other factors to consider apart from the number of nearest neighbours?
 A: First of all, we should differentiate the effect of increased dimensionality $d$ and increased coordination $z$ (number of neighbors). While it is usually going to  be the case that in higher dimension $z$ is larger, in general one can construct high dimensional lattices with low coordination. For example, the square lattice in 2D has $z=4$, the cubic lattice in 3D has $z=6$, but the diamond lattice in 3D has $z=4$ and the Laves graph in 3D has $z=3$ (it is a diamond lattice with some bonds "deleted").
When it comes to criticality, dimensionality is extremely important. For example, the Mermin-Wagner theorem states that systems with continuous symmetries (e.g. Heisenberg models) cannot have symmetry breaking phase transitions in two dimensions, because the fluctuations in the symmetry broken phase generically destroy the ordered structure (On the other hand, the XY model does have a phase transition, but it is topological rather than symmetry-breaking, the famous Kosterlitz-Thouless transition). The Ising model can have a phase transition in 2D because the symmetry ($\mathbb{Z}_2$) is discrete, not continuous, but in fact it cannot have a phase transition in 1D.
Coordination number $z$ does of course have an effect on the transition temperature, but this is usually not necessarily of much interest, because the precise critical temperature is one of the "non-universal" properties of phase transitions, and will generically depend on any number of the microscopic details of the model.
As a simple example, consider an Ising model in dimension $d$ (e.g. on a hypercubic lattice), with Hamiltonian
$$H = -J \sum_{\langle ij \rangle} s_i s_j$$
where $J$ is an energy scale and $s_i = \pm 1$, where $i,j$ label the sites of the lattice and the sum is over nearest-neigbor pairs. Note that we can write this as
$$H = - \sum_i s_i \left(\frac{J}{2}\sum_{j@i} s_j\right)$$
where the second sum is over all neighbors of $i$. Clearly the second sum contains exactly $z$ terms.
Now we perform a standard mean field approximation, by replacing the quantity in parentheses by
$$ h_{\mathrm{mf}} \equiv J z \langle s \rangle/2,$$
where
$$\langle s \rangle = \left\langle \frac{1}{N} \sum_i s_i \right\rangle = \langle s_i \rangle$$
is the average magnetization per spin. The mean field approximation assumes that the magnetization $\langle s \rangle$ is non-zero, i.e. the global $\mathbb{Z}_2$ symmetry sending $s_i \to -s_i$ is spontaneously broken, and that the magnetization is uniform (the same at every site).  The mean field Hamiltonian is then simply
$$H_{\mathrm{mf}} = - \sum_i s_i h_{\mathrm{mf}} $$
All of the spins are now decoupled from each other, and their energy is simply determined by their orientation relative to the mean field of their neighbors (i.e. compare the energy of a magnetic dipole in an external field, $U = - \mathbf{m} \cdot \mathbf{h}$). We can now compute $\langle s_i \rangle$ self-consistently:
$$\langle s_i \rangle = \frac{\sum_{s_i=\pm 1} s_i e^{\beta h_{\mathrm{mf}} s_i}}{\sum_{s_i = \pm 1} e^{\beta h_{\mathrm{mf}} s_i}} = \tanh(\beta h_{\mathrm{mf}})$$
i.e.
$$\langle s \rangle = \tanh(Jz \langle s \rangle/2T)$$
For $T>T_c$, this equation a single solution $\langle s \rangle = 0$. For $T < T_c$, there appear two solutions with $\langle s \rangle \neq 0$, as seen below.

At $T=T_c$ the slope on each side of this equation matches, i.e. using $\tanh(x) \approx x$
$$T_c = Jz/2$$
so increasing the coordination does indeed affect $T_c$. Higher coordination means higher $T_c$, as you argued. The mean field approximation is obviously very crude but it yields practically exact results in high dimensions ($d\geq 4$ for continuous symmetries). Coordination has a role to play in stabilizing the ordered phase, and one could argue that the higher the coordination, the more neighbors there are to stabilize each spin, and thus the mean field approximation would be reasonable. But as I stated at the start, dimensionality is more important, and you can have low coordination in high dimension, so coordination number is not that important for critical behavior and is generically a non-universal characteristic.
