values at the baseline redshift $z_S=1.2$, which is that of $\zeta =1$ ($\Lambda$CDM). The $M_{200}$ halo mass is used in the plot in order to get consistent comparisons with concentration and bias, where the same mass definition are adopted. We have checked our measured halo mass functions for $\Lambda$CDM simulation, which agree well with the fitting function proposed by [@Sheth2001].

In order to assess the difference, the Poisson error bars are only added for $\zeta=1$. Because of the same initial condition been used, the error bars for $\zeta\ne1$ have the similar forms. In Fig. \[fig:massfunction\], the lower panel shows the ratio of the measured mass functions of MG models to the one in $\Lambda$CDM. The Poisson error bars for $\zeta=1$ are added to assess the maximum limit. Since $D(z_\zeta,\zeta)=D(z_S,\zeta=1)$, the linear power spectrum and hence $\sigma(M,z)$ are identical. If the halo mass function is indeed universal and described by Eq. (\[eqn:ms\]) and Eq. (\[eqn:gf\]), the mass function in the MG model should be the same as that in $\Lambda$CDM. However, the result in Fig. \[fig:massfunction\] shows the difference of mass functions between $\zeta=0.8$ and $\zeta=1.0$($\Lambda$CDM) is about $20\%$ with the halo mass range from $10^{11.5}\msun$ to $10^{14.0}\msun$. The difference between 