  We study the concentration of a degree-$d$ polynomial of the $N$ spins of a
general Ising model, in the regime where single-site Glauber dynamics is
contracting. For $d=1$, Gaussian concentration was shown by Marton (1996) and
Samson (2000) as a special case of concentration for convex Lipschitz
functions, and extended to a variety of related settings by e.g., Chazottes et
al. (2007) and Kontorovich and Ramanan (2008). For $d=2$, exponential
concentration was shown by Marton (2003) on lattices. We treat a general fixed
degree $d$ with $O(1)$ coefficients, and show that the polynomial has variance
$O(N^d)$ and, after rescaling it by $N^{-d/2}$, its tail probabilities decay as
$\exp(- c\, r^{2/d})$ for deviations of $r \geq C \log N$.
