Cosmic Degeneracies III: Nbody Simulations of Interacting Dark Energy with NonGaussian Initial Conditions
Abstract
We perform for the first time Nbody simulations of Interacting Dark Energy assuming nonGaussian initial conditions, with the aim of investigating possible degeneracies of these two theoretically independent phenomena in different observational probes. We focus on the largescale matter distribution, as well as on the statistical and structural properties of collapsed halos and cosmic voids. On very large scales, we show that it is possible to choose the Interaction and nonGaussian parameters such that their effects on the halo power spectrum cancel, and the power spectrum is indistinguishable from a model. On small scales, measurements of the nonlinear matter power spectrum, halomatter bias, halo and subhalo mass function and cosmic void number function validate the degeneracy determined on large scales. However, the internal structural properties of halos and cosmic voids, namely halo concentrationmass relation and void density profile, are very different from those measured in the model, thereby breaking the degeneracy. In practice, the values of required to cancel the effect of interaction are already ruled by observations. Our results show in principle that the combination of large and smallscale probes is needed to constrain Interacting Dark Energy and Primordial nonGaussianity separately.
keywords:
dark matter – dark energy – cosmology: theory, largescale structure – galaxies: formation1 Introduction
According to the most recent measurements of cosmic microwave background (CMB) anisotropies performed by the Planck satellite mission (Ade et al., 2016a), the standard CDM cosmological model is still extremely successful in reproducing different observational datasets. This in turn favours the more economic cosmological constant as an explanation of the latetime cosmic acceleration over alternative and more complex Dark Energy (DE) or Modified Gravity (MG) models. Nonetheless, theoretical problems in understanding the energy scale and the time evolution of (known as the finetuning and coincidence problems, respectively, see e.g. Weinberg, 1989; Padilla, 2015) as well as recent observational tensions between CMB cosmological constraints and those inferred from independent probes in the local Universe (see e.g. Heymans et al., 2013; Hildebrandt et al., 2017; Simpson et al., 2016; Vikhlinin et al., 2009; Ade et al., 2016e) motivate the investigation of such alternative and more complex scenarios.
In particular, various possible realisations of Interacting Dark Energy models (hereafter IDE, see e.g. Wetterich, 1995; Amendola, 2000; Pettorino & Baccigalupi, 2008; Amendola et al., 2008; Baldi, 2011a, 2012c; Pourtsidou et al., 2013) based on a direct energymomentum exchange between a DE scalar field and the CDM particle sector, have attracted significant interest and for small values of the interaction strength appear still consistent with current CMB constraints (Salvatelli et al., 2013; Costa et al., 2014; Salvatelli et al., 2014; Ade et al., 2016b).
Similarly, measurements of higherorder statistics of the CMB anisotropies are consistent with a nearly Gaussian distribution of the primordial curvature perturbations by providing very tight constraints on the Primordial NonGaussianity (hereafter PNG) parameters: and for the local and equilateral configurations, respectively (Ade et al., 2016c). As some level of nonGaussianity in the primordial density distribution is a common and clean prediction of basically all models of inflation (see e.g. Maldacena, 2003) – i.e. the hypothetical mechanism driving the early exponential expansion of the Universe – measurements of PNG are considered as a smoking gun to discriminate between various inflationary models (Bartolo et al., 2004; Giannantonio et al., 2014).
Despite the tight constraints on the PNG amplitude from Planck data, LargeScale Structure (LSS) observations in the late universe coming from the next generation of widefield galaxy redshift surveys could outperform these constraints. More specifically, recent measurements of galaxy clustering and of the integrated SachsWolfe (ISW) effect already provide constraints of (Ross et al., 2013; Giannantonio et al., 2014; Leistedt et al., 2014), while future redshift galaxy surveys like Euclid (Laureijs et al., 2011) and SKA (Maartens et al., 2015; Camera et al., 2015) are forecast to outperform the CMB in constraining PNG, especially via the multitracer method (Alonso & Ferreira, 2015; Fonseca et al., 2015).
This is possible due to the various observational signatures that PNG imprints on LSS at late times, namely on the abundance of massive objects (which can be either enhanced or suppressed for positive and negative values of the PNG amplitude, respectively), on the bias between galaxies and the underlying matter distribution (that becomes scaledependent on large scales in the presence of some PNG) and on the 3point correlation function of galaxies that encodes the shape of PNG (for more details, see e.g. Desjacques & Seljak, 2010; Liguori et al., 2010; Desjacques et al., 2018).
Recent studies on the effects of IDE models on structure formation (Baldi et al., 2010; Moresco et al., 2014; Hashim et al., 2014; Duniya et al., 2015; Cui et al., 2012) have shown that similar features may arise also in these models. In particular, Hashim et al. (2014) showed that the observational signatures of IDE and PNG on the largescale galaxy power spectrum can mimic each other. This is due to the fact that some models of IDE introduce a scale dependence in the matter density contrast on very large scales, mimicking PNG scaledependent halo bias. Also on nonlinear scales, numerical simulations of IDE (see e.g. Baldi & Pettorino, 2011; Baldi, 2012b, a; Cui et al., 2012) and of PNG (Grossi et al., 2007; Pillepich et al., 2010; Wagner et al., 2010; LoVerde & Smith, 2011) scenarios showed that IDE enhances the abundance of massive halos in a similar way to PNG with a positive amplitude.
This degenerate behaviour between PNG and IDE indicates that separate observational constraints on the PNG amplitude and the IDE interaction rate could be misinterpreted or possibly that their joint effects could become indistinguishable from the standard CDM reference model. This represents the main motivation for the present work, where we will present for the first time a joint numerical analysis of nonlinear structures forming from PNG initial conditions through an IDE cosmological evolution. Our main goal is to test whether such degeneracy holds for all observables at all scales and if not to identify specific statistics that clearly disentangle the two phenomena.
This paper is organized as follows: in Sec. 2 we introduce IDE and PNG extensions to the standard scenario. In Sec. 3 we use the linear halo power spectrum as an observational probe to test the IDE–PNG degeneracy on large scales. In Sec. 4 we test the IDE–PNG degeneracy on nonlinear scales by running a set of NBody simulations for all models under consideration. In Sec. 5 we present all results for the nonlinear matter power spectrum, halomatter bias, halo mass function, subhalo mass function, halo concentrationmass relation, void number density and void density profiles. Finally, conclusions are summarised in Sec. 6.
2 NonStandard Cosmological Models
In this section, we present the two nonstandard extensions to the fiducial model that we will consider in this work. The first extension is based on the assumption of a nonGravitational interaction between CDM particles and a dynamical DE scalar field. The other extension relays on a nonGaussian distribution of the primordial density field as generically predicted by inflationary models.
2.1 Interacting Dark Energy
Various models of IDE have been proposed in the literature over the past two decades (see e.g. Amendola, 2000, 2004; Koyama et al., 2009; Baldi, 2011a; Clemson et al., 2012). In this paper we consider the most widely studied example of such models based on a quintessence dynamical scalar field playing the role of the DE, subject to a selfinteraction potential and to a direct interaction with the CDM fluid via energymomentum exchange (Bertolami & Martins, 2000; Amendola, 2000). The background evolution of such cosmological scenarios is governed by the KleinGordon equation for the scalar field: {ceqn}
(1) 
and by the continuity equations of the different components that contribute to the total energy density of the universe:
(2)  
(3)  
(4) 
as well as by the Friedmann constraint {ceqn}
(5) 
where and are the energy density of CDM, baryons and radiation, respectively. An overdot represents a derivative with respect to the cosmological time . The Hubble function is defined as where is the scale factor and . The parameter represents the energy density of the DE fluid defined as . The righthand side source terms in Eqs. (1) and (2) represent the interaction parameter between CDM particles and DE that is proportional to the CDM energy density through the dimensionless constant that sets the strength of the coupling. The sign of the term determines the direction of the energymomentum exchange between the two interacting components. In order to fulfil Bianchi identities and not violate total energymomentum conservation, the source terms in Eqs. (1) and (2) should be equal and have opposite sign.
By integrating the CDM conservation equation (2) one gets the time evolution of the CDM density as: {ceqn}
(6) 
which shows a basic property of IDE models: matter density is not separately conserved as the energy exchange results in a timedependent CDM particle mass. In this work, we consider the exponential form for the selfinteraction potential (Wetterich, 1988; Lucchin & Matarrese, 1985), {ceqn}
(7) 
where and are constants.
In the Newtonian gauge, the perturbed metric (assuming flatness and vanishing anisotropic stress) is given by {ceqn}
(8) 
where is the gravitational potential. The Poisson equation^{1}^{1}1We also ignore baryons for simplicity. reads as (Hashim et al., 2014): {ceqn}
(9) 
where are the comoving density contrasts and are the peculiar velocity potentials of CDM and of the scalar field , respectively. The coupling term in Eq. (9) is proportional to the velocity potentials and introduces a scaledependence to the matter growth factor on large scales. Since are gaugeinvariant, the resulting signal is an explicit coupling effect and not a false gauge effect.
The perturbed conservation equations are then given by (Hashim et al., 2014)
(10)  
where and indicates CDM and scalar field respectively, is the soundspeed of the th species (which is vanishing for CDM while for DE perturbations ), is the total equation of state, is the total peculiar velocity potential and is the momentum transfer potential given by (Koyama et al., 2009) {ceqn}
(12) 
The source term on the right hand side of Eq. (LABEL:CDM_pertrubEq) is given by
(13)  
These equations fully specify the evolution of the linear gaugeinvariant perturbations of the coupled system, we refer the interested reader to Hashim et al. (2014) for a more complete derivation of these equations.
As we will be interested in the evolution of the system at small scales and beyond the linear regime, we also recall (see e.g. Amendola, 2004) that in the Newtonian limit the evolution equation for CDM density perturbations, Eqs. (10) and (LABEL:CDM_pertrubEq) imply: {ceqn}
(14) 
since comoving and Newtonian density contrasts are equal, i.e. , and we ignore derivatives of scalar field perturbations. The coupling terms in Eq. (14) are: , which represents an extra friction arising as a consequence of momentum conservation, and , which is responsible for the fifth force acting on CDM perturbations.
2.2 Primordial NonGaussianity
Local type nonGaussianity in the primordial curvature perturbations, that maximizes the bispectrum in the squeezed shape, is parametrized by {ceqn}
(15) 
where is the Gaussian gravitational field and is the PNG parameter. Singlefield inflation models predict a very small value of (Maldacena, 2003), but multifield models can generate large nonGaussianity in squeezed configurations (Moroi & Takahashi, 2001; Lyth & Wands, 2002).
On large scales, PNG enhances the large peaks of matter perturbations (Matarrese & Verde, 2008; LoVerde et al., 2008; Matarrese et al., 2000). This introduces a scaledependent signal in the bias between the virial collapsed objects at high peaks and the underlying traced matter. By measuring the cross halomatter power spectrum in NBody simulations with localtype nonGaussian initial conditions, many authors have confirmed that the largescale bias is scaledependent (see e.g. Dalal et al., 2008; Pillepich et al., 2010): {ceqn}
(16) 
where is the matter autopower spectrum, is the Gaussian bias and {ceqn}
(17) 
with being the critical overdensity for halo collapse, the transfer function and the linear dark matter growth factor which is normalised to in the matter dominated era. On very large scales, and so . Since we only consider local type PNG in the current analysis, for simplicity we drop the loc superscript from our notation.
Since IDE introduces a scale dependence in the matter growth factor and nonnegligible DE perturbations in the Poisson equation (9), the scaledependent PNG bias for IDE models becomes {ceqn}
(18) 
where the effect of IDE appears in the scale dependence of and in the factor {ceqn}
(19) 
where is the DE growth factor and depends on the coupling parameter though the background equations (1) and (2). We can notice that on very large scales, behaves as .
3 Linear Halo Power Spectrum
Parameter  Value 

0.703  
0.0451  
0.2711  
0.729  
2.42  
0.966 
In this section, we will illustrate the degeneracy between IDE and PNG by computing the halo power spectrum on linear scales for both models and for their combination.
The halo power spectrum is given in general by {ceqn}
(20) 
In order to compute this we first numerically solve Eqs. (10) and (LABEL:CDM_pertrubEq) for the growth factors , and then calculate the matter power spectrum using (Ade et al., 2016d): {ceqn}
(21) 
where is the spectral index, is the spectral amplitude and is the pivot scale. We use (Lewis et al., 2000) to compute the transfer function . We then apply the bias relation, Eq. (18), to the matter power spectrum as given in Eq. (20)^{2}^{2}2For the Gaussian bias, we use the ansatz .. We adopt the cosmological parameters given in Table 1.
In computing the growth rate of CDM density perturbations we consider both the case where largescale perturbations in the DE scalar field are properly taken into account () and the the case where such perturbations are artificially set to zero (). The latter case, while being not fully consistent, allows us to match the approximations adopted in the numerical treatment that we will discuss below and to obtain a more direct correspondence between the PNG and IDE parameters that are expected to provide a strong degeneracy in the nonlinear regime under such approximations. In Fig. 1, we show the ratio of the linear halo power spectrum to the fiducial model at for the cases given in Table 2.
DE  

\@slowromancapi@  0.05  0.0  
\@slowromancapii@  0.05  0.0  
\@slowromancapiii@  0.0  151.51  
\@slowromancapiv@  0.0  166.66  
\@slowromancapv@  0.05  151.51  
\@slowromancapvi@  0.05  166.66 
These values of are obtained by minimizing the residual for the combined model, i.e. they correspond to the values of maximum degeneracy for a DECDM coupling parameter for the cases and . Therefore, for these combinations of parameters, as clearly seen in Fig. 1, IDE and PNG are strongly degenerate with each other, in the sense that their combination is indistinguishable from the fiducial CDM case^{3}^{3}3We chose and because the same degeneracy does not apply for negative in the Newtonian approximation, since the coupling enters also as a term in Eq. (14). This means that for and , we do not expect a degeneracy in the nonlinear regime. .
Although these derived values of are at least one order of magnitude larger than currently allowed by observational constraints, we will continue to use these values as a toy example of the IDEPNG degeneracy. Realistic values for with the standard assumption of a scaleindependent amplitude of PNG would have too weak effects on nonlinear structure formation to significantly influence the observational features for any nonnegligible coupling parameter . On the other hand, scaledependent PNG (see e.g. RenauxPetel, 2015), where evolves with wavenumber , may still provide an effective at scales relevant for nonlinear structure formation, while remaining consistent with current bounds around the Planck pivot scale Mpc.
In order to model the mimicking degeneracy relation between and that is illustrated in Fig. 1, we repeat the procedure of minimizing the residual for a wide range of the parameters and , for both perturbed and nonperturbed DE cases. We find that relation {ceqn}
(22) 
where and are constants, provides a good fit, with exponent . This is shown in Fig. 2, where the numerical results are overplotted with the fitting function Eq. (22) for perturbed and nonperturbed DE cases. Note that the degeneracy slope increases if we assume nonperturbed Dark Energy.
4 Degeneracy on nonlinear Scales
It is well known that IDE and PNG separately imprint characteristic features in the nonlinear regime of structure formation, which can be tested through different observational probes. For example, IDE affects the highmass tail of the halo mass function (HMF) by enhancing the abundance of halos (Cui et al., 2012), while PNG impacts the number of massive CDM halos, suppressing (increasing) it for negative (positive) (see e.g. Grossi et al., 2009; Wagner et al., 2010). It is therefore plausible that some form of degeneracy may appear also at these nonlinear scales, and in particular that the combination of IDE with a negative value of for PNG may result in a HMF hardly distinguishable from the reference CDM case at all masses.
IDE also shows distinctive features on other observational probes, including higherorder correlation functions and nonlinear bias, in a similar way to PNG (Moresco et al., 2014; Desjacques et al., 2009; Wagner & Verde, 2012). IDE further affects the structural properties of CDM halos and voids (Pollina et al., 2017, 2016; Giocoli et al., 2013; Baldi, 2014, 2011b), and PNG is also expected to show significant effects on these probes (Neyrinck & Yang, 2013; Abel et al., 2012; Sutter et al., 2014).
This implies that the mimicking degeneracy which we have found at linear scales for the halo power spectrum may persist (fully or partly) in some smallscale nonlinear observables, while it may be broken by others. In the following, we test the linear degeneracy relation, defined in Eq. (22), on nonlinear scales by analysing a suite of cosmological Nbody simulations that include IDE and PNG, both separately and in a combined way. To this end, we will consider various nonlinear probes, starting from the nonlinear matter power spectrum and the halomatter bias to the statistical and structural properties of CDM halos and voids.
4.1 NBody Simulations
In order to consistently account for the effects of IDE in the nonlinear regime, we made use of a modified version of the parallel TreePM NBody code GADGET (Springel, 2005) that incorporates all the specific features of the coupling between DE and CDM, i.e. modified background expansion, CDM particle mass time variation, the extra friction and the fifth force acting on CDM particles (see Baldi et al., 2010, for a detailed description of the modified Nbody algorithm). The simulations follow the evolution of CDM particles within a periodic cosmological box of Gpc per side, for all the cosmological parameters given in Table 1, with a mass resolution at of and softening length . Our numerical implementation of IDE assumes that DE perturbations are negligible in relation to structure formation processes compared to the dominant effects of background evolution, extra friction and fifth force. This is a valid approximation on subhorizon scales; it becomes less accurate at scales comparable with the cosmic horizon, but this is beyond the fundamental mode of our 1 boxes. For this reason, we have chosen to consider the same approximation (i.e. ) to select our combination of values for the parameters and , so as to ensure consistency between the degeneracy relation displayed in Fig. 2 and the outcomes of our Nbody simulations at small scales.
In order to generate the initial conditions for NBody simulation of all models considered in this paper, we slightly modified the publicly available code 2LPTic (Scoccimarro et al., 2012). The algorithm implements nonGaussian initial conditions with external Hubble and growth functions consistent with IDE modifications.
The nonGaussian initial conditions are generated for localtype PNG with an extra nonGaussian term according to Eq. (15), where is a random realization of a Gaussian field with the primordial power spectrum . Then, the linear density field is obtained from the nonGaussian potential through the Poisson equation: {ceqn}
(23) 
where the transfer function is computed using (Lewis et al., 2000) for the fiducial cosmology. We assume the transfer function is not affected by the latetime interaction (Baldi et al., 2010; Baldi, 2012b). The growth function for all models is normalized at to directly compare the impact of IDE on the structure growth in the period between and the present time. For PNG, we set as the value corresponding to the interaction rate on linear scales (with ), as determined by Eq. (22).
Particle positions are then displaced from a homogeneous glass distribution (Baugh et al., 1995) using the Zel’dovich approximation (Zeldovich, 1970) according to the displacement field at the initial redshift . In order to compute particle initial velocities, we used the relation , where the growth rate function is derived for each model by solving Eqs. (14) for the growth function. For the IDE–PNG combined model, we apply the growth function of IDE after transforming the initial Gaussian potential to the nonGaussian form according to Eq. (15). In order to minimize the sampling variance, we used the same initial random seed for all the simulations.
5 Results
In this section, we present the main results of our numerical simulations of IDE, PNG and the combined IDE–PNG extensions rescaled with respect to the fiducial model. We focus mainly on the nonlinear matter power spectrum, the halomatter bias and the statistical and structural properties of CDM halos and voids.
5.1 The nonlinear matter power spectrum
We computed the nonlinear matter power spectrum for each simulation by calculating the density field using a CloudinCell mass assignment on a cubic grid with the same resolution as the Particle Mesh grid used for the integration of the Nbody system (i.e. ). According to this procedure, the nonlinear matter power spectrum is determined up to the Nyquist scale, . We truncate the resulting power spectrum at the mode where the shot noise is below of the measured power. From the simulated power spectra, we can estimate the effects of IDE, PNG and, for the first time, the joint effects of IDE and PNG, on linear and nonlinear scales at different redshifts.
In Fig. 3 we display the ratio of the nonlinear matter power spectrum for IDE, PNG and their combination, to that of the standard cosmological model, at and . The plots show the following features.
IDE with Gaussian initial conditions (dashed green curve with solid diamonds) – shows the expected scaledependent power enhancement at nonlinear ranges due to the combined effects of the fifth force and of the extra friction associated with the DECDM interaction. Also, since we normalize the power spectrum at the redshift of the CMB, the normalization at linear scales () is increased by about relative to the standard model (i.e. ), due to the higher linear growth rate in the IDE case (Baldi, 2011a); this is consistent with previous works (see e.g. Baldi, 2011b). At higher redshift (right panel) the nonlinear power spectrum enhancement due to the IDE fifth force is slightly reduced, with the peak ratio shifted towards smaller scales.
NonGaussian initial conditions in (dashed blue curve with solid squares) – shows the expected suppression of power at small scales, relative to the standard case. The deviation is larger at higher redshifts and the minimum shifts towards smaller scales, in agreement with predictions from the halo model presented in Fedeli & Moscardini (2010).
IDE and PNG combined (dashed red curve with solid triangles) – the ratio no longer shows any significant scale dependence down to ranges corresponding to the location of the peak/minimum in the ratio for the two separate models. The difference in the power normalization at linear scales associated with the enhanced growth rate in IDE remains unchanged. This seems to indicate a mimicking degeneracy between IDE and PNG in the matter power spectrum on nonlinear scales while there is no degeneracy on linear scales. Remarkably, the figure shows that there is a nonlinear mimicking degeneracy for the same combination of parameters that produce mimicking degeneracy in the halo power spectrum at much larger scales, as described by Eq. (22) and Fig. 2.
In the figures we have also overplotted, for comparison, a black dotted curve representing a simple superposition of the two effects, i.e. the PNGonly (blue squares) deviation times the IDEonly (green diamonds) deviation. The very good agreement of this simple prediction with the actual power measured from the combined IDEPNG simulation seems to indicate that the two effects acting on structure formation are decoupled – which suggests that full combined Nbody simulations may be unnecessary in order to compute the combined power spectrum for other combinations of and .
5.2 HaloMatter Bias
Following the standard hierarchical clustering scenario of structure formation, halos and galaxies are biased tracers of the underlying matter distribution. In this section, we compute the linear bias between halos and the underlying dark matter density field, as the ratio between the haloCDM cross power spectrum and the auto power spectrum in Fourier space: {ceqn}
(24) 
(We suppress the dependence for simplicity.) This bias estimator is used to avoid shotnoise (Hamaus et al., 2010; Smith et al., 2007; Baldauf et al., 2010; Baldauf et al., 2013), and we follow the approach of VillaescusaNavarro et al. (2014) for the computation of the two power spectra.
In Fig. 4, we show the ratio of the halomatter bias for the IDE, PNG and IDE–PNG models, relative to the fiducial model. As expected, PNG introduces a clear scaledependence at large scales. On the contrary, the bias in the IDE model appears to have a slightly lower normalization than though retaining the same evolution with scale as the standard scenario. This different behaviour is most visible at higher redshifts, as shown in the right panel of Fig. 4, where the scaledependence of the PNG simulation is stronger. On nonlinear scales, both PNG and IDE show a maximum deviation relative to the reference model but in opposite directions, with the amplitudes of the peak/ minimum increasing and their position moving towards smaller scales at higher redshifts. These outcomes are all consistent with the previous literature (Matarrese & Verde, 2008; Desjacques et al., 2009; Moresco et al., 2014; Marulli et al., 2012) and qualitatively show how the halo bias is affected at similar scales for both IDE and PNG.
For the combined IDE–PNG scenario, we find that at the halo bias retains some scaledependence on large scales, i.e Mpc, while it is nearly scaleindependent on scales Mpc Mpc. Furthermore, it retains the lower normalization that characterizes the IDE model at all scales. This combination of the two effects is more clear at higher redshift, where we can clearly identify two distinct regions for scaledependent (Mpc) and scaleindependent (Mpc) deviations from the reference model. Also in this case, the simple superposition of the two separate effects, very accurately reproduces the behaviour of the combined IDE–PNG simulation, thereby suggesting that the two phenomena act on the biasing of collapsed structures independently.
Similar to the nonlinear matter power spectrum, the halomatter bias satisfies the degeneracy relation, Eq. (22), on nonlinear scales, while this is broken at larger scales. We argue that this may be due to the fact that our NBody implementation of IDE as discussed above (see also Baldi, 2011b), does not account for the contribution of largescale DE perturbations (i.e. it assumes the approximation ) which would result in a scaledependent growth function. Therefore, including the effects of largescale DE perturbations should boost in a scaledependent way the IDE linear power spectrum and consequently the halomatter bias on large scales. This would recover the result of a mimicking degeneracy at all scales that was obtained from linear perturbation theory (Sec. 3). A proper verification of this conjecture would require major modifications to our Nbody codes, that go beyond the scope of the present paper, and we defer an extensive study on this subject to future works.
5.3 Statistical and structural properties of CDM halos
In this section we test the degeneracy relation in the statistical and structural properties of CDM halos.
5.3.1 The Halo Mass Function
We identified collapsed halos in our simulations following a standard procedure, amounting to a first identification of particle groups by means of a FriendsofFriends (FoF) algorithm with linking length , where indicates the mean interparticle separation. On top of these FoF halos we run the SUBFIND algorithm (Springel et al., 2001) in order to identify gravitationally bound substructures present within each group. The latter procedure allows to assign to each FoF group the virial mass of its primary substructures, defined as the mass of a spherical region with its centre on the particle with the halo’s minimum potential enclosing a mean overdensity equal to times the critical density of the universe.
Given these halo catalogues, we computed the halo mass function for IDE, PNG and the combined IDE–PNG models by binning the halo masses into 13 logarithmically equallyspaced mass bins over the mass range . The lower mass bound is set by the minimum halo mass resolved in the fiducial model, composed of at least particles.
In Fig. 5 we show the ratio of the cumulative HMF to the model for IDE, PNG and the combined IDE–PNG models. As expected, IDE enhances the abundance of large mass halos with respect to the standard case, while PNG shows on the contrary a suppression of the abundance of halos in the highmass tail, consistent with previous results (Cui et al., 2012; Wagner et al., 2010).
The combined IDE–PNG model shows some level of degeneracy with the standard cosmology at , with the combined mass function being slightly lower than in the pure IDE case. The degeneracy becomes more clear, given the larger amplitude of the individual effects, at higher redshifts (), where the IDE and PNG deviations from the reference model reach about at the largest masses with abundance suppression by only in the combined case. Furthermore, the exponential dependence on halo mass of the deviation with respect to is also significantly weakened in the combined model. Nonetheless, as a mimicking degeneracy is never fully attained, the halo mass function seems not to follow the degeneracy relation of Eq. 22, thereby providing a possible way to disentangle these phenomena.
5.3.2 The subhalo mass function
As a further statistic of structure properties at small scales, we computed – for all our simulated cosmologies – the subhalo mass function, defined as the number of subhalos of mass within a main halo of virial mass . In Fig. 6, we display the ratio of the subhalo mass function with respect to the measurements in the simulation, as a function of the mass ratio . In order to avoid resolution effects, we consider only subhalos hosted by clustersize halos, i.e. systems with . We underline to the reader that the measured subhalo counts in the model are characterized by the typical slope of approximately consistent with different previous findings Gao et al. (2004); Giocoli et al. (2010); Despali & Vegetti (2017).
As can be seen from the figure, IDE suppresses the abundance of substructures over the whole range of subhalo fractional mass, even though the effect is small (about 35%). On the contrary, PNG enhances the abundance of subhalos up to about (for the highest values of the subhalo fractional mass) over the same mass range. The combined IDE–PNG case shows again a quite clear degeneracy, with a suppression never exceeding , marginally consistent with the Poissonian error range of the model. The simple superposition of IDE and PNG models is in reasonable agreement with the combined IDE–PNG simulation. These results underline that the degeneracy relation seems to remain valid also at the level of CDM halo substructures.
5.3.3 Halo Concentration
Finally, we conclude our investigation of the combined effects of IDE and PNG on structural properties of collapsed halos by computing the average halo concentration as a function of halo mass, which is usually known as the concentrationmass relation (Zhao et al., 2009; Giocoli et al., 2012). In order to compute the concentrations for the halos identified in our simulations, we adopt the NFW formula used in Springel et al. (2008): {ceqn}
(25) 
where is the characteristic overdensity, is the halo concentration, is the maximum circular velocity of the halo attained at radius . In Fig. 7, the ratio of the concentrationmass relation of IDE, PNG and the combined IDE–PNG models relative to is presented at (left panel) and (right panel).
As expected, IDE halos are found to be less concentrated with respect to the fiducial case, in agreement with results given in (Baldi, 2011b). Similarly, PNG with also suppresses halo concentrations (the opposite would occur for a positive ). Therefore, for the first time we encounter an observational probe showing deviations from CDM pointing in the same direction for IDE and our negative PNG scenarios.
The combined IDE–PNG simulation, accordingly, shows an even stronger suppression of the concentrationmass relation relative to the model than the two individual models separately. The effects are less pronounced at higher redshifts, while the trends and the relative ordering of the various models is preserved. Superposition of the individual effects of IDE and PNG seems to agree well with the combined simulation. This however indicates that the degeneracy is broken for the CDM halo concentrationmass relation, which might then represent another direct way to disentangle the models, when combined with another more degenerate probe. It is also reasonable to emphasize that this effect is relatively small; only future wide field observational campaigns – like the future ESAmission Euclid (Laureijs et al., 2011) – will be able to collect the large number of galaxy groups and clusters (Sartoris et al., 2016) necessary for these tests.
5.4 Statistical and structural properties of cosmic voids
In this section, we move our focus to underdense regions of the universe by testing whether cosmic voids also follow the degeneracy relation. In order to identify cosmic voids in our set of simulations, we employ the publicly available void finder VIDE (Sutter et al., 2015), which is based on the ZOBOV algorithm (Neyrinck, 2008). The cosmic void identification is mainly done by means of a Voronoi tessellation scheme that associates a polyhedrical cell to each particle tracing the CDM density field. Subsequently, cell volumes are compared in order to identify local density minima, i.e. cells with a larger Voronoi volume than all their surrounding cells. A hierarchy of identified voids is then obtained via the watershed transform algorithm (Platen et al., 2007), by joining Voronoi cells around a local density minimum. In our analysis, we consider only voids with a central density that is below the density of the universe by and a lower density contrast limit , corresponding to a probability of voids arising from Poisson noise below (Neyrinck, 2008).
5.4.1 Void number function
As a first statistics for cosmic voids, we study their abundance as a function of the void effective radius , defined as the radius of a sphere centred on the most underdense particle of a void and having the same volume as the Voronoi volume of the void: {ceqn}
(26) 
In Fig. 8 we show the ratio of the void number functions relative to the ones in the cosmology for all models under consideration, as a function of the effective radius at . From the figure we see that IDE suppresses the number of cosmic voids with effective radius Mpc by about relative to the case, and correspondingly enhances the abundance of larger voids by up to . The trend is qualitatively similar, though quantitatively weaker, for PNG, in agreement with previous results of Kamionkowski et al. (2009).
The combined IDE–PNG simulation shows suppression of the void number function for radii Mpc, similar to the IDE case and barely enhances the void abundance at Mpc relative to the case, so that it is indistinguishable within the Poisson error range at these radii. As we did for all previous observables, we also compute the simple superposition of the two effects, by taking the product of the two separate deviations with respect to the reference case. For the first time in our analysis, we see that such a superposition fails to reproduce the results of the combined simulation at large void effective radii: this follows from comparing the black dotted curve, representing the analytical superposition, with the blue squares, showing the combined simulation in Fig. 8.
This suggests that the two phenomena interplay in some way in shaping the growth of large cosmic voids, and cannot be considered as fully independent in this regime. In any case, we notice that the degeneracy is fulfilled by the abundance of cosmic voids with large effective radii (Mpc), while it does not seem to apply at smaller void radii.
5.4.2 Void density profiles
To further check the degeneracy on cosmic void structural properties, we computed the average void density profiles for two different bins of void radius, namely Mpc and Mpc. We do this by stacking individual density profiles of randomly selected voids, for each radius bin, corresponding among the different cosmological simulations. We display the ratio of the resulting void mean density profiles in Fig. 9 for all considered models, relative to at . The grey area represents the confidence limit, computed by means of a bootstrap resampling technique.
Again, we compare the observational signature of the individual IDE and PNG models with their combination. As can be seen from the plot, cosmic voids in the IDE case tend to have a lower inner density than their counterparts. This indicates that cosmic voids are emptier in the IDE case, fully consistent with previous results (see e.g. Pollina et al., 2016). Correspondingly, the compensating overdensity around the effective radius is found to be more prominent than in . On the other hand, PNG shows a negligible effect on cosmic void density profiles. It is then not surprising that the combined IDE–PNG model also shows lower density profiles in the central regions of the voids. This result also shows that cosmic voids do not seem to follow the same degeneracy relation that applies for most of the observables related to properties of the overdense regions of the universe.
6 Discussion and Conclusions
The concept of observational degeneracy in cosmology arises in several different forms: (1) Parameter Degeneracy represents the existence of large error correlations between different model parameters for specific measurements (Efstathiou & Bond, 1999; Crooks et al., 2003; Tereno et al., 2005; Howlett et al., 2012); (2) Dark Degeneracy reflects the fact that gravitational experiments measure the energymomentum tensor of the total dark sector and splitting into Dark Energy and Dark Matter is arbitrary (Kunz, 2009; Aviles & CervantesCota, 2011); (3) Mimicking Degeneracy occurs when cosmological models different from the standard mimic some of its specific features, like background expansion and the growth of matter perturbations (Fay et al., 2007; Setare & Mohammadipour, 2013; Fay, 2016).
Cosmic degeneracy of IDE has been investigated in the literature (Clemson et al., 2012; Väliviita & Palmgren, 2015), including the partial mimicking degeneracy of IDE and MG (Wei & Zhang, 2008; Koyama et al., 2009; Wei et al., 2013). A mimicking degeneracy between PNG in the power spectrum in the Newtionian approximation, and the correct general relativistic power spectrum with Gaussian initial conditions, has been shown by Bruni et al. (2012); Jeong et al. (2012). Also, parameter degeneracy has been investigated in the nonGaussian halo bias by Carbone et al. (2010). Moreover, Abramo & Bertacca (2017) investigated the degeneracy of largescale velocity effects on galaxy clustering with the (local) nonGaussianity parameter , by simulating galaxy surveys and combining the clustering of different types of tracers of largescale structure. They studied how largescale velocity contributions could be mistaken for the signatures of primordial nonGaussianity (see also Raccanelli et al., 2014, 2018).
In this paper – as part of a Cosmic Degeneracies paper series (Baldi et al., 2014; Baldi & VillaescusaNavarro, 2018) – we have considered the mimicking degeneracy between IDE and PNG that was first shown in linear perturbation theory by Hashim et al. (2014). Since IDE can mimic PNG, the possibility exists that we can choose IDE and PNG parameters such that the two effects cancel, i.e., produce standard behaviour. We confirmed this mimicking degeneracy in the halo power spectrum on very large scales, i.e. , based on purely analytical calculations in the linear regime. We then fitted the degeneracy relation with a power law, (depicted in Fig. 2), by minimizing the residual of the halo power spectrum for the combined IDE–PNG model with respect to the mimicked model.
To further investigate and validate the degeneracy, Eq. (22), at nonlinear scales, we employed a suite of specifically designed NBody simulations including the effects of IDE and PNG, both separately and combined with each other, and we extracted from these simulations a set of standard statistics, namely:

The nonlinear matter power spectrum, for which we observed that the mimicking degeneracy persists, remarkably, on nonlinear scales in the sense that the scaledependent deviation with respect the reference scenario characterising the two separate models at nonlinear scales disappears in the combined simulation even though the difference in the linear power normalisation due to the enhanced growth rate in IDE is not removed;

The halo matter bias, for which we find similarly to the nonlinear power spectrum, that the scaledependence imprinted by the two different models at nonlinear scales is also strongly suppressed in the combined simulation while on linear scales such scaledependent feature is retained and so breaks the observed degeneracy;

The halo mass function, which also shows some level of degeneracy though not satisfying Eq. (22) for the degenerate values thus allowing us to disentangle the observed degeneracy;

The subhalo mass function, also showing mimicking degeneracy over the whole subhalo mass range availabe in our simulations;

The halo concentrationmass relation, which we found to be the first observable to explicitly break the degeneracy as both PNG and IDE have qualitatively the same impact on halo concentrations, namely to suppress concentrations at a given mass with respect to the reference CDM scenario;

The void number function showing mimicking degeneracy for large voids (Mpc) while the degeneracy is broken for smaller void radii;

and The void density profiles for which, similarly to the case of the concentrationmass relation, the mimicking degeneracy is also not observed at all as both individual models predict a lower inner density of cosmic voids compared to CDM.
Therefore, we conclude that measurements of CDM halo and cosmic void internal structural properties, namely halo concentrationmass relation and void density profile would allow us basically to break the degeneracy when combined to any of the other probes that we investigated in this work.
In principle, this degeneracy creates difficulties in identifying the simultaneous presence of IDE and PNG, and in accurately constraining them separately. However, in practice, the degeneracy only arises for values of that are ruled out by current constraints. Nevertheless, our investigation has shown which nonlinear probes could be most useful for improving constraints on IDE and PNG.
Acknowledgements
MH, CG and MB acknowledge support from the Italian Ministry for Education, University and Research (MIUR) through the SIR individual grant SIMCODE, project number RBSI14P4IH. CG acknowledges support from the Italian Ministry of Foreign Affairs and International Cooperation, Directorate General for Country Promotion. The simulations described in this work have been performed on the Sciama High Performance Computing (HPC) cluster at ICG, Portsmouth University. DB acknowledges partial financial support by ASI Grant No. 201624H.0. During the preparation of this work DB was also supported by the Deutsche Forschungsgemeinschaft through the Transregio 33, The Dark Universe and Unidad de Excelencia “María de Maeztu". RM is supported by the South African SKA Project and the National Research Foundation of South Africa (Grant No. 75415), and by the UK Science & Technology Facilities Council (Grant No. ST/N000668/1).
References
 Abel et al. (2012) Abel T., Hahn O., Kaehler R., 2012, Mon. Not. Roy. Astron. Soc., 427, 61
 Abramo & Bertacca (2017) Abramo L. R., Bertacca D., 2017, Phys. Rev., D96, 123535
 Ade et al. (2016a) Ade P. A. R., et al., 2016a, Astron. Astrophys., 594, A13
 Ade et al. (2016b) Ade P. A. R., et al., 2016b, Astron. Astrophys., 594, A14
 Ade et al. (2016c) Ade P. A. R., et al., 2016c, Astron. Astrophys., 594, A17
 Ade et al. (2016d) Ade P. A. R., et al., 2016d, Astron. Astrophys., 594, A20
 Ade et al. (2016e) Ade P. A. R., et al., 2016e, Astron. Astrophys., 594, A24
 Alonso & Ferreira (2015) Alonso D., Ferreira P. G., 2015, Phys. Rev., D92, 063525
 Amendola (2000) Amendola L., 2000, Phys. Rev., D62, 043511
 Amendola (2004) Amendola L., 2004, Phys. Rev., D69, 103524
 Amendola et al. (2008) Amendola L., Baldi M., Wetterich C., 2008, Phys. Rev., D78, 023015
 Aviles & CervantesCota (2011) Aviles A., CervantesCota J. L., 2011, Phys. Rev., D84, 083515
 Baldauf et al. (2010) Baldauf T., Smith R. E., Seljak U., Mandelbaum R., 2010, Phys. Rev., D81, 063531
 Baldauf et al. (2013) Baldauf T., Seljak U., Smith R. E., Hamaus N., Desjacques V., 2013, Phys. Rev., D88, 083507
 Baldi (2011a) Baldi M., 2011a, Mon. Not. Roy. Astron. Soc., 411, 1077
 Baldi (2011b) Baldi M., 2011b, Mon. Not. Roy. Astron. Soc., 414, 116
 Baldi (2012a) Baldi M., 2012a, Mon. Not. Roy. Astron. Soc., 420, 430
 Baldi (2012b) Baldi M., 2012b, Mon. Not. Roy. Astron. Soc., 422, 1028
 Baldi (2012c) Baldi M., 2012c, Annalen Phys., 524, 602
 Baldi (2014) Baldi M., 2014, Phys. Dark Univ., 3, 4
 Baldi & Pettorino (2011) Baldi M., Pettorino V., 2011, Mon. Not. Roy. Astron. Soc., 412, L1
 Baldi & VillaescusaNavarro (2018) Baldi M., VillaescusaNavarro F., 2018, Mon. Not. Roy. Astron. Soc., 473, 3226
 Baldi et al. (2010) Baldi M., Pettorino V., Robbers G., Springel V., 2010, Mon. Not. Roy. Astron. Soc., 403, 1684
 Baldi et al. (2014) Baldi M., VillaescusaNavarro F., Viel M., Puchwein E., Springel V., Moscardini L., 2014, Mon. Not. Roy. Astron. Soc., 440, 75
 Bartolo et al. (2004) Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Phys. Rept., 402, 103
 Baugh et al. (1995) Baugh C. M., Gaztanaga E., Efstathiou G., 1995, Mon. Not. Roy. Astron. Soc., 274, 1049
 Bertolami & Martins (2000) Bertolami O., Martins P. J., 2000, Phys. Rev., D61, 064007
 Bruni et al. (2012) Bruni M., Crittenden R., Koyama K., Maartens R., Pitrou C., Wands D., 2012, Phys. Rev., D85, 041301
 Camera et al. (2015) Camera S., Santos M. G., Maartens R., 2015, Mon. Not. Roy. Astron. Soc., 448, 1035
 Carbone et al. (2010) Carbone C., Mena O., Verde L., 2010, JCAP, 1007, 020
 Clemson et al. (2012) Clemson T., Koyama K., Zhao G.B., Maartens R., Valiviita J., 2012, Phys. Rev., D85, 043007
 Costa et al. (2014) Costa A. A., Xu X.D., Wang B., Ferreira E. G. M., Abdalla E., 2014, Phys. Rev., D89, 103531
 Crooks et al. (2003) Crooks J. L., Dunn J. O., Frampton P. H., Norton H. R., Takahashi T., 2003, Astropart. Phys., 20, 361
 Cui et al. (2012) Cui W., Baldi M., Borgani S., 2012, Mon. Not. Roy. Astron. Soc., 424, 993
 Dalal et al. (2008) Dalal N., Dore O., Huterer D., Shirokov A., 2008, Phys. Rev., D77, 123514
 Desjacques & Seljak (2010) Desjacques V., Seljak U., 2010, Class. Quant. Grav., 27, 124011
 Desjacques et al. (2009) Desjacques V., Seljak U., Iliev I., 2009, Mon. Not. Roy. Astron. Soc., 396, 85
 Desjacques et al. (2018) Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rept., 733, 1
 Despali & Vegetti (2017) Despali G., Vegetti S., 2017, MNRAS, 469, 1997
 Duniya et al. (2015) Duniya D. G. A., Bertacca D., Maartens R., 2015, Phys. Rev., D91, 063530
 Efstathiou & Bond (1999) Efstathiou G., Bond J. R., 1999, Mon. Not. Roy. Astron. Soc., 304, 75
 Fay (2016) Fay S., 2016, Mon. Not. Roy. Astron. Soc., 460, 1863
 Fay et al. (2007) Fay S., Nesseris S., Perivolaropoulos L., 2007, Phys. Rev., D76, 063504
 Fedeli & Moscardini (2010) Fedeli C., Moscardini L., 2010, Mon. Not. Roy. Astron. Soc., 405, 681
 Fonseca et al. (2015) Fonseca J., Camera S., Santos M., Maartens R., 2015, Astrophys. J., 812, L22
 Gao et al. (2004) Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS, 355, 819
 Giannantonio et al. (2014) Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014, Phys. Rev., D89, 023511
 Giocoli et al. (2010) Giocoli C., Tormen G., Sheth R. K., van den Bosch F. C., 2010, MNRAS, 404, 502
 Giocoli et al. (2012) Giocoli C., Tormen G., Sheth R. K., 2012, Mon. Not. Roy. Astron. Soc., 422, 185
 Giocoli et al. (2013) Giocoli C., Marulli F., Baldi M., Moscardini L., Metcalf R. B., 2013, Mon. Not. Roy. Astron. Soc., 434, 2982
 Grossi et al. (2007) Grossi M., Dolag K., Branchini E., Matarrese S., Moscardini L., 2007, Mon. Not. Roy. Astron. Soc., 382, 1261
 Grossi et al. (2009) Grossi M., Verde L., Carbone C., Dolag K., Branchini E., Iannuzzi F., Matarrese S., Moscardini L., 2009, Mon. Not. Roy. Astron. Soc., 398, 321
 Hamaus et al. (2010) Hamaus N., Seljak U., Desjacques V., Smith R. E., Baldauf T., 2010, Phys. Rev., D82, 043515
 Hashim et al. (2014) Hashim M., Bertacca D., Maartens R., 2014, Phys. Rev., D90, 103518
 Heymans et al. (2013) Heymans C., et al., 2013, Mon. Not. Roy. Astron. Soc., 432, 2433
 Hildebrandt et al. (2017) Hildebrandt H., et al., 2017, Mon. Not. Roy. Astron. Soc., 465, 1454
 Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, JCAP, 1204, 027
 Jeong et al. (2012) Jeong D., Schmidt F., Hirata C. M., 2012, Phys. Rev., D85, 023504
 Kamionkowski et al. (2009) Kamionkowski M., Verde L., Jimenez R., 2009, JCAP, 0901, 010
 Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
 Koyama et al. (2009) Koyama K., Maartens R., Song Y.S., 2009, JCAP, 0910, 017
 Kunz (2009) Kunz M., 2009, Phys. Rev., D80, 123001
 Laureijs et al. (2011) Laureijs R., et al., 2011
 Leistedt et al. (2014) Leistedt B., Peiris H. V., Roth N., 2014, Phys. Rev. Lett., 113, 221301
 Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
 Liguori et al. (2010) Liguori M., Sefusatti E., Fergusson J. R., Shellard E. P. S., 2010, Adv. Astron., 2010, 980523
 LoVerde & Smith (2011) LoVerde M., Smith K. M., 2011, JCAP, 1108, 003
 LoVerde et al. (2008) LoVerde M., Miller A., Shandera S., Verde L., 2008, JCAP, 0804, 014
 Lucchin & Matarrese (1985) Lucchin F., Matarrese S., 1985, Phys. Rev., D32, 1316
 Lyth & Wands (2002) Lyth D. H., Wands D., 2002, Phys. Lett., B524, 5
 Maartens et al. (2015) Maartens R., Abdalla F. B., Jarvis M., Santos M. G., 2015, PoS, AASKA14, 016
 Maldacena (2003) Maldacena J. M., 2003, JHEP, 05, 013
 Marulli et al. (2012) Marulli F., Baldi M., Moscardini L., 2012, Mon. Not. Roy. Astron. Soc., 420, 2377
 Matarrese & Verde (2008) Matarrese S., Verde L., 2008, Astrophys. J., 677, L77
 Matarrese et al. (2000) Matarrese S., Verde L., Jimenez R., 2000, Astrophys. J., 541, 10
 Moresco et al. (2014) Moresco M., Marulli F., Baldi M., Moscardini L., Cimatti A., 2014, Mon. Not. Roy. Astron. Soc., 443, 2874
 Moroi & Takahashi (2001) Moroi T., Takahashi T., 2001, Phys. Lett., B522, 215
 Neyrinck (2008) Neyrinck M. C., 2008, Mon. Not. Roy. Astron. Soc., 386, 2101
 Neyrinck & Yang (2013) Neyrinck M. C., Yang L. F., 2013, Mon. Not. Roy. Astron. Soc., 433, 1628
 Padilla (2015) Padilla A., 2015
 Pettorino & Baccigalupi (2008) Pettorino V., Baccigalupi C., 2008, Phys. Rev., D77, 103003
 Pillepich et al. (2010) Pillepich A., Porciani C., Hahn O., 2010, Mon. Not. Roy. Astron. Soc., 402, 191
 Platen et al. (2007) Platen E., van de Weygaert R., Jones B. J. T., 2007, Mon. Not. Roy. Astron. Soc., 380, 551
 Pollina et al. (2016) Pollina G., Baldi M., Marulli F., Moscardini L., 2016, Mon. Not. Roy. Astron. Soc., 455, 3075
 Pollina et al. (2017) Pollina G., Hamaus N., Dolag K., Weller J., Baldi M., Moscardini L., 2017, Mon. Not. Roy. Astron. Soc., 469, 787
 Pourtsidou et al. (2013) Pourtsidou A., Skordis C., Copeland E. J., 2013, Phys. Rev., D88, 083505
 Raccanelli et al. (2014) Raccanelli A., Bertacca D., Doré O., Maartens R., 2014, JCAP, 1408, 022
 Raccanelli et al. (2018) Raccanelli A., Bertacca D., Jeong D., Neyrinck M. C., Szalay A. S., 2018, Phys. Dark Univ., 19, 109
 RenauxPetel (2015) RenauxPetel S., 2015, Comptes Rendus Physique, 16, 969
 Ross et al. (2013) Ross A. J., et al., 2013, Mon. Not. Roy. Astron. Soc., 428, 1116
 Salvatelli et al. (2013) Salvatelli V., Marchini A., LopezHonorez L., Mena O., 2013, Phys. Rev., D88, 023531
 Salvatelli et al. (2014) Salvatelli V., Said N., Bruni M., Melchiorri A., Wands D., 2014, Phys. Rev. Lett., 113, 181301
 Sartoris et al. (2016) Sartoris B., et al., 2016, Mon. Not. Roy. Astron. Soc., 459, 1764
 Scoccimarro et al. (2012) Scoccimarro R., Hui L., Manera M., Chan K. C., 2012, Phys. Rev., D85, 083002
 Setare & Mohammadipour (2013) Setare M. R., Mohammadipour N., 2013, JCAP, 1301, 015
 Simpson et al. (2016) Simpson F., et al., 2016, Phys. Rev., D93, 023525
 Smith et al. (2007) Smith R. E., Scoccimarro R., Sheth R. K., 2007, Phys. Rev., D75, 063512
 Springel (2005) Springel V., 2005, Mon. Not. Roy. Astron. Soc., 364, 1105
 Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, Mon. Not. Roy. Astron. Soc., 328, 726
 Springel et al. (2008) Springel V., et al., 2008, Mon. Not. Roy. Astron. Soc., 391, 1685
 Sutter et al. (2014) Sutter P. M., Elahi P., Falck B., Onions J., Hamaus N., Knebe A., Srisawat C., Schneider A., 2014, Mon. Not. Roy. Astron. Soc., 445, 1235
 Sutter et al. (2015) Sutter P. M., et al., 2015, Astronomy and Computing, 9, 1
 Tereno et al. (2005) Tereno I., Dore O., Van Waerbeke L., Mellier Y., 2005, Astron. Astrophys., 429, 383
 Väliviita & Palmgren (2015) Väliviita J., Palmgren E., 2015, JCAP, 1507, 015
 Vikhlinin et al. (2009) Vikhlinin A., et al., 2009, Astrophys. J., 692, 1060
 VillaescusaNavarro et al. (2014) VillaescusaNavarro F., Marulli F., Viel M., Branchini E., Castorina E., Sefusatti E., Saito S., 2014, JCAP, 1403, 011
 Wagner & Verde (2012) Wagner C., Verde L., 2012, JCAP, 1203, 002
 Wagner et al. (2010) Wagner C., Verde L., Boubekeur L., 2010, JCAP, 1010, 022
 Wei & Zhang (2008) Wei H., Zhang S. N., 2008, Phys. Rev., D78, 023011
 Wei et al. (2013) Wei H., Liu J., Chen Z.C., Yan X.P., 2013, Phys. Rev., D88, 043510
 Weinberg (1989) Weinberg S., 1989, Rev. Mod. Phys., 61, 1
 Wetterich (1988) Wetterich C., 1988, Nucl. Phys., B302, 668
 Wetterich (1995) Wetterich C., 1995, Astron. Astrophys., 301, 321
 Zeldovich (1970) Zeldovich Ya. B., 1970, Astron. Astrophys., 5, 84
 Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Boerner G., 2009, Astrophys. J., 707, 354