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 ≥ C log N .