Cosmic variance limited Baryon Acoustic Oscillations from the DEUSFUR CDM simulation
Abstract
We investigate the nonlinear evolution of Baryon Acoustic Oscillations (BAO) in the lowredshift matter power spectrum from the DEUSFUR CDM model simulation. This is the first cosmological Nbody simulation encompassing the full observable cosmic volume, thus allowing cosmic variance limited predictions at BAO scales. We control the effect of numerical systematic errors using a series of large volume highresolution simulations. The combined analysis allows us to measure the matter power spectrum between and to over the entire BAO range, , in bins of size . We define the BAO with respect to a nonlinearly evolved wigglefree spectrum and determine the characteristics of the BAO without recurring to extrapolation from global fitting functions. We quantify the effects of nonlinearities on the position and amplitude of the BAO extrema, and the coupling to the broadband slope of the power spectrum. We use these estimates to test nonlinear predictions from semianalytical models. Quite remarkably from the analysis of the redshift evolution of BAO we find that the second dip and third peak remains unaltered by nonlinear effects. Furthermore, we find that the square of the damping factor and the shift of the position of BAO extrema scale to good approximation as the square of the growth factor, in agreement with expectations from perturbation theory. This confirms the idea that, besides cosmic distances, an accurate measurement of BAO at different redshifts can directly probe the growth of cosmic structures.
keywords:
Cosmology; Nbody Simulations; LargeScale Structures;1 Introduction
The propagation of primeval acoustic waves in the coupled photonbaryon plasma before recombination (Sakharov, 1965; Silk, 1968; Peebles & Yu, 1970) generates a distinct pattern of temperature and polarization anisotropies in the Cosmic Microwave Background (CMB). These have been measured with unprecedented accuracy by observations of the Wilkinson Microwave Anisotropy Probe (WMAP) (Spergel et al., 2003, 2007; Komatsu et al., 2009) and more recently by the Planck satellite (Planck Collaboration XVI, 2013).
A similar imprint is present in the late time distribution of largescale structures in the form of an oscillatory pattern in the matter power spectrum, the socalled Baryon Acoustic Oscillations (BAO). These have been detected using measurements of twopoint galaxy correlation function from the Sloan Digital Sky Survey (SDSS) (Eisenstein et al., 2005; Huetsi, 2006) and the galaxy power spectrum from the 2degree Field (2dF) survey (Cole et al., 2005).
CMB observations have accurately determined the distance travelled by acoustic waves at decoupling (i.e. the sound horizon). Thus, the detection of BAO in the galaxy distribution can be used as a standard ruler to estimate the angular diameter distance as well as the Hubble rate at redshifts probed by galaxy surveys (see e.g. Eisenstein, Hu & Tegmark, 1998; Blake & Glazebrook, 2003; Seo & Eisenstein, 2003). A new generation of galaxy surveys has been designed to measure the BAO at different redshifts to few percent error from which it will be possible to infer stringent bounds on the cosmological parameters (see e.g. Ivezic et al., 2008; Amendola et al., 2012; Dawson et al., 2013). This demands for equally accurate cosmological model predictions.
In the linear regime of cosmic structure formation BAO appear as a series of damped oscillations superimposed to the broadband shape of the matter power spectrum. These can be computed to the desired level of accuracy using linear perturbation theory (Eisenstein & Hu, 1998; Lewis, Challinor & Lasenby, 2000). However, at low redshifts the onset of the nonlinear clustering of matter induces deviations from the linear theory which are much harder to predict. Such nonlinearities degrade the BAO pattern by altering the amplitude, shape and position of the extrema as well as the broadband slope of the power spectrum. If not accounted these effects can alter the cosmological parameter inference (see e.g. Angulo et al., 2008).
Nonlinear effects have been investigated using numerical Nbody simulations in numerous studies (e.g. Seo & Eisenstein, 2005; Seo et al., 2008; Angulo et al., 2008; Jeong & Komatsu, 2009; Carlson et al., 2009; Nishimichi et al., 2009; Seo et al., 2010; Orban & Weinberg, 2011). Alternatively, semianalytical approaches have been developed to compute the nonlinear modifications of the matter power spectrum at intermediate scales (Nishimichi et al., 2007; Crocce & Scoccimarro, 2008; Padmanabhan & White, 2009; Sherwin & Zaldarriaga, 2012; Taruya et al., 2012). The comparison with Nbody simulations has shown these to provide accurate predictions on BAO scales below hMpc (see e.g. Carlson et al., 2009; Nishimichi et al., 2009).
However, BAO occur on a relatively large scale ( hMpc) with a shape that depends on smallscale lagrangian displacements hMpc (Eisenstein et al., 2007a). Hence, the main limitation of numerical studies is due to sample variance errors on the one hand and the lack of large dynamical range on the other hand. For instance, Seo et al. (2010) have used a series of 5000 simulations with particles with box size hGpc from Takahashi et al. (2009) corresponding to an effective volume of hGpc. This allows them to drastically reduce sample variance errors on the matter power spectrum, nevertheless these are still nearly twice as large as cosmic variance uncertainty, while BAO remain under sampled in Fourier space. Because of this, Nbody studies of the effects of nonlinearities on the BAO have been mainly estimated from functional fits to the numerical power spectrum. Furthermore, averaging over a large ensemble of simulations does not resolve the problem of numerical systematics that needs to be carefully addressed especially when statistical errors are reduced to the percent level. These can be particularly important in the case of low mass and spatial resolution runs such as those from Takahashi et al. (2009). To date, Nbody simulation studies of the BAO accurate to percent level are still missing.
The work presented here aims to fill this gap. To this purpose we use a large volume Nbody simulation of a flat CDM cosmology from the Dark Energy Universe Simulation  Full Universe Runs (DEUSFUR) (Alimi et al., 2012). This is one of the three DEUSFUR simulations of Dark Energy models with particles and box length of hGpc. The size of the simulation box allows us to derive cosmic variance limited predictions for the matter power spectrum. However, in order to control the effect of numerical systematics we perform a convergence study using a set of Nbody simulations for which we have varied several numerical simulation parameters.
The combined analysis allows us for the first time to determine the effect of nonlinearities on all BAO peaks and troughs between redshift and to accuracy level (and in bins of size ). It is in this redshift interval that the power spectrum is most difficult to predict. On the other hand it is the probed range of current and future observational missions (to give an idea the median redshift of the planned Euclid survey is ).
We use the DEUSFUR power spectrum to test the validity of semianalytical models on the BAO scale. In addition, we precisely characterize the nonlinear modifications of the BAO extrema without resorting to curvefitting procedures. We show that the shift of the position of the BAO extrema remains a few percent effect at all redshift. In contrast, the damping of the BAO amplitude is the dominant modification induced by the nonlinear clustering of matter. This evolves as function of redshift proportionally to the linear growth rate, thus suggesting that measurements of the amplitude of BAO extrema at different redshifts can be used as an independent cosmological probe carrying information on the linear growth of structures.
The paper is organized as follows. In Section 2 we summarize the characteristics of the Nbody simulations used in our analysis, while in Section 3 we present a detailed study of the statistical and systematic errors affecting the measurement of the matter power spectrum. In Section 4 we define the BAO spectrum and discuss the comparison against semianalytical models. In Section 5 we quantify the nonlinear effects on the position and amplitude of the BAO extrema and the coupling to the broadband slope of the power spectrum. Finally, we discuss our conclusion in Section 6.
2 Nbody simulations
2.1 Simulation set
We use the DEUSFUR simulation of a flat CDM cosmology bestfit to the WMAP7yr data (Spergel et al., 2007). This has been realized using the application AMADEUS – A Multipurpose Application for Dark Energy Universe Simulation – expressly developed for the DEUSFUR project (Alimi et al., 2012). Gaussian initial conditions using the Zel’dovich approximation have been generated with an optimized version of the code MPGRAFIC (Prunet et al., 2008). The Nbody run has been performed with a version of RAMSES (Teyssier, 2002) that has been specifically improved to run on a large number of cores (). This is an Adaptive Mesh Refinement (AMR) code with a treebased structure that allows for recursive refinements above a userdefined density threshold. The Nbody solver evolves particles using a ParticleMesh (PM) method, while the Poisson equation is solved with a multigrid technique (Guillet & Teyssier, 2011). A detailed description of the algorithms, optimization schemes and the computing challenges involved with the realization of DEUSFUR is given in Alimi et al. (2012), while the characteristics of the three DEUSFUR simulations of Dark Energy models will be presented in a forthcoming paper.
The box length of the DEUSFUR simulations is set to hGpc, thus enclosing the horizon diameter of the three simulated cosmologies (hGpc in CDM).
These simulations employ particles and a coarse grid of resolution elements with levels of refinement, reaching a formal mass resolution of hM and a spatial resolution of hkpc (corresponding roughly to the mass and size of the Milky Way).
The initial redshift is set to , hence sufficiently high to avoid transient effects which occur in the case of a late start of the initial conditions using the Zel’dovich approximation (Scoccimarro, 1998; Crocce et al., 2006). The refinement threshold is set to such as to limit the total amount of AMR cells generated during the run (and the associated memory usage) without affecting the accuracy of the nonlinear power spectrum calculation. The final count of AMR cells is still of trillions. The model parameters of the simulated CDMW7 are quoted in Table 1, while the characteristics of the DEUSFUR simulation are listed in Table 2.
Model  h  

CDMW7  0.2573  0.72  0.801  0.963  0.04356 
L (hMpc)  n  n  n  m  n  m (hM)  (hkpc)  z  C  P 

21000  8192  6  14  40  106  0.2  CAMB 
L  n  m  z  C 

(hMpc)  
0.2  
0.2  
0.2  
0.2  
0.2  
0.2  
0.08  
0.5  

0.2  
0.2  
0.2  
0.2  

0.5  
0.5  
0.5  
0.5  
0.5  

0.2  
0.2 
We generate a set of Nbody simulations which we use to perform convergence tests and control the effect of numerical systematics on the DEUSFUR matter power spectrum. These consist of smaller volume simulations typically with particles in which we have varied the simulation box length, the refinement threshold, the starting redshift, the generation of initial conditions, the integration timestep and the mass resolution. The characteristics of these simulations are listed in Table 3. This benchmark is similar to that used in Reed et al. (2013), though in our case we have opted for higherresolution largervolume simulations.
In addition to this testbed of simulations, we have performed a Nbody run with an initial wigglefree linear matter power spectrum from Eisenstein & Hu (1998). We use the resulting spectrum as reference wigglefree power spectrum to define the BAO. The initial spectra for all other simulations has been computed using the CAMB code (Lewis, Challinor & Lasenby, 2000).
2.2 Power spectrum estimator
We compute the matter power spectrum from the Nbody runs with the code POWERGRID (Prunet et al., 2008) which we have optimized to handle a FFT grid with trillion elements. The FFT grid size is two times thinner than the coarse grid of the dynamical solver. The resulting power spectrum is corrected for the effect of smoothing due to the CellInCell (CIC) method. We do not perform any correction for the shot noise, which turns out to be completely negligible. Similarly, we do not correct for aliasing, since, by varying the size of the CIC grid we find that aliasing effects are negligible below half the Nyquist frequency of the CIC grid. This sets the computation of the power spectrum in the interval and , where is the simulation box length and is the spatial resolution at the coarse grid level.
3 Matter power spectra: statistical and systematic errors
We now present the DEUSFUR CDMW7 power spectrum and discuss the evaluation of statistical and systematic errors.
Figure 1 illustrates the evolution of the matter power spectrum between redshift to over a range of scales spanning nearly four order of magnitudes, from h Mpc to h Mpc.
One remarkable feature of the plotted spectra is the low level of noise. Because of the large size of the simulation box the noise is even on scales near the peak of the CDM power spectrum at h Mpc. We may also notice the excess of power at small scales ( h Mpc) due to the onset of the nonlinear clustering regime. At this alters by the prediction of the linear theory up to very large scales h Mpc. The imprint of the BAO is clearly distinguishable in all plotted spectra, extending over a decade from hMpc to hMpc with a peaktopeak amplitude of .
Numerical effects on the matter power spectrum have already been investigated in several works (see e.g. O’Shea et al., 2005; Heitmann et al., 2005, 2008; Joyce et al., 2009; Heitmann et al., 2010; Smith et al., 2012; Reed et al., 2013). Such studies have shown that the nature and amplitude of numerical systematics may depend on the specifics of the Nbody solver. Hence, it is critical for us to investigate numerical source of errors which may be specific to the solver we have used.
In principle, we could simply focus on the power spectrum normalized by the spectrum of the initial conditions. This would allow us to suppress the amplitude of statistical fluctuations on the large scale, but it will not reduce sample variance errors nor discount for possible numerical systematics, especially at small scales where nonlinearities induce phase correlations in Fourier space.
The error analysis that we present here extends previous studies to the case of very large volume simulations (hGpc) with large number of particles (). Our goal is to evaluate the error budget on the DEUSFUR matter power spectrum and select the range of scales where the power spectrum is determined to accuracy.
3.1 Statistical errors
The finite size of cosmological simulations introduces two types of statistical errors. Firstly, it reduces the number of accessible modes causing sample variance errors which dominate the error budget when probing scales near the size of the simulation box. Secondly, modes which are larger than the simulation box length are absent (Bagla & Ray, 2005) and since the gravitational coupling to these modes is missing, this results in a lower amplitude of the power spectrum (Heitmann et al., 2010). Both problems can be handled either by averaging the spectra of a large number of different realizations (e.g. Takahashi et al., 2008, 2009) or setting the simulation box length to be as large as possible. However, since observations are limited by the size of the cosmological horizon, statistical errors cannot be reduced to less than cosmic variance. Hence, by setting the box length of the DEUSFUR simulations to the diameter of the observable universe we are guaranteed to derive cosmic variance limited (i.e. minimal sample variance) predictions for the matter power spectrum.
In Figure 2 we plot the power spectrum at from the auxiliary simulations with varying from 1.3 hGpc to 10.5 hGpc relative to the DEUSFUR case. As expected, the statistical fluctuations on the power spectrum increase as the simulation box length decreases.
The rootmeansquare fluctuation of the initial power spectrum is given by (Jeong & Komatsu, 2009):
(1) 
where is the number of modes in a shell of size and is the mean particle number density. The propagation of this error into the nonlinear regime is nontrivial (Ngan et al., 2012), however we think that Eq. (1) still provides a reasonable approximation on the BAO scales. The error bars on the lefthand side of Figure 2 give a visual comparison of the expected statistical error at h Mpc for the different simulation box lengths. We can see that even a large volume simulation with box length of hGpc still leads to statistical errors of order of , whereas for DEUSFUR these are within .
We select the interval of interest by truncating the DEUSFUR power spectrum over modes where the statistical error given by Eq. (1) exceeds the level. This sets a lower bound h Mpc. Notice that for all simulations the measured fluctuations are indeed within level.
The finite size of the simulation box length is also relevant in determining the accuracy of the mode sampling, which is critical to accurately locate the position of the BAO extrema. In fact, the wavenumber binsize is given by , that sets the size of the “error bars” along the axis. These are illustrated in Figure 2 at h Mpc for different simulation box lengths. The arrows mark the wavenumber values where equals to of for a given simulation box length.
As already mentioned, running a large ensemble of smallbox simulations (as in the case of Takahashi et al., 2009) is effective to reduce statistical errors, however it does not improve the ksampling of the spectrum which is (usually) governed by the simulation box length. Common methods to reduce the uncertainty on the wavenumber value and thus the location of the BAO consist in using functional fitting procedures. For instance, one can use a functional form (e.g. polynomial expansions) to bestfit the measured power spectrum and subsequently infer the position of the BAO from the resulting bestfit function. Alternatively, one can bin the power spectrum or convolve it with a filter function. In any case, these procedures do not provide better information than that fixed by the size of (i.e. the simulation box length). Therefore the tradeoff of these methods is a loss of accuracy in the determination of the location of the BAO peaks and dips, which is key to the determination of cosmological distances. In fact, the position of the BAO as inferred from the bestfitting function may be sensitive to the specific choice of the functional form used to fit the power spectrum.
In the case of the DEUSFUR simulations this issue does not arise since the choice of the box length guarantees the ideal sampling.
3.2 Systematic errors
The determination of numerical systematic effects on the matter power spectrum is critical to reach the accuracy required by future surveys. Here, we assess their impact by comparing the DEUSFUR power spectrum with that obtained from the testbed simulations in which we have varied the mass resolution, the refinement strategy, the starting redshift, the generation of the initial conditions and the integration timestep. The results of this comparison are shown in Figure 3.
Mass resolution
The mass resolution is a key attribute of Nbody simulations. To date only a few studies have investigated the dependence of the matter power spectrum on the mass resolution of Nbody simulations (see e.g. Joyce et al., 2009; Heitmann et al., 2010). A brute force approach consisting in running simulations with fixed volume and increasingly large number of particles is optimal to evaluate the amplitude and mode dependence of this effect. It is in this way that Heitmann et al. (2010) have found that the power spectrum falls by at h Mpc for a mass resolution of hM.
Here, we perform a similar analysis by comparing the power spectra at from the testbed of different mass resolution simulations. The result of the comparison is shown in the topleft panel of Figure 3, where we plot the relative difference of the power spectrum (normalized to the linear prediction) of a given simulation relative to that of a reference case with mass resolution hM (black solid line) characterized by hMpc and particles. The simulations with mass resolution hM, similar to that of DEUSFUR have hMpc with particles (green solid line) and hMpc with (green dashed line) respectively. We also consider a simulation with a higherresolution with respect to the reference one with hM and characterized by hMpc with particles (blue dashed line).
We can see that the lower the mass resolution the larger the deviation of the power spectra at highk. In particular the reference simulation underestimates the power spectrum with respect to the higherresolution by up to h Mpc, while the simulations with the DEUSFUR resolution underestimates the power spectrum by which is compatible with the results by Heitmann et al. (2010).
This effect is a generic feature of codes used for large Nbody simulations which rely on a PM method to solve the largescale dynamics. The fact that the amplitude of this effect increases with does not mean that smallscale structures such as halos are not well resolved (some halos have 10000 particles), rather that most of the particles are still in underdense regions where the force calculation is not refined. Knebe et al. (2001) have performed Zel’dovich wave tests and shown that in void regions cells per particle are necessary to accurately follow the Zel’dovich wave, while inside virialized structures at least particles per cell are required to suppress the Poisson noise. However, such a refinement strategy is too expensive for large volume simulations such as DEUSFUR, since it would result in an extremely large number of AMR cells and consequently an excessive memory usage beyond computational capabilities. Nevertheless, one can correct for such an error by combining information from the higher resolution simulations.
As it can be seen from the plot in the topleft panel of Figure 3 the mass resolution effect on the matter power spectrum is a very smooth function of . It does not alter the BAO structure, but only affects the broadband shape. Hence, we can infer a precise estimate of the mass resolution effect by taking the ratio of the matter power spectrum from the simulation with hMpc box length and particles (that has the same mass resolution as DEUSFUR) to that of the simulation with same box length and particles (that has a higher resolution with hM). The correcting function for each redshift is then obtained by fitting this ratio with a polynomial as function of . As we can see in the topleft panel of Figure 3 the power spectrum of the simulation with particles and hMpc box length is within of the higher resolution simulation when divided by (red solid line) up to h Mpc. This sets the upper bound on the wavenumber interval of interest. Thus, the mass resolution effect is corrected in the DEUSFUR power spectrum by considering . We have tested that improving the mass resolution by another factor gives the same spectra at level. Therefore, we are confident that the estimated ratio accurately corrects the mass resolution effect on the DEUSFUR power spectrum as well.
The end result of this error analysis is the selection of the interval of interest h Mpc where the dominant systematic uncertainty due to mass resolution and the statistical errors are both within level. This is an unprecedented achievement that, as we will discuss in the next Section, allows us to precisely evaluate the effects of nonlinearities on the BAO.
Refinement strategy
The AMR algorithm allows for an accurate calculation of the force over a large dynamical range. The AMRtree dominates the memory usage, since it is nearly proportional to the number of grid cells. Thus, it can be a limiting factor when running large volume highresolution simulations. This has forced us to optimize the refinement strategy for DEUSFUR simulations, compromising between the accuracy on the matter power spectrum (and the mass function) and the memory usage. As an example, a simulation with hGpc box length and 8.6 billion particles at the end of the run has a total number of cells that varies between billions for a refinement threshold (number of particles per cell) and billions for a refinement threshold of 25. We have opted for a refinement threshold of 14 (corresponding to a total of 30 billion cells for the testbed simulation). As illustrated in the topright panel of Figure 3 the effect of the refinement does not alter the matter power spectrum by more than over the interval of interest and consequently it can be neglected.
Initial conditions
The generation of initial conditions and the choice of the starting redshift of the simulations are also potential sources of systematic errors. The effect of transients on the matter power spectrum (Scoccimarro, 1998) has been studied in detail by Crocce et al. (2006). These authors have shown that starting the initial conditions using the Zel’dovich approximation (ZA) at redshift leads to underestimated power spectra at the level (near h Mpc). This can be corrected by generating initial conditions using second order lagrangian perturbation theory (2LPT). On the other hand, Reed et al. (2013) have shown that using either ZA or 2LPT with a starting redshift leads to numerical errors of on the mass function and on the power spectrum (near h Mpc), whereas a starting redshift between and with 2LPT gives a converged mass function at the few percents level for halos larger than particles and a converged power spectrum at the percent level.
In the bottomleft panel of Figure 3 we plot the relative difference of the matter power spectrum for different initial condition generators and starting redshifts. We can see that using ZA or 2LPT at makes very little difference (). We have also explored a wide range of starting redshifts varying from to and found that differences increase only up to level between and h Mpc. Interestingly, for higher starting redshift () our PMAMR Nbody solver leads to a lower power spectrum on the large scales ( h Mpc). As a result of these tests we have set the DEUSFUR starting redshift to .
Time step
The integration timestep is another simulation parameter that may affect the accuracy of the matter power spectrum. RAMSES uses an adaptive timestep method where the timestep of the level is divided by a factor 2 compared to the time step of the level (Teyssier, 2002). The timestep at the coarse level is given by the minimum between (where is the Hubble constant), (where is the local freefall time) and (where is the spatial resolution, and is the maximum velocity of dark matter particles). At low redshift, the latter condition is the more stringent. The default value of the Courantlike factor is . For comparison we have run testbed simulations with such “large” timestep divided by a factor of () and (). In the bottomright panel of Figure 3 we can clearly see that the effect of these different timesteps on the matter power spectrum remains negligible. In the case of the DEUSFUR simulations we have set the Courantlike factor to since this allows us to increase the time resolution of the lightcone data.
4 NonLinear Evolution of Baryon Acoustic Oscillations
4.1 BAO Spectrum
The oscillatory pattern that characterizes the imprint of BAO on the matter power spectrum is usually defined relative to a smooth (wigglefree) function, , that accounts for the broadband slope of the underlying spectrum. Clearly, the choice of this function becomes critical if we aim to achieve accuracy on the BAO. This is because the relative amplitude of BAO is of order , hence even a change as small as in the definition of translates into a variation in the BAO amplitude.
Several approaches have been considered in the literature to define . The simplest choice is to assume a smooth version of the linear power spectrum. For instance Eisenstein & Hu (1998) provides a formula for a wigglefree power spectrum. This is an unphysical one that nonetheless accounts for the broadband slope of the linear matter power spectrum. However, its use is limited by the fact that the nonlinear regime increasingly boosts the amplitude of the spectrum at larger , thus altering the broadband slope of the spectrum with respect to the linear prediction in a way that erases the BAO signature.
Alternatively, one can define a smooth function directly from observations. For instance, Percival et al. (2007) use a ninenode cubic spline to fit the broadband slope of the measured power spectrum, having excluded BAO data points from the fit. Seo & Eisenstein (2005) use a similar approach, but instead of a cubic interpolation they assume polynomial fitting functions. Such methods, although free of cosmological assumptions, may introduce spurious effects in the analysis of the BAO. This is because the freedom in the choice of the fitting interval or the specific form of the fitting functions can potentially absorb part of the BAO signal. As we will argue this is inevitable to happen for measurements that aims accuracy due to a subtle coupling between the BAO and the broadband slope of the spectrum.
Given these limitations one may be tempted to simply use the linear Cold Dark Matter (CDM) power spectrum rather than the matter one (CDM + baryons). However, as shown by Eisenstein & Hu (1998), this has a different broadband slope which makes harder to highlight the BAO pattern.
In order to better account for the broadband slope Crocce & Scoccimarro (2008) define a smooth function by nonlinearly evolving an initial wigglefree spectrum using Renormalized Perturbation Theory. Here, we follow their approach and compute a smooth nonlinear matter power spectrum from a CDMW7 simulation (see Table 1) with particles and hMpc box (see Table 3) for which initial conditions have been generated using a wigglefree spectrum from (Eisenstein & Hu, 1998). Upon doing so we have ensured that the effective onedimensional amplitude of the largescale velocity flow is the same as that inferred from the linear power spectrum computed with CAMB. The spectrum obtained from this simulation is then fitted with a polynomial of order 8 which defines our smooth power spectrum, (see Appendix A). Here, the fact that we use a fit does not affect the determination of the BAO since we are directly fitting a smooth function. This is also the reason as to why can be inferred without the need of running an extremely large volume simulation.
The wiggleonly spectrum of BAO is then given by
(2) 
where is the matter power spectrum. Such definition solely depends on the choice of the initial smooth power spectrum that we have used to generate the initial conditions of the wigglefree Nbody simulation. Nevertheless, the conclusions depend only weakly on such choice since the initial wigglefree spectrum given by Eisenstein & Hu (1998) formula interpolates through the peaks and dips of the BAO at the initial redshift. In principle, one could opt to directly subtract the smooth spectrum computed from the Nbody simulation rather than using its polynomial fit. This has the advantage of cancelling out statistical fluctuations at large sales. However, this may also introduce additional uncertainties at smaller scales due to the fact that nonlinear effects cause mode couplings and thus the cancellation can lead to spurious effects.
We plot in Figure 4 the relative difference of the BAO spectrum with respect to the wigglefree DEUSFUR spectrum at for DEUSFUR CDMW7 (red solid line) and in the case of a simulation with hMpc box length (black solid line). This illustrates the advantage of studying BAO with DEUSFUR, since disposing of cosmic variance limited measurements allows us to finely resolve the BAO structure. Four peaks and troughs are clearly distinguishable, with the amplitude of the first oscillation (peaktotrough) , while that of the fourth one is of only halfpercent.
4.2 BAO vs SemiAnalytical Model Predictions
It is beyond the scope of this work to test against DEUSFUR all existing semianalytical models that aim to predict the nonlinear corrections to the linear matter power spectrum. Instead, we focus on two such models that are exemplary of substantially different approaches, one based on a perturbative calculation and the other on a combination of fit to simulations and the halo model. The former can provide accurate predictions on quasilinear scales, while failing deep in the nonlinear regime (see e.g. Carlson et al., 2009). The latter, while more accurate at small scales, is less at intermediate ones (0.10.3 h Mpc). In the first case we consider the two loop regularized multipoint propagator method (RegPT) for which we compute the nonlinear matter power spectrum using the code RegPT (Taruya et al., 2012), while in the other case we consider the widely used fitting prescription Halofit (Smith et al., 2003).
The results of the comparison are summarized in Figure 5. In the top panels we plot the relative difference of the matter power spectrum with respect to the DEUSFUR smooth spectrum at (left panel) and (right panel) respectively. The different lines correspond to the linear prediction (black solid line), RegPT (blue solid line), Halofit (green solid line), DEUSFUR (red solid line) and DEUSFUR fit (orange solid line, see Appendix A). First, it is worth noticing that nonlinear modifications of the matter power spectrum occur over the entire BAO range. At (topleft panel) the discrepancy between the linear theory and DEUSFUR is of on the first trough and on the first peak, while at (top right panel) this reduces to and respectively, but still above cosmic variance errors. We may also notice that RegPT and Halofit are in agreement with DEUSFUR to within on the first trough at , while larger deviations occur at smaller scales. In particular, Halofit deviates at level on the first peak and underestimates the damping of the BAO up to beyond the second peak, while it recovers to good approximation the position of the extrema. RegPT shows deviation of order up to the second peak, but afterward it rapidly diverges. At the nonlinear damping of the BAO is smaller, consequently the fourth peak has become more prominent, while even the fifth trough has appeared. Deviations from the linear prediction are now shifted to larger , with the linear amplitude being off by beyond the third trough. The maximal deviation of Halofit and RegPT is , with RegPT diverging beyond the fourth peak.
The discrepancy of the semianalytic models with respect to DEUSFUR is largely due to the broadband slope of the nonlinear power spectrum. In the bottom panels of Figure 5 we plot the BAO spectrum obtained by subtracting the smooth spectrum predicted by each model. In this case we also plot the prediction from the leading term of Renormalized Perturbation Theory (lRPT) given by
(3) 
with the linear power spectrum, the smooth (wigglefree) linear power spectrum and
(4) 
is the effective onedimensional amplitude of largescale velocity flows (Crocce & Scoccimarro, 2008). This expression coincides with that from the peakbackground split to lower order as well as the resummed Lagrangian Perturbation Theory, as such it has motivated several fitting formulae of the BAO.
We notice that Halofit now has deviations of order of on dips and on peaks at h Mpc both at and respectively.
RegPT shows a remarkable agreement with DEUSFUR with differences up to the third trough, which is well beyond the supposed range of validity of RegPT. Slightly larger deviations occurs only beyond the third trough at . The case of the leading RPT term (light blue line) also provides a good description to the DEUSFUR spectrum only a few percent worse than RegPT.
Next, we will perform a detailed quantitative analysis of the nonlinear effects on the BAO extrema.
5 Featuring BAO pattern
The nonlinear clustering of matter shifts the position of the BAO peaks and dips, damps their amplitude and alters the broadband slope of the power spectrum. These are not independent effects, since the damping and the broadband slope both contribute to shifting the position of the BAO extrema relative to the linear case. Nevertheless such decomposition provides a phenomenologically meaningful way of quantifying the effect of nonlinearities, testing the accuracy of semianalytical model predictions or comparing them against observational measurements (to this purpose we provide in Appendix A a fitting formula of the DEUSFUR spectrum).
Nonlinearities are usually considered to be a nuisance which degrades the cosmic distance information encoded in the BAO. We will show that an accurate modelling of these nonlinear features can provide additional cosmological information as they probe the linear growth rate of cosmic structures.
The nonlinear shift of the position of the BAO extrema has been quantified using perturbation theory in Nishimichi et al. (2007). Here, we perform a detailed numerical analysis using the BAO spectrum from DEUSFUR. Our aim is to evaluate in absolute terms the effect of nonlinearities on the characteristics of BAO extrema, hence differently from the previous Section we normalize the BAO spectrum to the amplitude of the DEUSFUR spectrum at on a linear scale, , such as not to alter the extrema.
As already stressed, differently from previous studies (see e.g. Seo & Eisenstein, 2005; Angulo et al., 2008; Seo et al., 2008, 2010), we will characterize each BAO extremum directly from the measured Nbody spectrum, rather than a bestfitting function. We detect the BAO peaks and dips performing a local secondorder leastsquare polynomial fit (SavitskyGolay filter) of width 100 data point which takes into account DEUSFUR error bars. In Figure 6 we illustrate our decomposition of the BAO extrema in terms of the horizontal position, peaktotrough amplitude and peaktotrough average.
In order to consistently compare with the semianalytical model predictions we have binned the corresponding power spectra as the DEUSFUR case and applied the same procedure to estimate the location and amplitude of the BAO extrema.
In Figure 7 we plot a zoom on the first four troughs (left panels) and peaks (right panels) from DEUSFUR (red line), linear theory (black line), RegPT (blue line) and Halofit (green line) respectively. Crosses mark the location of the detected extrema. In Figure 8 we plot the same zoom on the extrema having subtracted the smooth power spectrum predicted by each model.
5.1 Shift of BAO extrema
The physical origin of this shift is discussed in several papers (Eisenstein et al., 2007a; Nishimichi et al., 2007; Crocce & Scoccimarro, 2008; Padmanabhan & White, 2009; Sherwin & Zaldarriaga, 2012). An accurate determination of the nonlinear shift of BAO extrema is crucial to infer unbiased cosmic distance measurements. As shown in (Angulo et al., 2008) an error of on the peak position at leads to a error on a Dark Energy equation of state parameter (assuming the other cosmological parameters to be known).
In the topleft panel of Figure 9 we plot the shift of the BAO peaks (even points) and dips (odd points) at defined as the ratio between the location of the extrema from the linear theory (black line), RegPT (blue line), Halofit (green line) relative to DEUSFUR (red line). We can see that the linear theory misestimates the location of the extrema at beyond the first dip. In contrast, RegPT deviates at more than only beyond the third dip as consequence of the drop of the broadband slope. Halofit performs better with deviations up to the fourth dip.
Subtracting the wigglefree spectrum predicted by each of the models greatly reduces the discrepancies. This can be seen in the topright panel of Figure 9. In the case of RegPT the deviations are well within up to the fourth dip, similarly for Halofit and the leading RPT term. Even the linear theory is within up to the third peak.
We now focus on the redshift dependence of the shift of the position of each BAO extrema in DEUSFUR with respect to that in the linearly evolved initial power spectrum. This is plotted for different extrema in the left panel of Figure 10. We can see that the amplitude of the shift is at all redshifts. It is worth noticing that the second dip and the third peak remains unaltered by the nonlinear effects with shifts which are . Given the fact that the latest BAO measurements already test the shift of the BAO at a few percent level (Anderson et al., 2012; Padmanabhan et al., 2012; Xu et al., 2002), this result suggests that upcoming cosmic distance measurements from BAO may derive unbiased estimates by specifically focusing on these two extrema. We also remark that the redshift dependence is well described by the square of the linear growth function (dashed line) as expected from perturbation theory (Crocce & Scoccimarro, 2008; Padmanabhan & White, 2009; Seo et al., 2010). Hence, future surveys may have enough resolution to detect such a trend and independently probe the linear growth rate.
5.2 Damping of BAO amplitude
The suppression of BAO is the most visible nonlinear effect which is direct consequence of the relative displacement (510 hMpc) of pair of particles separated by hMpc (Eisenstein et al., 2007a). Because of this, a measurement of the BAO damping factor can in principle probe the lagrangian displacement field and provide constraints on the growth rate of cosmic structures (see e.g. Nomura, Yamamoto & Nishimichi, 2008; Nomura et al., 2009).
We estimate the damping of the BAO by measuring the absolute difference of the amplitude of two consecutive extrema (i.e. first trough to first peak, first peak to second trough,…). We plot the predicted amplitudes for the different models relative to DEUSFUR case in the middleleft panel of Figure 9. We can see that both the linear theory and Halofit overestimate the amplitude of the extrema with deviations exceeding at h Mpc. In the right panel we show the same amplitudes estimated after subtracting the wigglefree spectrum predicted by each model. In such a case we can clearly see that the linear calculation and Halofit exponentially overestimate the amplitude of the BAO extrema. In the case of RegPT the discrepancy amounts to up to h Mpc and diverges afterward, whereas it reduces to no more than up to h Mpc after subtracting the smooth component and increases up to on the extrema at the highest .
We now estimate the damping factor in DEUSFUR BAO spectrum relative to the linear prediction. This is defined as
(5) 
where and are the peak to trough amplitude in DEUSFUR and the linear theory (initial condition) respectively, and are the wavenumbers of the peak and trough respectively. The above formula converges to first order to the case of an exponential term with constant damping rate .
In the right panel of Figure 10 we plot as function of redshift for the BAO extrema in DEUSFUR. First, we can see that the damping does not behave as a constant term , rather it increases as function of from one extrema to the next. Secondly, it decreases as function of redshift as expected from perturbation theory. However, the exact value for each pair of extrema is not given by as expected from Eq. (3). Because of these features a global fit to the BAO would lead to a different estimation of the damping factor depending on the choice of the fitting range or the relative size of experimental errors. On the other hand, we can see that the redshift evolution of the damping rate approximately scales as (dashed line). Thus, it should be possible to infer the linear growth rate from the measurement of BAO amplitude at two different redshifts.
5.3 BAO coupling to broadband slope
The inaccuracy in reproducing the broadband slope of the nonlinear matter power spectrum is primarily responsible for the discrepancies between the semianalytic model predictions and the DEUSFUR spectrum. That is why subtracting the smooth wigglefree spectrum predicted by each of the models provides a much better description. We quantify this effect by computing the average amplitude of neighbouring pairs of extrema. We plot the prediction for different models at relative to DEUSFUR in the bottomleft panel of Figure 9. We can see that the linear theory overestimates the average amplitude up to at h Mpc, while it underestimates at larger . Halofit and RegPT are in agreement with DEUSFUR at level up to h Mpc. At larger wavenumbers Halofit remains within , while the prediction from RegPT sharply falls to more than above h Mpc. Subtracting the wigglefree spectrum of each model absorbs most of the broadband slope effect as shown in the bottomright panel of Figure 9. Discrepancies are now , however it is worth noticing a mode dependent trends which differ from one model to another. This indicates the presence of a residual coupling between the BAO and the continuum which complicate the modelling of the nonlinear effects on the BAO. That is why defining a wigglefree spectrum through a fit to the observed spectrum may introduce spurious effects at the percent level.
6 Conclusion
We have presented a study of the BAO in the realspace matter power spectrum from the DEUSFUR simulation of CDMW7 cosmology. The simulation box covers the size of the full observable volume of the simulated cosmology thus enabling a cosmic variance limited estimation of the matter power spectrum. Using a testbed of large volume highresolution simulations we are able to control numerical systematic errors to within . This allows us to investigate for the first time the full BAO range () at low redshift (from to ) with less than systematic and statistical errors.
Previous works dedicated to assessing the effect of nonlinearities on the BAO have either relied on perturbation theory (which eventually breaks down at large wavenumbers) or used Nbody simulations of smaller effective volumes without assessing the effect of numerical systematic errors. Furthermore, differently from previous studies, the very narrow sampling in (with bins of size ) accessible through DEUSFUR has allowed us to determine the characteristics of the BAO extrema by directly measuring of the peaks and dips in the BAO spectrum. In this case the choice of a wigglefree spectrum necessary to define the BAO is of particularly importance. Here, rather than using fitting formulae that may introduce spurious effects and absorb part of the BAO signal, we have computed a smooth spectrum from a Nbody simulation with an initial wigglefree spectrum.
Such characteristics make the BAO from DEUSFUR BAO an optimal reference for testing semianalytical model predictions and calibrate the amplitude of nonlinear effects expected in a standard CDM cosmology. To this purpose we provide in Appendix A polynomial fitting formulae of the BAO spectrum valid in the redshift interval , based on a decomposition of the nonlinear effects in term of a shift of the BAO extrema, a damping factor and a broadband slope.
We have tested two widely used models, Halofit (Smith et al., 2003) and RegPT (Taruya et al., 2012). We find that Halofit overestimates the amplitude of the BAO by a large factor (, and for the three first peaks respectively). In contrast RegPT accurately reproduce the BAO spectrum up to h Mpc beyond which discrepancies with respect to DEUSFUR diverge. This is mostly due to inaccurately predicting the broadband slope of the nonlinear power spectrum. Thus, subtracting a smooth spectrum as predicted using RegPT leads to much smaller deviations.
We have detected a small nontrivial coupling between BAO and the broadband slope. Furthermore, we have quantified the shift and damping of each BAO extrema and determined their evolution as function of redshift. We find that at a given redshift the second and third trough are not affected by nonlinearities and thus can provide unbiased measurements of the cosmic distances. The redshift dependence of the shift is proportional to , hence measurements of the shift of the BAO extrema at different redshifts may test the linear growth factor. However, given the relatively small amplitude of the shift () such measurements may turn to be very challenging even with future surveys. In contrast, the damping of the BAO amplitude has much larger impact and varies with redshift proportionally to . Hence, by measuring the ratio of peaktotrough amplitude at two different redshifts, it should be possible to constrain the growth factor and infer additional cosmological constraints besides those encoded in the cosmic distance to the BAO.
The nonlinear effects on the BAO spectrum vary with the underlying cosmological model. In a future work we will investigate such dependency extending this analysis to the two other DEUSFUR simulations of nonstandard Dark Energy models.
Acknowledgements
We are particularly grateful to the team of engineers at TGCC for technical support. We are very thankful to Romain Teyssier (RAMSES code), Simon Prunet & Christophe Pichon (MPGRAFIC and POWERGRID) for providing their original applications. We are also thankful to Ravi Sheth, Martin Crocce, Francisco Prada, ChiaHsun Chuang, Atsushi Taruya and Fabrice Roy for fruitful discussions. This work was granted access to HPC resources of TGCC, CCRT and IDRIS through allocations made by GENCI (Grand Équipement National de Calcul Intensif) in the framework of the “Grand Challenges” DEUS and DEUSS. The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007 2013 Grant Agreement no. 279954). We acknowledge support from the DIM ACAV of the Region Ile de France.
Appendix A Fitting formula for low redshift Baryonic Acoustic Oscillations
Here, we provide a polynomial fitting formula to the DEUSFUR CDMW7 BAO spectrum valid in the range to h Mpc in the redshift interval . Our template is inspired by that used in Seo et al. (2010) and formulated such as to preserve our decomposition of nonlinear effects in terms of shift, damping and continuum. When using functional fitting formula the physical interpretation of this quantities may depend on the exact form of the assumed template, hence the comparison with equivalent quantities derived from other fitting functions should be carefully performed. For this reason, in our analysis of the BAO we have preferred to measure the shift, damping and the continuum directly on the BAO spectrum. Nevertheless a fitting formula is still useful for model testing or comparison with noisy data. An implementation in various languages (fortran 90, C, IDL, etc.) can be found at http://www.deusconsortium.org/overview/results/bao/.
Our template reads as
(6)  
with the shift, the continuum and the damping given by
(7)  
(8)  
(9) 
and the redshift evolution of the coefficient is given as function of the scale factor in matrix form
(10)  
Statistical errors are given by the number of independent modes in a given bin as discussed in the text, these are added in quadrature to a systematic error which we have estimated to be .
The quality of this fitting formula is shown as orange line in the BAO plots and reproduces well the amplitude, shift and average of the DEUSFUR BAO extrema.
For sake of completeness we also provide here the polynomial fitting formula of at :
We have implemented an extension of this fit to higher redshift assuming the same form as above. A code for this fitting function is also available at:http://www.deusconsortium.org/overview/results/bao/.
References
 Alimi, J.M., Bouillot, V., Rasera, Y., et al., 2012, IEEE Computer Soc. Press, CA, USA, SC2012, Article No 73.
 Amendola, L., et al., 2012, arXiv:1206.1225
 Anderson, L. et al., 2012, MNRAS, 427, 3435
 Angulo, R. E., Baugh, C. M., Frenk, C. S., & Lacey, C. G., 2008, MNRAS, 383, 755
 Bagla, J. S., & Ray, S., 2005, MNRAS, 358, 1076
 Blake, C., Glazebrook, K, 2003, ApJ, 594, 665
 Carlson, J., White, M., & Padmanabhan, N., 2009, Phys.Rev.D, 80, 043531
 Cole, S. et al., MNRAS, 2005, 362, 505
 Crocce, M., Pueblas, S., & Scoccimarro, R., 2006, MNRAS, 373, 369
 Crocce, M., & Scoccimarro, R., 2008, Phys.Rev.D, 77, 023533
 Dawson, K.,S., et al., 2013, AJ, 145, 10
 Eisenstein, D. J., & Hu, W., 1998, ApJ, 496, 605
 Eisenstein, D. J., Hu, W., Tegmark, M., 1998, ApJ, 504, L57
 Eisenstein, D. J., et al., 2005, ApJ, 633, 560
 Eisenstein, D. J., Seo, H.J., & White, M., 2007a, ApJ, 664, 660
 Eisenstein, D. J., Seo, H.J., Sirko, E., Spergel, D.N., 2007b, ApJ, 664, 675
 Guillet, T., & Teyssier, R. 2011, Journal of Computational Physics, 230, 4756
 Heitmann, K., Ricker, P. M., Warren, M. S., & Habib, S., 2005, ApJs, 160, 28
 Heitmann, K., et al. 2008, Computational Science and Discovery, 1, 015003
 Heitmann, K., White, M., Wagner, C., Habib, S., & Higdon, D., 2010, ApJ, 715, 104
 Huetsi, G., A&A, 2006, 449, 891
 Ivezic, Z., et al., 2008, arXiv:0805.2366
 Jeong, D., & Komatsu, E., 2009, Apj, 691, 569
 Joyce, M., Marcos, B., & Baertschiger, T., 2009, MNRAS, 394, 751
 Knebe, A., Green, A., & Binney, J., 2001, MNRAS, 325, 845
 Komatsu, E., et al., 2009, ApJs, 180, 330
 Lewis, A., Challinor, A., & Lasenby, A., 2000, ApJ, 538, 473
 Nomura, H., Yamamoto, K., Nishimichi, T., 2008, JCAP, 10, 031
 Nomura, H., Yamamoto, K., Hutsi, Gert, Nishimichi, T., 2009, PRD, 79, 063512
 Nishimichi, T., Ohmuro, H., Nakamichi, M., et al., 2007, Publications of the Astronomical Society of Japan, 59, 1049
 Nishimichi, T., Shirata, A., Taruya, A., et al., 2009, Publications of the Astronomical Society of Japan, 61, 321
 Ngan, W., HarnoisDéraps, J., Pen, U.L., McDonald, P., & MacDonald, I., 2012, MNRAS, 419, 2949
 Orban, C., & Weinberg, D. H., 2011, Phys. Rev. D, 84, 063501
 O’Shea, B. W., Nagamine, K., Springel, V., Hernquist, L., & Norman, M. L., 2005, ApJs, 160, 1
 Padmanabhan, N., & White, M., 2009, Phys. Rev.D, 80, 063508
 Padmanabhan, N., et al., 2012, MNRAS, 427, 2132
 Peebles, P. J. E., Yu, I.T., 1970, ApJ, 162, 815
 Percival, W. J., Nichol, R. C., Eisenstein, D. J., et al., 2007, ApJ, 657, 51
 Planck Collaboration, 2013, arXiv:1303.5076
 Prunet, S., Pichon, C., Aubert, D., Pogosyan, D., Teyssier, R., & Gottloeber, S., 2008, ApJs, 178, 179
 Reed, D. S., Smith, R. E., Potter, D., et al., 2013, MNRAS, 431, 1866
 Sakharov, A.D., 1965, JETP, 49, 345
 Scoccimarro, R., 1998, MNRAS, 299, 1097
 Seo, H.J., & Eisenstein, D. J., 2003, ApJ, 598, 720
 Seo, H.J., & Eisenstein, D. J., 2005, ApJ, 633, 575
 Seo, H.J., Siegel, E. R., Eisenstein, D. J., & White, M., 2008, ApJ, 686, 13
 Seo, H.J., Eckel, J., Eisenstein, D. J., et al., 2010, ApJ, 720, 1650
 Sherwin, B. D., & Zaldarriaga, M., 2012, Phys.Rev.D, 85, 103523
 Silk, J., 1968, ApJ, 151, 459
 Smith, R. E., et al., 2003, MNRAS, 341, 1311
 Smith, R. E., Reed, D. S., Potter, D., et al., 2012, arXiv:1211.6434
 Spergel, D. N., et al., 2003, ApJs, 148, 175
 Spergel, D. N., et al., 2007, ApJs, 170, 377
 Takahashi, R., Yoshida, N., Matsubara, T., et al., 2008, MNRAS, 389, 1675
 Takahashi, R., Yoshida, N., Takada, M., et al., 2009, ApJ, 700, 479
 Taruya, A., Bernardeau, F., Nishimichi, T., & Codis, S., 2012, Phys.Rev.D, 86, 103528
 Teyssier, R., 2002, A&A, 385, 337
 Xu, X., et al., 2012, MNRAS, 427, 2146