We critically examine how well the evolution of large-scale density perturbations is followed in cosmological N-body simulations. We first run a large volume simulation and perform a mode-by-mode analysis in three-dimensional Fourier space. We show that the growth of large-scale fluctuations significantly deviates from linear-theory predictions. The deviations are caused by non-linear coupling with a small number of modes at largest scales owing to finiteness of the simulation volume. We then develop an analytic model based on second-order perturbation theory to quantify the effect. Our model accurately reproduces the simulation results. For a single realization, the second-order effect appears typically as \u27zig-zag\u27 patterns around the linear-theory prediction, which imprints artificial \u27oscillations\u27 that lie on the real baryon acoustic oscillations. Although an ensemble average of a number of realizations approaches the linear-theory prediction, the dispersions of the realizations remain large even for a large simulation volume of several hundred megaparsecs on a side. For the standard Λ cold dark matter (ΛCDM) model, the deviations from linear growth rate are as large as 10 per cent for a simulation volume with L = 500 h-1 Mpc and for a bin width in wavenumber of Δk = 0.005 h Mpc-1, which are comparable to the intrinsic variance of Gaussian random realizations. We find that the dispersions scales as α L-3/4 Δ-1/2 and the mean dispersion amplitude can be made smaller than a per cent only if we use a very large volume of L \u3e 2h-1 Gpc. The finite box size effect needs to be appropriately taken into account when interpreting results from large-scale structure simulations for future dark energy surveys using baryon acoustic oscillations
