The First Constraint from SDSS Galaxy–Galaxy Weak Lensing Measurements on Interacting Dark Energy Models

The current standard cosmological model, ΛCDM model, is widely accepted in explaining various astronomical observations (Begeman et al. 1991; Persic et al. 1996; Perlmutter et al. 1998, 1999; Riess et al. 1998; Chemin et al. 2011; Planck Collaboration et al. 2016). However, recently some observational tensions have been reported if the universe is described by the ΛCDM model. It was found that there is about 3σ mismatch for the Hubble constant inferred from the cosmic microwave background (CMB) measurements and from the direct local observations if the ΛCDM model is assumed (Riess et al. 2011, 2016). Besides, the Baryon Oscillation Spectroscopic Survey experiment showed that there is a 2.5σ deviation from the ΛCDM model in the measurement of the Hubble constant and angular distance at an average redshift z = 2.34 (Delubac et al. 2015). Furthermore, a "substantial discordance" at the level of 2.3σ was obtained between the weak lensing data taken from a 450 − deg² observing field of the Kilo Degree Survey and the Planck 2015 CMB data (Hildebrandt et al. 2017) if the ΛCDM model is supposed. Besides observational challenges, the ΛCDM model faces serious theoretical problems, such as the cosmological constant problem (Weinberg 1989) and the coincidence problem (Zlatev et al. 1999). This motivates us to find some more viable models to describe our universe.

In the framework of Einstein gravity, nearly 95% of the universe content is composed of dark matter (DM) and dark energy (DE). From the field theory point of view, it is a specific assumption that DM and DE live independently in the universe. More naturally we can consider some interactions between these two biggest components. The interaction between DM and DE has been discussed extensively in the literature, for a recent review please see Wang et al. (2016) and references therein. It is interesting to find that appropriate interaction between dark sectors can relieve discordances in observations as previously inferred from the ΛCDM model (Costa et al. 2017; Ferreira et al. 2017; An et al. 2018). Moreover, the coincidence problem can be alleviated if there is a proper interaction between DM and DE (He et al. 2011).

The influence of interacting dark energy (IDE) models on the background dynamics and the linear perturbation evolutions in the universe has been studied extensively, see the review of Wang et al. (2016) and the references therein. In the nonlinear regime, N-body simulations are essential to understanding the structure formation and evolution. A preliminary attempt on the N-body simulation by considering quintessence DE interacting with DM was proposed in Baldi et al. (2010) and Baldi (2011a), where the initial condition in the simulation was naively taken from the ΛCDM model and the DE perturbation was not consistently computed at different scales and redshifts. For general phenomenological IDE models, self-consistent N-body simulations are still lacking.

In this Letter, we devise a novel cosmological N-body simulation pipeline for cosmological models beyond ΛCDM. We consider self-consistent initial conditions for IDE models and include DE distributions from directly solving perturbation equations. We do not limit the DE in the quintessence region and consider general DE fluid phenomenologically interacting with DM. We apply our simulation pipeline to four types of IDE models (Wang et al. 2016; Costa et al. 2017) and try to explore the physics in the structure formation when there are interactions between dark sectors. With this self-consistent and effective pipeline, we open a new window to precisely study the nonlinear structure formation in IDE models at low redshifts. This enables us to employ a new probe, weak lensing, to put further constraints on IDE models. For the first time, we use the galaxy–galaxy weak lensing measurements from the SDSS data (Luo et al. 2017, 2018) to compare with our simulation results. We find that the improvement of the constraint for the interaction strength in some IDE models can reach up to 1250%. This shows the power of our cosmological N-body simulation pipeline in studying IDE models. With this tool, we can refine IDE models allowed by linear constraints (Costa et al. 2017). Finally, we are able to look for smoking gun signatures of interactions between the dark sectors, in the simulations constrained by current observations.

The interaction between dark sectors is well motivated from field theory and is widely discussed in the literature, see the recent review of Wang et al. (2016). With the interaction between dark sectors, the background continuity equations of DM and DE obey

$$\begin{eqnarray}&&\dot{{\rho }_{c}}+3H{\rho }_{c}=Q,\dot{{\rho }_{d}}+3H(1+{w}_{d}){\rho }_{d}=-Q.\end{eqnarray} \tag{ 1 }$$

Here we focus on the commonly assumed phenomenological interaction form $Q=3{\xi }_{1}H{\rho }_{c}+3{\xi }_{2}H{\rho }_{d}$, where ρ_c is the DM density, ρ_d is the DE density, and the dot denotes the derivative with respect to the conformal time, H is the Hubble parameter and ${w}_{d}={p}_{d}/{\rho }_{d}$ is the equation of state for DE. We do not limit DE to be a quintessence field with w_d > −1 (Baldi 2011b), but instead allow w_d to be a free value either in the quintessence or the phantom regions. ξ₁ and ξ₂ indicate the strength of interactions.

The linear evolutions of density and velocity perturbations for DM and DE were described in He et al. (2009a, 2009b); Zhang et al. (2018). In the subhorizon approximation, from linear level equations we can obtain the Poisson equation in the real space

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Psi }}=-\displaystyle \frac{3}{2}{H}^{2}[{{\rm{\Omega }}}_{c}{{\rm{\Delta }}}_{c}+(1-{{\rm{\Omega }}}_{c}){{\rm{\Delta }}}_{d}],\end{eqnarray} \tag{ 2 }$$

where ${{\rm{\Delta }}}_{d}$ (Δ_c) is the density perturbation of DE (DM), and Ω_c is the background density ratio of DM. It is clear that with the interaction, the gravitational potential is modified. The corresponding Euler equation in the real space reads

$$\begin{eqnarray}&&{\rm{\nabla }}\dot{{v}_{c}}+\left[H+3H\left({\xi }_{1}+\displaystyle \frac{{\xi }_{2}}{r}\right)\right]{\rm{\nabla }}{v}_{c}={{\rm{\nabla }}}^{2}{\rm{\Psi }},\end{eqnarray} \tag{ 3 }$$

where $r={\rho }_{c}/{\rho }_{d}$. The coupling between dark sectors introduces an additional acceleration on DM particles at each time step in the simulation. In the following we will concentrate our discussions on phenomenological IDE models listed in Table 1 (Costa et al. 2017), which are natural Taylor expansions of the interaction kernel Q into energy densities ρ_c and ρ_d. With Planck 2015, SNe Ia, baryon acoustic oscillations, and the Hubble constant observation, tight constraints on ξ₁ for Model III and ξ for Model IV were obtained in Costa et al. (2017). However, for Models I and II, the obtained constraints on the strengths of couplings are loose (Costa et al. 2017). This is expected because in Models I and II the interaction is proportional to the energy density of DE, which was subdominant when CMB was produced. It is more reasonable to expect that the observations at low redshifts shall provide tighter constraints on Models I and II, especially the small-scale structure information. For this purpose we resort to using N-body simulations to make an accurate analysis.

Table 1.  Phenomenological IDE Models

Model	Q	w
I	$3{\xi }_{2}H{\rho }_{d}$	$-1\lt {w}_{d}\lt -1/3$
II	$3{\xi }_{2}H{\rho }_{d}$	${w}_{d}\lt -1$
III	$3{\xi }_{1}H{\rho }_{c}$	${w}_{d}\lt -1$
IV	$3\xi H({\rho }_{c}+{\rho }_{d})$	${w}_{d}\lt -1$

Download table as:  ASCIITypeset image

Because IDE models are different from the ΛCDM model in every relevant equation, it is naive to count on empirical fits to the ΛCDM model, e.g., the halofit, to understand the physics in the nonlinear structure formation. We require a new N-body simulation pipeline to understand the structure developed in IDE models. There are four modifications we have considered in devising the new pipeline compared to the standard ΛCDM model. First, the preinitial condition is generated by the Capacity Constrained Voronoi Tessellation (CCVT) method (Liao 2018), instead of the classically used glass or grid. This makes sure that our preinitial condition is free of the Poisson equation, generating a geometrical equilibrium state of particle distribution. We have tested this using CCVT; grid or glass makes negligible differences for the preinitial condition in simulations. The choice of CCVT is mainly because of self-consistency consideration, rather than the accuracy consideration. Second, the initial matter power spectrum is generated by our modified CAMB (Lewis & Bridle 2002; Costa et al. 2017) with the coupling between dark sectors, which is different from the ΛCDM model. Third, the perturbations of the particle distribution are calculated by using 2LPTic (Crocce et al. 2006), which is properly modified to be consistent with our models. Fourth, the N-body simulation code is also heavily modified for consistency.

Instead of treating the DE perturbation as a constant excess of gravity at all scales and redshifts (Baldi et al. 2010; Baldi 2011a), we include the DE perturbation self-consistently as a function of scale (k) and redshift (z) by solving perturbation equations from the modified CAMB (Costa et al. 2017). We also modify the N-body simulation code Gadget2 (Springel 2005) into ME-Gadget. Technical details can be found in Zhang et al. (2018). We find that ME-Gadget is as efficient as the original Gadget2 code, and the testing results are consistent with Baldi (2011a) by using their models. Our convergence test results also show that our code can reach 5% accuracy as k is approaching the Nyquist limit for the nonlinear matter power spectrum at z = 0. We would like to emphasize that our N-body simulation pipeline is fully self-consistent, accurate, and efficient enough for general phenomenological IDE models. The simulation parameters we use are shown in Table 2, which were constrained from the combination of Planck 2015, SNe Ia, baryon acoustic oscillations, and the Hubble constant observation data sets (PBSH in short hereafter; Costa et al. 2017). We use a comoving box size of 400h⁻¹ Mpc and 256³ particles in our computations for the matter power spectrum. A comoving box size of 400h⁻¹ Mpc and 512³ particles is used for the measurements of galaxy–galaxy lensing signals in the simulations.

Table 2.  Cosmological Parameters

Parameter	IDE_I	IDE_II	IDE_III	IDE_IV	ΛCDM
${{\rm{\Omega }}}_{b}{h}^{2}$	0.02223	0.02224	0.02228	0.02228	0.02225
${{\rm{\Omega }}}_{c}{h}^{2}$	0.0792	0.1351	0.1216	0.1218	0.1198
$100{\theta }_{{MC}}$	1.043	1.04	1.041	1.041	1.04077
τ	0.08204	0.081	0.07728	0.07709	0.079
$\mathrm{ln}({10}^{10}{A}_{s})$	3.099	3.097	3.088	3.087	3.094
ns	0.9645	0.9643	0.9624	0.9624	0.9645
w	−0.9191	−1.088	−1.104	−1.105	−1
ξ1	⋯	⋯	0.0007127	0.000735	⋯
ξ2	−0.1107	0.05219	⋯	0.000735	⋯
H0	68.18	68.35	68.91	68.88	67.27
Ωm	0.2204	0.3384	0.3045	0.3053	0.3156

Download table as:  ASCIITypeset image

The matter power spectrum is used to quantify the large-scale structures. The linear evolution of the matter power spectrum can be simply calculated by linear growth theory. People usually use halofit (Takahashi et al. 2012) to estimate the nonlinear matter power spectrum at low redshifts. However, because halofit is an empirical fit to ΛCDM N-body simulations, it is not appropriate to use it to describe IDE models. To clarify, we compare the measured matter power spectrum from our N-body simulations with the prediction of halofit in Figure 1. The matter power spectrum is computed using the ComputePk code (L'Huillier 2014). We find that halofit can marginally be used to describe the nonlinear matter power spectra at z = 0 for IDE_III and IDE_IV models, although it is not exactly consistent with that from N-body simulations. This is because the strengths of the interaction in these two models are quite small (∼0.0007) so that the deviations from the ΛCDM model are negligible. However, for models IDE_I and IDE_II, it is clear that halofit cannot give the true matter power spectrum, especially at small scales, because their interactions are relatively large ($| {\xi }_{2}| \gt 0.05$), which causes large deviations from the standard ΛCDM model. The empirical fit to ΛCDM is no longer appropriate in these cases to describe the nonlinear structure, and the appropriate N-body simulations pipeline is called for.

Zoom In Zoom Out Reset image size

The galaxy images are distorted by the foreground gravitational potential, which is known as gravitational lensing. Assuming an isotropic distribution of both galaxy shape and orientation, any nonzero residual can be considered as such an effect, a.k.a. tangential shear γ_T. In galaxy–galaxy lensing, the signal is usually interpreted as the combination of γ_T and the geometry of a lensing system, referring to the critical density ${{\rm{\Sigma }}}_{\mathrm{crit}}({z}_{l},{z}_{s})=\displaystyle \frac{{c}^{2}}{4\pi G}\tfrac{{D}_{s}}{{D}_{{ls}}{D}_{l}}$, where ${z}_{l},{z}_{s}$ denote the redshifts of the lens and the source, D_l, D_s, and D_ls are the angular diameter distances of the lens, source galaxy, and the difference between them. The signal measured from galaxy–galaxy lensing actually reflects the differential change of 2D surface density-excess surface density (ESD),

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Sigma }}(R)={\rm{\Sigma }}(\leqslant R)-{\rm{\Sigma }}(R)={\gamma }_{t}{{\rm{\Sigma }}}_{\mathrm{crit}}({z}_{l},{z}_{s}),\end{eqnarray} \tag{ 4 }$$

here Σ( ≤ R) is the average surface density inside the projected distance R and Σ(R) is the surface density at the projected distance R. This signal is multiplied by a factor, i.e., a boost factor, to correct for the contamination by galaxies associated with lens galaxy. This factor is calculated following Mandelbaum et al. (2005)

$$\begin{eqnarray}&&B(R)=\displaystyle \frac{n(R)}{{n}_{\mathrm{rand}}(R)},\end{eqnarray} \tag{ 5 }$$

where n(R) and n_rand(R) are the background numbers around lens sample at radius R.

We use the shear catalog from Luo et al. (2017), which is based on the SDSS DR7 image data. For groups of galaxies, we employ the catalog from Yang et al. (2007) to identify the lens systems. Following the galaxy–galaxy lensing measurement procedure in Luo et al. (2018), we select the most luminous 3660 galaxy groups in the group catalog from redshift 0.01–0.2 as the lens. To estimate the abundance, we calculate the comoving volume of the SDSS DR7 north cap between redshifts 0.01 and 0.2, labeled as V_com, with completeness considered. The completeness of each galaxy is given by the NYU-VAGC (Blanton et al. 2005). The number of halos with same abundance in the simulation is then estimated as $3660\times {400}^{3}/{V}_{\mathrm{com}}\approx 1771$. We stack the tangential shear of these 3660 lens systems to measure the ESD. Taking the halos of the same abundance as that in the observation in our N-body simulations, we select the most massive 1771 halos and stack their particles to measure the ESD. The systematics introduced by photometric redshift is about 2.5%–2.7% estimated based on Equation (24) in Mandelbaum et al. (2005). About 5% of the galaxies are satellites according to Cacciato et al. (2009), the contribution to the ESD is roughly about 4% at $300\,{h}^{-1}\,\mathrm{kpc}$ to 10% contribution to the maximum at $1000\,{h}^{-1}\,\mathrm{kpc}$. These uncertainties and bias are neglected for the analysis below because they are too small to affect our final results.

The weak lensing measurements and simulation predictions of each cosmological model are shown in Figure 2. By comparing the ESD curves from the 256³ and 512³ simulations, we find that the 256³ results converge to the 512³ ones with a level of <5% at r > 600 kpc for all three cosmological models (LCDM, IDE_I, and IDE_II). Therefore, we expect that the ESD curves from the 512³ simulations shown in Figure 2 should have a convergence level of <5% at r > 300 kpc for all models. We find that the measured data points are systematically lower than the prediction from the ΛCDM model shown by the black dashed line, which is mainly due to the Eddington bias (Luo et al. 2018). The Eddington bias comes from the incorrect estimation of the halo mass using the galaxy luminosity or other indicators. The incorrect estimation will mistakenly identify lower mass halos as higher mass halos, thus contaminate the ESD signal. We corrected the Eddington bias by assuming a 0.3 dex scatter in mass–luminosity relation following Luo et al. (2018), shown as the solid lines. We have tested that the Eddington bias introduced in Luo et al. (2018) is similar for IDE models by using the halo catalogs from our simulations. The shaded area represents the dispersion due to finite width of redshift bin. The groups of galaxies we selected locate at different redshifts, central at z = 0.15 (range 0.01 < z < 0.2). Thus, the uncertainty due to the redshift difference was also taken into account in our analysis. We estimate the redshift bin by measuring the ESD signal from simulation snapshots at z = 0.1 and z = 0.2 separately. The solid lines show the central value of the shaded area with the same colors. Even with such conservative treatments, it is still quite clear that IDE_I and IDE_II are not favored by the SDSS galaxy–galaxy weak lensing data, even though these two IDE models are well constrained by PBSH. Therefore, tight constraints from comparing our simulations with observational galaxy–galaxy lensing signals are expected.

Zoom In Zoom Out Reset image size

We estimate the constraints from galaxy–galaxy lensing signals by assuming that the ESD signal deviation from the ΛCDM model in logarithmic space is linearly proportional to the interaction strength. We have tested that the above assumption is not significantly affected by the choice of logarithmic space or linear space. We have also tested that such an assumption is reasonably accurate using multiple simulations with different parameters. The likelihood is constructed as

$$\begin{eqnarray}&&L=\exp \left\{-\displaystyle \frac{1}{2}\displaystyle \sum _{i}\displaystyle \frac{{\left[{\rm{\Delta }}{\rm{\Sigma }}{({R}_{i})}_{\mathrm{sim}}-{\rm{\Delta }}{\rm{\Sigma }}{({R}_{i})}_{\mathrm{obs}}\right]}^{2}}{{\sigma }_{z}^{2}+{\sigma }_{\mathrm{obs}}^{2}}\right\}.\end{eqnarray} \tag{ 6 }$$

Here R_i denotes the measured five data points, σ_z = 0.288 times the width of the shaded area, representing the uncertainty due to the finite width of the redshift bin, and σ_obs is the error estimated from the lensing signal. We show the likelihood from our comparison in Figure 3. Comparing to the linear constraints given by Costa et al. (2017) shown in dashed lines, the constraints from our SDSS galaxy–galaxy weak lensing (SDSS-WL) are clearly tighter for models IDE_I and IDE_II. The improvements of constraints for models IDE_III and IDE_IV are negligible and we do not show them here. The joint likelihood of PBSH and SDSS-WL is about 1250% tighter than PBSH alone for IDE_I. The best-fitted ξ₂ for IDE_I is ξ₂ = 0, and the best-fitted ξ₂ for IDE_II becomes ξ₂ = 0.012.

Zoom In Zoom Out Reset image size

We have successfully devised a self-consistent N-body simulation pipeline to examine the influence of the interaction between dark sectors in structure formation at low redshifts. This formalism is appropriate for general IDE models and efficient in examining the signature of the interaction. With this tool in hand, we do not need to blindly count on halofit, which is an empirical fit to the ΛCDM model, to disclose nonlinear structures.

Considering that interactions in IDE_I and IDE_II models are proportional to the energy density of DE, which is subdominant at high redshifts, it is natural to find that the constraints of these interactions from PBSH are loose. With the self-consistent N-body simulation pipeline, however, we can examine these two models more carefully by using the nonlinear low-redshift observations, such as the SDSS galaxy–galaxy weak lensing. It is interesting to find that our first try of the pipeline can obtain up to 1250% improvement in the interaction strength constraint for the IDE_I model. Combining PBSH and SDSS galaxy–galaxy weak lensing measurements, we find the constraint of the interaction strength ${\xi }_{2}={0}_{-0.0105}^{+0.0}$ for the IDE_I model. The ΛCDM model is still favored. For the IDE_II model, combining PBSH and SDSS galaxy–galaxy weak lensing data sets, we obtain ${\xi }_{2}={0.0120}_{-0.012}^{+0.0166}$, which is also improved significantly by including the nonlinear structure information. For IDE models III and IV, the galaxy–galaxy lensing constraints by employing N-body simulations do not improve much of the constraints if we compare with the linear PBSH results. It is interesting that our pipeline is effective in disclosing physics in the structure formation when there is coupling between dark sectors and it can also help to refine IDE models. We would like to address that the likelihood and improvement from SDSS-WL is only a rough estimation about the parameter rather than a complete constraint. It will be useful in guiding our future study.

By combining the linear and nonlinear scale information, we can not only constrain the interaction strength between dark sectors with much higher precision, but also perform the simulations constrained by the current observations. In such simulations, we can look for smoking gun signatures of the DM DE interactions. These signatures can be directly tested by observations in the future.

J.Z. acknowledges the support from China Postdoctoral Science Foundation 2018M632097. W.L. acknowledges the support from NSFC 11503064 and Shanghai Jiao Tong University and University of Michigan Joint Foundation (AF0720054). The work of B.W. was partially supported by NNSFC. This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), Japan. We also thank Surhud More from IUCAA, India, and Feng Shi from KASI, Korea, for useful discussions.