  Current constraints on f(R) gravity from the large-scale structure are at the
verge of penetrating into a region where the modified forces become nonlinearly
suppressed. For a consistent treatment of observables at these scales, we study
cluster quantities produced in chameleon and linearized Hu-Sawicki f(R) gravity
dark matter N-body simulations. We find that the standard Navarro-Frenk-White
halo density profile and the radial power law for the pseudo phase-space
density provide equally good fits for f(R) clusters as they do in the Newtonian
scenario. We give qualitative arguments for why this should be the case. For
practical applications, we derive analytic relations, e.g., for the f(R) scalar
field, the gravitational potential, and the velocity dispersion as seen within
the virialized clusters. These functions are based on three degrees of freedom
fitted to simulations, i.e., the characteristic density, scale, and velocity
dispersion. We further analyze predictions for these fitting parameters from
the gravitational collapse and the Jeans equation, which are found to agree
well with the simulations. Our analytic results can be used to consistently
constrain chameleon f(R) gravity with future observations on virialized cluster
scales without the necessity of running a large number of simulations.
