# Constraints on Dynamical Dark Energy Models from the Abundance of Massive Galaxies at High Redshifts

N. Menci, A. Grazian, M. Castellano, P. Santini, E. Giallongo, A. Lamastra, F. Fortuni, A. Fontana, E. Merlin, T. Wang, et al.

► To cite this version:

N. Menci, A. Grazian, M. Castellano, P. Santini, E. Giallongo, et al., Constraints on Dynamical Dark Energy Models from the Abundance of Massive Galaxies at High Redshifts. The Astrophysical Journal, 2020, 900 (2), pp.108. 10.3847/1538-4357/aba9d2, hal-02987500

HAL Id: hal-02987500

https://hal.science/hal-02987500v1

Submitted on 9 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Constraints on Dynamical Dark Energy Models from the Abundance of Massive Galaxies at High Redshifts

N. Menci<fnref n="1" />, A. Grazian<fnref n="2" />, M. Castellano<fnref n="1" />, P. Santini<fnref n="1" />, E. Giallongo<fnref n="1" />, A. Lamastra<fnref n="1" />, F. Fortuni<fnref n="1" />, A. Fontana<fnref n="1" />, E. Merlin<fnref n="1" />, T. Wang<fnref n="3" />, D. Elbaz<fnref n="4" />, and N. G. Sanchez<fnref n="5" />

## Abstract

We compare the maximal abundance of massive systems predicted in different dynamical dark energy (DDE) models at high redshifts $z \approx 4-7$ with the measured abundance of the most massive galaxies observed to be already in place at such redshifts. The aim is to derive constraints for the evolution of the dark energy equation-of-state parameter $w$ that are complementary to existing probes. We adopt the standard parameterization for the DDE evolution in terms of the local value $w_0$ and of the lookback time derivative of the equation of state $w_a$. We derive constraints on combinations of ($w_0, w_a$) in the different DDE models by using three independent probes: (i) the observed stellar mass function of massive objects at $z \ge 6$ derived from the CANDELS survey; (ii) the estimated volume density of massive halos derived from the observation of massive, star-forming galaxies detected in the submillimeter range at $z \approx 4$; (iii) the rareness of the most massive system (with an estimated gas mass exceeding $3 \times 10^{11} M_{\odot}$) observed to be in place at $z \approx 7$, a far-infrared-luminous object recently detected in the South Pole Telescope survey. Finally, we show that the combination of our results from the three above probes excludes a sizable fraction of the DDE parameter space $w_a \gtrsim -3/4 - (w_0 + 3/2)$ presently allowed (or even favored) by existing probes.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Cosmological parameters (339); Dark energy (351)*

## 1. Introduction

The current theory of structure formation envisages all cosmic structures to form from the collapse and the growth of initially tiny density perturbations of the dark matter (DM) density field in a universe characterized by an accelerated expansion. Such an acceleration indicates that the dominant component (with density parameter $\Omega_{\Lambda} \approx 0.7$) of the cosmic fluid must be composed of some form of dark energy (DE), with equation-of-state parameter $w \equiv p/\rho \le -1/3$. Although the nature of such a component remains unknown, the simplest model assumes DE to be connected with the vacuum energy, the so-called cosmological constant, with equation-of-state parameter $w = -1$. When coupled with the assumption that DM is composed of nonrelativistic particles at decoupling, such a scenario leads to the $\Lambda$CDM standard cosmological model.

While measurements of the Cosmic Microwave Background (CMB) have provided a first, strong confirmation of such a scenario, tensions have recently emerged (mostly between measurements related to the early and late universe): these include the inadequacy of the model to provide a perfect fit to the Planck CMB temperature and polarization angular spectra (see, e.g., Addison et al. 2016), the discrepancy between the combined values of the power spectrum normalization $\sigma_8$ and matter density parameter $\Omega_M$ derived by Planck with respect to those derived from cosmic shear surveys such as CFHTLenS (Heymans et al. 2012) and KiDS-450 (Hildebrandt et al. 2017), and—most of all—the tension (at a more than $3\sigma$ confidence level) in the Hubble constant $H_0$ between the values derived from Planck ($h \approx 0.67$, in units km s$^{-1}$ Mpc$^{-1}$) and those

obtained from local luminosity distance measurements ($h \approx 0.74$, see Riess 2019 and references therein).

Such tensions have stimulated an extended effort toward the investigation of more complex cosmological models. One of the simplest physical alternatives is comprised of a DE with a time-dependent equation of state (dynamical dark energy, DDE). In fact, this kind of scenario is a possible solution to the abovementioned tension between the values of the Hubble constant derived from local indicators and from the CMB (Di Valentino et al. 2017; Pan et al. 2019). In addition, a constant equation of state is not expected in physically motivated scenarios in which DE originates from a scalar “quintessence” field $\phi$ (Peebles & Ratra 1988; Caldwell et al. 1998; Sahni & Wang 2000; Copeland et al. 2006; Frieman et al. 2008) evolving in a potential $V(\phi)$. In fact, its pressure $p_{\phi} = \dot{\phi} - V(\phi)$ and energy density $\rho_{\phi} = \dot{\phi} + V(\phi)$ lead to a DE with a time-evolving equation-of-state parameter $w \equiv p_{\phi}/\rho_{\phi}$. Parameterizing the evolution of $w$ with the expansion factor $a$ as $w(a) = w_0 + w_a(1-a)$; (Chevallier & Polarski 2001; Linder 2003), the dynamics of such models can be related to different combinations of ($w_0, w_a$); (see Caldwell & Linder 2005; Barger et al. 2006; Linder 2006). For example, “thawing” models (Scherrer & Sen 2008; Chiba 2009; Gupta et al. 2015) are characterized by the growth of $w$ with time starting from $w = -1$ in the early universe, while “freezing” models (Chiba 2006; Scherrer 2006; Sahlen et al. 2007) predict a decrease in $w$ with $a$, approaching a cosmological constant value $w = -1$.

Present observational constraints on the DDE parameter space $w_0-w_a$ (see, e.g., Zhai et al. 2017) provide contrasting

<footnote n="1">¹INAF—Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio, Italy</footnote>

<footnote n="2">²INAF—Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy</footnote>

<footnote n="3">³National Astronomical Observatory of Japan, Mitaka, Tokyo, Japan</footnote>

<footnote n="4">⁴AIM, CEA, CNRS, Universitéé Paris-Saclay, Universitéé Paris Diderot, Sorbonne Paris Citeé, Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵LERMA, CNRS UMR 8112, 61, Observatoire de Paris PSL, Sorbonne Universités, UPMC Univ. Paris 6, 61 Avenue de l'Observatoire, F-75014 Paris, France

*Received 2020 April 17; revised 2020 June 19; accepted 2020 July 24; published 2020 September 7*</footnote>
results. On the one hand, results coming from the CMB power spectrum (Ade et al. 2016) and weak lensing tomography (Amara & Refregier 2007; Massey et al. 2007; see Refregier 2003 for a review)—although they leave unconstrained a relatively large volume of the DDE parameter space—disfavor combinations of $w_0-w_a$ with $w_0 \lesssim -1.5$ and positive values of $w_a \gtrsim 0.7$ (especially when combined with constraints from baryonic acoustic oscillations and from supernovae (SNe) Ia; see Ade et al. 2016; Di Valentino et al. 2017; Scolnic et al. 2018 and references therein), as well as most of the combinations with $w_0 \gtrsim -1$ and $w_a \gtrsim 0.5$. On the other hand, the recent determination of the Hubble diagram of quasars in the range $0.5 \le z \le 5.5$ (Risaliti & Lusso 2019) favors large values of $w_a \gtrsim 0$ with negative values of $w_0 < -1$, which deviates from the $\Lambda$CDM model with a statistical significance of $4\sigma$. This measurement is based on quasar distances estimated from the ratio between their X-ray and ultraviolet emission. Although it strongly relies on the assumed invariance with redshift of the X-ray-to-ultraviolet ratio, the excellent agreement with the Hubble diagram derived from SNe Ia in the overlapping range $0.5 \le z \le 1.5$ and the evidence of non-volving UV and X-ray spectral properties strongly argue for the reliability of such results.

In such a context, a further, independent probe of the nature and evolution of DE is the evolution of the galaxy population over cosmic time. In fact, the inverse dependence of the amplitude of initial density perturbations on the mass scale (measured from fluctuations in the CMB, see, e.g., Tegmark & Zaldarriaga 2002, 2009; Aghanim et al. 2019) implies that the formation of galactic DM halos proceeds bottom-up. Although the physics of the assembly of baryons into the DM halos is a complex issue, a solid consequence of the above scenario is that—in any specific adopted cosmological model—large-mass DM halos must become progressively rarer with increasing redshift. Thus, viable cosmological models must allow for an evolution of the initial density perturbations that is fast enough to match the abundance of massive galaxies observed to be in place early on in history of the universe. Indeed, several observations concerning massive galaxies at high redshifts are already challenging the canonical $\Lambda$CDM cosmological model. For example, works by several authors (see, e.g., Hildebrandt et al. 2009; Lee et al. 2012; Caputi et al. 2015; Finkelstein et al. 2015; Merlin et al. 2019; see also Wang et al. 2019) have increased the tension between the expected evolution of the DM halo mass function and the observed galaxy luminosity and mass functions at $z \gtrsim 4$. While the present understanding of the baryonic processes leading to gas condensation and to star formation in DM halos struggles to describe the rapid evolution of the star formation that is needed to match the observed abundance of massive galaxies (see Steinhardt et al. 2016), an enhanced efficiency in converting baryons into stars at high redshifts could still allow for consistency between the $\Lambda$CDM predictions and the observed number density of luminous, massive galaxies (Behroozi et al. 2013; Behroozi & Silk 2015, 2018; Finkelstein et al. 2015; Sun & Furlanetto 2016; Moster et al. 2018).

The above degeneracy between baryonic effects and cosmology in determining the expected abundance of luminous galaxies can be bypassed by noticing that the ratio of galaxy baryonic components (stellar mass or gas mass) to DM halo mass has an absolute maximum at the cosmic baryon fraction $f_b \approx 0.16$ (Aghanim et al. 2019). In fact, the observed

abundance of galaxies with large mass in the baryonic component $M_b$ places a lower limit on the abundance of DM haloes with masses $M \gtrsim M_b/f_b \approx 6.3 M_b$. Such a constraint can be used to rule out cosmological models that do not allow for a sufficiently rapid growth of galactic DM halos. In fact, an observed abundance $\phi_{\text{obs}}(M_b, z)$ would rule out any cosmological models that predict a number density of DM haloes $\phi(M \gtrsim M_b/f_b, z) \le \phi_{\text{obs}}(M_b, z)$, independently of the details of the complex baryon physics involved in the galaxy formation process.

In this paper we apply such a probe to cosmological models based on DDE. We compare the maximal abundance of massive galaxies predicted in different DDE models at high redshifts with the measured abundance of the most massive systems observed to be already in place at the same redshifts.

The paper is organized as follows:

1. In Section 2 we present the method adopted to derive the halo mass function in different DDE models, and we describe how we compute the basic quantities that we will compare with observational data.

2. Such a comparison is performed in Section 3 for three different observations concerning the abundance of massive objects at high redshift. We first compare the observed stellar mass function of massive objects at $z \gtrsim 6$ as observed by the CANDELS survey with the halo mass function predicted in different DDE models (Section 3.1), deriving exclusion plots in the DDE parameter space $w_0-w_a$. In Section 3.2, we perform a comparison with the estimated volume density of massive halos as derived from the observation of massive, star-forming galaxies detected in the submillimeter at $z \gtrsim 4$ (Wang et al. 2019), which are expected to reside in the most massive DM haloes at their redshift. Finally, in Section 3.3, we consider the most massive object (with an estimated gas mass exceeding $3 \times 10^{11} M_\odot$) in place at $z \approx 7$ recently detected in the South Pole Telescope (SPT) survey, and we derive (adopting the most conservative assumptions for the baryon-to-DM ratio) the probability for such a massive object to be present in the area covered by SPT for different DDE models for different combinations ($w_0, w_a$). The constraints derived from the combination of the above observations are shown in Section 3.4.

3. The final Section 4 is devoted to a discussion and conclusions.

## 2. Method

To compute the expected abundance of DM haloes in different DDE we adopt the canonical Press and Schechter approach, which relates the number of halos of mass $M$ at redshift $z$ to the overdense regions of the linear density field with density contrast $\delta$ over a background average density $\bar{\rho}$. Successful models predict the halo mass function $\phi$ (i.e., the number of virialized DM haloes with mass in the range $M - M + dM$ per unit volume) to take the form

$$ \phi(M) = \frac{dN}{dM} = \frac{\bar{\rho}}{M^2} \frac{d \ln \nu}{d \ln M} f(\nu). \quad (1) $$

Here, $\nu = \delta_c/\sigma(M, z)$, where $\delta_c$ corresponds to the critical linear overdensity (equal to 1.69 in the $\Lambda$CDM scenario) and $\sigma(M, z)$ is the variance of the linear density field smoothed on
the scale $R = [3M/4\pi\bar{\rho}]^{1/3}$, and evolving with time according to the linear growth factor $D(z)$ of density perturbations. The function $f(\nu)$ is universal to the changes in redshift and cosmology. In the original Press-Schechter theory (Press & Schechter 1974) and in excursion set theory (Bond et al. 1991), the function $f$ takes the form $f_{\text{PS}} = (\sqrt{2/\pi})^{1/2}\nu \exp(-\nu^2/2)$ appropriate for spherical collapse. More recent approaches (Sheth & Tormen 1999; Jenkins et al. 2001; Warren et al. 2006; Tinker et al. 2008) provided more accurate forms that have been extensively tested against N-body simulations. Here, we adopt the form given by Sheth & Tormen (1999):

$$f(\nu) = 2A\left(\frac{1}{\nu'^{2q}} + 1\right)\frac{\nu'^2}{2\pi}e^{-\nu'^2/2} \quad (2)$$

with $\nu' = \sqrt{a}\nu$, $a = 0.71$, and $q = 0.3$. The normalization factor (ensuring that the integral of $f(\nu)$ gives unity) is $A = 0.32$. Corresponding theoretical advances (see, e.g., Sheth et al. 2001; Maggiore & Riotto 2010; Corasaniti & Achitouv 2011a; Achitouv & Corasaniti 2012) have revealed the physical meaning of the coefficients $a$ and $q$ in the framework of the excursion set theory (the normalization $A$ is not an independent parameter because its value is obtained from the requirement that the integral of $f(\nu)$ gives unity). This enables one to map the computation of the mass function in terms of a first-passage process of a random walk across a barrier. While a constant barrier corresponds to spherical collapse, the above authors showed that a drifting, diffusive barrier as a function of $\sigma$ accounts not only for the nonspherical form of the collapse, but also for more complex aspects of the underlying dynamics. The drift and diffusion coefficient of the barrier determine the coefficients $q$ and $a$ in Equation (2), which are then related to the collapse process of halos.

The universality of the above mass function and of the coefficients $a$ and $q$ for different cosmologies has been studied by various authors, starting with Sheth et al. (2001). These authors have tested the above expression against N-body simulations for a variety of CDM cosmologies including the cases of a critical universe ($\Omega_M = 1$, $\Omega_\Lambda = 0$, $h = 0.5$), an open universe ($\Omega_M = 0$, $3$, $\Omega_\Lambda = 0$, $h = 0.7$), and the $\Lambda$CDM case ($\Omega_M = 0.3$, $\Omega_\Lambda = 0$, $7$, $h = 0.7$). More recently, Despali et al. (2016) have tested the above mass function against the SBARBINE set of N-body simulations for a variety of combinations of $\Omega_M$ (ranging from 0.2 to 0.4) and $\Omega_\Lambda$ (ranging from 0.6 to 0.8), finding an excellent agreement with simulations when the same definition of halos (based on the virial value) is adopted for the same, fixed set of parameters. These authors concluded that—with the proper definition of a halo—the Sheth & Tormen mass function is universal as a function of redshift and cosmology. In all explored cosmologies, the deviations from the mass function that we adopt are shown to be within 5%. For nonstandard cosmologies, Achitouv et al. (2014) explored a Ratra-Peebles quintessence model of DE in addition to the standard $\Lambda$CDM model. They found that the abundance of dark matter halos described by a drifting diffusive barrier matches the results of N-body simulations to within 5% for all explored cosmologies. Most importantly, the parameters that define the diffusion of the barrier (determining the coefficients $q$ and $a$) change by less than 5% when passing from $\Lambda$CDM to the quintessence cosmology for the large masses $M \ge 10^{11} M_\odot$ relevant to this

paper. Motivated by the above results, in the rest of the paper we retain the same form of the mass function for all DDE cosmologies. We have verified that allowing for a 5% uncertainty in the mass function in Equation (2) does not change our results appreciably.

As for the threshold $\delta_c$, we notice that in principle this depends weakly on cosmology. Here we shall adopt the conservative value $\delta_c = 1.65$ for all DDE models. This constitutes a lower bound for the possible values taken in different DDE cosmologies (Mainini et al. 2003; Pace et al. 2010), thus maximizing the predicted abundance of massive DM halos.

The mass function $dN/dM$ (Equation (1)) allows us to compute the expected number of galaxies in a given range of mass and redshift over a given fraction of the sky, $f_{\text{sky}}$, as

$$N = f_{\text{sky}} \int dz \frac{dV}{dz} \int dM \frac{dN}{dM} \quad (3)$$

where the $z$ and $M$ integrals are over the region of the $(M, z)$ plane being considered.

The above expression depends on the assumed power spectrum of perturbations (determining the dependence of $\sigma$ on the mass $M$) and on cosmology, which affects the volume element $dV/dz$ and the growth factor of perturbations $D(z)$.

For the first, we adopt the CDM form (Bardeen et al. 1986), which has long been known to provide an excellent match to a wide set of observational data (see Tegmark & Zaldarriaga 2002, 2009; Hlozek et al. 2012). In this paper, we adopt a top-hat filter in real space to relate the variance $\sigma(M)$ to the power spectrum $P(k)$. As explained in detail in other papers (see, e.g., Benson 2012; Schneider et al. 2013) for CDM power spectra, which are essentially pseudo power laws with a slowly varying slope, the shape of $\sigma(M)$ varies little with the choice of filter function. In any case, even in cosmologies with a small scale cut-off (such as those of Warm Dark Matter models), the variance $\sigma(M)$ depends on the filter choice only at small masses $M \lesssim 10^9 M_\odot$, while in this paper we focus large masses $M \ge 10^{10} M_\odot$, for which the variance is in practice independent of the filter choice (see Benson et al. 2013).

The dependence on cosmology (and in particular on the DE equation of state) is our main focus here. In the present paper, we follow the approach adopted in Lamastra et al. (2012), to which we refer the reader for further details. Here, we summarize the key points.

We assume a spatially flat, homogeneous, and isotropic universe filled by nonrelativistic matter plus a dark energy component. For our analysis, we use the Chevallier-Polarski-Linder parameterization (Chevallier & Polarski 2001; Linder 2003) to describe the evolution in terms of the scale factor $a$ (normalized to unity at the present cosmic time):

$$w(a) = w_0 + w_a(1 - a) = w_0 + w_a \frac{z}{1 + z} \quad (4)$$

where the parameter $w_0$ represents the value of $w$ at the present epoch, while $w_a$ corresponds to its look-back time variation $w_a = -dw/d\alpha$. In the above parameterization, the standard $\Lambda$CDM cosmology corresponds to $w_0 = -1$ and $w_a = 0$. Using this parameterization, the cosmic expansion described by $H = \dot{a}/a$ is given by

$$E(z) \equiv H/H_0 = [\Omega_M a^{-3} + \Omega_\Lambda a^{-3(1+w_0+w_a)}e^{3w_a(a-1)}]^{1/2}. \quad (5)$$
The above equation also yields the line-of-sight comoving distance corresponding to a distant object at redshift z in any DE model:

$$ \chi(z) = \frac{c}{H_0} \int_{0}^{z} \frac{dz'}{E(z')} \qquad (6) $$

This enters into the expression for the luminosity and angular distances, and for the volume element $dV/dz$ (see, e.g., Weinberg 1972). In the following, we will indicate as $V_{w_0,w_a}$ the cosmic volume computed for DDE cosmologies, while $V_\Lambda$ is the same quantity computed for the case $w_0 = -1$ and $w_a = 0$ (the cosmological constant).

As for the growth factor, its expression in the $\Lambda$CDM case is given by Carrol et al. (1992) in the form:

$$ \delta(a) = \frac{5\Omega_M}{2a} \frac{da}{d\tau} \int_{0}^{a} \left( \frac{da'}{d\tau} \right)^{-3} da' \qquad (7) $$

where $\tau = H_0 t$. For the DDE models, we use the parameterization to the solution given in Linder (2005):

$$ \frac{\delta(a)}{a} = \exp\left(\int_{0}^{a} [\Omega(a)^{\gamma} - 1] d\ln a\right) \qquad (8) $$

where $\Omega(a) = \Omega_M a^{-3}/(H(a)/H_0)^2$, and $\gamma$ is the growth index given by the fitting formula (Linder 2005):

$$ \begin{aligned} \gamma &= 0.55 + 0.05(1 + w(z = 1))w(z = 1) \ge -1 \\ \gamma &= 0.55 + 0.02(1 + w(z = 1))w(z = 1) < -1. \end{aligned} \qquad (9) $$

This parameterization reproduces the behavior of the growth factor to within 0.1%–0.5% accuracy for a wide variety of dark energy cosmologies (Linder 2005; Linder & Cahn 2007) and allows for a rapid scanning of the parameter space of DDE models. We normalize the amplitude of perturbations $\delta$ in terms of $\sigma_8$, the local ($z = 0$) variance of the density field smoothed over regions of $8 h^{-1}$ Mpc. Present cosmological constraints based on Planck data, baryonic acoustic oscillations, and SNe Ia yield $\sigma_8 = 0.8$ for the $\Lambda$CDM cosmology (Planck Collaboration 2016); such a value varies by $\approx 2\%$ when different combinations ($w_0, w_a$) are assumed (Planck Collaboration 2016; Di Valentino et al. 2017; Mehrabi 2018).

Both the distance relation (Equation (6)) and the growth factors (Equation (8)) of DDE models deviate mildly from the cosmological constant case when the equation of state is negatively evolving with redshift $w_a < 0$, while models with $w_a > 0$ yield growth factors and cosmic times lower than those predicted in the $\Lambda$CDM case. This is because in the $w_a > 0$ models, the influence of DE at early times is strong even at high redshift (see Equation (5)), yielding shorter ages and implying a delay in the growth of DM perturbations compared to the $\Lambda$CDM case (see Figure 1 in Lamastra et al. 2012).

The impact of the above effects on the predicted abundance of halos (Equation (1)) is illustrated in Figure 1. We show the DM mass function derived from Equation (1) in selected DDE cases with four different values of $w_a$ (0, 0.5, 0.9, 1.1) and fixed $w_0 = -1$. For a prompt comparison, we have also shown the mass distribution of galaxies (with $2\sigma$ error bars) corresponding to the stellar mass function of Grazian et al. (2015), assuming that all baryons are converted into stars, i.e., $M_*/M = f_b$, to convert the observed stellar mass $M_*$ to the halo DM mass $M$ (see Section 3.1 for a detailed comparison).

<figure><img src="image_3.png" /><figcaption>Figure 1. The halo mass function at $z = 6$ predicted by different DDE models with $w_0 = -1$ and $w_a = 0, 0.5, 0.9$, and 1.1 (from left to right: light green, green, dark green, and red, respectively). The black squares correspond to the stellar mass function measured by Grazian et al. (2015) at $z = 6$; for the sake of simplicity, in this plot a conversion factor $M_*/M = f_b$ has been assumed to assign a DM mass to the stellar mass of each measured data point. We have shaded the region below each curve to stress that the DM halo mass function is an upper limit to the stellar mass function. In any assumed DDE scenario, the stellar mass function must entirely lie within the corresponding shaded region.</figcaption></figure>

While the mass function of observed galaxies is expected to be lower than the halo mass function due to the inefficient conversion of baryons into observable stellar mass, luminous galaxies cannot outnumber their host DM halos. Thus, the DM mass functions in the figure should be considered as upper limits for the mass distribution of observed galaxies.

The rapid, exponential decline of the DM mass function at large masses (Equations (1) and (2)) results in a large sensitivity of the predicted abundance of massive DM halos to the growth factor $D(z)$ of the different considered DDE models. As a result, in DDE models with $w_0 = -1$ and $w_a \ge 1.1$, the maximal abundance of massive DM halos is too low (deviation larger than $2\sigma$) to account for the observed number density of massive galaxies.

The example above shows that selecting extremely massive objects at high redshift is essential to provide constraints on DDE models. Thus, in the following sections, we will compare DDE predictions with different observations concerning the most massive objects already in place at high redshifts.

## 3. Results

Here, we compare the abundance of DM halos predicted by DDE models with different observations. When computing the mass function (Equation (1)), we assume a matter density parameter $\Omega_M = 0.31$ and a baryon density parameter $\Omega_b = 0.045$, corresponding to the values that provide the best fit to CMB data when $w_0$ and $w_a$ are allowed to vary (see Di Valentino et al. 2017); similar results are obtained if we convolve our predictions with a Gaussian uncertainty distribution centered on the above values and with variances $\sigma_{\Omega_M} = 0.02$ and $\sigma_{\Omega_b} = 0.005$, respectively. For the Hubble constant we take the value $h = 0.7$, although the best-fit values to the CMB (in combination with other probes) vary in the range $0.67 \le h \le 0.74$ when $w_0$ and $w_a$ are allowed to vary
(Di Valentino et al. 2017). In fact, the final constraints we obtain on the $w_0-w_a$ plane are weakly dependent on $H_0$ for the massive systems redshift range considered $z \approx 4-7$: for any given $w_0$, varying $H_0$ in the above range changes our constraints on $w_a$ by less than 2%.

### 3.1. The Stellar Mass Function at $z=6$ from CANDELS

We first compare with the observed stellar mass distribution of massive, distant galaxies. Since stellar mass is a time-integrated quantity, it is less sensitive to the details of the star formation history and can be more easily related to the DM mass of the host halo. However, an extended wavelength coverage is essential for estimating stellar masses from spectral energy distribution (SED) fitting, while measuring the abundance of massive, rare galaxies requires a combination of survey volume and depth. The CANDELS project (Grogin et al. 2011; Koekemoer et al. 2011) takes advantage of the optical/near-infrared/mid-infrared imaging provided by the Hubble Space Telescope (HST), Spitzer, and the Very Large Telescope, and provides an ideal data set to base such a measurement. Here, we use the high-redshift ($z = 5.5-6.5$) mass function derived by Grazian et al. (2015), who used a spectral-fitting technique to derive stellar masses for a galaxy sample with high-quality photometric redshifts based on the CANDELS-UDS, GOODS-South, and HUDF fields. The high redshifts we are considering ensure that at the largest masses $M_* \approx 10^{11} M_\odot$ probed by observations, the mass functions predicted by the different DDE models are in the full exponential regime and are steep enough to compare with the observed number density discriminant for the different DDE models (see Section 2 and Figure 1) with a high ($2\sigma$) confidence level, as we show below. To take into proper account the uncertainties related to the stellar mass measurements (that depend on age, dust extinction, metallicity, and star formation history), and to photometric redshifts, star formation histories, cosmic variance, and Poissonian statistical fluctuations, we have run a Monte Carlo simulation specific to the adopted data set. This allows us to derive, for any chosen stellar mass bin, the whole probability distribution functions $p(\phi_{\text{obs}})$; (PDF hereafter) of measuring an abundance $\phi_{\text{obs}}$.

We associate the stellar mass $M_*$ to the host halo DM mass $M$ using the relation $M_* = Ff_b M$, where $F$ describes the efficiency of baryon conversion into stars. While $F=1$ corresponds to the complete conversion, the different processes (gradual gas cooling, gas ejection, and stellar feedback) that take place in galaxies limit $F$ to lower values. In fact, the standard conversion for $\Lambda$CDM derived from abundance matching techniques (see, e.g., Behroozi & Silk 2018) yields $F \approx 0.25$. Such a value cannot be safely considered as a baseline for generic DDE models, since the stellar masses derived from abundance matching assume a $\Lambda$CDM halo mass function. However, we can study the effectiveness of baryon conversion into stars using hydrodynamical $N$-body simulations, since the physics of such a conversion is expected to weakly depend on the background cosmology. To this aim, we analyzed the public release of three simulations: the suite of Illustris simulations (Springel 2010; Genel et al. 2014; Vogelsberger et al. 2014a, 2014b), its updated version (IllustrisTNG, Weinberger et al. 2017; Pillepich et al. 2018), and the EAGLE simulation (Schaye et al. 2015). For Illustris TNG, we considered the highest resolution version of the largest and medium volume realizations (TNG300 and

TNG100). We computed the conversion efficiency $F = (M_*/M)/f_b$ from the ratio between the DM mass $M$ of each subhalo in the simulations and the stellar content $M_*$ associated with the considered subhalo, finding that $F = 0.5$ constitutes an effective, conservative upper limit for such a quantity, since (in all simulations) no massive ($M_* \ge 3 \times 10^{10} M_\odot$) galaxies have been found with $F \ge 0.5$ at $z = 5.5-6.5$.

Then, we consider a grid of DDE models characterized by different combinations ($w_0, w_a$). For each combination ($w_0, w_a$), we first correct the observed abundances $\phi_{\text{obs}}$ with the volume factor $f_{\text{vol}} = V_\Lambda/V_{w_0,w_a}$ (computed in the redshift range $z = 5.5-6.5$) to account for the fact that the mass function given in Grazian et al. (2015) has been derived assuming a $\Lambda$CDM cosmology. Analogously, we must take into account the fact that the stellar masses measured by Grazian et al. (2015) have been inferred from luminosities that were converted from observed fluxes assuming a $\Lambda$CDM cosmology. Thus, for each combination ($w_0, w_a$), we must correct the masses $M_*$ measured by Grazian et al. (2015) by a factor $f_{\text{lum}} = D_{L,w_0,w_1}^2/D_{L,\Lambda}^2$, where $D_{L,w_0,w_1}^2$ is the luminosity distance computed (at the considered redshift $z=6$) for a considered ($w_0, w_a$) combination, and $D_{L,\Lambda}^2$ is its value in the $\Lambda$CDM case.

We focus on the largest stellar mass bin (centered on $M_* = 8 \times 10^{10} M_\odot$, assuming a Salpeter initial mass function) analyzed by Grazian et al. (2015). For each combination ($w_0, w_a$), we compare the volume-corrected, observed abundance of galaxies $\tilde{\phi} = \phi_{\text{obs}} f_{\text{vol}}$ with stellar mass $M_* = 8 f_{\text{lum}} 10^{10} M_\odot$ at $z=6$ with the predicted number density $\phi_{w_0,w_a}(M)$; (Equation (1)) of DM halos with DM masses larger than $M = M_*/(F f_b)$ for the considered ($w_0, w_a$) combination. The confidence for the exclusion $P_{\text{excl}}$ of each considered DDE model is obtained from the PDF as the probability that the measured abundance is larger than number density predicted by the model, i.e., $P_{\text{excl}}(w_0, w_a) = \int_{\phi_{w_0,w_a}}^{\infty} p(\tilde{\phi}) d\tilde{\phi}$.

We show in Figure 2 the region of $w_0-w_a$ excluded at the $2\sigma$ confidence level (i.e., $P_{\text{excl}} \ge 0.95$) for $F=1$, $F=0.5$ and $F=0.25$. The exclusion region is overplotted on the regions allowed by CMB and weak lensing observations, and to the region derived by the combination of the same data with the Hubble diagram of supernovae and quasars (Risaliti & Lusso 2019). Our probe significantly restricts the region in DDE parameter space allowed by other methods. In particular, we exclude an appreciable part of the region favored by the distant quasar method.

We stress that our method allows for significant improvements when more extended databases will be available in the future. To stress this point, we show in Figure 3 the constraints that would be obtained from the Grazian et al. (2015) stellar mass function if we decrease by half the dispersion in the PDF around the average value. This (approximately) simulates the effects of the larger statistics that would be obtained analyzing the full CANDELS data set (five fields). In fact, the spread $\Delta\phi_{\text{obs}}$ in the measured stellar mass function at large masses $M_* \approx 10^{11} M_\odot$ is dominated by Poisson fluctuations that account for 70% of $\log \Delta\phi_{\text{obs}}$ (while cosmic variance and the uncertainties related to the assumed star formation histories approximately account for 15% and 8%, respectively). In this case, most of the region allowed by distant quasars would be excluded, even in the conservative case $F=0.5$.

Up to this point, we have focused our analysis on the stellar mass function by Grazian et al. (2015). This is because
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Exclusion regions (2σ confidence level) in the $w_0-w_a$ plane derived from the observed CANDELS stellar mass function at $z=6$ (Grazian et al. 2015). The brown, red, and orange regions correspond to the assumption of $F=1$, $F=0.5$, and $F=0.25$, respectively (see text). Our exclusion region is compared with the 2σ and 3σ contours allowed by CMB+weak lensing (green regions) and by the combination of the same data with the Hubble diagram of supernovae and quasars (blue region), derived from Figure 4 of Risaliti & Lusso (2019). The black dot corresponds to the ΛCDM case ($w_0 = -1, w_a = 0$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Figure 2, but assuming the error bars of the stellar mass function are reduced by half, to simulate the inclusion of the full set of CANDELS fields. This improvement would exclude most of the region allowed by distant quasars.</figcaption></figure>

it couples large stellar mass coverage at high redshifts (stellar masses extending to $M_* \approx 10^{11} M_\odot$) with a detailed analysis of the uncertainties, including those related to the assumption of different star formation histories, metallicities, ages, dust extinction, photometric redshifts, and cosmic variance. Indeed, for this specific measurement, we could exploit our previous work to perform a full computation of the whole PDF through Monte Carlo simulations that account for the above uncertainties.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Exclusion regions (2σ confidence level) in the $w_0-w_a$ plane derived from the observed CANDELS GOODS-South field stellar mass function at $z=6$ by Duncan et al. (2014).</figcaption></figure>

Of course, our method can be applied to other estimates of the stellar mass function at high redshifts. Although several observational works exist in the literature (see, e.g., Duncan et al. 2014; Song et al. 2016; Stefanon et al. 2017; Bhatawdekar et al. 2019), to our aim it is essential to compare with observed stellar mass functions that cover high redshifts ($z \gtrsim 6$) and large masses $M_* \gtrsim 5 \times 10^{10}$, where the upper limit provided by the halo mass function can be violated in some DDE models (see Figure 1). For example, the extremely deep mass functions measured by Bhatawdekar et al. (2019) for the Frontier Field galaxies do not provide any constraint to DDE models, since at high redshifts ($z \ge 6$) they reach masses $M_* \approx 5 \times 10^9 M_\odot$; these probe DM halo masses below $M \le 10^{11} M_\odot$, where the upper limit provided by the halo mass functions is essentially consistent with all DDE models (see Figure 1). The same argument applies to the measurements by Song et al. (2016), who analyze the CANDELS GOODS-South field to probe the distribution of stellar masses up to $M_* \lesssim 1.5 \times 10^{10} M_\odot$ at $z=6$ based on the stellar mass-to-light conversion. Stefanon et al. (2017) also use a stellar mass-to-light conversion to probe somewhat larger stellar masses extending up to $M_* \lesssim 3 \times 10^{10} M_\odot$ at $z \approx 6$, but the error bar of the most massive bin is too large to allow for effective discrimination among different DDE models.

A better constraining measurement for our scope has been performed by Duncan et al. (2014), who analyzed the CANDELS GOODS-South field. Assuming an observationally based stellar mass-to-light conversion, the above authors measured the stellar mass function at $z \approx 6$ up to large stellar masses $M_* \approx 10^{11} M_\odot$. Since we cannot compute the full PDF for this observational analysis, the 2σ confidence level of exclusion has been derived doubling the error bars presented in Duncan et al. (2014). To account for the uncertainty related to cosmic variance (not provided by the above authors for the large masses considered here) we have added in quadrature the value provided by Song et al. (2016; for the same GOODS-South field) for different stellar masses at various redshifts. The result is presented in Figure 4, and shows constraints on DDE
<figure><img src="image_7.png" /><figcaption>**Figure 5.** For each combination $w_0-w_a$, we show the value of the baryon efficiency $F$ (color coded as shown in the vertical bar) needed to match the stellar mass function at $z=4$ for $10.5 \le M_*/M_\odot \le 11$.</figcaption></figure>

models that are very similar to the ones derived using the measurement by Grazian et al. (2015).

Finally, we stress that the constraints presented above are likely to be extremely conservative. Although baryon conversion efficiencies larger than $F=0.5$ are in principle possible in different DDE models, values closer to the $\Lambda$CDM value $F \approx 0.2$ (derived from abundance matching; see, e.g., Behroozi et al. 2013; Moster et al. 2018) are much more probable. To address this point in a closer detail, we have computed the values of $F$ that would be needed to match (i.e., to lay within the $1\sigma$ error bar) the massive end ($10.75 \le M_*/M_\odot \le 11$) of the observed stellar mass function at $z=4$. Adopting the stellar mass function at $z=4$ by Grazian et al. (2015), we derive—for each DDE model and for the above range of $M_*$—the baryon efficiency $F$ shown in Figure 5. For all relevant $w_0-w_a$ combinations, we obtain $F < 0.5$, while for small values of $w_0 \lesssim -1.5$, we obtain $F \lesssim 0.25$. (Similar values of $F$ are obtained when comparing with the data by Ilbert et al. 2013; Duncan et al. 2014; Stefanon et al. 2015; Davidzon et al. 2017. Even lower values of $F$ are obtained comparing with the data by Song et al. 2016.) Notice that when the $\Lambda$CDM model is considered (the black point in Figure 5) our analysis yields $F = 0.18$, corresponding to a stellar-to-DM mass ratio $M_*/M = Ff_b = 0.025$, in excellent agreement with the value estimated by Behroozi et al. (2013; see also Figure 9 in Song et al. 2016) from abundance matching in the $\Lambda$CDM case.

Under the assumption of constant $F$ between redshift $z=4$ and $z=6$, we can then use the values of $F$ shown in Figure 5 to recompute the constraints on DDE models from the comparison with the CANDELS stellar mass function at $z=6$. The result, shown in Figure 6, shows that in this case we obtain even stronger constraints, since the values of $F$ for the different DDE models are now within the range $0.2 \lesssim F \lesssim 0.25$ for a wider set of combinations $w_0-w_a$. This shows that—if the baryon fraction $F$ has a slow evolution in DDE models between $z=4$ and $z=6$—the assumptions used to derive the constraints in Figures 2–4 are indeed conservative.

## 3.2. Massive Galaxies Detected in the Submillimeter at $z \approx 4$

The above population of galaxies (identified in rest-frame optical and ultraviolet) is known to under represent the most

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Figure 2, but assuming the baryon efficiency $F$ from Figure 5.</figcaption></figure>

massive galaxies, which have rich dust content and/or old stellar populations. These are, however, detectable at submillimeter wavelengths. Recently, Wang et al. (2019) performed detailed submillimeter (870 µm) observations at the Atacama Large Millimeter/submillimeter Array (ALMA) of a sample of Spitzer/Infrared Array Camera (IRAC)—bright galaxies. They detected 39 star-forming objects at $z > 3$, which are unseen in even the deepest near-infrared (H-band) imaging with the HST (“H-dropouts”), which proved to be massive galaxies with stellar mass extending up to $M_* \approx 3 \times 10^{11} M_\odot$, with a median mass $M_* \approx 4 \times 10^{10} M_\odot$.

For such objects, we follow a procedure similar to that explained in the previous section. We compute the number density of galaxies with stellar masses in the bin $10.25 \le \log(M_*/M_\odot) \le 10.75$ (dominating the statistics of observed objects) at redshifts $z = 4.5–5.5$, and derive the corresponding $2\sigma$ lower limit $\phi_{\text{low}}(M_*) = 1.8 \times 10^{-5} \text{ Mpc}^{-3}$. To relate the observed stellar mass $M_*$ to the DM mass $M$, we adopt the highly conservative assumption $M = M_*/f_b$. We then compute the number density $\phi_{w_0,w_a}$ of a DM halo mass corresponding to the observed $M_*$ for different combinations ($w_0, w_a$), and compared it with the observed $2\sigma$ lower limit $\phi_{\text{low}}$. For each combination ($w_0, w_a$), observed number densities and stellar masses (measured assuming a $\Lambda$CDM cosmology) have been rescaled with the factors $f_{\text{vol}}$ and $f_{\text{lum}}$ (see Section 2). The comparison allows us to exclude (at a $2\sigma$ confidence level) the combinations ($w_0, w_a$) for which $\phi_{w_0,w_a} < \phi_{\text{low}}$. The result is shown as a brown exclusion region in Figure 7.

Of course, the above approach is very conservative, since we assume that the whole baryonic mass is in stars, and that the baryon mass of DM haloes is related to the DM mass through the universal baryon fraction (no loss of baryons). In fact, the very fact that the objects are characterized by a high star formation rate $\gtrsim 200 M_* \text{ yr}^{-1}$ indicates that a sizable fraction of baryons is in the form of gas. Properly accounting for such a gas fraction would yield larger values of $M$ associated with the observed $M_*$ and, hence, tighter constraints. Although we have attempted to estimate the gas mass for the ALMA-detected H-dropout galaxies from the submillimeter continuum using
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Exclusion regions (2σ confidence level) in the $w_0-w_a$ plane (see text) derived from the observed abundance $\phi_{\text{obs}}$ of luminous submillimeter galaxies at $z = 4.5-5.5$ (Wang et al. 2019). The brown region corresponds to the assumption that the observed stellar masses $M_* = M/f_b$ are related to the DM mass through the baryon fraction $f_b$. The red region corresponds to adoption of the DM mass derived from the measured cross-correlation function of H-dropouts (see text).</figcaption></figure>

the dust mass and the dust-to-gas ratio, the inferred gas masses are affected by large uncertainties (they span a range between $5 \times 10^9 M_{\odot}$ and $5 \times 10^{10} M_{\odot}$) related to photometric redshifts (uncertainties are particularly critical given the steep shape of the spectrum in the far-infrared), the adoption of a single and simplified graybody at the average temperature, and the adoption of the mass-metallicity relation at $z = 3.5$ for all sources.

To bypass the uncertainty related to the gas fraction, and to derive more realistic constraints for DDE models, we analyze the clustering properties of the H-dropouts. We base our analysis on the procedure adopted by Wang et al. (2019), who estimated the two-point angular cross-correlation function $\omega(\theta)$ of H-dropouts with all CANDELS galaxies in the redshift range $3.5 \le z \le 5.5$. Assuming a power-law form for the cross correlation $\omega(\theta) = A_{\omega} \theta^{-\beta} - \text{IC}$ (with $\beta = 0.8$ and IC the integral constraint in Equation (4) of Wang et al. 2019), the above authors derived the amplitude $A_{\omega}$. This is related to the correlation length $r_0$ by the Limber equation (Croom & Shanks 1999; He et al. 2018):

$$r_0 = \left[ A_{\omega} \frac{c}{H_0 Q} \frac{\int N_{\text{H}}(z) dz \int N_{\text{G}}(z) dz}{\int N_{\text{H}}(z) N_{\text{G}}(z) \chi^{1-\gamma}(z) E(z) dz} \right]^{1/\gamma} \quad (10)$$

where $\gamma = \beta + 1$, the constant $Q = \Gamma(1/2)\Gamma(\gamma/2 - 1/2)/\Gamma(\gamma/s)$ is a combination of $\Gamma$ functions, $\chi(z)$ and $E(z)$ are given in Equations (5) and (6), and $N_{\text{H}}(z)$ and $N_{\text{G}}(z)$ are the redshift distributions of H-dropouts and CANDELS galaxies. The correlation length $r_0$ is then converted to galaxy bias through the relation (Peebles 1993)

$$b = \frac{72}{(3 - \gamma)(4 - \gamma)(6 - \gamma)2^{\gamma}\sigma_{8}(z)} \left[ \frac{r_{0}}{8 h^{-1} \text{Mpc}} \right]^{\gamma}, \quad (11)$$

that we assume to hold in all DDE models. Here, $\sigma_8(z)$ is the amplitude of the dark matter fluctuation on the scale of $8 h^{-1}$ Mpc. The DM mass is then derived from the relation $b = 1 + [\nu(M, z) - 1]/\delta_c$ (Mo & White 2002). For the standard $\Lambda$CDM case, the above procedure yields $M = 10^{13 \pm 0.3} M_{\odot}$ for the average DM mass (Wang et al. 2019).

For our comparison with DDE predictions, we cannot take the above DM mass at face value, since it has been derived assuming a $\Lambda$CDM cosmology. In fact, for generic DDE cosmologies, the above value will (weakly) change due to two factors: (i) the Limber equation (Equation (10)) relating the observed $A_{\omega}$ to $r_0$ depends on cosmology through the functions $E(z)$ and $\chi(z)$; (Section 2); and (ii) the different growth factor (Section 2, Equation (7) and below) affects the quantities $\nu(M, z)$ and $\sigma_8(z)$ entering the computation of the average mass $M$ (Equation (11) and below). Thus, we compute the maximal effect of cosmology on the value of $M$ derived by Wang et al. (2019) when our grid of values for the combinations ($w_0, w_a$) is considered. Assuming the same measured angular cross-correlation amplitude $A_{\omega}$, we consider the effect of different cosmologies on the derived 3D correlation length (Equation (10)) and on the bias factor (Equation (11)). We find that $M = 10^{13} M_{\odot}$ constitutes a $2\sigma$ lower limit for the value of the DM mass derived from cross correlation for any DDE model we consider. We then conservatively compute the DDE number density of objects with such a DM mass and compare it with the observed number density of H-dropouts with stellar mass $M_* = 10^{10.5} M_{\odot}$ (the average stellar mass of the sample). The resulting exclusion region in the ($w_0, w_a$) plane is shown in red in Figure 7.

### 3.3. SPT0311-58 at $z = 6.9$

The most massive system detected at $z \ge 6$ is a far-infrared-luminous object at redshift $z = 6.9$ originally identified in the 2500 deg$^2$ SPT survey (Marrone et al. 2018). Observation in the optical with the HST, infrared observations with the Spitzer Space Telescope, Gemini Optical/IR imaging, and spectroscopy subsequently allowed for a characterization of this source. High-resolution imaging revealed this source (denoted SPT0311-58) to be a pair of extremely massive star-forming galaxies, with the larger galaxy (SPT0311-58W) forming stars at a rate of 2900 $M_{\odot}$ yr$^{-1}$. An elongated faint object seen at optical and near-infrared wavelengths is consistent with a nearly edge-on spiral galaxy at $z \approx 1.4$ acting as a gravitational lens for the source, with an estimated magnification $\mu = 2$.

Measurements of the far-infrared continuum with ALMA led to an estimate of a huge H$_2$ gas mass ranging from $M_{H_2} \approx (7.6 \pm 2) \times 10^{10} M_{\odot}$ (based on the CO luminosity converted with a standard value $\alpha_{\text{CO}} = 1 \text{ km s}^{-1} \text{ pc}^2$) to $M_{H_2} \approx (3.1 \pm 1.9) \times 10^{11} M_{\odot}$ (estimated from a radiative transfer model in Strandet et al. 2017).

In the following, we shall adopt the latter value as a baseline, since it is based on a detailed fit with a radiative transfer model, built ad hoc to study the properties of interstellar medium of this object. As explained in detail in Strandet et al. (2017 and references therein), the model provides a fit to the far-infrared continuum sampled in seven broad bands, three CO lines and the [C I] line observed in this object by ALMA, Herschel, SPT, APEX, and ATCA. As mentioned in Marrone et al. (2018), such an estimate is more accurate than the one obtained by converting the CO line. In fact, the latter is derived using a
standard average ULIRG-like $\alpha_{CO} = 1$ km s$^{-1}$ pc$^2$. However, such a quantity strongly depends in specific properties of the considered object like star formation and metallicity (see Bolatto et al. 2013), so that a large scatter around the average value in different objects is expected (as, in fact, is observed; see, e.g., Weiß et al. 2007; Papadopoulos et al. 2012). Indeed, the radiative transfer model applied to SPT0311-58 yields $\alpha_{CO} \approx 4$ km s$^{-1}$ pc$^2$. The authors explain that the difference with respect to the typical ULIRG-like factor is due to the much higher density of the interstellar medium in this object.

To estimate the DM mass associated with this object, we cannot follow the procedure adopted in Marrone et al. (2018), since they derive the gas-to-DM conversion factor from abundance matching techniques (see, e.g., Behroozi & Silk 2018) that cannot be safely considered as a baseline for generic DDE models because they based on the $\Lambda$CDM halo mass function.

Therefore, to estimate a conversion fraction from the observed H$_2$ mass to the DM mass, we first adopted the conservative assumption that the total baryonic mass $M_b = M_\ast + M_{gas}$ (here, $M_{gas}$ is the total gas mass) is related to the DM mass through the baryon fraction $M = (M_\ast + M_{gas})/f_b$. Although no stellar light is convincingly seen from SPT0311-58W (probably due to the large extinction), a lower limit on the stellar content can be inferred from existing measurements of the molecular gas fraction $f_{H_2} = M_{H_2}/(M_\ast + M_{H_2})$. Measurements of high-z star-forming galaxies (ranging from relatively quiescent BzK galaxies to dusty starbursts), suggest $f_{H_2} = 0.2-0.8$ (e.g., Daddi et al. 2010; Tacconi et al. 2010, 2013; Geach et al. 2011; Magdis et al. 2012; Combes et al. 2013; see Casey et al. 2014 for a review). However, all theoretical models (Benson 2012; Daveé et al. 2012; Fu et al. 2012; Lagos et al. 2012; Popping et al. 2014; see also Gabor & Bournaud 2013; Ginolfi et al. 2019) predict typically smaller values in the range $f_{H_2} \lesssim 0.5$. One possible solution to this mismatch has been offered by Narayanan et al. (2012), who suggested that the canonical conversion factor CO-H$_2$ was too large for the most extreme systems at high redshift, and that the correct observed gas fractions are in the range $f_{H_2} = 0.1-0.4$. Similar conclusions were drawn by Tacconi et al. (2013), who suggested that the tension between galaxy gas fractions measured in observations and those in simulated galaxies may be due to the incomplete sampling of galaxies.

Even assuming that H$_2$ constitutes 80% of the gas mass (i.e., $f_g \equiv M_{H_2}/M_{gas} = 0.8$) at high redshifts (an upper limit according to Lagos et al. 2011, 2014), the estimated baryonic mass $M_b = M_{gas} + M_\ast = M_{H_2}(f_{H_2} + f_g - f_{H_2}f_g)/f_{H_2}f_g$ takes the value $M_b = 1.4 M_{H_2}$ if we adopt the most conservative estimate $f_{H_2} = 0.8$, and $M_b = 2.75 M_{H_2}$ if we adopt the estimate $f_{H_2} = 0.4$ suggested by theoretical models and by the effects suggested by Narayanan et al. (2012) or by Tacconi et al. (2013). Although the latter works refer to galaxies at $z \le 4$, the consideration of the above range of uncertainty is the best we can do with the present data and theoretical predictions. This leads to the association of a DM mass to the observed $M_{H_2}$ $\bar{M} = M_b/f_b = 2 \times 10^{12} M_\odot$ in the most conservative case, and to $\bar{M} = M_b/f_b \approx 6 \times 10^{12} M_\odot$ in the other case; we will consider both values in the following analysis. An even larger DM mass would be consistent with the observations if the object lost the majority of its molecular gas content.

To estimate the rareness of such a system in all the considered DDE cosmologies, we compute the Poisson

probability of finding such a massive object within the volume probed by the SPT survey, for different combinations ($w_0, w_a$). Following the method in Harrison & Hotchkiss (2013), as done in Marrone et al. (2018) for the $\Lambda$CDM cosmology, we first compute from Equation (3) the number $N(M, z)$ of systems with mass $M$ and higher at redshift $z$ and higher expected in the sky area $f_{sky}$ covered by the SPT survey, for a grid of values of $M$ and $z$. Then, we compute such a number $N(\bar{M}, \bar{z})$ for the values $\bar{M}$ and $\bar{z}$ associated with the observed systems (i.e., $\bar{z} = 6.9$ and $\bar{M} = (2-6) \times 10^{12} M_\odot$ as discussed above). Finally, we consider the number $N_{rare}$ defined as $N(M, z)$ computed only for the masses $M$ and redshifts $z$ for which $N(M, z) \ge N(\bar{M}, \bar{z})$, as discussed in Harrison & Hotchkiss (2013). The Poisson probability of observing at least one system with both greater mass and redshift than the one which has been observed is

$$R_{>\bar{M},>\bar{z}} = 1 - \exp(-N_{rare}). \quad (12)$$

The above probability depends on the region of the $M-z$ plane to which the SPT survey is sensitive (which provides the lower limit for the integration over redshift and mass in Equation (3)), and on $f_{sky}$. Following Marrone et al. (2018), we assume that the survey is complete for $z \ge 1.5$ and for $M \ge 10^{11} M_\odot$, a conservative assumption, as discussed in detail by the above authors. To take into account the uncertainties in the measured value of $M_{H_2}$ (in turn affecting the corresponding DM mass $M$), we follow Harrison & Hotchkiss (2013) and convolve Equation (12) with the probability of measuring a given $M_{H_2}$, assuming a Gaussian distribution around the central value $M_{H_2} = 3.1 \times 10^{11} M_\odot$ with standard deviation $1.9 \times 10^{11} M_\odot$. As for the effective fraction of the sky $f_{sky} = \Omega_{sky}/(41253 \text{ deg}^2)$ in Equation (3), the total area corresponding to the SPT survey is $\Omega_{sky} = 2500 \text{ deg}^2$. However, Marrone et al. (2018) noticed that the effective survey area is potentially much smaller. In fact, most of the objects in the survey are strongly lensed, indicating that a source must be gravitationally lensed to exceed the 20 mJy threshold for inclusion in follow-up observations. Given the uncertainties related to a proper accounting for such an affect, we show our results for both the total area ($\Omega_{sky} = 2500 \text{ deg}^2$) and for an effective area reduced by $1/10$ ($\Omega_{sky} = 250 \text{ deg}^2$) to illustrate the effect of such an uncertainty (Marrone et al. considered an even more extreme case, $\Omega_{sky} = 25 \text{ deg}^2$).

For each combination ($w_0, w_a$), we compute the expected number of systems like SPT0311-58 detectable in the SPT survey. Then we associate a rareness to the resulting predicted number after Equation (10), and we compute the associated exclusion regions in the $w_0-w_a$ plane. The result ($2\sigma$ confidence level) is shown in Figure 8 for the case $\Omega_{sky} = 2500 \text{ deg}^2$, for the two considered values $\bar{M} = 2 \times 10^{11} M_\odot$ (red region) and $\bar{M} = 6 \times 10^{11} M_\odot$ (orange region). We also show as a dashed line the bound of the exclusion region that would be obtained for $\bar{M} = 0.5 \times 10^{11} M_\odot$. This would correspond to a gas mass $M_{H_2} = (7.6 \pm 1.9) \times 10^{10} M_\odot$ based on the CO luminosity, converted with a standard value $\alpha_{CO} = 1$ km s$^{-1}$ pc$^2$ with $f_{H_2} = 0.8$.

In the case $\bar{M} = 6 \times 10^{11} M_\odot$, corresponding to the assumption of $f_{H_2} = 0.4$ for the H$_2$ gas fraction, a major portion of the $w_0-w_a$ plane is excluded, although the $\Lambda$CDM case ($w_0 = -1, w_a = 0$) remains allowed. The excluded region
<figure><img src="image_10.png" /></figure>

**Figure 8.** Exclusion regions (2σ confidence level) in the $w_0-w_a$ plane (see text) for two different inferred DM masses of SPT0311-58: $2 \times 10^{12} M_\odot$ (red area) and $6 \times 10^{12} M_\odot$ (yellow area). In both cases, the full SPT survey area $\Omega_{\text{sky}} = 2500 \text{ deg}^2$ has been assumed. The dashed line represents the bound of the exclusion region that would be obtained for $\bar{M} = 0.5 \times 10^{11} M_\odot$, corresponding to the gas mass derived from the CO measurement with conversion factor $\alpha_{\text{CO}} = 1 \text{ km s}^{-1} \text{ pc}^2$ (see text).

includes both the larger $w_a$ cases allowed by the quasar method (blue region) and the cases $w_0 \ge -0.6$ allowed by the CMB+ weak lensing results, which shows the potential impact of our results. Even tighter constraints are obtained for the case $\Omega_{\text{sky}} = 250 \text{ deg}^2$, shown in Figure 9.

### 3.4. Combining the Different Probes

In the previous Sections 3.1–3.3, we have shown the potential of different observables as constraints on DDE models, and we discussed how the effectiveness of each probe relies on how much the observed baryon-to-DM mass ratio is suppressed with respect to the baryon fraction limit. While future observations will allow for a more precise determination of the gas and stellar mass fractions (see the discussion in Section 4 below), strong constraints can be derived—even under the most conservative assumptions—when all the probes presented in Sections 3.1–3.3 are combined. In fact, the probabilities for each combination ($w_0, w_a$) to be consistent with each of the considered observations are independent. Therefore, we can derive a combined constraint by multiplying the probabilities of being consistent with each probe. The resulting exclusion region is shown in Figure 10, adopting—for each probe in Sections 3.1–3.3—the most conservative assumption for the relation between the observed baryonic component and the DM mass $M$. For the comparison with the CANDELS field, we assume that the observed stellar mass is $M_* = 0.5 f_b M$ (i.e., $F = 0.5$, see Section 3.1). For the comparison with the abundance of submillimeter galaxies (Section 3.1), we assume that the observed stellar mass is related to $M$ by the baryonic fraction limit. As for the rareness of SPT0311-58, we take the conservative values for the gas mass fraction, leading to a DM mass estimate $M = 2 \times 10^{12} M_\odot$ (see Section 3.3), and we consider the whole survey area ($\Omega_{\text{sky}} = 2500 \text{ deg}^2$).

<figure><img src="image_11.png" /></figure>

**Figure 9.** Same as Figure 8, but assuming an affective SPT area $\Omega_{\text{sky}} = 250 \text{ deg}^2$.

<figure><img src="image_12.png" /></figure>

**Figure 10.** Exclusion regions (2σ confidence level) in the $w_0-w_a$ plane derived from combination of the different probes. For each observable, the most conservative case has been considered: for the CANDELS field, we have assumed $F = 0.5$, for SPT0311-58 we have taken a DM mass $M = 2 \times 10^{12} M_\odot$, and for submillimeter galaxies, we have converted stellar masses to DM mass assuming $M = M_*/f_b$. The dashed line shows the analytical approximation for the boundary of the excluded region $w_a = -3/4 - (w_0 + 3/2)$.

Inspection of Figure 10 shows that a major fraction of the parameter space favored by distant quasars, combined with the CMB and weak lensing, is excluded at the 2σ confidence level, independent of the details of the assumed baryon physics.

## 4. Conclusions and Discussion

We have computed the abundance of massive systems predicted using different DDE models at high redshifts $z \approx 4-7$. Such predictions have been compared with different observational probes: the bright end of the stellar mass function at $z \ge 6$, the space density of luminous submillimeter galaxies
at $z = 4-5$, and the rareness of the extreme hyperluminous infrared galaxy SPT0311-58 at $z \approx 7$.

We have derived exclusion regions in the parameter space $w_0-w_a$ of DDE models from each of the above probes. Adopting the most conservative assumptions for the ratio between the observed baryonic component and the DM mass, we have combined the above results to derive conservative, robust constraints for the parameter space of DDE models that do not depend on the details of the baryon physics involved in galaxy formation. In addition, our results do not depend on the nature of the DM component when the present limit on the mass of DM particle candidates $m_X \gtrsim 3$ keV (see, e.g., Viel et al. 2013; Menci et al. 2016) is taken into account. In fact, for DM particle masses in the keV range (Warm Dark Matter), the associated power spectrum (Bode et al. 2001; Destri et al. 2013) on the mass scales investigated in this work $M \gtrsim 10^{10} M_\odot$ is identical to the CDM form assumed here, and our results are unchanged.

### 4.1. Implications of Our Results

1. When the most conservative values concerning the baryon-to-DM mass are assumed, our combined results allow us to rule out DDE models with

$$w_a \geq -3/4 - (w_0 + 3/2)$$

as displayed in Figure 10, thus excluding a major fraction of the parameter space favored by the quasar distances (Risaliti & Lusso 2019), including the best-fit combination $w_0 \approx -0.8$ and $w_a = -1.5$ obtained with such a probe.

2. Our results leave open the possibility that the present tension in the value of $H_0$ between the values derived from Planck and those obtained from local luminosity distance measurements may be solved in DDE models, since the combinations ($w_0, w_a$) that allow for the reconciliation of the different observations include values outside our exclusion region (see Di Valentino et al. 2017; Zhao et al. 2017).

3. On the other hand, our results almost entirely rule out the quintessence models in which initially $w > -1$ and $w$ decreases as the scalar rolls down the potential (“cooling” models), which occupy most of the region $w_0 > -1$, $w_a > 0$ (see Barger et al. 2006). These typically arise in models of dynamical supersymmetry breaking (Binetruy 1999; Masiero et al. 2000) and supergravity (Brax & Martin 1999; Copeland et al. 2000), including the “freezing” models in Caldwell & Linder (2005) in which the potential has a minimum at $\phi = \infty$.

4. For “phantom” models with $w_0 < -1$ (see Caldwell 2002), our constraint $w_a \geq -3/4 - (w_0 + 3/2)$ excludes a major portion of the parameter space that corresponds to models for which the equation of state crossed the phantom divide line $w = -1$ from a higher value.

### 4.2. Improving Constraints with Improved Measurements

For each of the observables we considered, our constraints can be greatly tightened when improved, reliable measurements of the actual baryon fraction in galaxies, and of the relative weight of each baryonic component, are available. For example, a stellar-to-halo mass ratio $M_{\text{star}}/M = 0.25 f_b$ (a value

favored by present hydrodynamical N-body simulations) would greatly tighten the constraints from the stellar mass function, allowing us to rule out all models with $w_a \geq 1$ that are presently allowed by the distant quasar method.

Also, the constraints from the abundance of submillimeter galaxies at high redshifts could be greatly tightened when the gas mass of H-dropouts are reliably measured. Spectroscopic follow-up of H-dropout galaxies with future facilities (e.g., the James Webb Space Telescope, JWST) will add a valuable improvement to the present analysis. Future measurements on the H₂ gas fraction at high redshift will also allow for the reduction of the present gap with respect to the theoretical expectations $f_{\text{H}_2} \approx 0.1-0.4$.

Such improved measurements will probably need future facilities. For example, while a more accurate estimate of the gas-to-stellar mass fraction for the SPT0311-58 pair could in principle be inferred from their stellar mass, the latter is currently poorly constrained: their rest-frame optical SED is only sampled by two (IRAC CH1 and CH2 for the Eastern source) and four (F125W, F160W, IRAC CH1, and CH2 for the Western source) photometric points, resulting in a $1\sigma$ uncertainty on the inferred stellar mass spanning a factor of 15-20. In the near future, JWST will easily improve the accuracy of the stellar mass measurement of the SPT0311-58 pair by providing a much more detailed characterization of the rest-frame optical and near-infrared SED of these galaxies.

### 4.3. Statistics

Increasing the statistics of high-redshift massive objects will also greatly tighten present constraints (as shown by the comparison between Figures 2 and 3). Large surveys from space with the Euclid (Laureijs et al. 2011) and the Wide Field Infrared Survey Telescope (WFIRST, Spergel et al. 2015) satellites will increase the number of massive, high-z galaxies by orders of magnitude with respect to current HST samples. The Euclid surveys will cover 15,000 deg² at $H \leq 24$ mag depth, and 40 deg² at $H \leq 26$, while the WFIRST High Latitude Survey will observe 2200 deg² at $H \leq 26.7$.

As a reference, the CANDELS GOODS-South sample comprises only one source with $M_{\text{star}} \approx 10^{11} M_\odot$ at $z \geq 6$ for $H \leq 24$, and seven such sources at $H \leq 26.7$ in an area $\approx 0.05 \text{ deg}^2$. The statistical uncertainty on the stellar mass function will thus be reduced by a factor 30-300 by the aforementioned surveys, extending also to higher masses than those probed today. Unfortunately, systematic uncertainties will then dominate the error budget, mostly because the observed H band samples the rest-frame UV at $z \geq 6$, resulting in a potentially biased and incomplete selection of massive sources. In addition, the lack of information in the optical rest-frame adds significant uncertainties to the physical parameters estimated from SED fitting. This problem will be overcome by JWST observations with the Mid-Infrared Instrument at 5.6-25 $\mu$m, albeit on a much smaller area than the Euclid and WFIRST observations.

Despite the lack of any plan for large mid-infrared surveys from space, the combination of H-selected samples from future cosmological surveys and improved characterization of high-z objects on smaller areas thanks to JWST, will lead to tighter constraints on the high-mass end of the stellar mass function at
z ≥ 6, and thus on the parameter space (w₀, wₐ) of DDE models.

We acknowledge support from INAF under PRIN SKA/ CTA FORECaST and PRIN SKA-CTA-INAF ASTRI/CTA Data Challenge. N.G.S. acknowledges CNRS for Emeritus Director of Research contract in LERMA-Observatoire de Paris-PSL-Sorbonne U. We thank the referee for helpful and constructive comments that helped to improve the paper.

## ORCID iDs

N. Menci https://orcid.org/0000-0002-4096-2680

A. Grazian https://orcid.org/0000-0002-5688-0663

M. Castellano https://orcid.org/0000-0001-9875-8263

P. Santini https://orcid.org/0000-0002-9334-8705

E. Giallongo https://orcid.org/0000-0003-0734-1273

A. Lamstra https://orcid.org/0000-0003-2403-913X

A. Fontana https://orcid.org/0000-0003-3820-2823

E. Merlin https://orcid.org/0000-0001-6870-8900

## References

Achitouv, I., Wagner, C., Weller, J., et al. 2014, JCAP, 2014, 077

Achitouv, I. E., & Corasaniti, P. S. 2012, JCAP, 02, 002

Addison, G. E., Huang, Y., Watts, D. J., et al. 2016, ApJ, 818, 132

Ade, P. A. R., Agnahim, N., Arnaud, M., et al. 2016, A&A, 594, A14

Aghanim, N., Douspis, M., Hurier, G., et al. 2019, arXiv:1907.12875

Amara, A., & Refregier, A. 2007, MNRAS, 381, 1018

Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

Barger, V., Guarnaccia, E., & Marfatia, D. 2006, PhLB, 635, 61

Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32

Behroozi, P. S., & Silk, J. 2018, MNRAS, 477, 5382

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Benson, A. J. 2012, NewA, 17, 175

Benson, A. J., Farahi, A., Cole, S., et al. 2013, MNRAS, 428, 1774

Bhatawdekar, R., Conselice, C. J., Margalef-Bentabol, B., & Duncan, K. 2019, MNRAS, 486, 3805

Binetruy, P. 1999, PhRvD, 60, 063502

Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440

Brax, P., & Martin, J. 1999, PhLB, 468, 40

Caldwell, R. R. 2002, PhLB, 545, 23

Caldwell, R. R., Dave, R., & Steinhardt, P. J. 1998, PhRvL, 80, 1582

Caldwell, R. R., & Linder, E. V. 2005, PhRvL, 95, 141301

Caputi, K. I., Ilbert, O., Laigle, C., et al. 2015, ApJ, 810, 73

Carrol, S. M., Press, W. H., & Turner, E. L. 1992, ARA&A, 30, 499

Casey, C. M., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Chevallier, M., & Polarski, D. 2001, IMPD, 10, 213

Chiba, T. 2006, PhRvD, 73, 063501

Chiba, T. 2009, PhRvD, 79, 083517

Combes, F., Garciéa-Burillo, S., Braine, J., et al. 2013, A&A, 550, A41

Copeland, E. J., Nunes, N. J., & Rosati, F. 2000, PhRvD, 62, 123503

Copeland, E. J., Sahni, M., & Tsujikawa, S. 2006, IMPD, 15, 1753

Corasaniti, P. S., & Achitouv, I. E. 2011a, PhRvD, 84, 023009

Corasaniti, P. S., & Achitouv, I. E. 2011b, PhRvL, 106, 241302

Croom, S. M., & Shanks, T. 1999, MNRAS, 303, 411

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Daveé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Despali, G., Giocoli, R. E., Angulo, R. E., et al. 2016, MNRAS, 456, 2486

Destri, C., de Vega, P., & Sanchez, N. G. 2013, PhRvD, 88, 3512

Di Valentino, E., Melchiorri, A., Linder, E. V., & Silk, J. 2017, PhRvD, 96, 023523

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015, ApJ, 814, 95

Frieman, J. A., Turner, S., & Huterer, D. 2008, ARA&A, 46, 385

Fu, J., Kauffmann, G., Li, C., & Guo, Q. 2012, MNRAS, 424, 2701

Gabor, J. M., & Bournaud, F. 2013, MNRAS, 434, 606

Geach, J. E., Smail, I., Moran, S. M., et al. 2011, ApJL, 730, L19

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

Ginolfi, M., Schneider, R., & Valiante, R. 2019, MNRAS, 483, 1256

Grazian, A., Fontana, A., Santini, p., et al. 2015, A&A, 575, 96

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gupta, G., Rangarajan, R., & Sen, A. A. 2015, PhRvD, 92, 123003

Harrison, I., & Hotchkiss, S. 2013, JCAP, 7, 022

He, W., Akiyama, M., Bosch, J., et al. 2018, PASJ, 70, S33, (SP1)

Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146

Hildebrandt, H., Pielorz, J., Erben, T., et al. 2009, A&A, 498, 725

Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 471, 4412

Hlozek, R., Dunkley, J., Addison, G., et al. 2012, ApJ, 749, 90

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Lagos, C. d. P., Baugh, C. M., Lacey, C. G., et al. 2011, MNRAS, 418, 1649

Lagos, C. d. P., Baugh, C. M., Zwaan, M. A., et al. 2014, MNRAS, 440, 920

Lagos, C. d. P., Bayet, E., Baugh, C. M., et al. 2012, MNRAS, 426, 2142

Lamstra, A., Menci, N., Fiore, F., Di Porto, C., & Amendola, L. 2012, MNRAS, 420, 2429

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Lee, K.-S., Ferguson, H. C., Wiklund, T., et al. 2012, ApJ, 752, 66

Linder, E. V. 2003, PhRvL, 90, 091301

Linder, E. V. 2005, PhRvD, 72, 043529

Linder, E. V. 2006, PhRvD, 73, 063010

Linder, E. V., & Cahn, R. N. 2007, APh, 28, 481

Magdis, G. E., Daddi, E., Sargent, M., et al. 2012, ApJL, 758, L9

Maggiore, M., & Riotto, A. 2010, ApJ, 717, 515

Mainini, R., Maccio, A. V., Bonometto, S. A., & Klypin, A. 2003, ApJ, 599, 24

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Natur, 553, 51

Masiero, A., Pietroni, M., & Rosati, F. 2000, PhRvD, 61, 023504

Massey, R., Rhodes, J., Leauthard, A., et al. 2007, ApJS, 172, 239

Mehrabи, M. 2018, PhRvD, 97, 083529

Menci, N., Grazian, A., Castellano, M., & Sanchez, N. G. 2016, ApJL, 825, L1

Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309

Mo, H. J., & White, S. D. M. 2002, MNRAS, 336, 112

Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822

Narayanan, D., Bothwell, M., & Daveé, R. 2012, MNRAS, 426, 1178

Pace, F., Waizmann, J.-C., & Bartelmann, M. 2010, MNRAS, 406, 1865

Pan, S., Yang, W., Di Valentino, E., Saridakis, E. N., & Chakraborty, S. 2019, PhRv, 100, 103520

Papadopoulos, P. P., van der Werf, P. P., Xilouris, E. M., et al. 2012, MNRAS, 426, 2601

Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton, NJ: Princeton Univ. Press)

Peebles, P. J. E., & Ratra, B. 1988, ApJL, 325, L17

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Refregier, A. 2003, ARA&A, 41, 645

Riess, A. G. 2019, NatRP, 2, 10

Risaliti, G., & Lusso, E. 2019, NatAs, 3, 272

Sahlen, M., Liddle, A. R., & Parkinson, D. 2007, PhRvD, 75, 023502

Sahni, V., & Wang, L. M. 2000, PhRvD, 62, 103517

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Scherrer, R. J. 2006, PhRvD, 73, 043502

Scherrer, R. J., & Sen, A. A. 2008, PhRvD, 77, 08351515

Schneider, A., Smith, R. E., & Reed, D. 2013, MNRAS, 433, 1573

Scolnic, D. M., Jones, D. O., & Rest, A. 2018, ApJ, 859, 101

Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Springel, V. 2010, MNRAS, 401, 791

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2017, ApJ, 843, 36

Stefanon, M., Marchesini, D., Muzzin, A., et al. 2015, ApJ, 803, 11

Steinhardt, C. L., Capak, P., Masters, D., & Speagle, J. S. 2016, ApJ, 824, 21

Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJL, 842, L15

Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Natur, 463, 781

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tegmark, M., & Zaldarriaga, M. 2002, PhRvD, 66, 103508
Tegmark, M., & Zaldarriaga, M. 2009, PhRvD, 79, 083530

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. A. 2013, PhRvD,
88, 3502

Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, MNRAS, 444, 1518

Vogelsberger, M., Genel, S., Springel, V., et al. 2014b, Natur, 509, 177

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Natur, 572, 211

Warren, M. S., Abazajian, K., Holz, D. E., & Teodoro, L. 2006, ApJ, 646, 881

Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955

Weinberg, S. 1972, Gravitation and Cosmolgy: Principles and Applications of the General Theory of Relativity (New York: Wiley)

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Zhai, Z., Blanton, M., Slosar, A., & Tinker, J. 2017, MNRAS, 850, 183

Zhao, Gong-Bo, Raveri, Marco, Pogosian, Levon, et al. 2017, NatAs, 1, 627