  Numerical simulations in modified gravity have commonly been performed under
the quasi-static approximation -- that is, by neglecting the effect of time
derivatives in the equation of motion of the scalar field that governs the
fifth force in a given modified gravity theory. To test the validity of this
approximation, we analyse the case of $f(R)$ gravity beyond this quasi-static
limit, by considering effects, if any, these terms have in the matter and
velocity divergence cosmic fields. To this end, we use the adaptive mesh
refinement code ECOSMOG to study three variants ($|f_{R}|= 10^{-4}[$F4$],
10^{-5}[$F5$]$ and $10^{-6}[$F6$]$) of the Hu-Sawicki $f(R)$ gravity model,
each of which refers to a different magnitude for the scalar field that
generates the fifth force. We find that for F4 and F5, which show stronger
deviations from standard gravity, a low-resolution simulation is enough to
conclude that time derivatives make a negligible contribution to the matter
distribution. The F6 model shows a larger deviation from the quasi-static
approximation, but one that diminishes when re-simulated at higher resolution.
We therefore come to the conclusion that the quasi-static approximation is
valid for the most practical applications in $f(R)$ cosmologies.
