  We compute the density and velocity power spectra at next-to-next-to-leading
order taking into account the effect of time- and scale-dependent growth of
massive neutrino perturbations as well as the departure from
Einstein--de-Sitter (EdS) dynamics at late times non-linearly. We determine the
impact of these effects by comparing to the commonly adopted approximate
treatment where they are not included. For the bare cold dark matter
(CDM)+baryon spectrum, we find percent deviations for $k\gtrsim
0.17h~\mathrm{Mpc}^{-1}$, mainly due to the departure from EdS. For the
velocity and cross power spectrum the main difference arises due to time- and
scale-dependence in presence of massive neutrinos yielding percent deviation
above $k\simeq 0.08, 0.13, 0.16h~\mathrm{Mpc}^{-1}$ for $\sum m_{\nu} = 0.4,
0.2, 0.1~\mathrm{eV}$, respectively. We use an effective field theory (EFT)
framework at two-loop valid for wavenumbers $k \gg k_{\mathrm{FS}}$, where
$k_{\mathrm{FS}}$ is the neutrino free-streaming scale. Comparing to Quijote
N-body simulations, we find that for the CDM+baryon density power spectrum the
effect of neutrino perturbations and exact time-dependent dynamics at late
times can be accounted for by a shift in the one-loop EFT counterterm,
$\Delta\bar{\gamma}_1 \simeq - 0.2~\mathrm{Mpc}^2/h^2$. We find percent
agreement between the perturbative and N-body results up to $k\lesssim
0.12h~\mathrm{Mpc}^{-1}$ and $k\lesssim 0.16h~\mathrm{Mpc}^{-1}$ at one- and
two-loop order, respectively, for all considered neutrino masses $\sum m_{\nu}
\leq 0.4~\mathrm{eV}$.
