  We study the prediction of the dark matter power spectrum at two-loop order
in the Effective Field Theory of Large Scale Structures (EFTofLSS) using high
precision numerical simulations. In our universe, short distance non-linear
perturbations, not under perturbative control, affect long distance
fluctuations through an effective stress tensor that needs to be parametrized
in terms of counterterms that are functions of the long distance fluctuating
fields. We find that at two-loop order it is necessary to include three
counterterms: a linear term in the over density, $\delta$, a quadratic term,
$\delta^2$, and a higher derivative term, $\partial^2\delta$. After the
inclusion of these three terms, the EFTofLSS at two-loop order matches
simulation data up to $k\simeq 0.34 \,h\, {\rm Mpc}^{-1}$ at redshift $z=0$, up
to $k\simeq 0.55\,h\, {\rm Mpc}^{-1}$ at $z=1$, and up to $k\simeq 1.1\,h\,
{\rm Mpc}^{-1}$ at $z=2$. At these wavenumbers, the cosmic variance of the
simulation is at least as small as $10^{-3}$, providing a high precision
comparison between theory and data. The actual reach of the theory is affected
by theoretical uncertainties associated to not having included higher order
terms in perturbation theory, for which we provide an estimate, and by
potentially overfitting the data, which we also try to address. Since in the
EFTofLSS the coupling constants associated with the counterterms are unknown
functions of time, we show how a simple parametrization gives a sensible
description of their time-dependence. Overall, the $k$-reach of the EFTofLSS is
much larger than previous analytical techniques, showing that the amount of
cosmological information amenable to high-precision analytical control might be
much larger than previously believed.
