Impact of cosmological signatures in two-point statistics beyond the linear regime

We are interested in the cosmological consequences, and particularly imprints on the matter power spectrum, generated by an extra energy density $\rho_{ex}$, beyond the standard $\Lambda$CDM, for $a\leq a_{c}$, which dilutes rapidly after a transition takes place at the scale factor $a_{c}$, and with a mode $k_{c}=a_{c}H(a_{c})$ crossing the horizon at such a time. For definiteness, we will assume that such component, henceforth called Rapid Diluted Energy Density (RDED) de la Macorra & et al. (2020), tracks the leading background energy density, in either radiation, matter or DE domination epochs and we study these different cases.

Let us consider a $\Lambda$CDM cosmology with and additional energy density $\rho_{ex}$, and cosmic abundance

where $H^{2}_{sm}$ and $H^{2}_{smx}$ are the Hubble parameters for the model with and without the RDED component, respectively.We consider the case for which $\Omega_{ex}$ is a constant for $a\leq a_{c}$, i.e. $\rho_{ex}$ tracks the leading background energy density with $w_{ex}=w_{sm}$ at the time of the transition ($w_{sm}=1/3$ in radiation and $w_{sm}=0$ in matter domination) and a rapid dilution of $\rho_{ex}$ takes place at $a_{c}$. At this time the EoS $w_{ex}\equiv p_{ex}/\rho_{ex}$ suffers a transition from $w_{c}\equiv w_{ex}(a_{c})=w_{sm}$ for $a\leq a_{c}$ to $w_{f}\equiv w(a_{f})>w_{c}$ for $a_{f}>a_{c}$. The value of $\Delta w\equiv w_{f}-w_{c}>0$ and the width $\Delta a\equiv(a_{f}-a_{c})/a_{c}$ set the steepness of the transition of the energy density of $\rho_{ex}(a_{c})$, i.e. how fast it dilutes. For $a>a_{c}$ we have

with $\Delta w\equiv w_{f}-w_{c}$ and for $w_{c}=1/3,w_{f}=1$ we have $\Delta w=2/3$ the amount of extra energy density $\Omega_{ex}$ will decrease. Notice that we would expect a massive particle to go from being relativistic for $T/m\gg 1$ with $w_{c}=1/3$ to non-relativistic $T/m\ll 1$ with $w_{f}=0$ giving a negative $\Delta w=w_{f}-w_{c}=-1/3$ and an increasing $\Omega_{ex}$.

This property is not accomplished by an RDED component, hence it is beyond the standard model. However, it can be implemented using a scalar field $\phi$, with energy density $\rho_{\phi}=E_{k}+V$ and pressure $p_{\phi}=E_{k}-V$, where $E_{k}$ is the kinetic term and $V$ the potential energy (for a review see Copeland et al. 2006). The evolution of the homogeneous scalar field $\phi(t)$ is given by the Klein Gordon equation

Widely used scalar potentials are either exponential or inverse power laws Steinhardt et al. (1999); de la Macorra & Stephan-Otto (2002), e.g.

with $n$ a dimensionless constant while $\Lambda,\alpha^{-1}$ have mass dimension. These potentials $V$ have different behaviours depending on the value of the parameters and the initial conditions of the scalar field. The energy density $\rho_{\phi}$ can track the background evolution (named tracker fields Steinhardt et al. (1999); Armendariz-Picon et al. (2001)) either in radiation or matter domination with $w_{i}=0$ or $w_{i}=1/3$ and later leap to $w_{f}=1$ (for $V\gg E_{k}$) for a long period of time, diluting $\rho_{\phi}$ and generating a rapid dilution $\rho_{ex}$ situation.

For example, for inverse power low (IPL) scalar field potential $V=\Lambda^{4+n}\phi^{-n}$ de la Macorra & Stephan-Otto (2002) the evolution of $\phi$ goes through an epoch of kinetic-term domination with $w=1$, before transitioning to potential-term domination region with $w=-1$ to finally grow to a value $w\approx-0.9$ close to present time for $n\sim 1$ Chevallier & Polarski (2001); Linder (2003) and de la Macorra & Stephan-Otto (2002); de la Macorra & Piccinelli (2000) while BDE model has $n=2/3$ de la Macorra & Almaraz (2018). On the other hand a value of $n>5$ has a tracking behaviour, i.e. the scalar field abundance evolves as the dominant background energy density, but these models have an EoS $w>-0.7$ and are ruled out by current cosmological observations Planck Collaboration (2018). Models with $n<2$ are viable as DE models.

An example of an IPL model is BDE, which is derived from particle physics where the fundamental particles are massless and contained in a gauge group similar to the strong QCD force, however due to the underlying dynamics they develop a IPL potential with $n=2/3$ at the transition with a scale factor $a_{c}\sim 10^{-6}$. In this case the EoS is $w_{c}=1/3$ at early times and at $a_{c}$ it leads to $w_{f}=1$ to a later evolve to $w\sim-1$ at $z\sim 1000$ and finally ending up at a value of $w\sim-0.93$ at present time. The complete evolution of $w(\phi)$ is a consequence of the dynamics of the equation of motion and is not set by hand. The BDE example satisfies the condition ($\Delta w=w_{f}-w_{c}>0$) of RDED twice. The first is during radiation domination at $a_{c}\sim 10^{-6}$ with a $k_{c}=0.94$ and the second takes place close to present time ($z\sim 0.28$). The transitions we have been studying here are clearly beyond the standard model of particle physics but are expected in phase transitions for dynamical scalar fields and are motivated to explain the nature of dark energy.

Matter perturbations are affected by RDED, leaving distinctive features in cosmological observables such as the matter perturbations. The amplitude of these perturbations is increased in RDED cosmology compared to a standard $\Lambda$CDM cosmology, and a bump feature is produced in the ratio of power spectrum $P_{smx}/P_{sm}$. The bump is generated for modes $k>k_{c}$, with $k_{c}\equiv a_{c}H(a_{c})$, entering the horizon at a scale factor $a<a_{c}$.

We will sketch how this bump is generated and we will study in detail the properties of these bumps in both the linear and non-linear approach. Let us explore here the linear regime of matter density fluctuations. In this case, the evolution of CDM perturbations obey the equation of motion:

where the sum runs over all the fluids with sound speed $c_{s,i}^{2}=\delta P_{i}/\delta\rho_{i}$ and the prime denotes derivatives with respect to conformal time, with $\mathcal{H}=a^{\prime}/a$. We see that the $\rho_{ex}$ has a twofold effect on matter perturbations by modifying first the expansion rate $\mathcal{H}$, and second by the adding of an extra source term proportional to $\Omega_{ex}\delta_{ex}$.

Initially, for modes outside the horizon, the matter perturbations do not evolve over time, leaving the ratio $Q_{m}\equiv\delta_{m}^{\textrm{smx}}/\delta_{m}^{\textrm{sm}}$ fixed. However, matter perturbations in the model with $\rho_{ex}$ are initially suppressed compared to a $\Lambda$CDM model because the initial amplitude of the matter perturbations, through the gravitational potential $\Psi$, depends on the fraction of relativistic particles $R_{\nu}=\rho_{eff}/(\rho_{eff}+\rho_{\gamma})$, with $\rho_{\gamma}$ the density of photons, $\rho_{eff}^{\textrm{sm}}=\rho_{\nu}$ and $\rho_{eff}^{\textrm{smx}}=\rho_{\nu}+\rho_{ex}$ giving an initial suppression factor $Q_{ini}=(1+(4/15)R_{\nu}^{\textrm{sm}})/(1+(4/15)R_{\nu}^{\textrm{smx}})$, which depends solely on $\rho_{ex}$.

Modes start evolving once they cross the horizon $a_{h}$, defined implicitly in terms of the Hubble radius by $k=a_{h}H(a_{h})$. However, there is a marked difference depending on whether they cross before $a_{c}$, with a corresponding mode $k_{c}\equiv a_{c}H(a_{c})$, or after the extra energy density $\rho_{ex}$ has diluted. Small modes $k>k_{c}$ are further suppressed with respect to $\Lambda$CDM because the crossing time in this case is delayed by the presence of the $\rho_{ex}$ . Comparing the same $k_{h}=(a_{h}H(a_{h}))\big{|}_{smx}=(a_{h}H(a_{h}))\big{|}_{sm}$ we have

and therefore small modes $k>k_{c}$ cross the horizon earlier in $\Lambda$CDM than in $\Lambda$CDMex. This is reflected in an early suppression in $\Delta\delta_{m}=\delta\rho_{m}^{smx}/\delta\rho_{m}^{sm}$, supplemented by a smaller amplitude at horizon crossing. However, after the initial suppression at horizon crossing, matter perturbations in the $smx$ model have a higher growing rate than in the $\Lambda$CDM model that not only compensates but also reverses the initial suppression (de la Macorra & et al., 2020).

In radiation domination the matter perturbations evolve as $\delta_{m}\propto\ln a$ with the growth function $f=d\ln\delta_{m}/d\ln a\propto 1/\delta_{m}$, and since $\delta_{m}(smx)<\delta_{m}(sm)$ and $f(smx)>f(sm)$ and $Q_{m}$ increases. This increase is boosted by the rapid dilution of $\rho_{ex}$, where the EoS leaps abruptly from $w_{i}=1/3$ to $w_{f}=1$ at $a_{c}$ affecting only the modes crossing the horizon before $a_{c}$, i.e. $k\geq k_{c}$. This is the characteristic signature of RDED and was presented in Almaraz & de la Macorra (2019) and Jaber-Bravo et al. (2020). However, to fully asses the growth of structure in these RDED cosmologies we must go beyond the linear regime, which is the main goal of this work.

To gain physical intuition on how the rapid diluted energy density affects matter-density fluctuations well inside the radiation dominated epoch we can analyze a simplified version of the equations where the gravitational potentials are subdominant: We can estimate this increase in radiation domination, where the amplitude $\delta_{m}=\delta\rho_{m}/\rho_{m}$ has a logarithmic growth

with $a_{h}$ and $H_{h}$ the scale factor and Hubble parameter at horizon crossing. We compare the evolution of the same mode $k^{smx}=k^{sm}$ crossing the horizon before $a_{c}$ and using equations (7) and (6) we find the ratio $\Delta\delta_{m}=\delta_{m}^{smx}/\delta_{m}^{sm}=(\delta_{mi}^{smx}/\delta_{mi}^{sm})([\ln(a/a_{h}^{smx})+1/2]/[\ln(a/a_{h}^{sm})+1/2])$. The evolution of $\Delta\delta_{m}$ after the RDED transition takes place at $a_{c}$ is given by de la Macorra & et al. (2020)

where we have taken into account in $\Lambda$CDMx the contribution of $\rho_{ex}$ in $H_{+}^{smx}(a_{c})$ for $a<a_{c}$ and $\rho_{ex}=0$ in $H_{-}^{smx}(a_{c})$ for $a>a_{c}$, giving the ratio $H_{-}^{smx}/H_{+}^{smx}=a_{h}^{\mathrm{smx}}/a_{h}^{\mathrm{sm}}=\sqrt{1-\Omega_{ex}}$ (c.f. equation (6)). The matter power spectrum therefore shows a bump for modes $\Delta\delta_{m}$ for $a\gg a_{c}$ for modes $k>k_{c}$ entering the horizon before $a_{c}$ when we compare RDED and $\Lambda$CDM cosmologies.

We ran $12$ $N$-body simulations using the cosmological simulation code gadget-2 (Springel, 2005), one each for the cosmologies detailed in Table 1. We assume a background $\Lambda$CDM cosmology with total matter density $\Omega_{m}=0.3$, baryon density $\Omega_{b}=0.05$, dark energy density $\Omega_{\Lambda}=0.7$, amplitude of the matter power spectrum $\sigma_{8}=0.8$, spectral index $n_{s}=0.96$, and dimensionless Hubble constant $h=0.7$. Initial conditions were generated at $z=99$ using ngenic,³³3https://www.h-its.org/2014/11/05/ngenic-code/which implements the Zeldovich approximation to displace particles from an initially Cartesian grid and assigns them velocities based on a ballistic trajectory. We chose to run simulations with different box sizes for the different $k_{T}$ values, to ensure that there was always a good sampling of modes around $k_{T}$. The box sizes of the simulations are $L_{box}=256$, $512$ and $1024\,h^{-1}\mathrm{Mpc}$, for $k_{T}=0.05$, $0.1$ and $1\,h\mathrm{Mpc}^{-1}$ respectively. Each simulation uses $512^{3}$ particles to approximate the density field, which is quite modest compared to modern simulation standards. One may worry that measurements from these simulations would be systematically biased as a result of the low mass resolution. However, we checked for convergence with respect to low-resolution $256^{3}$ particle simulations and found that our results for the power spectrum response were only significantly affected ($>1$ per cent) on scales smaller than half of the particle Nyquist frequency, although there were some noticeable, sub per-cent differences for scales as large as one tenth of the particle Nyquist frequency. This is different to a general suppression that affects the power spectrum (rather than the response) when using low-resolution simulations (Mead et al., 2015), and must be because some of the bias in power that arises when using a low particle number cancels when constructing a response; so the same (or a similar) bias must occur in the cosmologies of both the numerator and denominator that make up the response.

We look at predictions for the non-linear matter power spectrum of our cosmological models using the hmcode (Mead et al., 2015) model. hmcode is an augmented version of the traditional halo-model calculation for the non-linear power spectrum (e.g., Seljak, 2000; Peacock & Smith, 2000; Cooray & Sheth, 2002) where the augmentations account for defects and replace ingredients that are missing from the standard halo-model calculation. This improves the accuracy of the calculation from $\sim 30$ per cent to $\sim 5$ per cent for a wide range of cosmologies and redshifts. hmcode takes as input the linear power spectrum of the cosmology in question, and then uses some information about the background cosmological parameters and power-spectrum shape and amplitude in order to make its predictions. Although hmcode was not calibrated on the bump cosmologies investigated in this paper, the grounding of the model in physical reality means that we can hope that it will make reasonable predictions. hmcode models the power spectrum as a sum of (almost) linear theory and a one-halo term, which accounts for small-scale, deeply non-linear power under the assumption that all such power originates from the clustering within haloes. The bump cosmologies will therefore affect the hmcode predictions in two ways. The first, is trivially that the hmcode prediction at large scales is essentially linear theory and because linear theory contains the bump then so will hmcode. The second, is that within hmcode the one-halo power is determined by the halo mass function, which itself depends on the linear power spectrum via the variance in the power spectrum as a function of scale. The bump will therefore affect the halo mass function and we expect that it will boost the predicted numbers of haloes in certain mass ranges.

Before running hmcode, we can make the prediction that it should generally boost power in the one-halo term for a bump cosmology compared to a cosmological model that lacks a bump. The traditional halo model calculation has a problem in the transition region between the two- and one-halo terms ($k\sim 0.05\,h\mathrm{Mpc}^{-1}$ at $z=0$), and generically underestimates the true non-linear power spectrum in this region; this probably arises due to an improper treatment of non-linear halo bias (Smith et al., 2007). The solution to this problem in hmcode is a smoothing of the transition region. Based on this discussion, we could predict that the hmcode predictions for the bump cosmologies may be better in the deeply one-halo regime ($k>1\,h\mathrm{Mpc}^{-1}$), and that they may be less impressive in the transition region. In the linear region they should be perfect, given that the hmcode prediction is identical to linear theory at large scales. In the quasi-linear regime ($k\sim 0.1\,h\mathrm{Mpc}^{-1}$ at $z=0$) we would expect perturbation theory to perform better than hmcode since the latter lacks any formal grounding in analytical perturbation theory.

In this paper we consider two variations of perturbation theory, the first, ‘standard’ perturbation theory is constructed in Eulerian space, where we consider the evolution of the density field at fixed positions. We consider the phenomenological, but physically well motivated, redshift-space models of Scoccimarro (2004) and Taruya et al. (2010) to describe the redshift-space power spectrum, which improve the description of redshift-space distortion (RSD) effects on the power spectrum compared to Kaiser linear theory Kaiser (1987). In the following we shall refer to these models as Sc04 and TNS, respectively.

The anisotropic redshift-space clustering originates from the peculiar velocities $\mathbf{v}$ of matter, or more generally any tracer of it, such that an object located at a real space position $\mathbf{r}$ is observed to be located at an apparent redshift-space position $\mathbf{s}$. The relation between coordinates system is inferred via the Doppler effect to be $\mathbf{s}=\mathbf{r}+\hat{\mathbf{n}}v_{\parallel}(aH)^{-1}$, where $\hat{\mathbf{n}}$ is a the line-of-sight direction of the point-process sample, and $v_{\parallel}=\mathbf{v}\cdot\hat{\mathbf{n}}$. That is, we are using the plane-parallel approximation on which the observer is located at a distant position of the sample of objects over which we perform the statistics. The redshift-space power spectrum, $\langle|\delta^{s}(\mathbf{k})|^{2}\rangle$ is given by

where $\mathbf{r}=\mathbf{x}-\mathbf{x^{\prime}}$ and $\Delta v_{\parallel}=v_{\parallel}(\mathbf{x})-v_{\parallel}(\mathbf{x^{\prime}})$ and $\mu=\hat{\mathbf{n}}\cdot\hat{\mathbf{k}}$ is the angle between the wave vector and the line-of-sight direction.

The RSD correction at linear order, known as the Kaiser formula, is given by $\delta^{s}_{L}(\mathbf{k})=(1+f\mu^{2})\delta_{L}(\mathbf{k})$, where $f=d\log D_{+}(a)/d\log a(t)$ and $D_{+}(t)$ the linear growth function. The redshift-space power spectrum at linear order becomes

To move beyond linear theory it is common to define the dimensionless velocity divergence $\theta=-\nabla{\bf\cdot v}/(aHf)$, such that at linear order we have $\theta=\delta$.

The exponential oscillatory factor inside the correlator in equation (3.3) is due to virialized, non-coherent random motions of dark-matter particles along the line-of-sight direction, hence it is in essence non-perturbative. In (Scoccimarro, 2004) this factor is replaced by a phenomenological damping function that accounts for the velocity dispersion $\sigma^{2}_{v}=\langle\theta^{2}\rangle$. By rotational symmetry around the line-of-sight direction one obtains a simple prescription

where $P_{\delta\delta}$, $P_{\theta\theta}$ and $P_{\delta\theta}$ are respectively the non-linear matter density, velocity divergence, and density-velocity divergence power-spectra. In linear theory equation (3.3) reduces to the Kaiser power spectrum (equation 12) times the damping factor. Several other functional forms for this damping factor have been used in the literature. In this work, we opt for the most common – Gaussian damping. The velocity dispersion is physically motivated to be given by PT as

but due to its non-perturbative origin it is commonly replaced by a free parameter $\sigma^{2}_{\text{FoG}}$, especially for parameter estimation. However in this work, we regard $\sigma^{2}_{v}$ as the linear velocity dispersion, which is obtained from the above equation by replacing $P_{\theta\theta}$ by its linear value $P_{L}$. This prescription has shown to give good results when comparing theory to simulations (e.g., Taruya et al., 2010).

The TNS formalism, on the other hand, expands in cumulants the correlator in equation (3.3), and then replaces a residual exponential factor of the form $\exp\big{[}\langle e^{ik\mu\Delta v_{\parallel}/(aH)}\rangle_{c}\big{]}$, by a position independent, phenomenological damping factor that can be brought out of the integral. The standard formula for TNS is

where the new correction terms, $A(k,\mu)$ and $B(k,\mu)$, are given by

with the bispectrum

and

Hence, the TNS model partially accounts for the interaction between the Kaiser effect and the non-linear random motion of particles, reflected in the two extra functions $A$ and $B$. On the other hand, the Sc04 model factorizes the linear and Finger-of-God effects, so each can be treated separately.

The real-space matter power spectrum is given by $P_{\delta\delta}(k)=\langle|\delta(\mathbf{k})|^{2}\rangle$, which can be obtained from both of the above redshift-space formalisms by considering only the perpendicular to line-of-sight components (i.e., $\mu=0$). After that, one can simply use rotational symmetry to obtain the 1-loop SPT power spectrum

where $P_{22}$ and $P_{13}$ are the usual 1-loop corrections (e.g., Bernardeau et al., 2002). The linear power spectrum is the dominant term at large scales. However, at smaller scales, the loop corrections contribute with similar magnitude to the linear power, therefore significantly contributing to the total power spectrum.

In contrast to the Eulerian scheme, in a Lagrangian picture we follow the trajectories of individual particles with initial position $\mathbf{q}$ and current position $\mathbf{x}$. The Lagrangian coordinates $\mathbf{q}$ are related to the Eulerian coordinates $\mathbf{x}$ through the coordinates transformation

where $\Psi(\mathbf{q},t)$ is the Lagrangian displacement and $\Psi(\mathbf{q},t=t_{ini})=0$, meaning that at a sufficiently early time $t_{ini}$ both coordinates system coincide. Furthermore, mass conservation implies the relation between Lagrangian displacements and overdensities (Taylor & Hamilton, 1996)

from which the LPT correlation function is obtained as

where $\Delta^{i}=\Psi^{i}(\mathbf{q}_{2})-\Psi^{i}(\mathbf{q}_{1})$ are the Lagrangian displacement differences at two positions $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ separated by a distance $q=|\mathbf{q}_{2}-\mathbf{q}_{1}|$. The idea behind CLPT is to expand the correlator in eq. (23) in cumulants and thereafter to expand out of the exponential all but linear terms in the Lagrangian displacement, such that the k-integral can be performed analytically using standard Gaussian integration techniques (see Carlson et al. (2012)), by which we arrive at

with cumulants $A_{ij}=\langle\Delta_{i}\Delta_{j}\rangle_{c}$ and $W_{ijk}=\langle\Delta_{i}\Delta_{j}\Delta_{k}\rangle_{c}$, and tensors $G_{ij}=A^{-1}_{L\,ij}-g_{i}g_{j}$, $\Gamma_{ijk}=A^{-1}_{L\,ij}g_{k}+A^{-1}_{L\,jk}g_{i}+A^{-1}_{L\,ki}g_{j}-g_{i}g_{j}g_{k}$, and $g_{ij}=A_{L\,ij}^{-1}(r_{j}-q_{j})$. The label ‘$L$’ in the $A$ function denotes the linear piece and ‘loop’ the pure 1-loop piece, such that $A_{ij}=A^{L}_{ij}+A^{\text{loop}}_{ij}$. Notice that the ‘1’ in the squared brackets of the above equation corresponds to the Zeldovich correlation function and the other two terms yield the next-to-leading order, 1-loop contributions.

To solve numerically the CLPT correlation function of equation (3.4) we use the code mgpt (Aviles et al., 2018),⁴⁴4https://github.com/cosmoinin/MGPT which also accounts for the exact kernels for a $\Lambda$CDM background cosmology, instead of the most commonly used Einstein-de Sitter ($\Omega_{\mathrm{m}}=1$) kernels. Using the correct kernels can make a small, but noticeable, difference to the power spectrum – as much as $0.8$ per cent at quasi-linear scales at $z=0$.

The power spectra data are extracted from the simulations using a cloud-in-cell mass-assignment scheme, on a grid with resolution of $N_{grid}=512^{3}$ cells. They are binned in $100$ evenly log-spaced $k$-points over a range $[k_{min},k_{\mathrm{Ny}}]$, where $k_{\mathrm{Ny}}=N_{grid}\pi/L$ is the Nyquist frequency, and $k_{min}=2\pi/L$. $L$ is the size of the box, which is given in Table 1 for the different simulations. Traditionally power spectra measured in this way are considered to be accurate up to half of the Nyquist frequency, with the power at smaller scales potentially corrupted by aliasing from smaller scales. In our plots we show all measured scales, up to $k_{\mathrm{Ny}}$, and the reader should keep this in mind.

Figs. 1, 2, and 3 show the responses using the bump cosmology from equation (9) located at the scales $k_{T}=0.05$, $0.1$, and, $1\,h\mathrm{Mpc}^{-1}$, respectively. The non-linear power is computed using our different methods and then divided by their counterparts in $\Lambda$CDM to create the response function. This analysis compares how the bump cosmologies are modified by non-linearities within the different approaches.

As expected, at higher redshifts non-linear effects are smaller and the responses of all methods are very similar. In the $\sigma=1$ case of Fig. 1 as $z$ decreases we see two things in the simulated measurements: First, at the small-scale edge of the bump we see the power grow above the bump, an effect that is captured extremely well by SPT and less well by hmcode. Second, we see the generation of a second bump at much smaller scales, with a peak $k\sim 1\,h\mathrm{Mpc}^{-1}$. This second bump is well modelled by hmcode but not captured at all by SPT. The failure of SPT at these scales (which even predicts a decrease in the response at higher $z$) is not surprising given that these scales are far beyond its reach. In Fig. 1, corresponding to $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$, this non-linear ‘second bump’ peaks at $10$ times smaller scale, $k\approx 0.5\,h\mathrm{Mpc}^{-1}$, and reaches a maximum at $z=0$, peaking in amplitude at $\sim 10$, $2$, and $1$ per cent for the $\sigma=1$, $0.3$, and $0.1$ cases respectively. For $k_{T}=0.1$, shown in Fig. 2, this second non-linear bump is even clearer, peaking at $k\approx 1\,h\mathrm{Mpc}^{-1}$ for $z\simeq 0$, but at smaller scales for higher redshifts. For some of our cosmologies, this second bump is even larger than the primordial bump. For example, it is larger for $\sigma=1$ at both redshifts $z=0.0$ and $z=0.5$, contributing $\approx 18$ and $16$ per cent to the whole response respectively. We note that as the width of primordial bump decreases, the second, non-linear bump amplitude decreases. In all cases, the location and amplitude of the non-linear bump as seen in the simulation response is in remarkable agreement with the predictions from hmcode. For the cases of $k_{T}=0.05$ and $0.1\,h\mathrm{Mpc}^{-1}$, shown in Figs. 1 and 2, hmcode and the simulated data provide similar results at lower redshifts and in the non-linear regime. Conversely, 1-loop SPT gives results closer to those obtained from simulations at higher redshifts and in the mildly non-linear regime ($k\lesssim 0.2\,h\mathrm{Mpc}^{-1}$).

The non-linear bumps in the response functions are a consequence of the interaction between the primordial bump and the one-halo term, being highly enhanced for the wider bumps simply because these provide a greater enhancement of linear power. Physically, this can be understood within the halo model (and therefore within hmcode) as follows: the one-halo term is given by the integral of the halo mass function multiplied by the squared Fourier-space halo profile. The halo mass function itself is related to the standard-deviation in the density field when smoothed over the Lagrangian radius of the halo, $\sigma_{R}$. Our bump cosmology increases the power over a certain range of scales, such that $\sigma_{R}$ will also increase, and therefore so will the mass function. Hence, one effect of the linear bump is to accelerate halo formation relative to a cosmology with no bump. A different way to think about this is via a Press & Schechter (1974) type argument, where the increased amplitude of some modes, given here by our bump cosmology, is helping more small scale fluctuations to cross over the critical threshold to collapse. This occurs even when the bump is at very linear scales as long as the width is sufficiently large, because those long wavelength modes exist underneath the smaller-scale fluctuations enhancing the collapse to form the actual haloes. We remark that this is a highly non-linear effect, such that the second bump is not well captured by PT, where the main non-linear effect is the spreading and enhancement of the primordial bump.

In Fig. 3 we show the $k_{T}=1\,h\mathrm{Mpc}^{-1}$ bump cosmologies. The simulated data show that the primordial bump become erased by the non-linear evolution. Such effect is well captured by the hmcode– PT simply does not capture the non-linearities of the evolution of the bump since it is far out of the reach of its regime of validity, although at the highest redshifts it works moderately well. Interesting is the non presence of the second bump, because the primordial bump and the 1-halo term are located about the same scale.

Figs. 1, 2 and 3 demonstrate that perturbation theory provides an accurate model for the power spectrum response at large scales, with the accuracy extending to smaller scales at higher redshift. On the other hand, hmcode provides a reasonable (though less perfect) model for the response at the smaller scales, those typically associated with halo formation. It does not take a great leap of the imagination to consider combining these two approaches to provide an accurate model for the response that would be valid across a wider range of scales. This approach has been explored in the literature previously (e.g., Mohammed & Seljak, 2014; Seljak & Vlah, 2015; Philcox et al., 2020), but has never been applied to the bump cosmology specifically. It may be possible to directly add perturbation theory to the halo-model, perhaps replacing the two-halo term with some perturbative expansion, or else to use the response from perturbation theory at one extreme and hmcode at the other, interpolating between the two around some ‘non-linear’ wavenumber, perhaps defined via $\sigma(a/k_{\mathrm{nl}})=\alpha$ where $\sigma$ is the standard deviation in the density field when smoothed on a given physical scale and $a$ and $\alpha$ are fitted constants of order unity.

In this Section, we study the responses of the redshift-space monopole and quadrupole power spectra using the non-linear models Sc04 and TNS. For comparison we also use Kaiser linear theory. The multipole power spectra have been extracted from the $N$-body simulations using a triangular-shaped-cloud mass-assignment function, on a grid with $N_{grid}=512^{3}$ cells, using the n-bodykit package⁵⁵5https://nbodykit.readthedocs.ioof Hand et al. (2018). The response functions are defined for each multipole, such that $R_{\ell}(k)=P_{\ell}^{\text{bump}}(k)/P_{\ell}^{\text{$\Lambda$CDM}}(k)$.

In the top panel of Fig. 4 we show the response functions in the monopole $P_{0}$ at redshifts $z=0$ and $z=1$ for the bump cosmologies with $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$. We notice, as expected, that non-linear models perform closer to the simulated data than the linear theory in the quasi-linear regime, in the sense that the responses are closer to those from $N$-body simulations. However, for $k>0.2\,h\mathrm{Mpc}^{-1}$ both theoretical models make deviations from the simulated data. In the bottom panel of Fig. 4 we show the response functions for the bump cosmologies located at $k_{T}=0.1\,h\mathrm{Mpc}^{-1}$. The behaviour of the PT methods is qualitatively similar to the case of $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$, but the non-linear features of the bump are less well captured when comparing to the simulated data. This must be since the scale of the bump is located at the edge of the perturbative regime of validity. This has the consequence that the linear theory response provides a better match to the simulations for some non-linear scales, obviously the good performance of linear theory here is only a lucky coincidence. As in real space, the simulated data show the appearance of the non-linear second bump at $\approx 1\,h\mathrm{Mpc}^{-1}$, which is not present in the $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$ case, even though in real space this second non-linear bump was visible in both cases. We suggest that the reason for this is the damping along the line-of-sight direction of the redshift-space power spectrum, which ultimately comes from the highly oscillatory behavior of the correlator inside the integral of equation (3.3) at large $k$. This redshift-space effect damps the multipoles in all bump cosmologies, but since the second (real space) bump is larger for the case of $k_{T}=0.1\,h\mathrm{Mpc}^{-1}$, it can overcome the damping and it still shows up in the redshift-space responses.

We also show, in Fig. 5, the quadrupole redshift-space power spectra for $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$ and $k_{T}=0.1\,h\mathrm{Mpc}^{-1}$. Although the simulated quadrupole measurements are noisier than for the monopole, we note a similar trend to that predicted by the analytical approaches, most obviously for the case of $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$ at $z=1$. We also observe that the non-linear, second bumps do not appear in the quadrupole. We suggest that this is because this multipole gives maximum weight to the line-of-sight direction where the damping effect occurs, while the monopole gives equal weight to all directions.

In this sub-section we compare the measured real-space correlation function responses, defined as $\xi_{\text{bump}}(r)/\xi_{\text{$\Lambda$CDM}}(r)$, from both the simulated data and the analytical correlation function calculated according to the CLPT method. All simulated correlation functions have been measured by employing the n-bodykit code, using $60$ linearly spaced bins in the range $1$–$121\,h^{-1}\mathrm{Mpc}$ for the bump cosmologies at $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$, $30$ bins in the range $1$–$61\,h^{-1}\mathrm{Mpc}$ for models at $k_{T}=0.01$, and $15$ bins in the range $1$–$31\,h^{-1}\mathrm{Mpc}$ for models at $k_{T}=1.0\,h\mathrm{Mpc}^{-1}$.

Fig. 6 shows the response functions for the bump at $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$ for four different redshifts $z=2$, $1$, $0.5$, and $0$ and the bump cosmologies with widths $\sigma=1$, $0.3$, and $0.1$. We find a dip in the response around $85\,h^{-1}\mathrm{Mpc}$, which is due to a long wavelength modulation of the whole correlation function because of the localized bump feature in the corresponding power spectrum. By considering equation (9), the linear correlation function response becomes

with

where $j_{0}$ is the zero-order spherical Bessel function. Hence for scales $r\gg k_{T}^{-1}$ we have that $\tilde{F}(r)\rightarrow 0$ because of large cancellations provided by $j_{0}(kr)$; while as $r\rightarrow 0$, $\tilde{F}(r)\rightarrow(2\pi^{2})^{-1}\int dkk^{2}F(k)$, a constant that when computed turns out to be much smaller than the correlation function at those small scales; such that the linear response goes to unity at both large scales and small scales limits. On the other hand, $k^{2}F(k)$ reaches its maximum at scales $k_{max}>k_{T}$, the larger the width $\sigma$, the bigger are $k_{max}$ and the maximum values $k_{max}^{2}F(k_{max})$; however for larger widths, the function $k^{2}F(k)$ overlaps with more oscillations and $\tilde{F}(r)$ is suppressed. The inverse transformed bump feature, $\tilde{F}(r)$, turns out to be a sinc-like function in real space with amplitude proportional to the $\sigma$ and first trough around $r_{T}\lesssim(3/2)\pi/k_{T}$. Moreover, the aforementioned effects compete and this becomes more pronounced for the width $\sigma=0.3$ case. Coincidentally, for the case $k_{T}=0.05$, this trough is located at about the same scale as the dip, at about $87\,h^{-1}\mathrm{Mpc}$, in the correlation function, and both reinforce each other to produce the large dip observed in Fig. 6.

More generally, for the different choices of $k_{T}$ the changes in the correlation-function response will follow a similar pattern, but the wave-length modulations will be related to the characteristic scales: For larger $k_{T}$ and smaller $\sigma$, the oscillations in the correlation function have higher frequency and are damped at smaller scales. This can be seen in Figs. 7 and 8, for the cases $k_{T}=0.1$ and $1\,h\mathrm{Mpc}^{-1}$, respectively. The main qualitative difference compared to the $k_{T}=0.05\,h\mathrm{Mpc}^{-1}$ case is that the reinforcement with the BAO characteristic dip is not present in these cases.

The peaks and troughs are more pronounced in the linear theory since the correlation function scales simply with the scale-independent linear growth function. On the other hand for CLPT, even in the Zeldovich approximation, overdense regions are partially depleted, while underdense regions populated, due to the free-streaming of coherent matter flows over a scale settled by the Lagrangian displacements 1-dimensional variance $\sigma^{2}_{\Psi}=\int_{0}^{\infty}dkP_{L}(k)/(6\pi^{2})$.⁶⁶6This effect has the same origin as the smearing of the BAO peak observed in LPT and simulations (see, e.g. Tassev 2014). This effect is captured very well by our numerical results, and the damping in the responses become very similar for simulated data and CLPT. More remarkably for the cases of $k_{T}=1\,h\mathrm{Mpc}^{-1}$ with $\sigma=0.3$ and $0.1$, corresponding to middle and right panels of fig. 8, where the oscillations in the power spectrum are practically erased. The reason for this is, of course, that the peaks and troughs are more closer to each other and the free-streaming of particles can cover such distances. Indeed, particles travel, on average, a distance settled by the standard deviation of the displacement field $\sigma_{\Psi}(z)\sim 6\times D_{+}(z)\,h^{-1}\mathrm{Mpc}$ with $D_{+}(z)$ the linear growth function, so the process of particles moving out from more dense regions becomes very efficient.