The standard procedure to generate initial conditions in numerical simulations of structure formations 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, particularly for velocity perturbations 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 the standard formulation of non‐linear perturbation theory (PT) to include transients as non‐linear excitations of decaying modes caused by the initial conditions. These new non‐linear solutions interpolate between the initial conditions and the late‐time solutions given by the exact non‐linear dynamics. To quantify the magnitude of transients, we focus on higher order statistics of the density contrast δ and velocity divergence Θ, characterized by the Sp and Tp parameters. These describe the non‐Gaussianity of the probability distribution through its connected moments 〈δp〉c ≡ Sp〈δ2〉p−1, 〈Θp〉c ≡ Tp 〈Θ2〉p−1. We calculate Sp(a) and Tp(a) to leading order in PT with top‐hat smoothing as a function of the scale factor a. We find that the time‐scale of transients is determined, at a given order p, by the effective spectral index neff. The skewness factor S3 (T3) attains 10 per cent accuracy only after a ≈ 6 (a ≈ 15) for neff ≈ 0, whereas higher (lower) neff demands more (less) expansion away from the initial conditions. These requirements become much more stringent as p increases, always showing slower decay of transients for Tp than Sp. For models with density parameter Ω ≠ 1, the conditions above apply to the linear growth factor; thus an Ω = 0.3 open model requires roughly a factor of 2 larger expansion than a critical density model to reduce transients by the same amount. The predicted transients in Sp are in good agreement with numerical simulations. More accurate initial conditions can be achieved by using second‐order Lagrangian PT (2LPT), which reproduces growing modes up to second order and thus eliminates transients in the skewness parameters. We show that for p > 3 this scheme can reduce the required expansion by more than an order of magnitude compared to the ZA. Setting up 2LPT initial conditions requires only minimal, inexpensive changes to ZA codes. We suggest simple steps for its implementation.