  Inaccuracies in the initial conditions for cosmological N-body simulations
could easily be the largest source of systematic error in predicting the
non-linear large-scale structure. From the theory side, initial conditions are
usually provided by using low-order truncations of the displacement field from
Lagrangian perturbation theory, with the first and second-order approximations
being the most common ones. Here we investigate the improvement brought by
using initial conditions based on third-order Lagrangian perturbation theory
(3LPT). We show that with 3LPT, truncation errors are vastly suppressed,
thereby opening the portal to initializing simulations accurately as late as
z=12 (for the resolution we consider). We analyse the competing effects of
perturbative truncation and particle discreteness on various summary
statistics. Discreteness errors are essentially decaying modes and thus get
strongly amplified for earlier initialization times. We show that late starting
times with 3LPT provide the most accurate configuration, which we find to
coincide with the continuum fluid limit within 1 per cent for the power- and
bispectrum at z=0 up to the particle Nyquist wave number of our simulations (k
~ 3h/Mpc). In conclusion, to suppress non-fluid artefacts, we recommend
initializing simulations as late as possible with 3LPT. We make our 3LPT
initial condition generator publicly available.
