We use N-body simulations to study the nonlinear structure formation in
brane-induced gravity, developing a new method that requires alternate use of
Fast Fourier Transforms and relaxation. This enables us to compute the
nonlinear matter power spectrum and bispectrum, the halo mass function, and the
halo bias. From the simulation results, we confirm the expectations based on
analytic arguments that the Vainshtein mechanism does operate as anticipated,
with the density power spectrum approaching that of standard gravity within a
modified background evolution in the nonlinear regime. The transition is very
broad and there is no well defined Vainshtein scale, but roughly this
corresponds to k_*~ 2 at redshift z=1 and k_*~ 1 at z=0. We checked that while
extrinsic curvature fluctuations go nonlinear, and the dynamics of the
brane-bending mode C receives important nonlinear corrections, this mode does
get suppressed compared to density perturbations, effectively decoupling from
the standard gravity sector. At the same time, there is no violation of the
weak field limit for metric perturbations associated with C. We find good
agreement between our measurements and the predictions for the nonlinear power
spectrum presented in paper I, that rely on a renormalization of the linear
spectrum due to nonlinearities in the modified gravity sector. A similar
prediction for the mass function shows the right trends. Our simulations also
confirm the induced change in the bispectrum configuration dependence predicted
in paper I.Comment: 19 pages, 13 figures. v2: corrected typos, added more simulations,
better test of predictions in large mass regime. v3: minor changes, published
versio