  Improving the theoretical description of galaxy clustering on small scales is
an important challenge in cosmology, as it can considerably increase the
scientific return of forthcoming galaxy surveys -- e.g. tightening the bounds
on neutrino masses and deviations from general relativity. In this paper, we
propose and test a new model for the clustering of galaxies that is able to
accurately describe redshift-space distortions even down to small scales. This
model corresponds to a second-order perturbative Lagrangian bias expansion
which is advected to Eulerian space employing a displacement field extracted
from $N$-body simulations. Eulerian coordinates are then transformed into
redshift space by directly employing simulated velocity fields augmented with
nuisance parameters capturing various possible satellite fractions and
intra-halo small-scale velocities. We quantify the accuracy of our approach
against samples of physically-motivated mock galaxies selected according to
either Stellar Mass (SM) or Star Formation Rate (SFR) at multiple abundances
and at $z=0$ and $1$. We find our model describes the monopole, quadrupole, and
hexadecapole of the galaxy-power spectra down to scales of $k\approx 0.6
[h/$Mpc] within the accuracy of our simulations. This approach could pave the
way to significantly increase the amount of cosmological information to be
extracted from future galaxy surveys.
