good agreement between $b_{hm}$ and $b_{hh}$.

In the $\Lambda$CDM cosmological model, the linear halo bias can be expressed as a function of $\nu$, [@Sheth2001; @Mo1996; @Seljak2004; @Tinker2010], where $\nu=\delta_c/\sigma$ is the ratio of the critical over-density required for collapse to the rms density fluctuation.

Fig. \[fig:bias\] shows measured halo biases from N-body simulations for different $\zeta$ values at the baseline redshift $z_S=1.2$. The error bars is the standard deviation of the $5$ largest wavelength modes. The error bars for different $\zeta$ values have very similar amplitude and only $\zeta=1.0$ has been plotted in the figure. The lower panel shows the ratio of the measured halo bias of MG models to the one in $\Lambda$CDM. We find that the difference of halo bias for different $\zeta$ values becomes larger as the increase of halo mass. However, the deviations of halo bias differ by less than $10\%$ from GR ($\zeta=1$) for all the $\zeta\ne1$ values used in our simulation. This behavior is expected due to the halo bias is a function of $\nu$ in $\Lambda$CDM model. The simulations for different $\zeta$ values have the same linear power spectrum and growth factor, therefore, according to Eq. (\[eqn:gf\]) they have the same density fluctuation 