  The standard procedure to generate initial conditions (IC) in numerical
simulations is to use the Zel'dovich approximation (ZA). Although the ZA
correctly reproduces the linear growing modes of density and velocity
perturbations, non-linear growth is inaccurately represented because of the ZA
failure to conserve momentum. This implies that it takes time for the actual
dynamics to establish the correct statistical properties of density and
velocity fields. We extend perturbation theory (PT) to include transients as
non-linear excitations of decaying modes caused by the IC. We focus on
higher-order statistics of the density contrast and velocity divergence,
characterized by the S_p and T_p parameters. We find that the time-scale of
transients is determined, at a given order p, by the spectral index n. The
skewness factor S_3 (T_3) attains 10% accuracy only after a=6 (a=15) for n=0,
whereas higher (lower) n demands more (less) expansion away from the IC. These
requirements become much more stringent as p increases. An Omega=0.3 model
requires a factor of two larger expansion than an Omega=1 model to reduce
transients by the same amount. The predicted transients in S_p are in good
agreement with numerical simulations. More accurate IC can be achieved by using
2nd order Lagrangian PT (2LPT), which reproduces growing modes up to 2nd order
and thus eliminates transients in the skewness. We show that for p>3 this
reduces the required expansion by more than an order of magnitude compared to
the ZA. Setting up 2LPT IC only requires minimal, inexpensive changes to ZA
codes. We suggest simple steps for its implementation.
