Quintessential Isocurvature in Separate Universe Simulations

TL;DR

This paper investigates how quintessence isocurvature perturbations affect structure formation using separate universe simulations, revealing scale-dependent growth and bias modifications, and providing models to predict these effects.

Contribution

It introduces the first separate universe simulations for quintessence isocurvature cosmology, demonstrating their impact on power spectrum response and halo bias, and develops models for predicting these effects.

Findings

Growth rate depends on initial conditions and scale.

Isocurvature mode alters halo bias and power spectra.

Simulation responses can be linearly combined for mixed modes.

Abstract

In a universe with quintessence isocurvature, or perturbations in dark energy that are independent from the usual curvature perturbations, structure formation is changed qualitatively. The existence of two independent fields, curvature and isocurvature, causes the growth rate of matter perturbations to depend on their initial conditions. The quintessence perturbations cause their growth to depend on scale. We perform the first separate universe simulations for this cosmology. We demonstrate that the power spectrum response and the halo bias depend on scale and initial conditions and that the presence of the isocurvature mode changes the mapping from these quantities to the halo auto- and cross-power spectra, and the squeezed-limit bispectrum. We compare the bias to several models, finding reasonable agreement with both a power-spectrum-response model with one free parameter and a model that fits two independent bias parameters for curvature and isocurvature sourced fluctuations. We also verify that simulation responses to pure isocurvature and pure curvature modes can be linearly combined to reproduce responses with different ratios of isocurvature and curvature. This allows our results to be used to predict the halo power spectrum and stochasticity with arbitrary large-scale curvature and isocurvature power spectra. In an appendix, we study the generation of quintessence isocurvature during inflation and show that a modified kinetic term is typically required to produce observable isocurvature modes in a field with $w_Q\approx -1$.

Tables

Table 1. Table 1: Cosmological and N-body simulation parameters.

Parameter	Value
${\Omega }_Q$	0.7
${\Omega }_m$	0.3
${\Omega }_b$	0.05
$h$	0.7
$n_s$	0.968
$A_s$	2.137$\times {10}^{-9}$
$N_p$	${\left(512\right)}^3$
$L_W$	$500\mathrm{Mpc}/$$h$
$M_p$	$1.108\times {10}^{11}{M}_{\odot }$

Equations

$$\mathcal{L}$$

$$X$$

$$c_a^2\equiv \frac{{\dot{\bar{p}}}_{{}_Q}}{{\dot{\bar{\rho }}}_{{}_Q}}=w_{{}_Q},$$

$$Q^{\prime \prime }+\left(3c_{{}^Q}^2+\frac{H^{\prime }}{H}\right)Q^{\prime }+\frac{c_{{}^Q}^{2}V{,}_Q}{H^2}{\left(\frac{H^2{Q}^{\prime }{}^2}{2\Lambda }\right)}^{\frac{c_{{}^Q}^2-1}{2c_{{}^Q}^2}}=0,$$

$$Q^{\prime }$$

$$f(a)$$

$$c_a^2=c_{{}^Q}^2-\frac{1+c_{{}^Q}^2}{3}\frac{1}{f(a)}.$$

$$\displaystyle\delta\rho_{Q}\big{|}_{\rm{u}}$$

$$\displaystyle\delta p_{Q}\big{|}_{\rm{u}}$$

$$c_s^2{\delta }_Q\equiv \frac{\delta p_Q}{{\bar{\rho }}_Q}=c_{{}^Q}^2{\delta }_Q+3\left(c_{{}^Q}^2-c_a^2\right)\frac{aH}{k}{u}_Q.$$

$${\delta }_Q^{\prime }+3\left(c_s^2-w_{{}_Q}\right){\delta }_Q=$$

$$u_Q^{\prime }+(1-3w_{{}_Q})u_Q=$$

$${\delta }_m^{\prime \prime }+\left(2+\frac{H^{\prime }}{H}\right){\delta }_m^{\prime }=$$

$$+\frac{3}{2}\frac{H_0^2}{H^2}\frac{{\Omega }_Q}{a^{3(1+w_{{}_Q})}}\left(1+3c_s^2\right){\delta }_Q.$$

$${\delta }_i=\frac{\delta {\rho }_i}{{\bar{\rho }}_i}-1,$$

$${\delta }_Q(a_i,k)$$

$$u_Q(a_i,k)$$

$${\delta }_m(a_i,k)$$

$${\delta }_Q(a,k)$$

$$u_Q(a,k)$$

$${\delta }_m(a,k)$$

$$P_{mm}(k)=$$

$$+2T_m^{\mathcal{R}}(a,k)T_m^{\mathcal{I}}(a,k)P_{\mathcal{R}\mathcal{I}}(k)$$

$$+(T_m^{\mathcal{I}}(a,k){)}^2{P}_{\mathcal{II}}(k)$$

$$\equiv$$

$${\bar{\rho }}_{mW}(a;k)={\bar{\rho }}_m(a)\left(1+{\delta }_L(k,a)\right).$$

$$\frac{{\Omega }_m{H}_0^2}{a^3}(1+{\delta }_L)=\frac{{\Omega }_m{}_W{H}_{0W}^2}{a_W^3}.$$

$${\Omega }_m{}_W{H}_{0W}^2$$

$$a_W$$

$$H_W$$

$$\frac{d}{d\log a_W}$$

$${\delta }_{mW}(k_S,a;k)=D_W(a;k){\delta }_m(k_S,a_i),$$

$$\frac{d^2{D}_W}{d\log a_W{}^2}+\left(2+\frac{d\log H_W}{d\log a_W}\right)$$

$$-\frac{3}{2}$$

$$\displaystyle D_{W}(k,a)=D(a)\Big{(}1+\epsilon(a;k)\Big{)}\,.$$

$$D^{\prime \prime }+\left(2+\frac{H^{\prime }}{H}\right)D^{\prime }$$

$${ϵ}^{\prime \prime }+\left(2+2\frac{D^{\prime }}{D}+\frac{H^{\prime }}{H}\right){ϵ}^{\prime }$$

$$ϵ(a;k)=T_{ϵ}^{\mathcal{R}}(a;k)\mathcal{R}(k)+T_{ϵ}^{\mathcal{I}}(a;k)\mathcal{I}(k).$$

$$T_{ϵ}^{\mathcal{R}}(a_i;k)$$

Full text

††preprint: YITP-SB-18-43

I Introduction

The nonlinear regime of large-scale structure can provide a wealth of information about the initial conditions of our universe and its dynamical evolution. However, analytic calculations of nonlinear growth present formidable conceptual and technical challenges. The separate universe formalism has proven to be a useful tool for studying mode coupling in the evolution of large-scale structure, providing a conceptually simple framework that gives insight into the nonlinear regime. In the separate universe formalism, the nonlinearity of structure growth is studied through the effect that long-wavelength modes have on small-scale observables, such as the local power spectrum and local number densities of collapsed objects. Since the long-wavelength modes are far above the nonlinear scale, their evolution can be calculated using linear perturbation theory, while the small-scale observables, extracted from N-body simulations, can be determined deep into the nonlinear regime.

Separate universe techniques have been applied to studies of mode coupling in the production of primordial density perturbations during inflation Maldacena (2003); Creminelli and Zaldarriaga (2004), including isocurvature modes in multifield models Grin et al. (2011). The formalism was then developed for large-scale structure formation throughout matter and dark energy domination McDonald (2003); Sirko (2005); Wagner et al. (2015); Dai et al. (2015); Hu et al. (2016); Li et al. (2014a); Chiang et al. (2014); Manzotti et al. (2014); Baldauf et al. (2011, 2016); Lazeyras et al. (2016); Gnedin et al. (2011); Li et al. (2014b). This approach was applied to investigations of halo assembly bias Paranjape and Padmanabhan (2017), adiabatic perturbations in the presence of clustering quintessence Chiang et al. (2016), and the cosmological effects of massive neutrinos Chiang et al. (2017). These latter two studies found that scale dependence in structure formation can arise from scale dependent growth of the long-wavelength matter perturbations. Although, the effect was relatively small in both cases.

In this work, we apply these separate universe techniques to the case of clustering quintessence which sources isocurvature perturbations. Our goal is to study the clustering of matter on large scales in a broader context, in which the growth histories of long-wavelength matter perturbations can differ dramatically. For this purpose, we chose a model that makes the effects of history dependence large, allowing us to test our understanding of the scale dependence and time evolution of quantities such as the halo bias. The quintessence model is chosen as an academic example of a scenario with isocurvature perturbations that are important at late times and we do not require consistency of the model parameters with current data. We will, however, comment on the observational viability and physical imprint of this type of quintessence on the observed matter power spectrum, halo bias, and halo stochasticity.

Quintessence isocurvature is an extension of standard $\Lambda\rm{CDM}$ cosmology, in which a scalar field is both responsible for the dark energy content of the universe, and also gives rise to primordial entropy perturbations. These perturbations are additional degrees of freedom for the initial configuration of energy densities and pressures in the early universe. Their presence affects the evolution of structure growth and halo formation in a scale dependent way. The effect can be made arbitrarily large or small by tuning the ratio of primordial curvature and isocurvature.

The effect of quintessence isocurvature on cosmic structure enables us to study the history dependence of structure formation. Matter perturbations that are sourced by different amplitudes of primordial curvature and isocurvature will evolve differently in time. Additionally, the clustering property of the quintessence introduces a new physical scale, the quintessence Jeans scale, which causes scale dependent evolution for different modes of the matter density perturbation. Quintessential scale dependent growth was first studied in Chiang et al. (2016), using separate universe simulations with a model that contained purely adiabatic perturbations. In the current work we use the same model but with different parameters, allowing for isocurvature perturbations that are used to amplify the scale dependence in the evolution of the matter. This allows us to study structure formation in regions for which the large-scale matter and quintessence fields arrive at the same final state, but evolve through significantly different histories.

This paper is organized as follows. Section II presents the model of clustering quintessence, focusing on features of the model that allow for isocurvature perturbations that are conserved on superhorizon scales. Section III reviews the separate universe formalism and the numerical linear perturbation theory calculations that were used to fix the expansion histories of our separate universe simulations. Section IV summarizes how the simulations were performed, the parameter choices that were made, and how the results were analyzed. Section V shows the results for the matter power spectrum responses, comparing them to 1-loop perturbation theory calculations. The halo biases measured from the simulations are presented in Section VI. In Section VI we also verify that linear combinations of the individual responses measured from simulations with pure isocurvature and pure adiabatic long wavelength modes reproduce the net response for simulations with different fractions of initial isocurvature perturbations. This allows us to generate predictions for halo bias and the squeezed-limit bispectrum for arbitrary initial curvature and isocurvature power and cross-power spectra without the need for additional simulations. Section VII compares the simulation bias results to models for scale-dependent bias. In Section VIII we consider the long-wavelength power spectrum, clustering bias, and stochasticity arising from ensemble averaging over the long-wavelength primordial curvature and isocurvature modes. In the Appendix we discuss the viability of some mechanisms to generate quintessence isocurvature perturbations in the early universe.

II Isocurvature perturbations

We work with a K-essence Armendariz-Picon et al. (2000) type scalar field dark energy Lagrangian of the form,

[TABLE]

The quintessence field $Q$ is minimally coupled to the metric, but its kinetic term is the standard kinetic term $X$ raised to a power involving the model parameter $c_{{}^{Q}}^{2}$. The constant $\Lambda$ is not the cosmological constant here, it is another model parameter with the same dimensions as $X$, and its value is important for ensuring that $Q$ behaves as dark energy at late times. This dark energy model was considered in Garriga and Mukhanov (1999); Gordon and Hu (2004), and was shown to contain quintessence perturbations that are conserved outside the horizon. These dark energy perturbations act as a source of growing isocurvature modes in cold dark matter.

In the uniform field gauge (which is also the comoving gauge), where perturbations in $Q$ vanish, the sound speed is the constant $c_{Q}$. For $c_{Q}=1$, the kinetic term in Eq. (1) reduces to the canonical one. If this sound speed is chosen to be small compared to the speed of light, the dark energy perturbations cluster after horizon crossing.

The equation of state parameter, $w_{{}_{Q}}$, can be fixed by choosing an appropriate potential. However, if $w_{{}_{Q}}$ is chosen to be constant, then the adiabatic sound speed is

[TABLE]

where $\bar{\rho}_{Q}$ and $\bar{p}_{Q}$ are the energy density and pressure of the homogenous, background quintessence field. The dots denote derivatives with respect to time. As we will see in the following, gravitational coupling between the matter and quintessence perturbations makes it impossible to have conserved isocurvature on superhorizon scales under these conditions. The back reaction of matter onto quintessence will cause it to evolve even on superhorizon scales.

To avoid this, the equation of state can instead be fixed by choosing $\Lambda$ and $V(Q)$ so that the potential is large in comparison with the kinetic term. In this way, the equation of state is held at $w_{{}_{Q}}\simeq-1$ by the large value of the potential, while the adiabatic sound speed depends only on the slope of the potential.

The equation of motion for the homogeneous background field is

[TABLE]

where $H=H(a)$ is the Hubble rate and primes are derivatives with respect to $\log(a)$.

For simplicity we will assume the potential has a constant slope $V\!\!,_{Q}$, which allows us to find an exact solution,

[TABLE]

The integration constant $f(a_{i})$ fixes the value of $Q^{\prime}$ at some initial time $a_{i}$. Due to the factor of $a^{-3}H(a)$ in Eq. (6), and the small value of the sound speed, the integration constant is unimportant relative to the first term in $f(a)$ so we set $f(a_{i})$ to zero. In particular, during inflation the solution for $Q^{\prime}$ rapidly evolves to a constant that is independent of the initial condition. The slope of the potential is not required to be small, since the change in the field strength is determined by the constant $\Lambda$.

Integrating Eq. (5) gives the background evolution of the quintessence field. However, by assumption, the background field is dominated by a large integration constant $Q_{0}$. For a given potential, $V(Q_{0})$ is fixed to reproduce the present day energy density of the universe. The constant $\Lambda$ is then used to ensure that the time dependent term in $Q$ is always small compared to $Q_{0}$. This is different from the typical slow-roll assumptions, which enforce both that $w\simeq-1$, and that the fractional change in $w$ over a Hubble time is small. In our case, the fractional change in the equation of state is not small during eras when $H^{\prime}/H$ is not small, but the equation of state was so close to $w=-1$ initially that $Q$ still behaves as dark energy today.

The adiabatic speed of sound for this solution is

[TABLE]

At early times during matter domination, for $c_{{}^{Q}}^{2}\ll 1$, this gives $c_{a}^{2}\simeq-3/2$. Note, the adiabatic sound speed is only the sound speed of physical perturbations if they are adiabatic. For stability, the negative value of $c_{a}^{2}$ indicates that there must nonadiabatic stress, and therefore entropy perturbations.

The equations for linear perturbations are convenient to analyze in the synchronous gauge Bardeen (1980); Kodama and Sasaki (1984). Linear perturbations in the scalar components of the quintessence stress-energy tensor are gauge-dependent quantities. In the uniform $Q$ gauge and synchronous gauge, they are related by

[TABLE]

Here, $\delta\rho_{Q}$ is the linear perturbation to the quintessence energy density, $\delta p_{Q}$ is the linear perturbation to the its pressure, and $u_{Q}$ is the velocity potential of the quintessence field, all evaluated in the synchronous gauge. The quantities on the left are evaluated in the uniform field gauge, which is comoving with $Q$, so the velocity potential $u_{Q}|_{\rm{u}}$ vanishes.

From the above transformation, we can obtain the sound speed in the synchronous gauge, which is

[TABLE]

The equation of state and synchronous gauge sound speed can be used to eliminate the pressure and its linear perturbation from the quintessence continuity and Euler equations. The continuity equation for the pressureless matter can be used to eliminate the scalar metric degree of freedom. The resulting system of differential equations describing the linear perturbations of the combined matter-quintessence cosmic fluid is

[TABLE]

Energy density contrasts $\delta_{i}$ are defined

[TABLE]

with $i=m$ for matter or $Q$ for quintessence.

If $w_{{}_{Q}}\neq-1$, then there are no scaling solutions, $\delta_{i}\propto a^{\gamma_{i}}$ for some constant exponents $\gamma_{i}$. In this case, the quintessence perturbations evolve outside the horizon. Alternatively, if $w_{{}_{Q}}=-1$, then the matter does not appear in the quintessence continuity equation, and the condition for constant superhorizon quintessence becomes $c_{s}^{2}=-1$ for $k\ll aH$. In the limit where $c_{{}^{Q}}^{2}\ll 1$, the initial conditions for these perturbations during matter domination are

[TABLE]

where $f(a_{i})=2/9$ during matter domination. The matter perturbations are sourced by both primordial isocurvature fluctuations ($\mathcal{I}$) and curvature fluctuations ($\mathcal{R}$). The latter correspond to the homogeneous solution of the linear growth equation for matter.

Introducing transfer functions for the quintessence energy density contrast, velocity potential, and the separate curvature and isocurvature components of the matter density contrast

[TABLE]

The total linear matter power spectrum is then

[TABLE]

where $P_{\mathcal{R}\mathcal{R}}$, $P_{\mathcal{R}\mathcal{I}}$ and $P_{\mathcal{I}\mathcal{I}}$ are the auto and cross-power spectra for primordial curvature and isocurvature. In the last line we have defined $P_{mm}^{XY}$ as parts of the matter perturbation sourced by curvature and isocurvature spectra. On sub-Jeans scales, the total matter power spectrum is just the contribution from adiabatic perturbations, $P_{mm}^{\mathcal{R}\mathcal{R}}$. In this regime we will always choose $P_{mm}^{\mathcal{R}\mathcal{R}}(k)$ to be the usual $\Lambda$CDM matter power spectrum (our specific parameter choices will be given in Table 1. On larger scales, we leave the $P_{\mathcal{R}\mathcal{R}}$, $P_{\mathcal{R}\mathcal{I}}$ and $P_{\mathcal{I}\mathcal{I}}$ unspecified, but in Section VIII we will show how different assumptions about the power spectra impact observables.

The curvature-sourced transfer function depends on $k$ only through its initial condition. The isocurvature-sourced transfer function has $k$-independent initial conditions, but it has $k$-dependent evolution. This scale dependence ultimately originates in the velocity gradient in the quintessence continuity equation, which becomes important after horizon crossing. Scale-dependent growth occurs from the horizon down to the Jeans scale, below which the quintessence perturbations are pressure supported so they do not grow. Far below the Jeans scale, the matter is dominated by its adiabatic component with scale-independent growth.

Solutions to the linear growth equations are plotted in FIG. 1, with initial conditions fixed by setting the final $\delta_{m}{}_{0}$ to be the same value for all matter modes. To emphasize the effects of isocurvature we fix the ratio of primordial curvature to isocurvature to be $\mathcal{I}/\mathcal{R}=\pm 100$ outside the horizon. The comoving sound speed is set to $c_{Q}=0.1$, which fixes the quintessence Jeans scale today to be $k_{J}\simeq 10\ H_{0}$. Both the curvature and quintessence fluctuations are constant outside the horizon. After horizon crossing, the quintessence grows until crossing the Jeans scale. Below the Jeans scale, these perturbations oscillate and quickly decay away, so $k\simeq 100\ H_{0}$ is essentially a purely adiabatic mode.

During matter domination, the matter components sourced by curvature and isocurvature grow as $a$ and $a^{3}$ respectively. There are two scale-invariant regimes for matter growth, the isocurvature-dominated regime at large scales and the curvature-dominated regime at small scales. The ratio $\mathcal{I}/\mathcal{R}$ and the quintessence Jeans scale determine where the transition is between these two regimes. The larger $\mathcal{I}/\mathcal{R}$ becomes, the closer the transition is to the quintessence Jeans scale. The scale dependent growth is most extreme near these intermediate scales.

For the matter perturbations, the behavior depends on the relative sign of the primordial curvature and isocurvature fluctuations. Since the two terms in the initial condition for the matter have opposite signs, correlated curvature and isocurvature tend to decrease the matter perturbations below their adiabatic value, whereas they are increased by anticorrelated primordial fluctuations. The net effect is that matter perturbations grow rapidly at late times in the correlated case, while they approach a maximum value and then decrease when $\mathcal{R}$ and $\mathcal{I}$ are anticorrelated.

III Separate universe

The following summarizes previous work by Hu et al. Hu et al. (2016), in which separate universe techniques were developed for large-scale structure in the presence of a long-wavelength Jeans scale. In the separate universe formalism, we consider a region embedded in a long-wavelength matter perturbation $\delta_{L}$, which is approximately spatially homogeneous across the region. In this region, the long-wavelength perturbation appears as a shift in the average matter energy density,

[TABLE]

Here, the local quantities in the “windowed” region are denoted with a subscript $W$. The wave number of a long-wavelength perturbation under consideration is $k$. We restrict our analysis to windowed regions that are smaller than the quintessence Jeans scale so that the quintessence perturbations can be ignored, and the only fluctuations around $\bar{\rho}_{W}$ are fluctuations in matter. Wave numbers of modes smaller than the size of this region, and therefore smaller than the quintessence Jeans scale, will be denoted below as $k_{S}$.

Requiring that the local matter density evolves as $a_{W}^{-3}$ with respect to the local scale factor gives

[TABLE]

Since the long-wavelength perturbation is negligible at early times, the two cosmologies initially coincide. At leading order in $\delta_{L}$, this condition fixes the transformation from the global to the local cosmology,

[TABLE]

III.1 Linear perturbations

Consider a small-scale mode of the local matter perturbation $\delta_{mW}$. The evolution of this mode with respect to the local cosmology satisfies the same growth equation as adiabatic matter perturbations with respect to the global cosmology. For linear perturbations, we write $\delta_{mW}$ in terms of a local linear growth factor,

[TABLE]

where $k$ is the wave number of long-wavelength mode, and we are choosing $D_{W}(a_{i};k)=1$. If the initial condition for $\delta_{mW}$ at scale factor $a_{i}$ is set early enough, its dependence on the long-wavelength mode is negligible. The evolution of the linear growth factor is given by

[TABLE]

Small-scale perturbations are well below the Jeans scale, so at these scales the quintessence perturbations are negligible, which is why they do not show up on the right-hand side of the above expression. The local matter growth factor can be decomposed into a term equal to the global growth factor and a response term sourced by the long-wavelength mode,

[TABLE]

Using the transformations in Eqs. (26)–(28) to rewrite Eq. (30) in terms of the global cosmology gives

[TABLE]

The small-scale mode is only sourced by curvature, but the long-wavelength mode has both curvature and isocurvature contributions. We can define transfer functions for the response,

[TABLE]

During matter domination, the initial conditions are

[TABLE]

III.2 Power spectrum growth response

The local matter-matter power spectrum of a region within a long-wavelength mode will differ from the global power spectrum. The difference can be characterized by three contributions: change in comoving wavelengths, change in the definition of the average background energy density, and change in the local growth factor Takada and Hu (2013). These contributions are summarized by

[TABLE]

The first term refers to the dilation of comoving wavelengths, which is due to the local scale factor’s dependence on the long-wavelength mode. This term can be calculated

[TABLE]

The second term refers to the shift in the mean energy density, which changes the definition of energy density contrasts. This is simply

[TABLE]

The final term, referring to change in the growth factor due to the long-wavelength mode, can be estimated by taking the finite difference derivative between the local power in overdense and underdense separate universe regions,

[TABLE]

Note, this last term is the only contribution to the power spectrum response that depends dynamically on the long-wavelength mode. The other two contributions can be calculated without separate universe simulations. All of the references to the power spectrum response below refer only to the growth part of the response.

Since $R_{g}$ depends dynamically on the long-wavelength mode, it will depend on the evolutionary history of $\delta_{L}$ Ma (2007). Therefore $R_{g}$ depends on the particular ratio of adiabatic and isocurvature modes that comprise $\delta_{L}$. For a single realization, $\delta_{L}=T_{m}^{\mathcal{R}}\mathcal{R}+T_{m}^{\mathcal{I}}\mathcal{I}$, the power spectrum response is given by

[TABLE]

From Eq. (41) it is clear that the total growth response in Eq. (40) from a single realization of $\delta_{L}$ is in general a random variable, dependent on the realization of $\mathcal{I}$ and $\mathcal{R}$. This is in contrast to a cosmology in which a single degree of freedom determines $\delta_{L}$ and each mode of $\delta_{L}(k)$ has the same evolutionary history (e.g., for adiabatic perturbations in $\Lambda$CDM Wagner et al. (2015); Dai et al. (2015); Hu et al. (2016), quintessence Chiang et al. (2016), or neutrino-CDM cosmologies Chiang et al. (2017)).

In what follows it will be helpful to define individual growth responses for the pieces of $\delta_{L}$ generated by $\mathcal{I}$ and $\mathcal{R}$,

[TABLE]

Unlike the total growth response in Eq. (41), the individual responses to the adiabatic and isocurvature terms are not realization-dependent. The existence of two modes, $\mathcal{I}$ and $\mathcal{R}$, changes the usual separate universe relationships between $\epsilon$, $R_{g}$, and the squeezed-limit bispectrum. Determining the change in the power spectrum due to the presence of a long-wavelength mode from Eq. (37) requires knowledge of the local values of $\mathcal{R}$ and $\mathcal{I}$ that generated the long-wavelength mode. As we shall see, the squeezed-limit bispectrum will depend on $P_{\mathcal{R}\mathcal{R}}$, $P_{\mathcal{R}\mathcal{I}}$, and $P_{\mathcal{I}\mathcal{I}}$.

The power spectrum response for single realizations of $\delta_{L}$ is plotted in FIG. 2 as a function of the large-scale wave number. Far above the Jeans scale, the response is isocurvature dominated and scale invariant. Far below the Jeans scale, the response approaches the purely adiabatic value. Between the horizon crossing scale and Jeans scale the response is scale dependent. In the correlated case ($\mathcal{I}/\mathcal{R}>0$), the initial conditions for the curvature and isocurvature components of the matter modes have opposite sign. At a given redshift there is a scale at which the two components cancel, so the total matter perturbation vanishes, causing the response to diverge at this scale.

The squeezed-limit bispectrum can be determined by computing the correlation between the local power spectrum and $\delta_{L}$, $\langle P_{mmW}(a,k_{S}|\delta_{L})\delta_{L}\rangle$. Computing the ratio we find

[TABLE]

The coefficients of the growth responses depend on the primordial power spectra, $P_{\mathcal{R}\mathcal{R}}$, $P_{\mathcal{R}\mathcal{I}}$, and $P_{\mathcal{I}\mathcal{I}}$, which depend on the mechanism, inflationary or otherwise, that generates these fluctuations in the early universe. Here we make no attempt to find a model that would generate these primordial fluctuations in the early universe. However, the Appendix at the end of this paper considers the viability of certain models that could produce isocurvature modes that are significantly correlated or anticorrelated with the primordial curvature modes.

The last two terms in Eq. (45) can be interpreted as the growth part of the power spectrum response averaged over the long-wavelength mode in the following sense:

[TABLE]

For the purposes of illustrating the effects of averaging over long-wavelength modes, we make the assumption that both of the auto power spectra and the cross power spectrum of the primordial curvature and isocurvature fluctuations are scale invariant. FIG. 2 shows examples of the averaged power spectrum responses with the isocurvature amplitude taken to be 100 times larger than the curvature amplitude to emphasize its effect. Examples are shown with cross-correlation vanishing, and with the cross-correlation $\pm 90\%$ its bound with respect to the Cauchy-Schwartz inequality, $|P_{\mathcal{R}\mathcal{I}}|\leq\sqrt{P_{\mathcal{R}\mathcal{R}}P_{\mathcal{I}\mathcal{I}}}$.

Responses of small-scale observables to a long-wavelength mode that are described by the separate universe formalism are closely related to angle-averaged, equal-time cosmic consistency relations. The presence of quintessence isocurvature violates the assumptions used to derived the standard cosmic consistency relations Valageas (2014). The expression for the squeezed-limit bispectrum shown in Eq. (45) is a generalization of the bispectrum angle-averaged, equal-time cosmic consistency relation to the case of quintessence isocurvature.

IV Simulations

Separate universe simulations involve computing the evolution of an N-body system of particles within the context of a long-wavelength mode that is treated as a shift in the homogeneous background energy density. In our case, the particles are cold dark matter and baryons, which are treated equivalently because we are interested in scales much larger than the baryonic Jeans scale. The system includes $N_{p}=(512)^{3}$ nonrelativistic particles interacting only through Newtonian gravity on an Friedmann-Robertson-Walker expanding background. The comoving size of the box is fixed locally to length L_{W}=500\ \rm{Mpc}$$/h.

The expansion is characterized by the local Hubble rate $H_{W}$ as a function of local scale factor $a_{W}$, which differs from $\Lambda$CDM by a perturbatively small contribution due to the presence of the long-wavelength mode $\delta_{L}$. At a given wave number, $\delta_{L}$ is computed by numerically integrating Eqs. (15)–(17). Values of the local scale factor and Hubble rate are tabulated using Eqs. (26)–(28). The background expansion is fixed entirely by $\delta_{L}$ and the choice of global cosmology parameters $\Omega_{m}=0.3$, $\Omega_{Q}=0.7$, and $w_{{}_{Q}}=-1$.

The N-body simulations were run using the code Gadget2 Springel (2005), modified to read in tabulated values of the separate universe scale factor and Hubble rate for a given long-wavelength mode. Sets of 20 simulations were run for the long-wavelength modes $k/H_{0}$=1, 7, 10, 14, and 100 with $\mathcal{I}/\mathcal{R}=-100$, and $k/H_{0}$=1, 10, and 100 with $\mathcal{I}/\mathcal{R}=100$. These range from the isocurvature to the curvature dominated regime, capturing features of the scale dependent power spectrum response. Additional simulations were run for a pure isocurvature mode at $k/H_{0}=10$, and a pure curvature mode. With both overdense and underdense separate universe simulations for each mode, 400 separate universe simulation were run in total. An additional 20 were from for a larger global universe box of size $L=1000\ \rm{Mpc}/h$ with $(512)^{3}$ particles with standard adiabatic initial conditions. These were used to measure the clustering bias as a check on our separate universe results for the adiabatic response bias.

Initial conditions were generated as realizations of Gaussian random fields from the matter power spectrum at the initial simulation time $z_{i}=49$. Although the baryons behave indistinguishably from the dark matter in the N-body interactions, their presence affects the shape of the power spectrum through the baryonic acoustic oscillations. To account for this, the power spectrum is first calculated at redshift $z=0$ with the code CLASS Blas et al. (2011), using the cosmological parameters given in Table 1. The matter power spectrum is then rescaled back to the initial simulation time. At early times, the difference between local and global cosmology is negligible, so the primordial fluctuations are the same for both. However, the linear growth factors are not the same. To account for this, we rescale the power spectrum back using the global growth factor and forward to the initial simulation time using local scale factor

[TABLE]

For each long-wavelength mode, the same 20 random seeds were used to generate the initial conditions for the particles’ positions and velocities. These were corrected using a second order Lagrangian perturbation theory code to reduce transients Crocce et al. (2006). Simulation snapshots were taken at values of the global redshifts $z=1.0,\ 0.75,\ 0.5,\ 0.25,\ \mathrm{and}\ 0.0$. Since the simulations are run in local time, their snapshot output times are adjusted according to Eq. (26) to make sure overdense and underdense boxes are matched at the same the global time.

V Power spectrum response

For a given simulation snapshot, the particle positions were converted to a density field, which was Fourier transformed using FFTW Frigo and Johnson (2005). The density field was estimated by distributing the particles to sites on a $(1024)^{3}$ grid, using the cloud-in-cell method. The power spectrum was then calculated from the Fourier transformed density field, and the power spectrum response was estimated through the finite difference derivative in Eq. (40), from the overdense and underdense separate universe simulations. By resampling over the 20 random realizations, the bootstrap variance in the power spectrum response was calculated.

Effects of weakly nonlinear growth can be captured by calculating the power spectrum using perturbation theory. At 1-loop we have Jeong and Komatsu (2006)

[TABLE]

where the first term is the linear power spectrum and last two terms are proportional to $\left(D_{W}\right)^{4}$. The 1-loop correction to the power spectrum response is

[TABLE]

Here $R_{g}$ is the linear power spectrum response. The power spectrum responses measured from our simulations at $z=0.25$ are plotted in FIG. 3, along with a 1-loop calculation shown for each value of $k$ that was simulated. Above the nonlinear scale, the response measured from the simulations agrees well with the linear power spectrum response. Agreement with the 1-loop calculation continues until $k_{S}\simeq 6\times 10^{-2}\ \rm{Mpc}^{-1}$, which is approaching the nonlinear scale where perturbation theory is no longer valid.

VI Halo bias

The Lagrangian bias was determined for each simulation set by using the abundance matching method Li et al. (2016), which involves calculating the mass shift at a fixed cumulative number density between overdense and underdense separate universes. For a given set of simulations, a cumulative catalog of halo masses, from all 20 realization of the initial conditions, was constructed and then sorted in descending order. These sorted mass catalogs give an estimate of the cumulative number density of objects above some threshold mass,

[TABLE]

From our sorted halo catalogs, we construct the following lists:

[TABLE]

The $\pm$ superscripts refer to overdense and underdense separate universes. By fitting the above lists with splines, we obtained estimates of the functions $s(\log\!M)$, which is the mass shift, and $n(\log\!M)$, which is the halo mass function. The average Lagrangian bias above mass $M$ is then estimated

[TABLE]

Note, all biases reported and discussed below will be biases averaged above a mass threshold.

The halo finding program Rockstar Behroozi et al. (2013) was used to extract catalogs of bound objects and their spherical overdense masses. Halo masses are calculated in Rockstar by finding the outermost particle position from the halo’s center of mass at which the average spherical overdensity inside this radius is greater than some threshold,

[TABLE]

The density threshold used here is the virialization threshold

[TABLE]

Since Rockstar determines densities with respect to the local, separate universe comoving distances, the threshold has to be converted to these coordinates,

[TABLE]

That is, an overdense separate universe has an apparently lower density threshold to form collapsed objects.

At small masses, the discreteness of this halo mass determination introduces an artificial spread in the values $s_{i}$. We modified this mass calculation, based on the work Li et al. (2016), to include a contribution from particle $N+1$, assuming its mass to be uniformly distributed in the spherical shell between $r_{N+1}$ and $r_{N}$,

[TABLE]

Here $V_{N}$ is the volume inside the radial position of the $N^{th}$ particle from the halo’s center. Masses calculated this way are continuous, and produce less spread in the values of the mass shift at low mass.

We repeated this calculation of the Lagrangian bias, resampling with replacement over our sets of random initial conditions, and calculated the bootstrap variance. The biases at redshift $z=0$ are shown in FIG. 4, along with the ratio of the Lagrangian biases to the adiabatic ($\mathcal{I}=0$) Lagrangian bias. We will refer to the latter quantity as the relative bias.

The scale dependence is clearly demonstrated within the bootstrap variances from the simulations, especially for the anticorrelated case. The recovery of the adiabatic bias below the Jeans scale is demonstrated for the $k/H_{0}=100$ mode.

The relative biases are consistent with being mass independent for the masses we are sensitive to $(\sim 10^{13}$–$10^{15}\ \rm{M}_{\odot})$. Small variations in the relative bias as a function of mass is due to the spline fitting. At masses above $10^{15}\ \rm{M}_{\odot}$, halos are rare and the spline fitting becomes very sensitive to the choice of knots, so the bias is not well determined there.

As a test of our bias determinations, we compared them to measurements of the clustering bias from a simulation with box size $(1000\ \mathrm{Mpc})^{3}$, number of particles $(512)^{3}$, and a standard, global $\Lambda$CDM background cosmology with purely adiabatic perturbations. The clustering bias, $b_{{{}_{C}}}$, is measured from the matter-halo cross power spectrum on linear scales Dvali et al. (2004a); Desjacques et al. (2010); Assassi et al. (2014),

[TABLE]

Treating $b_{2}$ as a nuisance parameter, we fit the above function for value of $k_{S}$ up to $0.05\ \rm{Mpc}^{-1}$. We found the cluster bias to be in excellent agreement with our separate universe response biases for the adiabatic mode.

VI.1 Linearity of responses

Our simulations involve two independent components of the long-wavelength modes: the curvature and the quintessence-sourced isocurvature. We fix the initial relative amplitudes of these as part of our simulation parameters. However, since both components of the long-wavelength mode stay perturbatively small, it should be the case that responses can be decomposed into two terms, one for each component.111It is always, of course, a working assumption of the separate universe approach that responses computed from simulations with a particular realization of the long-wavelength mode can be combined to produce the response for a different realization, e.g. for a larger or smaller amplitude $\delta_{L}$ or a $\delta_{L}$ with a different density profile. We nevertheless want to verify this assumption for the more extreme examples of different $\delta_{L}$ considered in this paper.

Suppose we have a small-scale observable $\mathcal{O}_{W}$, which responds to the long-wavelength mode,

[TABLE]

In this way, we can take linear combinations of responses from simulations with purely adiabatic and purely isocurvature long-wavelength modes and obtain the response for any realization of $\delta_{L}$ with initial relative amplitude $\mathcal{I}/\mathcal{R}$.

The evolution of the total long-wavelength matter perturbation with both curvature-sourced and isocurvature-sourced components can differ significantly from the purely isocurvature-sourced and purely curvature-sourced perturbations. The evolution even becomes nonmonotonic in the case with $\mathcal{I}/\mathcal{R}=-100$. In order to test the linearity of responses, an additional set of simulations was run with long-wavelength mode that was purely isocurvature-sourced, at wave number $k/H_{0}=10$. The power spectrum response and biases obtained by linearly combining simulation results with the pure isocurvature and pure curvature long-wavelength modes are shown in FIG. 5.

Since the responses are estimated using finite difference derivatives with respect to the long-wavelength mode, we expect numerical errors on the order of $\delta_{L}\sim 1\%$. In addition to this, we have cosmic variance from the limited set of realizations for the initial conditions. The difference in the power spectrum response obtained from linear combination and the total long-wavelength mode simulations is $<2\%$ even into the nonlinear regime, which is consistent with our expected numerical accuracy.

For the biases, the error is dominated by shot noise and stochasticity. At masses between $10^{13}$–$10^{15}\ \rm{M}_{\odot}$, the two estimates of the Lagrangian bias are within $1.5\ \sigma$, so linearly combining purely isocurvature-sourced and purely curvature-sourced simulation results gives a good estimate of the small-scale observable responses. As a result, our simulation responses can be rescaled and linearly combined to estimate what the responses would be in the presence of different long-wavelength modes with particular realizations of $\mathcal{I}$ and $\mathcal{R}$, and ultimately, different primordial power spectra and cross spectra for $\mathcal{I}$ and $\mathcal{R}$ as discussed in III.

VII Bias models

The linearity of the biases allows us to identify two independent bias coefficients, $b_{\mathcal{R}}^{L}(a)$ and $b_{\mathcal{I}}^{L}(a,k)$, for purely curvature-sourced and isocurvature-sourced matter perturbations respectively. In general, these have different mass dependence, so their ratio depends on both mass and wave number. For a given realization of $\mathcal{I}/\mathcal{R}$, the total relative bias is

[TABLE]

To predict the total relative bias, we must model the relative bias for the purely isocurvature-sourced modes. This amounts to describing the scale dependence and time evolution of the pure isocurvature bias. We define the isocurvature relative bias,

[TABLE]

Different bias models correspond to different choices for the $k$-dependence and evolution of $\beta^{L}$.

Below we consider two classes of models. The first class attempts to reproduce the pure isocurvature relative bias using the transfer functions $T_{m}^{\mathcal{R}}$ and $T_{m}^{\mathcal{I}}$. This fixes the scale dependence and evolution of the isocurvature bias up to an overall multiplicative factor, which we take to be independent of scale and fit to the isocurvature dominated scales ($k/H_{0}\leq 1.$). The transfer function models generically involve at least one free parameter.

The second class of models assumes that the bias is proportional to another small-scale observable’s response to the long wavelength mode, such as the power spectrum response, or the response of the critical linear density for spherical collapse. The relative bias for this class of models is just the ratio between the pure isocurvature and pure curvature responses, so these models involve no free parameters.

VII.1 Transfer function models

Models of halo bias that are local in time assume that the fluctuations in halo number density are fixed by the configuration of the matter density field at a single redshift. The evolution of halo abundance is then an initial value problem, with initial data given by the matter density field’s configuration on a single time slice. This is, for instance, the point of view adopted in the excursion set theory approach to large scale structure. In this approach, the matter density field at early times is smoothed over a range of scales and compared to the critical density for spherical collapse linearly evolved back to the redshift of the initial data. From this point of view, it is natural to model the evolution and scale dependence of the halo bias using the transfer functions for the linear evolution of the different components of matter fluctuations.

For these models the Lagrangian biases have simple evolutions, inversely proportional to the transfer functions. We will consider two different possibilities for the evolution and scale dependence of $\beta(a,k)$. The first assumes passive evolution of halo abundance, so that the number density of halos at a given mass is conserved. Scale dependence in the bias then arises from scale dependence in the isocurvature transfer function at late times. The second model assumes that the individual biases are scale invariant at all times. Scale-dependent bias in this case arises from when the two bias terms are combined, through the scale dependence of the total matter perturbation.

The models we consider are based on the matter transfer functions $T_{m}^{\mathcal{R}}$ and $T_{m}^{\mathcal{I}}$, which is motivated by the linearity of the Lagrangian response bias. Another approach would be to include a bias term for the quintessence density contrast ($b_{Q}\delta_{Q}$), combined with a bias term for either the full matter perturbation, or just the curvature sourced part. This is similar to including a neutrino bias term for the case of scale dependent growth arising from massive neutrinos. This approach was found not to reproduce biases from separate universe simulations Chiang et al. (2017). Similarly, we find that the scale dependence from the quintessence transfer function, when combined with the matter transfer functions, does not reproduce the scale dependence of the simulation biases. Fits from these models have $\chi^{2}\sim 10^{3}$ (per degree of freedom), so we do not consider them below.

VII.1.1 Passive evolution model

A simple example of bias evolution, based on work by Hui and Parfrey Hui and Parfrey (2008); Parfrey et al. (2011), assumes that the number density contrast of halos is conserved,

[TABLE]

This passive evolution model (PE in what follows) is accurate so long as merger events between different halos are sufficiently rare. The evolution of the individual biases is given by

[TABLE]

Each of these individual bias terms are independent of the realization of the long-wavelength mode. The scale dependence from the initial condition in $T_{m}^{\mathcal{R}}$ cancels taking the ratio. While the initial condition for $T_{m}^{\mathcal{I}}$ does not depend on $k$, its subsequent evolution after horizon crossing does, and this is where the scale dependent bias originates in this model.

The relative bias for pure isocurvature is

[TABLE]

The free parameter $\beta^{L}_{PE}(M,a_{i})\equiv b^{L}_{\mathcal{I}}(M,a_{i})/b^{L}_{\mathcal{R}}(M,a_{i})$, can be fit by requiring this relative bias to match the simulation results on isocurvature dominated scales.

The total relative bias can be written

[TABLE]

Since the initial transfer function $T_{m}^{\mathcal{R}}(a_{i},k)$ grows as $k^{2}$, the correct small-scale relative bias $b^{E}/b^{E}_{\mathcal{R}}=1$ is recovered in the curvature dominated regime ($k\rightarrow\infty$).

The model fits at mass $M=10^{14}\ \rm{M}_{\odot}$ are shown in FIG. 6 at redshifts $z=0$ and $z=0.25$. The model does not reproduce the scale dependence of our simulation results. Forcing the model to agree with the isocurvature dominated results causes the relative bias to be dramatically over estimated for $\mathcal{I}/\mathcal{R}=-100$. For $\mathcal{I}/\mathcal{R}=100$, the relative bias at the Jeans scale is underestimated. Similar versions of this model were also shown to give poor reproduction of simulation results for adiabatic quintessence Chiang et al. (2016) and massive neutrinos Chiang et al. (2017).

VII.1.2 Constant $\mathcal{R}$-$\mathcal{I}$ bias model

Suppose that, instead of taking the evolution of $b^{L}_{\mathcal{I}}$ to be determined by the assumption of halo number density contrast conservation, we instead take it to be scale independent at all times. Then the total relative bias is given by

[TABLE]

We can again fit the parameter $\beta_{\mathcal{R}\mathcal{I}}^{L}(M,a)$ to match the isocurvature dominated scales at each redshift. We call this model the constant $\mathcal{R}$-$\mathcal{I}$ bias model (c$\mathcal{R}\mathcal{I}$). Although the individual components of the bias are assumed to be independent of scale, the total bias does depend on scale through the isocurvature transfer function. The component biases each depend on mass and redshift, and there is nothing requiring the relative bias to be independent of mass for this model.

Fits for this model are have $\chi^{2}\simeq 2$ per degree of freedom. For $\mathcal{I}/\mathcal{R}=100$, the Jeans scale relative bias is slightly over estimated, whereas the $\mathcal{I}/\mathcal{R}=-100$ biases are underestimated.

Overall, the transfer function models appear to capture the qualitative features of the relative biases. The scale dependence seems to be dominated by the shape of the function $\delta_{m}(a,k)^{-1}$, which is what one expects for passive halo evolution. Quantitatively, the c$\mathcal{R}\mathcal{I}$ model gives a much improved reproduction of the simulation biases compared with the PE model. The quality of these model fits is not strongly dependent on redshift. The relative success of the c$\mathcal{R}\mathcal{I}$ model indicates that the scale dependence of the isocurvature relative bias is fairly weak.

VII.2 Response models

If the number density of collapsed objects cannot be determined from the configuration of the density field at a single time, but instead relies on the cumulative growth history of the matter density, then the bias is nonlocal in time. In this case, we can try to characterize the local mass function as depending on some other small scale observable, which, in the separate universe context depends on the long-wavelength modes that are present in a given region,

[TABLE]

The response bias is then proportional to the growth response of $\mathcal{O}_{W}$,

[TABLE]

The first factor on the right side of the above equation is evaluated in the global universe, so it is only a function of mass, time, and the global value of the quantity $\mathcal{O}$. That is, this factor is independent of the long-wavelength mode. Then the relative bias is

[TABLE]

Notice that for this class of models, the relative bias is predicted to be independent of mass222If Eq. (73) were evaluated at the object formation time, rather than the observation time, then the relative bias would be mass-dependent since objects of different masses typically form at different times and are therefore sensitive to $R_{\mathcal{O}}(a,k)$ at different epochs.. Also, unlike the transfer function models, these models have no free parameters. The relative biases are completely determined by linear perturbation theory and the separate universe formalism.

We considered response models based on two possible choices for $\mathcal{O}$: the local power spectrum, and the local critical density for spherical collapse LoVerde (2014). These were shown in Chiang et al. (2016, 2017, 2018) to perform better than bias models based on transfer functions. We found that the results of the spherical collapse model match the results of the power spectrum response model, so we consider only the latter model in what follows.

Assume the local mass function in a region depends on the local power spectrum in a universal way, then the pure isocurvature relative bias is given

[TABLE]

This model is compared in FIG. 6 with relative biases from the simulations at $M=10^{14}\ \rm{M}_{\odot}$ for redshifts $z=0$ and $z=0.25$. We obtained $\chi^{2}=1.4$ per degree of freedom at the later redshift, and $\chi^{2}=2.6$ at the earlier redshift. It slightly underpredicts the Jeans scale bias in the $\mathcal{I}/\mathcal{R}=100$ simulations, and overpredicts the $\mathcal{I}/\mathcal{R}=-100$ biases.

At smaller values of $\mathcal{I}/\mathcal{R}$, the power spectrum response model does much better. The linearly combined pure isocurvature and pure curvature relative biases are shown in FIG. 7, for $k/H_{0}=10$. Agreement between the linear combination bias and the $R_{P}$ model is good for $|\mathcal{I}/\mathcal{R}|<100$.

VIII Observables

In real surveys or in large box simulations, the regions covered can be orders of magnitude above the quintessence Jeans scale considered here. The long-wavelength mode in these cases can also be sampled over. The effects of averaging over $\delta_{L}$ to obtain the squeeze limit bispectrum were already shown in FIG. 2, assuming scale invariant auto and cross primordial power spectra. Here will consider additional observables, including the large-scale matter power spectrum, the clustering bias, and the stochasticity.

Clustering bias and stochasticity are defined by taking correlations between the long-wavelength matter perturbations and the halo number of density fluctuations (e.g. Matsubara (1999); Hui and Parfrey (2008); Smith and LoVerde (2011)). Writing separate bias terms for the curvature and isocurvature modes

[TABLE]

where $\delta_{h}$ is the density contrast for the number halos and we have suppressed all arguments, including the $k$-dependence in $b_{\mathcal{I}}$, $\delta_{L}^{R}$, and $\delta_{L}^{\mathcal{I}}$, for notational clarity. We then have the following three power spectra for matter perturbations and halos:

[TABLE]

Note that for $b_{\mathcal{R}}\neq b_{\mathcal{I}}$ the above set of equations is nondegenerate so that equations can be inverted to solve for the individual contributions to the matter power spectrum sourced by $\mathcal{I}$ and $\mathcal{R}$, e.g. $P_{mm}^{\mathcal{R}\mathcal{R}}$, $P_{mm}^{\mathcal{I}\mathcal{R}}$, and $P_{mm}^{\mathcal{I}\mathcal{I}}$, from the observed $P_{mm}$, $P_{mh}$ and $P_{hh}$. The clustering bias is defined as

[TABLE]

The stochasticity is defined

[TABLE]

Here we will assume that all three primordial power spectra are scale invariant, and we take the amplitude of isocurvature fluctuations to equal the amplitude of scalar curvature fluctuations. Finally, we consider cases where the cross power spectrum is $65\%$ correlated or anticorrelated, and the case where it is totally uncorrelated.

Since we take the isocurvature amplitude to equal the curvature amplitude, we are well within the regime where the parameter-free power spectrum response model adequately describes the bias.

[TABLE]

where $X$ can be either $\mathcal{R}$ or $I$. The biases in this section are all Eulerian ($b^{E}=1+b^{L}$), although we will omit the superscript $E$ that indicate this.

Using the biases from adiabatic mode simulations, the mass and redshift dependent function can be fit and used to calculate the power spectrum response model biases at large-scales in the presence of isocurvature. With the bias determined at five redshifts, we fit this for mass $M=1\times 10^{14}\ \rm{M}_{\odot}$ as a quadratic function of scale factor.

These observables are plotted in FIG. 8 for redshifts ranging from $z=0$–$3$. The most significant feature predicted by this model is the growth of power at very large scales with $k<10^{-4}\ \rm{Mpc}^{-1}$. This is the isocurvature curvature dominated regime, above the quintessence Jeans scale. Current observations do not yet extend to scales with $k\lesssim 10^{-3}\ \rm{Mpc}^{-1}$ Tegmark and Zaldarriaga (2002); Hlozek et al. (2012).

The growth of the power spectrum at large, isocurvature dominated scales is due to that fact that $T_{m}^{\mathcal{I}}/T_{m}^{\mathcal{R}}\sim k^{-4}$. As long as $P_{\mathcal{I}\mathcal{I}}/P_{\mathcal{R}\mathcal{R}}$ does not scale as a higher power of $k$ than $k^{4}$, $P^{\mathcal{I}\mathcal{I}}_{mm}$ dominates the power spectrum on scales far above the Jeans scale ($k\ll k_{Jeans}\simeq 2\times 10^{-3}\ \rm{Mpc}^{-1}$), and the power spectrum will increase with decreasing wave number. Below the Jeans scale, $P^{\mathcal{R}\mathcal{R}}_{mm}$ dominates, so the typical $\Lambda$CDM matter power spectrum with adiabatic initial conditions is recovered on these scales.

Near the Jeans scale, the shape and evolution of the matter power spectrum depends on the curvature-isocurvature cross-correlation term $P^{\mathcal{R}\mathcal{I}}_{mm}$, as well as the ratio $P_{\mathcal{I}\mathcal{I}}/P_{\mathcal{R}\mathcal{R}}$. For positive cross-correlation, the total matter power spectrum has excess power compared with $P^{\mathcal{R}\mathcal{R}}_{mm}$ at the Jeans scale due to the peak in the isocurvature transfer function. The power spectrum drops below $P^{\mathcal{R}\mathcal{R}}_{mm}$ at larger scales due to the change in sign of the isocurvature transfer function. In the limit of $100\%$ correlation, there is a wave number above the Jeans scale at which the power spectrum vanishes. The total power spectrum increases as $P_{\mathcal{I}\mathcal{I}}$ starts to dominate near and above the horizon scale. Similarly, in the anticorrelated case, the total matter power spectrum at the Jeans scale is deficient compared with $P^{\mathcal{R}\mathcal{R}}_{mm}$. For the parameter choices in FIG. 8, the deviation between $P_{mm}$ and $P^{\mathcal{R}\mathcal{R}}_{mm}$ is about $1\%$ at the Jeans scale. For parameters comparable to the single realizations of our separate universe simulations, $A_{\mathcal{I}\mathcal{I}}=10^{4}A_{\mathcal{R}\mathcal{R}}$, the Jeans scale matter power spectrum differs from $P_{\mathcal{R}\mathcal{R}}$ by about $300\%$. However, the deviation is less than $1\%$ at scales below $k\gtrsim 8\times 10^{-3}\ \mathrm{Mpc^{-1}}$.

There is also clear scale dependence of the clustering bias above the Jeans scale. This scale dependence is sensitive to the cross-correlation between curvature and isocurvature fluctuations. The stochasticity is peaked at a scale above the Jeans scale. The location of this peak depends on the ratio of the curvature and isocurvature amplitudes. For larger isocurvature amplitudes, this scale approaches the quintessence Jeans scale. The height and width of the peaks in the stochasticity are sensitive to the cross-correlation. Over time, both the clustering bias and stochasticity evolve to approach unity, so the scale dependence is most prominent at earlier redshifts. However, the scale dependence shows up at larger scales at earlier times.

For completeness we note that the presence an isocurvature component in the matter perturbations will also have an effect on redshift space distortions. An object, such as a halo, has redshift space coordinate $\vec{s}$ defined by

[TABLE]

where $r$ is the radial distance to the object, $\vec{v}$ is the peculiar velocity of the object, and $\hat{r}$ is the line-of-site direction. Here we will consider the form of the redshift space matter power spectrum derived by Kaiser Kaiser (1987). The Jacobian of the transformation from real space to redshift space, and the linear matter continuity equation can be used to express the Fourier modes of the density contrast for halo numbers in redshift space,

[TABLE]

The coefficient $\mu(\vec{k})$ is the cosine of the angle between the wave vector and the line of sight. The superscripts $s$ and $r$ denote the density contrast in redshift space and real space respectively. Substituting the full matter density contrast with transfer functions and the halo biases, the halo power spectrum in redshift space, $P^{s}_{mm}$, is

[TABLE]

where we have defined $f_{\mathcal{R}}\equiv d\log T_{m}^{\mathcal{R}}/d\log a$, $f_{\mathcal{I}}\equiv d\log T_{m}^{\mathcal{I}}/d\log a$. Unlike in the case of purely adiabatic fluctuations, taking the ratio of the $\mu^{4}$ terms with the $\mu^{2}$ terms does not isolate the growth rate and bias from the primordial spectra. On the other hand, at each order in $\mu^{2}$ the terms in the redshift space power spectrum contain distinct weightings of $P_{mm}^{\mathcal{R}\mathcal{R}}$, $P_{mm}^{\mathcal{I}\mathcal{R}}$, and $P_{mm}^{\mathcal{I}\mathcal{I}}$ so that, even if $b_{\mathcal{R}}=b_{\mathcal{I}}$, the isocurvature and curvature power spectra could be solved for in terms of $f_{\mathcal{R}}$, $f_{\mathcal{I}}$, $b_{\mathcal{R}}$, and $b_{\mathcal{I}}$.

Since the quintessence isocurvature perturbations are only important at large scales and at late times, it is difficult to constrain the amplitude of the isocurvature auto power spectrum and the correlation between isocurvature and curvature fluctuations. For instance, an analysis of the 7-year WMAP data allowing for dark energy isocurvature leaves the dark energy sound speed and primordial isocurvature power spectra and cross-spectra virtually unconstrained Liu et al. (2011). A full analysis of this model with current CMB temperature and polarization data, along with large-scale structure data is needed to determine allowed ranges of the parameters $c_{Q}$, $P_{\mathcal{I}\mathcal{I}}$, and $P_{\mathcal{I}\mathcal{R}}$.

IX Conclusion

Quintessence isocurvature perturbations provide a scenario in which large-scale structure formation can be studied in a context that is more general than the canonical $\Lambda$CDM model. The scalar field dark energy perturbations introduce both scale dependent growth, through the presence of a quintessence Jeans scale, and also initial condition dependent growth history. As shown in FIG. 1, large-scale matter perturbations can have dramatically different, and even nonmonotonic evolution in this context, depending on scale and initial conditions.

Using separate universe techniques and N-body simulations, we have studied the responses of small-scale observables such as the local power spectrum (FIG. 3) and halo mass function (FIG. 4) to the presence of long-wavelength perturbations. Our methods have been validated by comparing the results of the power spectrum response to 1-loop calculations from perturbation theory and by comparing the response biases for the adiabatic case to the clustering bias in larger volume simulations. We have also verified, in FIG. 5, that the individual responses to long-wavelength curvature and isocurvature-sourced perturbations can be linearly combined to study the effects of varying the initial conditions. That is, the separate universe simulations can be used to make predictions for any initial large-scale curvature and isocurvature power spectra.

The linearity of small-scale observable responses guarantees that the evolution of the total bias depends on both scale and initial conditions. While scale-dependent bias has previously been found, originating from scale dependent growth due to neutrinos and quintessence without isocurvature, this is the first instance of bias evolution that depends on initial conditions of the long-wavelength matter modes.

In section VII, we tested a number of simple bias evolution models against our simulation results. These models break down into two classes: transfer function models and response models. While all of the models roughly capture the dependence on scale and initial conditions of the bias evolution, FIG. 6 demonstrates that two of the models considered reproduce the scale dependent features in the bias. One is the the transfer function model that assumed scale invariant bias factors for both the curvature-sourced and isocurvature-sourced components of the matter perturbations (the c$\mathcal{R}\mathcal{I}$ model). The success of this bias model is a consequence of the linearity of the response observables, and the weak scale dependence of the purely isocurvature-sourced halo bias. The other viable bias model was based on the power spectrum response ($R_{g}$). This model outperformed the c$\mathcal{R}\mathcal{I}$ model at low redshift, but was worse by redshift $z=0.25$. On the other hand, the response model is quite accurate for more modest amplitudes of isocurvature perturbations (see. FIG.7). Since this model has no free parameters, it can be used to make predictions without running additional simulations. The third model we considered, based on an assumption of passive halo evolution (PE), gave poor reproduction of the simulation results at all masses and redshifts.

The dependence of small-scale observables on growth history, through the wave number and initial conditions of a long-wavelength mode, indicates that structure formation is nonlocal in time. That is, it is not enough to know the statistics of the density field at a single redshift in order to predict the statistics of halos at another redshift. In this work we have provided a demonstration of history dependent structure formation in a model where the effect can be made arbitrarily large, in the sense that the ratio of isocurvature to curvature can be large. In principle, the relative strength of isocurvature fluctuations to curvature, as well as the quintessence sound speed can be constrained by large-scale structure observations and CMB measurements.

Acknowledgements.

DJ and ML are grateful for helpful correspondence and conversations with Niayesh Afshordi, Chi-Ting Chiang, Wayne Hu, Yin Li, and Ravi Sheth. Results in this paper were obtained using the high-performance computing system at the Institute for Advanced Computational Science at Stony Brook University. DJ is supported by Grants No. NSF PHY-1620628 and DOE DE-SC0017848. ML is supported by Grant No. DOE DE-SC0017848.

Appendix A On the generation of quintessence isocurvature perturbations

In this appendix we will make some basic observations about the possibility of generating quintessence isocurvature perturbations that are large enough, in comparison with adiabatic perturbations, to have observational impact at late times.

Dark energy isocurvature perturbations that are perfectly anticorrelated with the primordial curvature perturbations have been studied by Gordon and Hu (2004), among others, for the purposes of producing a low quadrupole amplitude in the CMB temperature power spectrum. The mechanism studied involves a delayed, postinflationary production of curvature perturbations sourced by the same field that sources the dark energy isocurvature perturbations. This is achieved, for instance, by making the quintessence field the field that modulates reheating Dvali et al. (2004b, a) or via curvatons Linde and Mukhanov (1997); Lyth and Wands (2002); Moroi and Takahashi (2001); Enqvist and Sloth (2002).

In Gordon and Hu (2004), the authors argue that inflationary generation of anticorrelated quintessence isocurvature perturbations generally makes the tensor-to-scalar ratio too large. However, in Gordon and Hu (2004) the scalar perturbations sourced by the inflaton were neglected. If included, they keep the tensor-to-scalar ratio small but will unavoidably dominate over the scalar curvature perturbations generated after inflation. This necessarily makes the isocurvature perturbations small and only weakly correlated to the adiabatic modes. On the other hand, we find it is generally possible to satisfy these criteria with a modified kinetic term, including the one used in the main text of this paper [see Eq. 1)].

Consider a scalar field $Q$ with a Lagrangian that is a general function of $X$ and $Q$,

[TABLE]

where

[TABLE]

The energy density and pressure are

[TABLE]

For a homogeneous background, we have

[TABLE]

which gives equation of motion for $Q$

[TABLE]

where the effective sound speed is given by

[TABLE]

It will also be convenient to define

[TABLE]

In absence of derivative couplings of the field to itself, $P_{,XQ}=0$, so we have $c_{X}^{2}=-1$.

A.1 $Q$ as dark energy

We define slow-roll parameters for $Q$ as

[TABLE]

which can be rewritten as

[TABLE]

In the limit $\epsilon_{Q}\ll 1$, we have $w_{Q}\approx-1$. The additional condition $\eta_{Q}\ll 1$ can be imposed to keep $\epsilon_{Q}\ll 1$ but is not necessary for our purposes. In particular, the example solution used in the main text of the paper has $\eta_{Q}\approx\frac{3}{2}(1+c_{Q}^{2})(1+w_{tot})\sim\mathcal{O}(1)$ during the matter and radiation eras while $\epsilon_{Q}\ll 1$ throughout.

In what follows it will be convenient to rewrite the continuity equation in terms of the sounds speeds in Eq. (92) and Eq. (93) and the slow-roll parameters in Eq. (96) and Eq. (97). This gives an expression for the change in energy due to change in the field strength as

[TABLE]

A.2 Inflationary generation of $\delta Q$

We assume that the field $Q$ is a spectator field during inflation and that mixing between the inflaton and $Q$ can be ignored. The field $Q$ then acquires perturbations $\delta Q$ that will have power spectrum Fujita et al. (2015),

[TABLE]

where the subscript $i$ indicates these quantities are to be evaluated during inflation. The $\delta Q$ will be our isocurvature modes. For the isocurvature modes to be correlated with $\mathcal{R}$, we assume at some time after inflation fluctuations in the energy density of $Q$ create adiabatic fluctuations via some mechanism with an efficiency factor $\gamma$,

[TABLE]

where $Q_{i}$ is the value of the field at the end of inflation. The subscript $Q$ in Eq. (100) indicates that these perturbations were inherited from $\delta Q$ after inflation.

In addition, there will be adiabatic perturbations generated during inflation, we label these perturbations $\mathcal{R}_{i}$. For simplicity we assume slow-roll inflation so that

[TABLE]

where $\phi$ is the inflaton and the power spectrum of $\delta\phi$ is

[TABLE]

so that

[TABLE]

where $\epsilon_{i}$ is the usual inflationary slow-roll parameter $\epsilon_{i}\equiv-H_{i}^{\prime}/H_{i}$.

The net curvature perturbation used in the body of this paper is the sum of the two components,

[TABLE]

and its power spectrum is

[TABLE]

where we have defined $\xi$ as the ratio of curvature perturbations generated by the inflaton, to curvature perturbations generated by $\delta Q$,

[TABLE]

The isocurvature perturbations used in this paper are

[TABLE]

where $c_{s}^{2}$ is the synchronous gauge sound speed, defined in Eq. (10) and in the last line we have used the fact that outside of the horizon, $c_{s}^{2}\simeq-1$.

The isocurvature auto- and cross-power spectra are then

[TABLE]

The ratio of isocurvature to adiabatic perturbations is

[TABLE]

and the correlation between them is given by

[TABLE]

So, to have substantial correlations between $\mathcal{R}$ and $\mathcal{I}$, we need \xi\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}1, and to have the isocurvature also be important ($\Delta_{\mathcal{R}\mathcal{R}}^{2}\sim\Delta_{\mathcal{I}\mathcal{I}}^{2}$) we need

[TABLE]

On the other hand, if we do not care about correlations between $\mathcal{I}$ and $\mathcal{R}$, $\xi$ can be large and we only need

[TABLE]

to have isocurvature effects today.

From Eq. (107), we can see that both correlated and uncorrelated cases require

[TABLE]

This can be rewritten in terms of $\Omega_{Q}$, and the slow-roll parameters, using Eq. (98) and the Friedmann equation, as

[TABLE]

where the subscript $i$ indicates that these quantities are evaluated during inflation, and quantities without this subscript are evaluated at late times.

For a standard kinetic term, Eq. (118) becomes

[TABLE]

which cannot be satisfied by a field that dominates the energy density today and behaves as dark energy.

On the other hand, with a nonstandard kinetic term it becomes easier to satisfy Eq. (118). For the Lagrangian in Eq. (1) we have

[TABLE]

For the solution we considered in this paper, $\left({X_{i}}/{X}\right)^{\frac{1-c_{Q}^{2}}{2c_{Q}^{2}}}\sim H/H_{i}$, which allows the left-hand side of Eq. (120) to be small even if, $\Omega_{Q}\sim 1$. The factor on the right-hand side involving $\eta_{Q}$ is given by the background solution,

[TABLE]

where $f$ is defined in Eq. (6). During matter domination, this gives $1+\eta_{Q}/(3(1+c_{Q}^{2}))=\frac{3}{2}$. This factor, combined with other factors of $f$ from the expression for $X$ evaluated on the background solution, combine to a numerical constant of order $\sim 10$ at most on the left-hand side of the above inequality.

Finally, note that the tensors generated during inflation are

[TABLE]

so the tensor-to-scalar ratio is

[TABLE]

So, if $\xi\gg 1$ (inflaton perturbations dominate) we have the usual expression for the tensor-to-scalar ratio, and as $\xi$ decreases, the tensor-to-scalar ratio gets even smaller. Note that even secondary production of gravitational waves by a spectator field like $Q$ cannot enhance the tensor-to-scalar ratio Fujita et al. (2015).