We investigate the properties of dark matter haloes and subhaloes in an $f(R)$ gravity model with $|f_{R0}|=10^{-6}$, using a very high-resolution N-body simulation. The model is a borderline between being cosmologically interesting and yet still con
sistent with current data. We find that the halo mass function in this model has a maximum 20% enhancement compared with the $Lambda$CDM predictions between $z=1$ and $z=0$. Because of the chameleon mechanism which screens the deviation from standard gravity in dense environments, haloes more massive than $10^{13}h^{-1}M_odot$ in this $f(R)$ model have very similar properties to haloes of similar mass in $Lambda$CDM, while less massive haloes, such as that of the Milky Way, can have steeper inner density profiles and higher velocity dispersions due to their weaker screening. The halo concentration is remarkably enhanced for low-mass haloes in this model due to a deepening of the total gravitational potential. Contrary to the naive expectation, the halo formation time $z_f$ is later for low-mass haloes in this model, a consequence of these haloes growing faster than their counterparts in $Lambda$CDM at late times and the definition of $z_f$. Subhaloes, especially those less massive than $10^{11}h^{-1}M_odot$, are substantially more abundant in this $f(R)$ model for host haloes less massive than $10^{13}h^{-1}M_odot$. We discuss the implications of these results for the Milky Way satellite abundance problem. Although the overall halo and subhalo properties in this borderline $f(R)$ model are close to their $Lambda$CDM predictions, our results suggest that studies of the Local Group and astrophysical systems, aided by high-resolution simulations, can be valuable for further tests of it.
