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 'zig-zag' patterns around the linear-theory prediction, which imprints artificial 'oscillations' 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 A cold dark matter (ACDM) 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/2 Δk -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 > 2 h -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.