halo abundance can be expressed as [@Jenkins2001] $$\label{eqn:ms}
\frac {{\rm d}n} {{\rm d}M}(M,z) = f(\sigma) \frac {\rho_0}{M}
\frac {{\rm d~ln}[\sigma^{-1}(M,z)]}{{\rm d}M},$$ where $\rho_0$ is the background matter density at redshift $z=0$, and $\sigma^2(M,z)$ is the variance of the linear matter power spectrum over a length $R$, $$\label{eqn:gf}
\sigma^2(M,z) = \frac{D^2(z)}{2\pi^2} \int_0^\infty k^2P(k)W^2(kR(M))dk,$$ where $W(x)=3[\sin(x)-x\cos(x)]/x^3$ is the Fourier transformation of the top-hat filter, $R(M)=(3M/4\pi \rho_0)^{1/3}$ is the smoothed radius with halo mass $M$, and $D(z)$ is the growth factor.

This definition of mass function has been widely examined against N-body simulations and useful fitting functions $f(\sigma)$ have been provided by several authors [@Jenkins2001; @Warren2006; @Tinker2008; @Bhatt2010]. Nevertheless, deviations from the universality have been reported. For example, @Bhatt2010 found through their N-body simulations that, for their $\omega$CDM cosmological models (where the dark energy equation of state parameter $\omega$ is constant in time, but $\omega\ne-1$), the universality of the mass function is systematically broken at a level of $10\%$. Deviations from the universality are also detected in interacting dark energy models @Cui2012. @Cui2012 reported that deviations can exceed $\sim10\%$ for most of the models in the high mass end.

Fig. \[fig:massfunction\] shows the measured halo mass functions from our simulations for different $\zeta$ 