High accuracy power spectra including baryonic physics in dynamical Dark Energy models

The next generation mass probes will obtain information on non--linear power spectra P(k,z) and their evolution, allowing us to investigate the nature of Dark Energy. To exploit such data we need high precision simulations, extending at least up to scales of k 10 h/Mpc, where the effects of baryons can no longer be neglected. In this paper, we present a series of large scale hydrodynamical simulations for LCDM and dynamical Dark Energy (dDE) models, in which the equation of state parameter is z-dependent. The simulations include gas cooling, star formation and Supernovae feedback. They closely approximate the observed star formation rate and the observationally derived star/Dark Matter mass ratio in collapsed systems. Baryon dynamics cause spectral shifts exceeding 1% at k>2-3 h/Mpc compared to pure n-body simulations in the LCDM simulations. This agrees with previous studies, although we find a smaller effect (~50%) on the power spectrum amplitude at higher k's. dDE exhibits similar behavior, even though the dDE simulations produce ~20% less stars than the analogous LCDM cosmologies. Finally, we show that the technique introduced in Casarini et al. to obtain spectra for any $w(z)$ cosmology from constant-w models at any redshift still holds when gas physics is taken into account. While this relieves the need to explore the entire functional space of dark energy state equations, we illustrate a severe risk that future data analysis could lead to misinterpretation of the DE state equation.


INTRODUCTION
There can be little doubt that Dark Energy (DE) is a necessary ingredient for any cosmological model approaching data (see, e.g., the data fits in WMAP5 release, Komatsu et al. 2009). In turn, the quest for the nature of DE is perhaps the main challenge in contemporary physics in spite of the fact that no available data conflict with ΛCDM models, in which the DE state parameter is w ≡ −1, as is true for a false vacuum. A false vacuum interpretation of DE, however, leads to an unnatural value for its density (fine tuning problem) and to a number of coincidences characterizing our epoch and the primeval fluctuation amplitude (coincidence problem).
Many efforts have therefore been devoted to finding vi-⋆ E-mail: casarini@mib.infn.it able alternatives. Among them, the possibility that DE is a scalar field, φ, self-interacting through a potential V (φ) (Wetterich 1989, Ratra & Peebles 1989, has been widely explored. In this case, both the ratio between kinetic and potential energies,φ 2 /2V (φ) and the DE state parameter w(z), depend on redshift, z. In this paradigm, observational data on w(z) would yield V (φ). We will refer to this class of models, as dynamical DE (dDE).
Future large tomographic shear surveys, such as the Euclid project (see, e.g., Refrégier et al. 2010), will measure matter density fluctuations and their evolution from the linear to the non-linear regime with unprecedented precision, approaching ≈ 1% (Huterer & Takada 2005). To translate these data into information on DE, we need to predict power spectra and the evolution of power spectra with similar precision for a wide set of w(z) laws. The resulting constraints on w(z) will be a substantial step toward understanding DE physics (see, e.g., Hu & Tegmark 1999, Dodelson & Zhang 2005, Manera & Mota 2005, Mota 2008, La Vacca et al. 2008, 2009, Kristiansen et al. 2010.
The evolution of the baryonic component of the Universe is highly uncertain and can affect predictions for lensing observables, manifesting as modified structure growth for a fixed cosmic distance scale , Hearin & Zentner 2009.
Collisionless N-body simulations have been extensively used to obtain non-linear power spectra (e.g. Smith et al. 2003). These simulations have ignored the effect of baryons, which comprise ∼ 15 % of the matter content of the observed Universe, (Komatsu et al. 2009). Although baryons make up a minor fraction of the Universal matter content and is distributed exactly the same as Cold Dark Matter (CDM) at the onset of the non-linear growth, the final distribution of baryons in haloes is significantly different than CDM, because of its separate later dynamical evolution. Jing et al. (2006) have shown that, in the spectra of ΛCDM models, the difference between N-body simulations and simulations including baryon physics exceeds 1% at k > 2-3 hMpc −1 , and increases to ∼ 10-20 % on smaller scales when k approaches 10 hMpc −1 . Rudd et al. (2008) and Guillet et al. (2010) recently confirmed the need for simulations including baryonic physics, even though their results differed from Jing et al. (2006).
The discrepant results originate from the different prescriptions used for star formation, SN feedback, AGN effects, and possibly, the numerical approach. To compare different w(z) laws, these uncertainties need to be under control. Using the same prescriptions for all models is essential to compare their spectra, as we do in this work. Their being tuned to observations is also a need, when aiming to fit future data. The prescription used here are tuned on available observations, such tuning being however carried on within the frame of ΛCDM cosmologies. This is to be borne in mind, in view of model dependent effects that we shall detect, e.g. for what concerns star production in haloes.
Dynamical Dark Energy (dDE) simulations, with a variable state parameter w(a) deduced from scalar field potentials admitting a tracker solution, have been performed since 2003 (e.g. Klypin et al. 2003, Linder & Jenkins 2003, Dolag et al. 2004, Macciò et al. 2004, Solevi et al. 2005. These have been compared with constant-w simulation outputs. Observables considered in these papers, only marginally included spectra. Recently Francis et al. (2007) have shown how spectral predictions for constant-w models at z = 0 can also be used to fit spectra of cosmologies with a state parameter given by a first degree polynomial: More precisely, the spectrum P (k, 0), in an assigned (A) model with state parameter w(z), can be approached by a suitable auxiliary constant-w model (W) if: (i) in A and W, Ω b,m,tot , h and σ8 are equal, (ii) the constant DE state parameter w of W is tuned so that W and A have the same comoving distance from the Last Scattering Band (LSB) and z = 0. In that case, spectral discrepancies stay < 1 %, up to k ∼ 2-3 h Mpc −1 (symbols here above have their usual meaning). Francis et al. (2007) also tried extending their results to higher redshifts, but discrepancies between A and W in-creased to several percent; not enough to exploit forthcoming data. The required precision was, however, recovered through a new technique introduced by Casarini, Macciò & Bonometto (2009, Paper I hereafter) and tested with Nbody simulations.
The final aim of this work is to extend the techniques use in Paper I to scales that only hydrodynamical simulation can test. This is a serious challenge for the description of physical processes included in hydrodynamical algorithms. We will devote a part of our work to understanding these challenges.
The focus of this paper is twofold: (i) Study the importance of the effects of baryons on the matter power spectrum on the small scales that will probed by future surveys, bearing in mind that the physical assumptions will have to be carefully calibrated to approach the unprecedented precision requirements. (ii) Extend the commonalities found in N -body simulations in Paper I among spectra of different models down to scales where baryon physics cannot be disregarded. Accordingly, we shall use the same baryon physics in a series of simulations to compare the spectra of different dDE models.
The paper is therefore organized as follows: In §2, we describe the approach used in Paper I. In §3, we motivate the choice of the models considered. §4 is devoted to describing our simulations and the techniques used to analyse them. In §5, we present our hydrodynamical simulations. §6 is subdivided in subsections that concern various aspects of the power spectra analysis. In §6 we test: (i) the relation between N-body and hydrodynamical spectra (ii) the effect of box size and numerical resolution (iii) the level at which model parameters can be discriminated; and, finally, (iv) the success of the extension of Paper I approach to the higher k range. A general discussion of the results closes our paper.

APPROACHING DYNAMICAL DE SPECTRA
Here is a short summary of the technique presented in Paper I, which improved the results of Francis et al. (2007) at z > 0. Given an assigned model A, we introduce an auxiliary model W(z), for each z. A and W(z) are required to share the values of both ω b,c,m = Ω b,c,m h 2 and σ8, at the selected z .
The first condition is satisfied at any redshift using the relationship between H and critical density ρcr where H is the Hubble parameter. Multiplying both sides of this equation by Ωm (or Ω b , Ωc) produces The r.h.s. of this equation scales as a −3 , independent of model, since ωm ∝ ΩmH 2 (or ω b , ωc). Therefore, once the A and W models share ω b,c,m at z = 0, they also share ω b,c,m at every other z. Accordingly, each W(z) model shares the same ω b,c,m values. Conversely, the evolution of σ8 depends on the state equation of DE and its value at z = 0 and can be obtained only from the constant DE state parameters w(z), which is likelihood ellipses on the wo-w ′ plane and the position of our models on that plane. The significant shift from WMAP5 to WMAP7 is not due to new CMB data, but mostly from the improved HST determination of H 0 . Dotted lines yield the best likelihood model. The selected model lies near the 1-σ contour for WMAP7, and within 2-σ's, for WMAP5 as well. Notice how much closer the W models (empty polygons) are than the assumed cosmology (filled triangle). This highlights a serious observational risk (see §7).
fixed by requiring that W(z) and A have equal distances between z and the LSB.
The choice of H0 (the Hubble parameter at z = 0) is unconstrained, so we set it equal to H0 in A to allow the simulated volumes to have the same side L in both Mpc and h −1 Mpc units. Notice that a simple-minded generalization of Francis et al. (2007) criterion to high z would require equal Ω b,c,m (z) and, thence, H(z), and creates serious problems of sample variance and model comparison, at most redshift values.
In Paper I, the mapping between W(z) and A spectra was tested, using N-body simulations, up to k ≃ 2-3 h Mpc −1 . The discrepancies were mostly in the per-mil range. The techniques used in Paper I worked better for at higher z values.

COSMOLOGICAL MODELS
To extend our results to hydrodynamical simulations, we use three cosmological models, each of which is consistent with both WMAP5 and WMAP7 data. Two models share most cosmological parameters: baryon and CDM density Ω b = 0.046, Ωc = 0.228, Hubble parameter (in units of 100 km/s/Mpc) h = 0.7, r.m.s matter fluctuation amplitude at z = 0 on 8 h −1 Mpc scales, σ8 = 0.81, and primeval spectral index n = 0.96. The models differ in the assumed dark energy equation of state, either w ≡ −1 for ΛCDM or using equation (1) with wo = −0.8, w ′ = −0.754 for dDE (also sometimes referred to as the A model). We compare the dDE with various auxiliary constant-w models W(z). We plan to extend the tests performed on dDE to other w(z) expressions in future work, but there is currently no reason to believe that our dDE model is peculiar. In Figure 1, the A model and the auxiliary models are set on the WMAP likelihood ellipses on the wo-w ′ plane.
In order to extend the model comparison, we also considered another ΛCDM model (ΛCDM1), with a greater CDM density Ωc = 0.254, which yields a matter density parameter Ωm = 0.3. Thus, ΛCDM1 is ∼ 2 σ's from ΛCDM using the WMAP5 data. ΛCDM is run in boxes of different sizes to test the effectiveness of the mass and spatial resolution used. A summary of models and boxes used in this paper is presented in Table 1. ΛCDM coincides with W(z = 0).

SIMULATIONS AND THEIR ANALYSIS
We run both N-body and hydrodynamical simulations (here below also dubbed DMO and DMG, representing "DMonly" and "DM & gas" respectively). For DMO simulations we use pkdgrav (Stadel 2001) that has been modified to enable variable w(z) as described in Paper I. For the hydrodynamical simulations, we use gasoline, a multi-stepping, parallel Treesmoothed particle hydrodynamic (SPH) Nbody code (Wadsley et al. 2004).
gasoline includes radiative and Compton cooling for a primordial mixture of hydrogen and helium. The star formation algorithm is based on a Jeans instability criteria (Katz 1992), but simplified so that gas particles satisfying constant density and temperature thresholds in convergent flows spawn star particles at a rate proportional to the local dynamical time (see Stinson et al. 2006). The star formation efficiency was set to 0.05 based on simulations of the Milky Way that satisfied the Kennicutt (1998) Schmidt Law. The code also includes supernova feedback as described by Stinson et al. (2006), and a UV background following Haardt & Madau (1996); see Governato et al. (2007) for a more detailed description of the code. We applied the same mod- ifications to gasoline that we implemented in pkdgrav in Paper I to handle the dDE models.
All of the simulations use 256 3 CDM particles (and 256 3 gas particles) and most of them are boxes 256 h −1 Mpc on a side. For the ΛCDM, dDE models, and the related auxiliary models, the CDM particles each have a mass of mch/M⊙ = 7.61 × 10 10 in the DMO simulations. In the DMG simulations, the CDM particles have a mass of mch/M⊙ = 6.33 × 10 10 and the gas particles have a mass of m b h/M⊙ = 1.28 × 10 10 . For the ΛCDM1 model, mch/M⊙ = 8.33 × 10 10 in the DMO and mch/M⊙ = 7.05 × 10 10 , m b h/M⊙ = 1.28 × 10 10 in the DMG.
The force resolution (softening) is 1/40 of the intraparticle separation. For L = 256 h −1 Mpc, this corresponds to a distance ǫ ≃ 25 h −1 kpc (wavenumber κ = 2π/ǫ ≃ 150 h Mpc −1 ). This softening provides spectral resolution up to k = 10hMpc −1 . Each simulation is started at zin = 24 using the same random realization of the density field. In order to have equal σ8 at the various redshift zi, where model spectra are compared, different models are started with (slightly) different σR at zin. (Here σR is the m.s. fluctuation over the scale R; e.g.: σ8 ≡ σ 8 h −1 Mpc .) The smallest resolved scale Rr yields the greatest relevant σR r . The selected zin is such that residual non-linear effects cannot induce discrepancies among models, over the scale Rr, exceeding 1:10 4 (Crocce et al. 2006). For the ΛCDM model (DMO and DMG) we also ran a smaller box L = 64 h −1 Mpc (ΛCDM-small), to examine the resolution effects.
One of the key analyses in this work are power spectra of the matter density field. The matter density field is found by interpolating the particle distribution onto a regular Ng × Ng × Ng grid using the Cloud-in-Cell algorithm. The power spectra are calculated using a FFT (Fast Fourier Transform) of the matter density field. For our analysis we use Ng = 2048 for the large boxes (L = 256h −1 Mpc) and Ng = 512 for ΛCDM -small, in order to cut the spectra at the same frequency.
Another key piece of analysis was finding haloes. We identified collapsed structures using the spherical overdensity (SO) algorithm. We used a time varying virial density contrast from the fitting formulae of Mainini et al. (2003). The halo catalogue includes all structures with more than 200 DM particles. For DMG simulations we estimated the halo masses by taking into account all particles inside the virial radius (see Macciò et al. 2008 for further details on our halo finding algorithm).

MASS FUNCTIONS AND STAR FORMATION
A significant test of our simulations are halo mass functions. Figure 2 compares the cumulative mass function N (> M ) in our ΛCDM simulations with the Sheth & Tormen (2002) predictions for z = 0 and z = 1. Figure 2 shows the haloes found in the DMO and DMG simulations for both 256 and 64 h −1 Mpc boxes. The consistency between the DMO and DMG simulations is readily apparent, as is the continuity of the mass function over 4 orders of magnitude. We make a more direct comparison with observations in   Figure  3 shows the average Mstar/M halo for each mass bin. These are compared with the predictions from the halo occupation model of Moster et al. (2010). Each simulation point is obtained from equal halo numbers and horizontal error bars yield the m.s. spread of log(M halo ). Even though the simulations are run at low resolution relative to star formation scales, the number of stars formed in the haloes generally follows the halo occupation distribution trend. It is possible that too many stars form in haloes in the highest M halo bin. This is likely due to overcooling that plagues hydrodynamic simulations and might be fixed by including AGN effects in the simulations.
Most important to the current study, Figure 3 shows a significant reduction of Mstar/M halo from the ΛCDM simulations to those run with dDE. Such reduction is significant and one should not be misled by the apparent consistency within error bars: they are large because of the spread of Mstar/M halo through individual haloes, but the reduction of Mstar/M halo is not randomly distributed, exhibiting just a slight increasing trend from smaller to greater haloes. We shall return on this point below.
Let us however remind that star formation in hydrodynamical simulations strongly depends on numerical (mainly mass) resolution (e.g. Mayer et al. 2008 and references therein). We therefore expect the result shown in Figure   3 to change if the same simulations are run at higher resolution (not only because we would resolve lower mass haloes). Naively one could expect to see a higher stellar production and, therefore, a significant deviation of simulated curves from observed ones. The point of our comparison with real data is not to show that our resolution allows intrinsically correct results, but that, with resolution adopted in this specific work, hydro simulations produce a reasonable amount of stars. Figure 4 shows the cosmic star formation history in the 256 h −1 Mpc boxes for both the ΛCDM and dDE simulations. The top panel shows the simulated star formation rate density (SFR,ρ * ) normalized to the observed z = 0 value from Hopkins & Beacon (2006). Due to the low resolution of the simulations, they underproduced stars by a factor of ≈ 1.5. We renormalized these values to show that the shape of the star formation history is similar to observations. The observational results are shown with constant ±20 % error bars that roughly approximate the 2-σ uncertainty of their results. Besides the slight deficiency in star formation, there are two key differences between the simulations and observed cosmic SFH. The Hopkins & Beacon (2006) observations peak at z ≃ 2.3 while the simulations peak around z ≃ 1.5. The onset of star formation is also later in the simulations than the observations. This discrepancy is again a result of low resolution in the 256 h −1 Mpc boxes. Christensen et al. (2010) showed that stars don't form in halos with less than 200 gas particles. Such haloes have a mass of ∼ 1.3 × 10 12 h −1 M⊙ in our simulations. Sufficient numbers of these haloes do not form in our simulations until z = 4 to form stars. The SFH of the 64 h −1 Mpc box more closely follows the observed SFH with a peak around z ≃ 2.5 and an early onset of star formation. The results in this case are also affected by significant sample variance due to the limited size of the simulated volume.
The bottom frame of Figure 4 shows the ratio betweeṅ ρ * the dDE and ΛCDM simulations, highlighting the reduced star formation found using the dDE model. The decrease is ∼ 20% and it is worth noting that it originates simply from a change of the DE equation of state. The models, in fact, have exactly the same age at z = 0 and identical cosmological parameters, apart from w. The ∼ 20 % discrepancy shows how different the evolutionary histories and concentration distributions are between the two cosmologies. The magnitude of the discrepancy is larger than anticipated and we will devote future work to determining the cause of such a significant difference.
Before concluding this section let us however return on the issue concerning resolution. The comparison with real data just shows that our hydro simulations give a reasonable star production rate, as function of halo mass, and a reasonable distribution of such rate, as a function of redshift. This ensures us that the effect of baryonic matter (mainly stars) on the total matter power spectrum should be realistic. This will be true especially on the spatial scales (R) that we will address in this work (for k ≈ 10 → R ≈ 300h −1 kpc), scales that are much larger than the typical sizes of galaxies.

SPECTRA
Now that we have shown that our hydrodynamic model is reasonable on the large scales we wish to study, we turn our attention to how the density power spectra depend on the DE model used. Figure 5 shows a comparison between DMO and DMG spectra at z = 0. Both the DM and total power spectra of the collisional case show a substantial enhancement especially at scales k > 2h Mpc −1 . This scale corresponds to the most massive structure in our simulation and reflects the fact that gas has condensed and cooled inside those structures driving a contraction of the haloes themselves (e.g. Gnedin et al. 2004). The gas shows a bias with respect to the other DMG spectra at large scales k < 8 − 10h Mpc −1 . This effect has been shown in previous work (Jing et al. 2006 and is related to gas moving from large scales to small scales due to cooling. The cooling is most obvious at k > 10h Mpc −1 , where the gas spectra has more power than DMO, reflecting the presence of dense clumps of cold gas at small scales, in the center of dark matter halos.

Hydro vs. N-body simulations
The differences between DMO and DMG can be better appreciated in Figure 6 where the ratio between the different spectra is shown as a function of redshift. It is interesting to note that the large scale bias of the gas spectrum already shown in Figure 5 is strongly redshift dependent and starts to develop only after z ≈ 1 and is related to gas cooling and consumption within dark matter haloes. Another interesting feature at z = 0 is the higher power in the DM spectrum with respect to the total one in the DMG simulation at intermediate scales (2 < k/(h −1 Mpc) < 10). This is due to a combination of two effects: the large bias in the gas spectrum (which reduces the total one) and the small contribution of the stellar component on these scales, that are larger than typical galactic scales. For k > 10h Mpc −1 the gap between the total and DM spectra vanishes thanks to the increased contribution of cold gas and stars. We will come back to discuss the effect of stars and resolution on DMG spectra in section 6.2 A further interesting result is the comparison between DMO or DMG spectra and the halofit expressions (Smith et al. 2003). Only the 256 h −1 Mpc box can be used for this comparison. The results are shown in Figure 7. The validity of halofit up to k < 2-3 h Mpc −1 is confirmed, although discrepancies O(6 %) are found in some spectral ranges. A clear conclusion is that, when aiming at 1 % precision, the halofit expressions need to be significantly improved. Figure 8 compares the power spectra of the large and small DMO boxes with the halofit expressions at z = 0. While the power spectrum from the 256 h −1 Mpc box agrees with halofit across the entire range of its validity, the spectrum from the 64 h −1 Mpc box exhibits a shortage of power in the spectral region where non-linearity begins to display its effects. This shortage of power is due to the lack of long wavelength components in the small box, and clearly shows that L ≈ 60-70 h −1 Mpc boxes are effected by significant sample variance. Thus, we cannot draw quantitative conclusions above k ≃ 0.5-1 h Mpc −1 from such small boxes. Conversely hydrodynamical simulation results could be dependent on numerical resolution, and the higher resolution possible in small boxes make them a useful tool to address resolution issues.  The large box is in reasonable agreement with halofit in its validity range (k < 3hMpc −1 ); while the small box shows a shortage of power in that range, due to the lack of long wavelength components. Figure 9 shows the fractional difference between DMO and DMG simulations in the small box. The gas component shows a slightly different behaviour than it did in Figure 6. In the small box, the bias at large scales is already present at z = 1.2 and grows substantially to z = 0, propagating its effects at smaller and smaller scales. This is a consequence of the enhanced SF in the small box that, as mentioned in §5, peaks at z = 2.4. Owing to the (smaller) discrepancies in the overall spectra, it seems natural to argue that a systematic power depletion in the gas spectrum occurs because of stars formation in small scale sites. This is confirmed by a comparison with Rudd et al. (2008) spectra, shown in their Figure  2, right-hand panels. In their case, the power depletion in the gas spectra is even stronger, so that the overall spectra are systematically above CDM. However, Rudd et al. (2008) state in their paper that their simulations are characterized by an excess in star formation by a factor ∼ 3. Because of the high star formation, they see the same depletion in their gas power spectra, only stronger. The enhanced star formation in the small box also has an effect on the ratio between total and DM spectrum in the DMG simulation. In the small box, the total spectrum is higher than the DM spectrum for k > 10h Mpc −1 thanks to the stellar component.

Box size and resolution effects
The effect more important than resolution is the intrinsic lack of power in the small box. Thus, the total spectrum is strongly reduced in the small box with respect to the large one. This difference is on the order of 50% at k = 10h Mpc −1 and even larger, up to 100% for k = 20h Mpc −1 (see also Levine & Gnedin 2006). It is interesting to note that our results on the total, gas and DM spectra behaviour in the large box are closest to the results of Jing et al. 2006 (who used a box of 100 h −1 Mpc ), who also found the DM spectrum to lie above the total one for intermediate values of k.
The two conclusions we draw from this comparison between the small and large box are: (i) L ∼ 60-70 h −1 Mpc boxes are effected by significant sample variance, frequently yielding potentially misleading results, and hardly allowing quantitative conclusions above k ≃ 0.5-1 h Mpc −1 ; (ii) in turn, boxes with a side L ∼ 250-260 h −1 Mpc yield results that will only require slight quantitative adjustments (possibly at higher resolution) for k > 6-7 h Mpc −1 .
We are thus left with no reasons why we should not trust spectral comparisons between different models, shown in the forthcoming Figures. If some (minor) bias is present, it should affect related spectra in a similar fashion.

dDE and spectral regularities
We now turn our attention from the effects of resolution and baryons on power spectra to the effects of using different cosmologies. Figure 10 shows the fractional spectral differences between ΛCDM and dDE (upper panel) and ΛCDM and ΛCDM 1 (lower panel). Two main effects are apparent: (i) the curves are much smoother for dDE than for ΛCDM 1. In particular, in the top panel, the discrepancy at z = 0 almost vanishes, for the total spectrum. This is because ΛCDM is the W(z = 0) model of dDE, so these findings are expected. In particular, the discrepancies up to ∼ 3 % at higher z affirm the findings of Francis et al. (2007). The discrepancies between ΛCDM and ΛCDM 1 (bottom panel) are far more exaggerated and are typical of models lying ∼ 2σ's away from WMAP5 data. This shows that spectral analysis that can detect 1 % differences will successfully discriminate among models within the 2-σ curves in Figure 1. (ii) The ΛCDM1 model re-affirms that DMO-DMG spectral discrepancies become most significant at k ≃ 2-3 hMpc −1 .
In Section 2 we reviewed the technique presented in Paper I, that we shall explicitly test here. We start from the assigned model A, which is our dDE model, and find the auxiliary models W(z) for at z = 0, 1 and 2 . They must share the values of ω b,c,m = Ω b,c,m h 2 and σ8 with A. Let us recall that the former request is easily fulfilled, as H 2 = (8πG/3)ρcr and, by multiplying both sides of this relation by Ωm (or Ω b , Ωc) we have that ωm ∝ ρm ∝ a −3 , independent of the model. Conversely, the evolution of σ8 depends on the state equation of DE and its values at z = 0 and z = 24 can only be worked out once we know the constant DE state parameters w(z) of the W(z) models for z = 0, 0.5, 1, 2 . The final requirement, causing the dependence on z of the constant w's is that w is tuned so that W(z) and A have equal distances between z and the LSB.
In Figures 11 and 12 we compare spectra at z = 0 and 0.5, 1, 2 . Notice the extreme expansion of the ordinate units in these plots. Besides the gas case, the whole ordinate range is within ±1 %, for k < 10h Mpc −1 .
The basic result is that regularities persist when gas dynamics is important, although residual discrepancies (of the order of a few permils) are greater in DMG simulations. The largest discrepancy is found in the gas spectra, even if the gas spectrum distortion at z = 0, attaining ∼ 5 %, apparently does not imply a significant distortion in the global spectrum. Altogether it is fair concluding that even the inclusion of hydrodynamics keeps the discrepancies between the auxiliary model and the true model in the permil range and this is one of of the major results of this work. 6.4 Impact of the results on next weak lensing data survey Matter density spectra P (k, z) provide a basic link between observational data and theoretical models of the nature of DE. We envisage comparing the models with data from tomographic shear surveys. Such surveys allow one to work out the angular power spectrum P κ ij (ℓ), for weak lensing convergence between the i-th and j-th tomographic beams cover- ing the redshift intervals ∆i and ∆j. It is then easy to show that P κ ij (ℓ) is directly related to P (k, z), being Here H(z) is the Hubble parameter at the redshift z, r(z, z ′ ) is the comoving angular diameter distance between z and z ′ , while Wi(z) = 1.5 Ωm(1 + z) ∆ i du ng(u) r(u, z)/r(0, z) , ng(z) being the comoving galaxy number distribution on redshift. Using eq. (4), tomographic shear survey data, for 100 < ℓ < 3000, will yield matter fluctuation spectra for 0.3 ∼< k /hMpc −1 <∼ 10. The essential point, here, is that such spectra arise from any possible gravitational source. Any link with assumptions on the mass/luminosity ratio is therefore broken. Furthermore, the expected accuracy of lensing surveys will allow to recover P (k, z) with a precision O(1%).
For gravitational lensing shear predictions, at first order, there seems to be no need to take into account the baryon component. Baryons are comprise a small fraction of the Universal mass fraction compared to DM, and their distribution can be only slightly different from CDM. But what makes baryons extremely important is the precision that shear measurements will reach (O(1%)) in the near future. This presents a new challenge for hydrodynamical simulations. Formerly, a slight inaccuracy in baryon physics was a second order problem compared with the high precision achieved in tracing the total matter distribution. While baryons are less massive than CDM, the way CDM reacts to the evolution of the baryon distribution makes baryons an essential ingredient.
Let us then outline a very crucial issue concerning the extrapolation of w(z) from tomographic cosmic shear data. This point, already made in Paper I, is strengthened by this work, taking into account baryon physics. The point is that the conversion between the observed evolution of w(z) and the true one of the underlying model is not a straightforward operation. Figure 13 shows the redshift evolution of the state parameter of the A model (solid line) compared with the values of the constant-w in the different auxiliary models (W) at different redshifts (solid squares). The figure outlines that cosmic shear measures, fitted assuming constant w in each z-bin (i.e. the solid squares) would infer an evolution of w with redshift (the dotted line) radically different from the true underlying model A (solid line).
The mapping between the w values best fitting data in the tomographic layer about z and the real evolution of w(z) should be done with extreme care as shown in Casarini (2010).

DISCUSSION AND CONCLUSIONS
In this work, we employed large scale hydrodynamical simulations, including gas cooling, star formation and Supernovae feedback. Our simulations perform quite well when compared to observational results. They produce a sizable stellar component in collapsed dark matter haloes and star formation histories comparable with observations. We are aware that these results strongly depend on the adopted (low) resolution, nevertheless the similarity to observations makes us confident that the effect of baryons on the matter power spectrum in our simulation is sufficiently realistic. That said, further work is needed to better address the impact of the physics of star formation on the final results and, most important, to provide better tests of sample variance. In particular, this work makes clear that no really quantitative predictions, unaffected by sample variance, can be drawn from simulations, in L ≈ 50h −1 Mpc boxes. Matching observed accuracies of O(1 %) will require boxes a few hundred Mpcs on a side. In principle a box with L = 256 h −1 Mpc, as the one we mostly used here, should be adequate to delve into haloes, providing spectra up to k ≃ 10 h Mpc −1 , with the required accuracy. However, the spectra obtained from our large box differed in a non negligible way from the ones obtained from a smaller box (L = 64h −1 Mpc). Although being fully compatible with the sample variance of the smaller box, this suggests further investigation will be required.
These residual uncertainties do not prevent us, though, from testing the efficiency of a method introduced in Paper I, that allows us to work out power spectra prediction for any DE model based on models with a constant equation of state (w). We found up to z = 2, that the auxiliary W models (with constant w) yield spectra consistent with our assigned A model (with a polynomial w(z)) within a few permils. The test was performed for a single A dDE model. We plan to extend such a test to more dDE models in forthcoming work. There is however no reason to expect that the point selected on the wo-w ′ plane owns peculiar properties.
The residual discrepancies between A and W models exhibit a specific trend, already visible in Paper I: discrepancies (of a few permils) are greater for smaller z values. Also DMG results follow the same rule, which is even even more pronounced for baryon spectra. Gas spectra discrepancies often exceed the permils but are systematically compensated by other components. The success of our technique was not unexpected, but discrepancies at such a tiny level, also for k > 2-3 h Mpc −1 , exceeded our expectations.
The requirement that two models have the same ω c,b and σ8(z) introduces basic similarity: they will have similar baryon acoustic oscillations and typical fluctuation amplitudes. Then, the requirement that w yield the same conformal time τ elapsed since recombination for assigned and auxiliary models apparently sets an equal timing for the gradual onset of the non-linear evolution on analogous structures. However, while τ is a fair time coordinate for structures slightly detached from the cosmic background, it no longer plays such a role in cosmic islands, which have abandoned the cosmic expansion and evolve on their own, forgetting the continuing increase of the scale factor a . Within these structures, space-time is substantially Minkowskian, and time evolution is set by the ordinary time t .
These are, however, two extreme situations. The scales which "recently" abandoned the overall expansion and entered a non-linear regime are in an intermediate position.
Our inspection mostly concerns such scales. Accordingly, the success of an equal-τ requirement is anything but granted.
Here we use hydrodynamical simulations to study a larger k domain (compared to Paper I), where the effects of baryons cannot be neglected anymore. By introducing non-gravitational effects we expect dynamics to be regulated by the ordinary time t, thus producing larger discrepancies at scales k > 2-3 h Mpc −1 . Figure 11 rejects this hypothesis by showing that the matter distribution is still ruled by the primeval clock. In particular, for k > 2-3 h Mpc −1 , the overall DMG spectral discrepancies, at z = 0, are perhaps even smaller than DMO. Curiously enough, the opposite seems true at some greater z's, suggesting that most discrepancies found are pure noise, rather than systematics. Clearly for larger k's, when we delve into regions that virialized long ago, A and W models unavoidably differ. However, while we restrict ourselves to the scales that forthcoming weak lensing surveys will study, the Paper I requirement still works unexpectedly well.
One ought to beware, however, that the success of our approach does not imply that constant-w and variable-w models are indiscernible. What makes them observationally different is the measured dependence of w on the redshift z.
An essential issue, however, is that the observed constant state parameter w(z) will NOT be the physical one of the DE component. Although we expect data to fit, at any z, the same density and Hubble parameter values, the observed w(z) will be the state parameter that a constant-w model must have, at such z, to grant the same distance between there and the LSB. The relation between the observed w and the intrinsic one, for the DE component, should be put under control; otherwise there is a severe risk that future data analysis lead to DE state equation misinterpretation.
A final point we wish to make is that halofit expressions, so useful for predicting spectra until now, have been tested and found compatible with their claimed accuracy. Unfortunately, accuracy at (∼ 6 %) no longer matches our needs. Furthermore, among models that approach the present data, halofit expressions only need to handle ΛCDM . If one aims to improve halofit expressions, there is no longer any need to cover the wide array of models that once were relevant, but now are far from observation.
Amelioration, however, is needed in two separate directions: (i) The spectral range covered by predictions must reach k ≃ 10 h Mpc −1 , which requires that gas be considered; Rudd et al. (2008) had actually attempted an extension in this direction, but limited themselves to ΛCDM cosmologies. (ii) DE state parameters different from constant w ≡ −1 should be included. Moreover, they should aim for a precision of O(1 %), which presents a severe challenge, even for a much more restricted model range. Nevertheless when halofit expressions of the required precision will be available, even only for constant-w models, thanks to the results of this work it will be straightforward to extend them to any w(z) law.