Non-linear density-velocity dynamics in $f(R)$ gravity from spherical
  collapse
  We investigate the joint density-velocity evolution in $f(R)$ gravity using
smooth, compensated spherical top-hats as a proxy for the non-linear regime.
Using the Hu-Sawicki model as a working example, we solve the coupled
continuity, Euler and Einstein equations using an iterative hybrid
Lagrangian-Eulerian scheme. The novel aspect of this scheme is that the metric
potentials are solved for analytically in the Eulerian frame. The evolution is
assumed to follow GR at very early epochs and switches to $f(R)$ at a
pre-determined epoch. Choosing the `switching epoch' too early is
computationally expensive because of high frequency oscillations; choosing it
too late potentially destroys consistency with $\Lambda$CDM. To make an
informed choice, we perform an eigenvalue analysis of the background model
which gives a ballpark estimate of the magnitude of oscillations. There are two
length scales in the problem: the comoving Compton wavelength of the associated
scalar field and the width of the top-hat. The evolution is determined by their
ratio. When the ratio is large, the evolution is scale-independent and the
density-velocity divergence relation (DVDR) is unique. When the ratio is small,
the evolution is very close to GR, except for the formation of a spike near the
top-hat edge, a feature which has been noted in earlier literature. We are able
to qualitatively explain this feature in terms of the analytic solution for the
metric potentials, in the absence of the chameleon mechanism. In the
intermediate regime, the evolution is profile-dependent and no unique DVDR
exists.
