Is it true that all correlation functions of any even number of spins in the ferromagnetic Ising model with nearest neighbors interaction are nonnegative in any spatial dimension?
In the one-dimensional case, it is easy to check all correlations directly. In the general case, the positivity of all even correlations in the ferromagnetic model looks obvious. Is there simple yet rigorous proof of this "obvious" fact?
Yes, it is a simple consequence of the first Griffiths (or GKS) inequality, which states that $$ \langle \sigma_A \rangle \geq 0 $$ for any finite set of vertices $A$. Above, I have used the standard notation $\sigma_A = \prod_{i\in A}\sigma_i$.
Griffiths' first inequality holds at any temperature and for any nonnegative magnetic field. Actually, it also holds for (ferromagnetic) interactions of arbitrary (possibly infinite) range. The proof is very easy (one simply expands the Boltzmann weight in a Taylor series and sum the resulting expression over the spins) and can be found, for instance, in Section 3.8.1 of this book (the version given there actually covers a substantially more general situation than described here).
In fact, one can even show that $$ \langle \sigma_A \rangle > 0 $$ at all finite temperatures, for any finite set $A$ containing an even number of vertices. This stronger version follows from the second Griffiths inequality stated below (see Exercise 3.25 in the book for the 2-point function and apply the second Griffiths inequality for general $A$ containing an even number of vertices).
Maybe one comment: it is crucial that the boundary condition is free, periodic or $+$ (or a mixture of $+$ and free, for instance). You cannot use a mixture of $+$ and $-$ boundary conditions (such as Dobrushin boundary conditions), as one can construct counterexamples to the claim in this case.
For pure $-$ boundary conditions, the result remains true when $h=0$ and $A$ contains an even number of vertices, since in that case the expectation coincides with the corresponding expectation under $+$ boundary condition by symmetry.
Just to be complete, the inequality is called the first Griffiths inequality, because there is a second one, which you might also find interesting: under the same assumptions, for any finite sets of vertices $A$ and $B$, $$ \langle \sigma_A\sigma_B \rangle \geq \langle \sigma_A \rangle \langle \sigma_B \rangle , $$ which shows that the random variables $\sigma_A$ and $\sigma_B$ are positively correlated (their covariance is nonnegative).
