We study the covariance matrix of the cluster mass function in cosmology. We adopt a two-line attack: first, we employ the counts-in-cells framework to derive an analytic expression for the covariance of the mass function. Secondly, we use a large ensemble of N-body simulations in the Λ cold dark matter framework to test this. Our theoretical results show that the covariance can be written as the sum of two terms: a Poisson term, which dominates in the limit of rare clusters; and a sample variance term, which dominates for more abundant clusters. Our expressions are analogous to those of Hu & Kravtsov for multiple cells and a single mass tracer. Calculating the covariance depends on: the mass function and bias of clusters, and the variance of mass fluctuations within the survey volume. The predictions show that there is a strong bin-to-bin covariance between measurements. In terms of the cross-correlation coefficient, we find r≳ 0.5 for haloes with M≲ 3 × 1014 h−1 M⊙ at z= 0. Comparison of these predictions with estimates from simulations shows excellent agreement. We use the Fisher matrix formalism to explore the cosmological information content of the counts. We compare the Poisson likelihood model, with the more realistic likelihood model of Lima & Hu, and all terms entering the Fisher matrices are evaluated using the simulations. We find that the Poisson approximation should only be used for the rarest objects, M≳ 5 × 1014 h−1 M⊙, otherwise the information content of a survey of size V∼ 13.5 h−3 Gpc3 would be overestimated, resulting in errors that are nearly two times smaller. As an auxiliary result, we show that the bias of clusters, obtained from the cluster-mass cross-variance, is linear on scales >50 h−1 Mpc, whereas that obtained from the auto-variance is non-linea
