The impact of dark energy on galaxy formation. What does the future of our Universe hold?

Abstract

We investigate the effect of the accelerated expansion of the Universe due to a cosmological constant, Λ, on the cosmic star formation rate. We utilize hydrodynamical simulations from the eagle suite, comparing a ΛCDM (cold dark matter) Universe to an Einstein–de Sitter model with Λ = 0. Despite the differences in the rate of growth of structure, we find that dark energy, at its observed value, has negligible impact on star formation in the Universe. We study these effects beyond the present day by allowing the simulations to run forward into the future (t > 13.8 Gyr). We show that the impact of Λ becomes significant only when the Universe has already produced most of its stellar mass, only decreasing the total comoving density of stars ever formed by |${\approx }15\hbox{ per cent}$|⁠. We develop a simple analytic model for the cosmic star formation rate that captures the suppression due to a cosmological constant. The main reason for the similarity between the models is that feedback from accreting black holes dramatically reduces the cosmic star formation at late times. Interestingly, simulations without feedback from accreting black holes predict an upturn in the cosmic star formation rate for t > 15 Gyr due to the rejuvenation of massive (>10¹¹ M_⊙) galaxies. We briefly discuss the implication of the weak dependence of the cosmic star formation on Λ in the context of the anthropic principle.

cosmology: theory, galaxies: evolution, galaxies: formation

1 INTRODUCTION

Precise observational data from the past two decades have allowed us to measure the cosmic history of star formation back to very early times (z ≈ 8). The star formation rate (SFR) density of the Universe peaked approximately 3.5 Gyr after the Big Bang (z ≈ 2), and declined exponentially thereafter (for a review, see Madau & Dickinson 2014).

Galaxy formation and evolution are a highly self-regulated process, in which galaxies tend to evolve towards a quasi-equilibrium state where the gas outflow rate balances the difference between the gas inflow rate and the rate at which gas is locked up in stars and black holes (BHs) (e.g. White & Frenk 1991; Finlator & Davé 2008; Schaye et al. 2010; Davé, Finlator & Oppenheimer 2012). Consequently, the cosmic SFR density is thought to be determined both by the formation and growth of dark matter haloes and by the regulation of the gas content in these haloes. The former depends solely on cosmology, whereas the latter depends on baryonic processes such as radiative cooling, stellar mass loss, and feedback from stars and accreting BHs.

Which of these factors is most responsible for the decline in cosmic star formation? It could be driven by the ‘freeze out’ of the growth of large-scale structure, caused by the onset of accelerating cosmic expansion. As galaxies are driven away from each other by the repulsive force of dark energy, accretion and merging are slowed and galaxies are gradually starved of the raw fuel for star formation. Or, it could be caused primarily by the onset of efficient stellar and BH feedback.

The discovery of the accelerating expansion of the Universe was a breakthrough achievement for modern cosmology (Riess et al. 1998; Perlmutter et al. 1999). However, the driving force behind the acceleration (generically known as dark energy) is still unknown. At present, all cosmological observations are consistent with a cosmological constant, or a form of energy whose density remains constant as the Universe expands. One such form of energy is vacuum energy: the energy of a quantum field in its ground state (zero particles). The present best-fitting cosmological model, known as the concordance model, or ΛCDM, includes both a cosmological constant Λ and cold (i.e. non-relativistic) dark matter. This model has been very successful in matching the observational data.

Nevertheless, the model raises a number of fundamental problems. Predictions from quantum field theory for the vacuum energy density overestimate the observed value of Λ by many orders of magnitude (for a review, see Weinberg 1989; Carroll 2001). In addition, the energy density of matter and the cosmological constant are within a factor of a few of each other at the present time, making our epoch unusual in the evolution of the Universe. This is known as the coincidence problem. These problems have motivated the search for alternative models of dark energy and modifications of gravity that might explain the acceleration of the Universe more naturally. For example, quintessence models propose that the density of matter and dark energy track each other. In many models, however, fine-tuning of the model parameters is still required to explain their observed similarity (see e.g. Weinberg 2000).

An alternative approach is therefore to explain the observed value of Λ on anthropic grounds. This has already been applied very promisingly to the coincidence problem. Since the coincidence concerns the time that we observe the universe, the nature and evolution of observers in the Universe are highly relevant. For example, Lineweaver & Egan (2007) argue that the production of planets in our Universe peaks when matter and dark energy are roughly coincident (see, however, Loeb, Batista & Sloan 2016).

For the cosmological constant and other fundamental parameters, anthropic reasoning requires a multiverse. Many models of inflation, such as eternal inflation, imply that the Universe as a whole is composed of a vast number of inflationary patches or sub-universes. Each sub-universe inherits a somewhat random set of physical constants and cosmic parameters from a wide range of possible values. Sub-universes in which the cosmological constant is large and positive will expand so rapidly that gravitational structures, such as galaxies, are unable to form (e.g. Weinberg 1987; Martel, Shapiro & Weinberg 1998; Efstathiou 1995; Sudoh et al. 2017). Large negative values will cause the universe to collapse rapidly, also preventing the formation of galaxies. Only sufficiently small values of Λ will lead to the formation of universes that are able to host observers. This argument eliminates extreme values of Λ. For example, Weinberg (2000) estimates an upper bound on a positive vacuum energy density to allow for the formation of galaxies of about 200 times the present mass density.

Refining Weinberg’s estimate requires us to more accurately explore the sensitivity of galaxy formation to the presence of Λ. Here, we use a suite of hydrodynamical simulations to take a first look at this problem by calculating the effect of the cosmological constant on galaxy and star formation in our Universe. Specifically, we compare the formation of galaxies in our Universe with a hypothetical universe that is indistinguishable from ours at early times but has no cosmological constant. Because Λ is negligibly small in the early universe, these two universes will evolve in nearly identical ways for the first ≈2 Gyr of cosmic time (when the dark energy density is less than 0.03 times the matter density). This means that the epochs of nucleosynthesis, recombination,¹ and reionization are indistinguishable.

In recent years, the accuracy of our understanding of galaxy formation has improved considerably, reaching the point at which it is possible to undertake this comparison meaningfully. The increased realism of simulated galaxies (in particular disc galaxies with more realistic sizes and masses) has been achieved due to the use of physically motivated subgrid models for feedback processes (e.g. Schaye et al. 2015; Dubois et al. 2016; Pillepich et al. 2017). One of the key ingredients that have allowed this progress is the inclusion of realistic models for the impact of feedback from the growth of super massive BHs (e.g. Bower et al. 2017). All successful models now demonstrate the need for active galactic nuclei (AGNs) as an additional source of feedback that suppresses the formation of stars in high-mass haloes (e.g. Benson et al. 2003; Bower et al. 2006; Croton et al. 2006; Crain et al. 2015; Pillepich et al. 2017). One of the aims of this paper is to compare the impact of the cosmological constant with that resulting from the inclusion of BHs in the simulation. In a previous study, van de Voort et al. (2011) found that by preventing gas from accreting on to the central galaxies in massive haloes, outflows driven by AGN play a crucial role in the decline of the cosmic SFR.

Different groups have used hydrodynamical simulations to study the effect of different dark energy or modified gravity models on cosmological, galactic, and sub-galactic scales (e.g. Puchwein, Baldi & Springel 2013; Penzo et al. 2014, 2016). Taking a different approach, in this paper we investigate the effect of the accelerated expansion of the Universe on galaxy formation by asking the following question:

How different would the Universe be if there had been no dark energy?

For our study, we use a suite of large hydrodynamical simulations from the Evolution and Assembly of GaLaxies and their Environment (eagle) project (Crain et al. 2015; Schaye et al. 2015). Using state-of-the-art subgrid models for radiative cooling, star formation, stellar mass loss, and feedback from stars and accreting BHs, the simulations have reproduced many properties of the observed galaxy population and the intergalactic medium both at the present day and at earlier epochs (e.g. Furlong et al. 2015, 2017; Lagos et al. 2015; Rahmati et al. 2015, 2016; Schaller et al. 2015; Trayford et al. 2015; Bahé et al. 2016; Rosas-Guevara et al. 2016; Segers et al. 2016). Given that the physics of the real Universe is reasonably well captured by the phenomenological sub-grid models implemented in the simulations, with the use of appropriate assumptions, we can run the simulations beyond the present time, and explore the consequences of our models for the future. Furthermore, the simulation re-scaling strategy developed here will be used in a companion paper (Barnes et al. 2018), considering a wider range of Λ values and determining the likelihood distribution of possible Λ values conditioning the existence of observers.

The layout of this paper is as follows. In Section 2, we develop a simple analytic model of the cosmic SFR that captures the suppression due to a cosmological constant. In Section 3, we briefly describe the simulations from which we derive our results and discuss our criteria for halo and galaxy definitions. In Section 3, we also describe our motivations to run our cosmological simulations into the future, and our assumptions in doing so. Section 4 provides a detailed discussion of our re-scaling strategy for the alternative cosmological models. In Section 5, we explore the dependence of the star formation history of the universe on the existence of a cosmological constant and the presence of BHs. We also explore their impact on other galaxy population properties, both up to the present time, and into the future. Finally, we summarize and discuss our results in Section 6.

2 A SIMPLE ANALYTIC MODEL FOR THE COSMIC SFR DENSITY

2.1 Comparing different cosmological models

The star formation history of the Universe is determined by the interplay of cosmic expansion and the time-scale at which cold gas can turn into stars. These processes happen on time-scales that differ by several orders of magnitude, but are coupled through the accretion rate of gas on to gravitationally bound haloes. The aim of our paper is to compare theoretical universes in which the star formation time-scales are the same, but the cosmological time-scales vary. We need, therefore, to be careful when comparing the different models, since the choice of coordinates that vary with cosmological parameters will obscure the similarities of the models. In particular, the expansion factor at the present day, a₀, is often treated as an arbitrary positive number, and it is common practice to set a₀ = 1. In this paper, we need to take a different approach since we want to compare the properties of different universes at the same cosmic time (measured in seconds or a multiple of key atomic transitions). Assuming a common inflationary origin, normalizing out a₀ is not appropriate, since the expansion factor at the present day (t₀ = 13.8 Gyr) would be different for each universe.

We still need to define a scale on which to measure the size of the universes we consider. Using a hat notation |$(\, \hat{}\,)$| to denote quantities in our observable Universe, we set |$\hat{a}_0 \equiv \hat{a} (t_0) =1$|⁠. We want to emphasize that the cosmological models that we consider all start from very similar initial conditions. It therefore makes sense to normalize them to the same value of the expansion factor at an early time, t₁. We therefore set |$a_1 \equiv a(t_1) = \hat{a}(t_1)$|⁠. We choose |$\hat{a}_1 = 1/(1+127)$|⁠, corresponding to a redshift of |$\hat{z}=127$| for a present-day observer in our Universe.² At this moment, the age of the universe is t₁ = 11.98 Myr. This applies to all of the universes we consider since the cosmological constant term has negligible impact on the expansion rate at such early times.

Although time (in seconds) is the fundamental coordinate that we use to compare universes, it is sometimes useful, for example when comparing to observational data, to express time in terms of the redshifts measured by a present-day observer in our Universe, |$\hat{z}$|⁠. We convert between cosmic time t (which is equivalent between universes) and |$\hat{a}$| by inverting the time-redshift relation for our Universe:

\begin{eqnarray*} \hat{z} = \frac{\hat{a}_0}{\hat{a}(t)} - 1. \end{eqnarray*}

(1)

It is important to note that |$\hat{z}$| is not the redshift that would be measured by an observer in an alternative universe.

In this paper, we will focus our comparison on two cosmological models, a standard ΛCDM universe as inferred by Planck Collaboration XVI (2014), and an Einstein deSitter (EdS) universe. Assuming both cosmological models have a common inflationary origin, the models can be normalized as follows:

• For the ΛCDM model (see Table 1), we set |$\hat{a}_0 = \hat{a}(t_0)=1$|⁠, where t₀ = 13.82 Gyr is the present-day age of the universe. At time t₁ = 11.98 Myr, |$\hat{a}_1 = \hat{a} (t_1) = 0.007813$|⁠. At this time, the expansion rate, as measured by the Hubble parameter is, |$\hat{H}_1 = \hat{H} (t_1)= 54,377$| km s⁻¹ Mpc⁻¹ =55.6 Gyr⁻¹.

• We require the EdS model to have the same early expansion history, i.e. a(t₁) = 0.007813 and H(t₁) = 55.6 Gyr⁻¹. In this universe, at the present day (i.e. t = t₀ = 13.82 Gyr) the universe has a size, a₀ = a(t₀) = 0.8589 and an expansion rate, H(t₀) = 0.0482 Gyr⁻¹ = 47.16 km s⁻¹ Mpc⁻¹.

Table 1.

The cosmological parameters for the eagle simulations used in this study. ΛCDM model refers to parameters inferred by Planck Collaboration (2014). EdS refers to an Einstein–de Sitter universe. |$\Omega_\mathrm{m}, \Omega_\Lambda , \Omega_\mathrm{b}$| are the average densities of matter, dark energy, and baryonic matter, respectively, in units of the critical density at redshift zero; H₀ is the Hubble constant, σ₈(t₁) is the square root of the linear variance of the matter distribution at the initial cosmic time of the simulations (t₁ = 11.98 Myr) when smoothed with a top-hat filter of radius 11.8 cMpc (8 h⁻¹ cMpc for a ΛCDM model), n_s is the scalar power-law index of the power spectrum of primordial adiabatic perturbations, and Y is the primordial abundance of helium. Values in bold show differences with respect to the ΛCDM values.

Cosmological Parameter	ΛCDM (Ref)	EdS
Ωm	0.307	1
|$\Omega_\Lambda$|	0.693	0
Ωb	0.048 25	0.157 17
h ≡ H0/(100 km s−1 Mpc−1)	0.6777	0.3754
σ8(t1)	0.0083	0.0083
ns	0.9611	0.9611
Y	0.248	0.248

Fig. 1 shows the cosmic scale factor as a function of time for the two cosmological models. As expected, at t = t₀, an EdS universe is smaller in size at the present day, as the cosmic expansion has not been accelerated by the effect of Λ. As the two universes evolve into the future, the size differences and relative expansion rates grow, e.g. at t = 20 Gyr, the scale factor for the ΛCDM models is |${\approx }25\hbox{ per cent}$| larger than for the EdS, and the expansion rate is |${\approx }50\hbox{ per cent}$| larger for our Universe.

Figure 1.

Cosmic scale factor as a function of time for two cosmological models. The model for the cosmological parameters for a standard ΛCDM universe as inferred by Planck Collaboration XVI (2014) is shown in blue. An Einstein–de Sitter universe is shown in orange. Note that by construction the scale factors are indistinguishable when the universes are less than 1 Gyr old. The power series approximation of equation (10) is shown with a dashed green line.

2.2 Cosmological expansion history as a function of time

In the standard model of cosmology for a homogeneous and isotropic universe, the geometry of space–time is determined by the matter-energy content of the universe through the Einstein field equations as described by the Friedmann–Lemaître–Robertson–Walker metric in terms of the scale factor a(t) and the curvature K, yielding the well-known Friedmann equation,

\begin{eqnarray*} \left(\frac{\dot{a}}{a}\right)^2 = H^2(t)= \frac{8 \pi G}{3} \rho - \frac{K c^2}{a^2} + \frac{\Lambda c^2}{3}, \end{eqnarray*}

(2)

where H(t) is the Hubble parameter. As the inflationary models predict that the Universe should be spatially flat, we only consider universes with no spatial curvature, i.e. K = 0.

The density of equation (2) includes the contribution of non-relativistic matter and radiation (ρ_m and ρ_r). The radiation content of the Universe dominated its global dynamics at very early times (a → 0), but its contribution is negligible thereafter. Ignoring ρ_r and using the energy density at an arbitrary time t₁, equation (2) can be written as

\begin{eqnarray*} \left(\frac{\dot{a}}{a}\right)^2 = \frac{8 \pi G}{3} \rho_{\mathrm{m},1} \left(\frac{a}{a_1}\right)^{-3} + \frac{\Lambda c^2}{3}, \end{eqnarray*}

(3)

where ρ_{m, 1} is the matter density of the universe at t = t₁, and a₁ = a(t₁). We choose t₁ such that it corresponds to a sufficiently early epoch, when the contribution of the cosmological constant term is negligible. As discussed in the previous section, at this time any universe closely approximates an EdS universe and we can assume that |$a_1 = \hat{a}_1$| and |$\rho_{\mathrm{m},1} =\hat{\rho }_{\mathrm{m},1} \Rightarrow \hat{\rho }_{\mathrm{m},0} (\hat{a}_0 / \hat{a}_1)^3 = \rho_{\mathrm{m},0} (a_0 / a_1)^3$|⁠. Then, equation (3) can be written as

\begin{eqnarray*} \left(\frac{\dot{a}}{a}\right)^2 = \frac{8 \pi G}{3} \hat{\rho }_{\mathrm{m},0} \left(\frac{a}{\hat{a}_0}\right)^{-3} + \frac{\Lambda c^2}{3}. \end{eqnarray*}

(4)

Note that in equation (4), the evolution of the scale factor for any arbitrary cosmology is written in terms of the matter density of our Universe at the present time |$\hat{\rho }_{\mathrm{m},0}$|⁠. We have left the factor of |$\hat{a}_0$| explicit in the equation, but it can be set to |$\hat{a}_0 = 1$|⁠, noting that a₀ ≠ 1 for any cosmological model different to our Universe.

The LHS of equation (4) has units of time⁻² and we will later find it useful to represent the RHS as the sum of two time-scales. The cosmological constant is often written as an energy component with energy density |$\rho_\Lambda = \Lambda c^2/8\pi G$|⁠, however, we can express this as a time-scale as follows:

\begin{eqnarray*} t_\Lambda = \sqrt{\frac{3}{\Lambda c^2}} = \frac{1}{H_0 \sqrt{\Omega_{\Lambda ,0}}}. \end{eqnarray*}

(5)

Similarly, the matter content of the Universe can be expressed as a time-scale,

\begin{eqnarray*} t_\mathrm{m} = \sqrt{\frac{3}{8\pi G\hat{\rho }_0}} = \frac{1}{\hat{H}_0 \sqrt{\hat{\Omega }_{\mathrm{m},0}}}. \end{eqnarray*}

(6)

Using the cosmological parameters for our Universe, |$t_\Lambda =\hat{t}_\Lambda =17.33$| Gyr and |$t_m = \hat{t}_m= 26.04$| Gyr. For an EdS universe, |$t_\Lambda \rightarrow \infty$|⁠.

Using this notation, equation (2) can be written as

\begin{eqnarray*} \left(\frac{\dot{a}}{a}\right)^2 = t_m^{-2} \left(\frac{a}{\hat{a}_0}\right)^{-3} + t_\Lambda ^{-2}, \end{eqnarray*}

(7)

which can be solved analytically to express the expansion factor as a function of time and the parameters t_m and |$t_\Lambda$|⁠:

\begin{eqnarray*} a(t) = \left[\frac{1}{2} {\rm e}^{-3{t}/{2t_\Lambda }} \left( {\rm e}^{3{t}/{t_\Lambda }} -1 \right) \left(\frac{t_\Lambda }{t_m}\right)\right]^{2/3} \end{eqnarray*}

(8)

In the limit |$t_\Lambda \rightarrow \infty$| this reduces to the familiar EdS solution,

\begin{eqnarray*} \lim_{t_\Lambda \rightarrow \infty } a(t) = \left[\frac{3}{2} \frac{t}{t_m}\right]^{2/3} \end{eqnarray*}

(9)

In order to explore the significance of the |$t/t_\Lambda$| term more clearly, we can expand equation (8) as a Taylor series:

\begin{eqnarray*} a(t) \approx \left[\frac{3}{2} \frac{t}{t_m}\right]^{2/3}\left(1 + \frac{1}{4} \left(\frac{t}{t_\Lambda }\right)^2 + \frac{1}{80} \left(\frac{t}{t_\Lambda }\right)^4 +\cdots \right) \end{eqnarray*}

(10)

The coefficients of the series decreased rapidly so that the first three terms provide a good approximation up to |$t = 2 t_\Lambda$| and beyond. Fig. 1 shows how well this power series approximation works.

2.3 The growth of density perturbations

In the standard model of cosmology, structures such as galaxies and clusters of galaxies are assumed to have grown from small initial density perturbations. Expressing the density, ρ, in terms of the density perturbation contrast against a density background,

\begin{eqnarray*} \rho (\mathbf {x},t) = \bar{\rho }(t)[1+\delta (\mathbf {x},t)], \end{eqnarray*}

(11)

the differential equation that governs the time dependence of the growth of linear perturbations in a pressureless fluid, such as dark matter, can be written as (for a review, see Peebles 1980; Mo, van den Bosch & White 2010),

\begin{eqnarray*} \frac{{\rm d}^2 \delta }{{\rm d}t^2} + 2 \frac{\dot{a}}{a} \frac{{\rm d}\delta }{{\rm d}t} - 4\pi G \bar{\rho } \delta = 0. \end{eqnarray*}

(12)

The growing mode of equation (12) can be written as

\begin{eqnarray*} \delta (t) = {D(t)} \delta (t_0), \end{eqnarray*}

(13)

where D(t) is the linear growth factor, which determines the normalization of the linear matter power spectrum relative to the initial density perturbation power spectrum, and is computed by the integral

\begin{eqnarray*} D(t) \propto \frac{\dot{a}}{a} \int_0^t \frac{{\rm d}t^\prime }{\dot{a}^2(t^\prime )}. \end{eqnarray*}

(14)

Using the hat notation as before, we normalize D(t) so that

• |$\hat{D}_+(t_0) = 1$|

• |$D(t_1) = \hat{D}(t_1)$|

In general, the growing mode can be obtained from equation (14) numerically. Fig. 2 shows the growth factor as a function of cosmic time for the two cosmological models. As expected, the figure shows that linear perturbations grow faster in an EdS universe compared to those in a ΛCDM universe.

Figure 2.

The linear growth factor for the ΛCDM and EdS cosmological models. The rates are all normalized such that |$\hat{D}(t_0) = 1$|⁠, for the ΛCDM model, and |$D(t_1) = \hat{D}(t_1)$|⁠. The bottom panel shows the growth factor at a given time, divided by the growth factor for ΛCDM. The presence of a cosmological constant suppresses the growth of structure in the ΛCDM model (blue) compared to that in the EdS model (orange). The power series approximation of equation (15) is shown with a dashed green line.

It is possible to gain more insight by integrating the power-series approximation for a(t) from equation (10). Expanding the solution again as a power series in |$(t/t_\Lambda )$|⁠, retaining the leading terms, yields,

\begin{eqnarray*} D(t) &=& \left[\frac{3}{2} \frac{t}{t_m}\right]^{2/3} \frac{2}{5} t_m^{2} K_D\nonumber\\ &&\times \,\left(1 - 0.1591 \left(\frac{t}{t_\Lambda }\right)^2 + 0.0366 \left(\frac{t}{t_\Lambda }\right)^4\right),\end{eqnarray*}

(15)

where K_D is a normalization constant. Requiring |$\hat{D}(t_0) = 1$| gives K_D = 4.70 × 10⁻³ Gyr⁻². Fig. 2 shows that equation (15) provides a good approximation up to |$t = t_\Lambda$|⁠.

This demonstrates that although the |$t_\Lambda$| term slows down the growth of perturbations, its effect is less than 10 per cent until |$t \sim t_\Lambda \left( 0.1 / 0.1591 \right)^{1/2} \approx 0.8 t_\Lambda$| corresponding to ≈13.8 Gyr (⁠|${\approx }\hat{t}_0$|⁠) in our Universe.

As we discuss in the following section, the quantity of fundamental interest for the accretion rate of dark matter haloes is the relative rate of growth of density perturbations, |$\frac{1}{D}\frac{{\rm d}D}{{\rm d}t}$|⁠. We show this for the numerical solution in Fig. 3. We can also compute the relative growth rate by differentiating the power-series approximation of equation (15). Retaining the lowest order terms, we find

\begin{eqnarray*} \frac{1}{D}\frac{{\rm d}D}{{\rm d}t} = \frac{2}{3t}\left(1 - 0.4773 \left(\frac{t}{t_\Lambda }\right)^2 + 0.1435\left(\frac{t}{t_\Lambda }\right)^4\right) \end{eqnarray*}

(16)

This expression does not depend on the constants t_m or K_D because we are focusing on the relative change in the growth factor. The impact of the cosmological constant term is relatively large, creating an |${\approx }50\hbox{ per cent}$| increase in growth rate for the EdS model compared to ΛCDM when |$t\approx t_\Lambda$|⁠.

Figure 3.

The relative rate of growth of density perturbations, |$\frac{1}{D}\frac{{\rm d}D}{{\rm d}t}$| for the ΛCDM and EdS cosmological models. The bottom panel shows the ratio, at a given time. The presence of a cosmological constant slows down the growth of structure in the ΛCDM model (blue) compared to that in the EdS model (orange). The power series approximation of equation (16) is shown with a dashed green line.

2.4 Impact on halo accretion rates

The growth rates of linear perturbations do not directly predict the growth rates of haloes; however, we can directly connect the two through the approach developed by Press & Schechter (Press & Schechter 1974; Bond et al. 1991; Bower 1991; Lacey & Cole 1993). Correa et al. (2015) showed that the accretion rates of haloes can be written as (see also Neistein, van den Bosch & Dekel 2006),

\begin{eqnarray*} \frac{1}{M_h} \frac{{\rm d}M_h}{{\rm d}t} = \sqrt{\frac{2}{\pi }} \frac{(\delta_c/D)}{S(M_h)^{1/2}\left(q^\gamma -1\right)^{1/2}} \frac{1}{D}\frac{{\rm d}D}{{\rm d}t}, \end{eqnarray*}

(17)

where M_h is the halo mass and S(M_h) is the variance of the density field on the length scale corresponding the halo mass. δ_c is a parameter that represents a threshold in the linearly extrapolated density field for halo collapse. The parameters, q and γ, are related to the shape of the power spectrum around the halo mass M_h. Approximating the scale dependence of the density field as a power law, |$S=S_0 M_h^{-\gamma}$|⁠, Correa et al. (2015) find S ≈ 3.98, γ ≈ 0.3, and q ≈ 3.16, giving [S(M_h)(q^γ − 1)]^−1/2 ≈ 0.78 for 10¹² M_⊙ haloes. These values depend only on the initial power spectrum (which we assume to be the same in all the universes we consider) and do not depend on the cosmological parameters. This formulation thus neatly separates the contribution of the power-spectrum shape from the cosmological parameters. We are therefore able to assume that q and γ are the same for all the universes that we consider, and focus on the dependence on D(t).

For the numerical values of the power-spectrum parameters around a halo mass of 10¹² M_⊙, equation (17) reduces to

\begin{eqnarray*} \frac{1}{M_h} \frac{{\rm d}M_h}{{\rm d}t} = 1.0456 \frac{1}{D^2} \frac{{\rm d}D}{{\rm d}t}. \end{eqnarray*}

(18)

This dependence can be understood as the combination of two factors. The first reflects the relative growth rate of density fluctuations |$\frac{1}{D}\frac{{\rm d}D}{{\rm d}t}$|⁠. The second factor of 1/D comes from the rarity of haloes, reflecting the higher growth rate of fluctuations in the tail of the density field distribution.

Further insight can be gained by using the series approximation. This gives

\begin{eqnarray*} \frac{1}{M_h} \frac{{\rm d}M_h}{{\rm d}t} = \frac{566.61}{\sqrt{S} \, t^{5/3}t_m^{4/3}} \left( 1 - 0.3182 \left(\frac{t}{t_\Lambda }\right)^2 + 0.0563 \left(\frac{t}{t_\Lambda }\right)^4 \right). \nonumber\\ \end{eqnarray*}

(19)

This explicitly shows how the presence of a cosmological constant modulates the halo growth rate. In our Universe, the impact of the cosmological constant term is relatively modest, however; at t = 13.8 Gyr, we expect the difference to be 20 per cent, growing to 40 per cent at t = 20 Gyr.

As an example, in Fig. 4 we show the accretion rate of haloes of M_h = 10¹² M_⊙, both numerically and using equation (19).

Figure 4.

The specific accretion rate of haloes of mass M_h = 10¹² M_⊙, |$\frac{1}{M_h} \frac{{\rm d}M_h}{{\rm d}t}$| for the ΛCDM and EdS cosmological models. The bottom panel shows the ratio at a given time. The presence of a cosmological constant slows down the specific accretion rates of haloes in the ΛCDM model (blue) compared to that in the EdS model (orange). The power series approximation of equation (19) is shown with a dashed green line.

2.5 Impact on the SFR of the Universe

In order to link the SFR of haloes of mass M_h to their accretion rate, as a first approximation, we assume a time-independent galaxy specific SFR to host halo specific mass accretion rate relation (e.g. Behroozi, Wechsler & Conroy 2013a; Tacchella, Trenti & Carollo 2013; Rodríguez-Puebla et al. 2016),

\begin{eqnarray*} \frac{\dot{M}_*/{M}_*}{\dot{M}_h/{M}_h} = \frac{\partial \mathrm{log} M_*}{\partial \mathrm{log} M_h} = \epsilon (M_h), \end{eqnarray*}

(20)

where the star formation efficiency ε, of haloes of mass M_h, is the slope of the stellar-halo mass relation. From this equation, the star formation as a function of halo mass can be written as

\begin{eqnarray*} \dot{M}_*(M_h) = \epsilon_*(M_h) \dot{M}_h, \end{eqnarray*}

(21)

where ε_*(M_h) ≔ ε(M_h) × (M_*/M_h) is completely defined by the stellar–halo mass relation. As there is no a priori knowledge of the functional form of ε_*(M_h), we use the abundance matching results from Behroozi, Wechsler & Conroy (2013b) to estimate ε_*(M_h). The efficiency ε_*(M_h) peaks at masses similar to Milky Way-sized haloes (∼10¹² M_⊙) and falls steeply for higher and lower masses. ε_*(M_h) can be well approximated by a broken power law as

\begin{eqnarray*} \epsilon_*({M_h})\propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\left(\frac{M_h}{10^{12}\mathrm{M}_{\odot }}\right)^1 & \text{if $M_h \le 10^{12}\mathrm{M}_{\odot } $} \\\\\left(\frac{M_h}{10^{12}\mathrm{M}_{\odot }}\right)^{-1} & \text{if $M_h \,\,{>}\,\,10^{12}\mathrm{M}_{\odot } $} \\\end{array}\right. \end{eqnarray*}

(22)

At low masses, SFR is suppressed because of the efficiency of feedback from star formation, at higher masses the cooling of the inflowing gas is suppressed by heating from BHs (White & Frenk 1991; Benson et al. 2003; Bower et al. 2006; Haas et al. 2013; Crain et al. 2015; Dubois et al. 2016; Bower et al. 2017).

In order to complete the analysis, we need to combine the specific halo mass accretion rate with an estimate of the halo abundance.

In the Press & Schechter analysis, the comoving abundance of haloes of mass M_h at time t is given by (Press & Schechter 1974; Bond et al. 1991; Bower 1991; Lacey & Cole 1993),

\begin{eqnarray*} \frac{{\rm d}n(M_h,t)}{{\rm d}M_h} = \frac{\hat{\rho }_0}{M_h^2} \frac{\delta_c \gamma }{\sqrt{2\pi } S^{1/2}} \frac{1}{D} \exp \left( - \frac{\delta_c^2}{2 S D^2}\right) \end{eqnarray*}

(23)

where we have assumed that the density power spectrum is a power law with exponent γ and written the comoving density of the Universe as |$\hat{\rho }_0$| following our convention. Note that we compute comoving densities. At the same cosmic time, the different expansion rates will result in different physical (proper) halo and SFR densities, simply because of the more rapid expansion of the ΛCDM cosmology.

The total cosmic SFR density is given by the integral of all star formation in all haloes,

\begin{eqnarray*} \dot{\rho }_*(t) &=& \int \dot{M}_*(M_h) \frac{{\rm d}n(M_h,t)}{{\rm d}M_h} {\rm d}M_h\nonumber\\ &=& \int \epsilon_*(M_h) \dot{M}_h \frac{{\rm d}n(M_h,t)}{{\rm d}M_h} {\rm d}M_h \end{eqnarray*}

(24)

Using the power series approximation equation (19) together with equations (22) and (23), the contribution to the cosmic SFR density from haloes of mass M_h (the integrand of equation (24)) is given by

\begin{eqnarray*} \frac{{\rm d}\dot{\rho_*}}{{\rm d}M_h} & =& \epsilon_*(M_h)\left[ \frac{1}{M_h} \frac{{\rm d}M_h}{{\rm d}t} \right] M_h \frac{{\rm d}n(M_h,t)}{{\rm d}M_h} \nonumber\\ &=& \epsilon_*(M_h) \frac{46230.9 \hat{\rho }_0}{M_h S t^{7/3} t_m^{8/3}} \left(1 - 0.1590 \left(\frac{t}{t_\Lambda }\right)^2 - 0.0056 \left(\frac{t}{t_\Lambda }\right)^4 \right) \nonumber\\ && \times \exp \left[ -\frac{232382}{S \, t^{4/3} t_m^{8/3}} \left(1 + 0.3182 \left(\frac{t}{t_\Lambda }\right)^2 + 0.0028 \left(\frac{t}{t_\Lambda }\right)^4\right)\right]. \nonumber \end{eqnarray*}

(25)

The cosmological constant term enters through both the multiplier and exponential terms, with a balance that depends on the halo mass through S (see equation 17). While smaller haloes are more abundant than large objects, a smaller fraction of the inflowing material is converted into stars. As a result, the SFR density is dominated by the largest haloes in which star formation is able to proceed without generating efficient BH feedback. The smaller haloes only contribute significantly at very early times, when the abundance of larger objects is strongly suppressed by the exponential term. We see therefore that the level of suppression expected for ≈10¹² M_⊙ haloes is representative of most of the SFR in the Universe.

The predictions for the contributions of different halo masses are shown in Fig. 5, together with the total expected cosmic SFR density, for the two cosmologies that we consider in this paper. We will compare this approximation in Section 5 to the results from the eagle simulations.

Figure 5.

The predicted SFR history of the Universe, and the expected influence of the cosmological constant using the simple model developed in Section 2.5. Coloured lines show the contributions from dark matter haloes of different masses (per dex), using the star formation efficiency described by equation (25). The total SFR for the ΛCDM universe calculated numerically is shown in blue. An Einstein–de Sitter universe is shown in orange. The integrated SFR calculated using the approximation of equations (24) and (25) is shown with a dashed green line. The bottom panel shows the ratio at a given time. The predicted suppression of SFR due to Λ at the present time is |${\approx }19\hbox{ per cent}$|⁠. At t ≈ 30 Gys the predicted SFR density for the EdS model is double than ΛCDM, and ≈6 times higher at t = 50 Gyr. The approximation of equations (24) and (25) ceases to work for t ≳ 25 Gyr.

3 THE EAGLE SIMULATIONS

The simple analytic model provides a basis for interpreting the results, but it is highly simplified. We therefore compare the analytic model to numerical hydrodynamic simulations based on the eagle project. The eagle simulation suite³ (Crain et al. 2015; Schaye et al. 2015) consists of a large number of cosmological hydrodynamical simulations that include different resolutions, simulated volumes, and physical models. These simulations use advanced smoothed particle hydrodynamics (SPH) and state-of-the-art subgrid models to capture the unresolved physics. The simulation suite was run with a modified version of the gadget-|${\scriptstyle 3}$| SPH code (last described by Springel 2005) and includes a full treatment of gravity and hydrodynamics. The calibration strategy is described in detail by Crain et al. (2015), who also presented additional simulations to demonstrate the effect of parameter variations.

The halo and galaxy catalogues for more than 10⁵ simulated galaxies of the main eagle simulations with integrated quantities describing the galaxies, such as stellar mass, SFRs, metallicities, and luminosities, are available in the eagle data base ⁴ (McAlpine et al. 2016). A complete description of the code and physical parameters used can be found in Schaye et al. (2015).

The eagle reference simulations used cosmological parameters measured by Planck Collaboration XVI (2014). In this paper we introduce three main eagle simulations that use the same calibrated sub-grid parameters as the reference model, but change the cosmological model by setting the cosmological constant to zero, and/or removing feedback from BHs. The values of the cosmological parameters used for the simulations are listed in Table 1. The values of other relevant parameters adopted by all simulations featured in this study are listed in Table 2. Together, these parameters determine the dynamic range and resolution that can be achieved by the simulations.

Table 2.

Box-size, number of particles, initial baryonic, and dark matter particle mass, comoving and Plummer-equivalent gravitational softening, inclusion of AGN feedback, cosmological model, and Hubble parameter for the eagle simulations used in this paper. Values in bold show differences with respect to the Ref simulation. The three bottom small box models were used for convergence tests.

Identifier	L	N	mgas	mDM	εcom	εprop	AGN	Cosmology	h
	(cMpc)		(M⊙)	(M⊙)	(ckpc)	(pkpc)
ΛCDM (Ref)	50	2 × 7523	1.81 × 106	9.70 × 106	2.66	0.70	Yes	Planck 14	0.6777
ΛCDM (No AGN)	50	2 × 7523	1.81 × 106	9.70 × 106	2.66	0.70	No	Planck 14	0.6777
EdS	50	2 × 7523	1.81 × 106	9.70 × 106	2.66	0.70	Yes	EdS	0.3754
EdS (No AGN)	50	2 × 7523	1.81 × 106	9.70 × 106	2.66	0.70	No	EdS	0.3754
Λ = 0 L12 h0_3754	12.50	2 × 1883	1.81 × 106	9.70 × 106	2.66	0.70	Yes	EdS	0.3754
Λ = 0 L12 h0_6777	8.43	2 × 1883	1.81 × 106	9.70 × 106	|$\mathbf {1.79}$|	0.70	Yes	EdS	0.6777
Λ = 0 L12 h0_4716	10.73	2 × 1883	1.81 × 106	9.70 × 106	|$\mathbf {2.28}$|	0.70	Yes	EdS	0.4716

Fig. 6 shows the projected gas density for the ΛCDM and EdS cosmological models both at the present day and into the future. At t = 13.8 Gyr, the general appearance of both models is similar, but over the next 6.8 Gyr, the effect of Λ becomes more significant slowing down the growth of structure.

Figure 6.

The evolution of the projected gas density for each eagle model centred on the most massive halo at the present time (t = 13.8 Gyr). The length of each image is 43 (proper) Mpc on a side, to highlight the difference on cosmic expansion. Left: ΛCDM universe. Right: EdS universe. Top: Cosmic time t = 13.8 Gyr. Bottom: Cosmic time t = 20.7 Gyr. The colour coding represents the (proper) surface gas density projected along the line of sight. At t = 13.8 Gyr, the general appearance of both models is similar, as the phases of the initial fluctuations are the same. Over the next 6.8 Gyr, the effect of Λ becomes more significant, slowing down the growth of structure compared to the EdS model.

3.1 Subgrid models

Processes that are not resolved by the simulations are implemented as subgrid physical models; they depend solely on local interstellar medium (ISM) properties. A full description of these subgrid models can be found in Schaye et al. (2015). In summary:

• Radiative cooling and photoheating are implemented element by element as in Wiersma, Schaye & Smith (2009a), including the 11 elements found to be important, namely, H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe. Hydrogen reionization is implemented by switching on the full Haardt & Madau (2001) background at the proper time corresponding to redshift z = 11.5 in our ΛCDM Universe.

• Star formation is implemented stochastically following the pressure-dependent Kennicutt–Schmidt relation as in Schaye & Dalla Vecchia (2008). Above a metallicity-dependent density threshold |$n^{\hat{a}}{\acute{\rm L}} \mathring{\rm U}_{\rm H}(Z)$|⁠, which is designed to track the transition from a warm atomic to an unresolved, cold molecular gas phase (Schaye 2004), gas particles have a probability of forming stars determined by their pressure.

• Time-dependent stellar mass loss due to winds from massive stars and AGB stars, core collapse supernovae and type Ia supernovae, is tracked following Wiersma et al. (2009b).

• Stellar feedback is treated stochastically, using the thermal injection method described in Dalla Vecchia & Schaye (2012).

• Seed BHs of mass M = 1.48 × 10⁵ M_⊙, are placed in haloes with a mass greater than 1.48 × 10¹⁰ M_⊙ and tracked following the methodology of Springel, Di Matteo & Hernquist (2005) and Booth & Schaye (2009). Accretion on to BHs follows a modified version of the Bondi–Hoyle accretion rate which takes into account the circularization and subsequent viscous transport of infalling material, limited by the Eddington rate as described by Rosas-Guevara et al. (2015).⁵ Additionally, BHs can grow by merging with other BHs as described in Schaye et al. (2015) and Salcido et al. (2016).

• Feedback from AGN is implemented following the stochastic heating scheme described by Schaye et al. (2015). Similar to the supernova feedback, a fraction of the accreted gas on to the BH is released as thermal energy with a fixed heating temperature into the surrounding gas following Booth & Schaye (2009).

For the eagle simulations, the subgrid parameters were calibrated to reproduce three properties of galaxies at redshift z = 0: the galaxy stellar mass function (GSMF), the galaxy size–stellar mass relation, and the BH mass-stellar mass relation.⁶ The calibration strategy is described in detail by Crain et al. (2015), who explore the effect of parameter variations.

3.2 Halo and galaxy definition

Haloes were identified running the 'Friends-of-Friends' (FoF) halo finder on the dark matter distribution, with a linking length equal to 0.2 times the mean inter-particle spacing. Galaxies were identified as self-bound over-densities within the FoF group using the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). A ‘central’ galaxy is the substructure with the largest mass within a halo. All other substructures within a halo are ‘satellite’ galaxies.

Comparing haloes from simulations with different cosmologies is not a well-defined task, as halo masses are usually defined in terms of quantities that depend on the specific cosmological parameters. Typically, this is done by growing a sphere outwards from the potential minimum of the dominant dark matter sub-halo out to a radius where the mean interior density equals a fixed multiple of the critical or mean density of the Universe, causing an artificial ‘pseudo-evolution’ of dark matter haloes by changing the radius of the halo (Diemer, More & Kravtsov 2013). Star formation, however, is governed by the amount of gas that enters these haloes and reaches their central regions. Wetzel & Nagai (2015) show that the growth of dark matter haloes is subject to this ‘pseudo-evolution’, whereas the accretion of gas is not. Because gas is able to cool radiatively, it decouples from dark matter, tracking the accretion rate near a radius of |$R_{200{\bar{\rho }}}$|⁠, the radius within which the mean density is 200 times the mean density of the universe, |$\bar{\rho}$|⁠. As we try to connect the accretion of dark matter haloes to star formation, we define halo masses as the total mass within |$R_{200{\bar{\rho }}}$|⁠,

\begin{eqnarray*} M_{200_{\bar{\rho }}} = 200 \frac{4\pi }{3} R^3_{200{\bar{\rho }}} \bar{\rho }. \end{eqnarray*}

(26)

Additionally, as |$\bar{\rho } = \Omega_{m}(t) \rho_{c}(t)$| is given in comoving coordinates, the mean density of the universe remains constant in time for each cosmological model.

Following Schaye et al. (2015) and Furlong et al. (2015), galaxy stellar masses are defined as the stellar mass associated with the subhalo within a 3D 30 proper kilo parsec (pkpc) radius, centred on the minimum of the subhalo’s centre of gravitational potential. This definition is equivalent to the total subhalo mass for low-mass objects, but excludes diffuse mass around very large subhaloes, which would contribute to the intracluster light.

3.3 Continuing the simulations into the future

As Λ continues driving the accelerated expansion of the universe, the linear growth of density perturbations, D(t), is suppressed (see equation 15). Further insight can be obtained if we analyse the evolution of the potential perturbations given by the perturbed Poisson equation for an expanding space,

\begin{eqnarray*} \nabla ^2 \Phi = 4 \pi G \bar{\rho }_{\mbox{p}} a^2 D\delta_0, \end{eqnarray*}

(27)

where the Laplace operator is with respect to comoving coordinates, and the mean density |$\bar{\rho }_{\mbox{p}}$| is given in proper coordinates. As |$\bar{\rho }_{\mbox{p}}$| evolves ∝ a⁻³, it follows that ∇²Φ ∝ D/a. Using equations (10) and (15), we can see that for an EdS universe, both D and a are ∝ t^2/3 and the potentials are expected to stop evolving (they are frozen in). On the other hand, the suppression of growth of density perturbations due to a cosmological constant causes a decay in the potentials as the universe expands. As shown in Fig. 2, according to linear theory, these two scenarios have comparable growth factors at the present time (⁠|${\approx }10 \hbox{ per cent}$| difference; see equation 15), but the difference becomes increasingly important in the future. Furthermore, star formation is expected to eventually exhaust the finite reservoir of cold gas in galaxies, shutting off the production of stars in the universe forever (e.g. Fukugita, Hogan & Peebles 1998; Loeb, Batista & Sloan 2016).

In order to study the impact of Λ in galaxy formation beyond the present day, and hence explore the uniqueness of the present epoch, and in order to determine the total mass of stars ever produced by the universe, we allow the simulations to run into the future, i.e. t > t₀ (e.g. Barnes et al. 2005; Loeb et al. 2016). The subgrid models for star formation, stellar mass loss, stellar feedback, BH seeding, and feedback from AGN were kept as described in Section 3.1 as the simulations ran into the future. On the other hand, as there is no information about the UV and X-ray background radiation from quasars and galaxies into the future, for simplicity, we assumed that the background radiation freezes out, i.e. we kept its value at t = t₀ constant into the future. We consider this to be a good simplification as the UV background only affects star formation in very low mass haloes, and hence does not affect the cosmic SFR at late times (e.g. Schaye et al. 2010).

4 SIMULATIONS RE-SCALING

In this section we describe our simulation re-scaling strategy. At early epochs, the universe was matter dominated, and so we can neglect the contribution of Λ. Hence, any universe with non-zero matter density, i.e. ρ_{m, 0} ≠ 0, will be close to an EdS universe at early epochs. Therefore, we can assume identical initial conditions for all cosmological models of interest here.

The initial conditions for the reference ΛCDM model were created in three steps. First, a particle load, representing an unperturbed homogeneous periodic universe, was produced. Secondly, a realization of a Gaussian random density field with the appropriate linear power spectrum was created over the periodic volume. Thirdly, the displacements and velocities, consistent with the pure growing mode of gravitational instability, were calculated from the Gaussian realization and applied to the particle load producing the initial conditions. The initial density perturbation power spectrum is commonly assumed to be a power law, i.e. |$P_i(k) \propto k^{n_s}$|⁠. From the Planck results (Planck Collaboration XVI 2014), the spectral index n_s, has a value of n_s = 0.9611. A transfer function with the cosmological parameters shown in Table 1 was generated using CAMB (version Jan_12; Lewis, Challinor & Lasenby 2000). The linear matter power spectrum was generated by multiplying the initial power spectrum by the square of the dark matter transfer function evaluated at the present day t = t₀, i.e. P(k, t) = P_i(k)T²(k)D²(t).⁷

The eagle version of gadget uses an internal system of units that includes both comoving coordinates and the dimensionless Hubble parameter, h. For the alternative cosmological models, we have the freedom to choose the present time t₀ for each simulation, and we re-scale all the initial condition such that they are identical in physical 'h-free' units at an early time t₁ = 11.98 Myr. Table 3 shows the parameters that have been re-scaled in the initial conditions.

Table 3.

Parameters re-scaled in the initial conditions. Hat notation indicates parameters for our Universe.

Parameter	Units	Re-scaling factor
Box size	cMpc h−1	|$(\hat{h}^{-1} h) \times (\hat{a}_{1} a_{1}^{-1})$|
Particle masses	M⊙ h−1	|$(\hat{h}^{-1} h)$|
Particle coordinates	cMpc h−1	|$(\hat{h}^{-1} h) \times (\hat{a}_{1} a_{1}^{-1})$|
Particle velocities	cMpc s−1	|$(\hat{a}_{1} a_{1}^{-1}) ^{1/2}$|

The same tables of radiative cooling and photoheating rates as a function of density and temperature were used for all cosmological models. The corresponding redshifts for the cooling tables were re-scaled such that they correspond to the same cosmic time for each cosmology. That is, using equation (8), we find the scale factor, a, for which the alternative cosmology satisfies

\begin{eqnarray*} {t} (\hat{a}) = t (a). \end{eqnarray*}

(28)

The average baryonic density Ω_b has been re-scaled in such a way that the baryon fraction (f_b = Ω_b/Ω_m) is equal in both cosmologies, i.e.

\begin{eqnarray*} \frac{\hat{\Omega }_\mathrm{b}}{\hat{\Omega }_\mathrm{m}} = \frac{\Omega_\mathrm{b}}{\Omega_\mathrm{m}}. \end{eqnarray*}

(29)

Table 4 shows additional parameters that have been re-scaled to be equivalent in h-free physical units. Finally, hydrogen and He ii reionization were also re-scaled in such a way that redshifts correspond to the same cosmic time.

Table 4.

Additional parameters re-scaled in the simulations. Hat notation indicates parameters for our Universe.

Parameter	Units	Re-scaling factor
Comoving softening	ckpc h−1	|$(\hat{h}^{-1} h) \times (\hat{a}_{1} a_{1}^{-1})$|
Max softening	pkpc h−1	|$(\hat{h}^{-1} h)$|
Seed BH mass	M⊙ h−1	|$(\hat{h}^{-1} h)$|
Min MFOF for new BH	M⊙ h−1	|$(\hat{h}^{-1} h)$|

In order to demonstrate that this re-scaling strategy works correctly, Fig. 7 shows the global SFR density for the three small box EdS simulations used for convergence (see Table 2). They each represent the same physical scenario, but choose a different proper time to be 'today', t₀. This has the effect of altering the values of the Hubble parameter h and the redshift of the initial conditions, so that the simulations begin at proper time t(a₁) = 11.98 Myr in all models. Despite the small size of the simulation boxes (hence the noisy curves), the figure shows consistent SFRs as a function of cosmic time for the three models. Therefore, our re-scaling strategy allows us to simulate any cosmological model, regardless of the value of h.

Figure 7.

Global SFR density for three EdS models scaled by the ratio of the initial scale factors for each model. The initial conditions for each model have been re-scaled such that the time at which we start the simulations remains unchanged, i.e. t(a₁) = 11.98 Myr.

5 RESULTS: THE EVOLUTION OF STAR FORMATION

5.1 The past history of the cosmic SFR

Fig. 8 shows the global SFR density as a function of cosmic time for our simulation models. For comparison, observations from Cucciati et al. (2012) [FUV], Bouwens et al. (2012) [UV], Robertson et al. (2013) [UV], and Burgarella et al. (2013) [FUV + FIR] are shown as well. Solid lines in the figure show the evolution of the (comoving) cosmic SFR density for the reference ΛCDM eagle run (blue) and for an EdS universe (orange). Dashed lines show simulation models without feedback from AGN. Dotted lines show the prediction for the cosmic SFR density using equation (25). We focus first on the evolution of the models up to the present age of the universe, t = t₀ = 13.8 Gyr.

Figure 8.

Global SFR densities. The eagle reference and the EdS models are shown in solid blue and orange lines, respectively. Observational data from Cucciati et al. (2012) [FUV], Bouwens et al. (2012) [UV], Robertson et al. (2013) [UV], and Burgarella et al. (2013) [FUV + FIR] are shown as symbols. The eagle reference simulation (solid blue) reproduces the shape of the observed SFR density remarkably well, with a small offset of 0.2 dex at t ≳ 2 Gyr (for clarity, all eagle models have been shifted by 0.2 dex). The analytical model of equations (24) and (25) is shown with dotted lines for both cosmologies. Power-law fittings for the SFR density for t > 8 Gyr, as per equation (30), for the ΛCDM and EdS models are shown in pink and red lines, respectively. The horizontal time axis has been plotted in a logarithmic scale for t ≤ 8 Gyr changing to a linear scale for t > 8 Gyr. The black vertical dotted line shows the transition from logarithmic to linear scale.

In linear time, the SFR rises very rapidly and most of the plot is dominated by the slow decline (for an example, see Furlong et al. 2015). Hence, in order to emphasize the growth and decline of the SFR, and to reproduce the familiar shape of the star formation history (Madau & Dickinson 2014), the horizontal time axis has been plotted in a logarithmic scale for t ≤ 8 Gyr. In order to explore the SFR in detail at the present epoch and into the future, the horizontal time axis changes to a linear scale for t > 8 Gyr. The black vertical dotted line shows the transition from logarithmic to linear scale. For reference, the redshifts, |$\hat{z}$|⁠, for an observer at t₀ in the ΛCDM universe, are given along the top axis. As discussed in detail in Furlong et al. (2015), the reference simulation (solid blue line) reproduces the shape of the observed SFR density remarkably well, with a small offset of 0.2 dex at t ≳ 2 Gyr. While the simulations agree reasonably well with the observational data at redshifts above 3, we caution that these measurements are more uncertain.

Remarkably, the shape of the cosmic SFR history is very similar for both the ΛCDM and EdS models: the SFR density peaks ≈3.5 Gyr after the Big Bang and declines slowly thereafter. The similarity of the universes prior to the peak is expected, since the Λ term in the Friedman equation is sub-dominant in both cases. At later times, however, we might naively have expected the decline to be more pronounced in the ΛCDM cosmology, since the growing importance of the Λ term slows the growth of density perturbations.

From Fig. 2, the linear growth factors of the two cosmological models differ by |${\approx }10\hbox{ per cent}$| at the present time, and so we might have expected a similar difference in the (comoving) cosmic SFR density (⁠|${\approx } 15 \hbox{ per cent}$|⁠, read from Fig. 8). This naive expectation is not borne out because of the complexity of the baryonic physics. Because of stellar and AGN feedback, haloes have an ample reservoir of cooling gas that is able to power further star formation, regardless of the change in the cosmic halo growth rate.

Our simulation demonstrates that the existence of Λ does play a small role in determining the (comoving) cosmic SFR density. However, these differences are minor. In order to put the differences into context, we compare with a pair of simulations in which the BH feedback is absent. These runs are shown as dashed lines in the plot. We focus here on the behaviour before t = 13.8 yr. As can be seen, the absence of AGN feedback has a dramatic effect on the shape of the cosmic star formation density (Schaye et al. 2010; van de Voort et al. 2011). Interestingly, however, while the normalization of the SFR density is considerably higher, the time of the peak is similar. BH feedback is not solely responsible for the decline in star formation after t ≈ 3.5 Gyr. This hints that the existence of the peak results from the interaction of the slowing growth rates of haloes (after the peak) and the star formation time-scale (set by the ISM physics) which limits the rate at which the galaxy can respond to convert in-falling material into stars (before the peak).

The SFR history predicted by the simple model developed in Section 2.5 (dotted curves) is in remarkable agreement with the observational data and predicts the relative difference of the cosmological simulations, both at the present time and into the future. We want to emphasize that the model is not a parametric fit to the data, but rather an analytical model derived from a simple relation of star formation to halo mass accretion.

5.2 The future of the cosmic SFR (t > 13.8 Gyr)

In order to explore whether the relative SFR densities will diverge as the impact of Λ becomes more pronounced, we ran the simulations for both cosmological models into the future (i.e. beyond a cosmic time of t = 13.8 Gyr). As the simulations run into the future, the small differences seen at t = 13.8 Gyr become larger, reading an |${\approx }40\hbox{ per cent}$| difference at t = 1.5 × t₀ = 20.7 Gyr, which is in agreement with the predictions shown in Fig. 5.

The decline in the SFR density can be approximated by a power law for both the reference ΛCDM and EdS models (red and pink dashed lines),

\begin{eqnarray*} \dot{\rho }_* = a(ct)^{-k}, \end{eqnarray*}

(30)

with the parameters a, c, and k given in Table 5. We used the reduced chi-square statistic for goodness-of-fit testing. We use this fitting function to extrapolate the results from the simulations further into the future.

Table 5.

Power-law parameter fitting for the median SFR shown in Fig. 8.

Model	a	c	k
	(M⊙ yr−1 cMpc−3)	(Gyr−1)
ΛCDM	5.42	0.81	2.77
EdS	3.99	0.94	2.41

We note here a striking feature of the universes with no BH feedback: the SFR increases again in the future, for both the ΛCDM and EdS cosmologies. In Section 5.4 we discuss how this effect originates from massive galaxies (M_* > 10¹¹ M_⊙) that rejuvenate in the future. As it is not clear from the simulation whether the SFR will continue to rise, or when it would start declining, we have not fitted any functional form to the 'No AGN' models.

5.3 The stellar mass density

To study the build-up of stellar mass, we present the growth in stellar mass density, ρ_*, across cosmic time in Fig. 9. The colour coding is the same as in Fig. 8. The lower panel shows the ratio of the stellar density compared to the reference ΛCDM model. Furlong et al. (2015) show that the reference eagle simulation is in good agreement with the observed growth of stellar mass across cosmic time. In contrast to Furlong et al. (2015), where the stellar mass density was obtained from aperture measurements to facilitate comparison with observations, we calculate ρ_* by integrating the SFR density from Fig. 8,

\begin{eqnarray*} \rho_* = \int_0^t \dot{\rho }_* {\rm d}t^\prime , \end{eqnarray*}

(31)

in order to provide an estimate of the total mass of stars produced by the universe.⁸ For the models with AGN feedback, we have extrapolated the power-law fit described in Section 5.1 far into the future, up to 10 trillion years, and considered this as the 'Total stellar mass density' of the universe. As suggested by the analytic model in equation (25), in universes with feedback from star formation and AGN, the cosmic SFR density is expected to continue decreasing into the future. At late times, the SFR becomes orders of magnitude lower than that of the peak at 3.5 Gyr. Fig. 9 shows that the total stellar mass density is dominated by the contribution from the peak in star formation and reaches a plateau at t ≈ t₀. Hence, the formal uncertainties in the extrapolation into the future are unimportant for the predicted total stellar mass density. The right-hand axis of Fig. 9 represents the percentage of the total stellar density compared to ΛCDM. For the reference ΛCDM model, the universe has already produced most (⁠|${\approx }88\hbox{ per cent}$|⁠) of its eventual stellar mass by the present day, adding up very little stellar mass into the future. Although there is no Λ to slow down the formation of cosmic structure, the EdS model closely resembles a ΛCDM cosmology with only |${\approx }15\hbox{ per cent}$| more stellar mass produced.

Figure 9.

The stellar mass density as a function of time in the eagle simulation models. The colour coding is the same as in Fig. 8. The right-hand axis represents the percentage of the total stellar density compared to the ΛCDM model. Both the ΛCDM and EdS models have already produced most of the stars in the universe by the present day, building up very little stellar mass into the future. The models without AGN feedback quickly deviate from the reference model, producing almost twice the mass in stars by the end of the simulation (≈20.7 Gyr). The figure shows that the effect of dark energy on the overall star formation is negligible.

As discussed in Section 3.3, a universe with Λ has a very different expansion history compared to one without dark energy. This produces a different growth of density perturbations, in particular into the future (see Figs 1 and 2). Nevertheless, as seen in Fig. 9, since both cosmologies have already produced most of the stars in the universe by the present day, when the contribution of Λ is becoming increasingly important, the effect of dark energy on the overall star formation is negligible.

In contrast, the models without feedback from BHs (dashed curves) quickly deviate from the reference model, starting from t ∼ 1 Gyr, producing almost twice the mass in stars by the end of the simulations, 20.7 Gyr after the Big Bang.

5.4 Other galaxy population properties

In this section we will compare the galaxy population properties of the two simulation models at the present time and into the future. In particular, we compare the GSMF and the specific star formation rate (SSFR) of galaxies. To compare with observational data, in each of the figures discussed below, the left-hand panel shows properties at t = 12.5 Gyr, equivalent to redshift |$\hat{z}=0.1$| for an observer at the present time in a ΛCDM universe. The right-hand panel shows the same property but at t = 1.5 × t₀ = 20.7 Gyr. To guide the eye, each property for the reference ΛCDM model at t = 12.5 Gyr is plotted with a black line in each plot at t = 1.5 × t₀ = 20.7 Gyr. Finally, we have also included the ratios of each quantity to the reference ΛCDM model at the corresponding time at the bottom of each panel.

5.4.1 The galaxy stellar mass function

The effect of Λ on the GSMF can be seen in Fig. 10. The colour coding is the same as in Fig. 8. For comparison, observational data from Baldry et al. (2012) and Li & White (2009) are shown as well. As discussed in Schaye et al. (2015), the observed GSMF at redshift 0.1 was used to infer the free parameters of the subgrid physics used in the simulation. The reference eagle model reproduces the shape of the observed GSMF reasonably well, with a slight under-abundance of galaxies at its knee.

Figure 10.

The GSMF at t = 12.5 Gyr, equivalent to redshift |$\hat{z}=0.1$| for an observer at the present time in a ΛCDM universe (left), and t = 1.5 × t₀ (right) for the eagle simulation models. The colour coding is the same as in Fig. 8. Observational data from Baldry et al. (2012) and Li & White (2009) are shown as symbols. The reference and 'No AGN' ΛCDM models at t = 12.5 Gyr are plotted in the right-hand panel for reference (solid and dashed black lines, respectively). The effect of dark energy on the GSMF is negligible, with very little evolution into the future. As expected, the models without AGN feedback predict a higher number density of massive galaxies (M_* > 10¹¹ M_⊙). This effect becomes more significant into the future, going from ≈0.7 to ≈1.1 dex.

Fig. 10 shows that the effect of Λ on the GSMF is negligible, with very little evolution into the future. The models without AGN feedback predict a higher number density of massive galaxies (M_* > 10¹¹ M_⊙) compared to the models with AGN feedback. This effect becomes more significant into the future, going from 0.7 to 1.5 dex. The origin of this difference is explored in the next section.

5.4.2 Specific SFRs

Galaxies can be broadly classified into largely distinct star-forming and passive populations according to their SSFR,

\begin{eqnarray*} \mathrm{SSFR} = \frac{\dot{M}_*}{M_*}. \end{eqnarray*}

(32)

For star-forming galaxies, there is a well-defined star-forming sequence, with SSFR observed to be approximately constant as a function of stellar mass (e.g. Noeske et al. 2007; Karim et al. 2011). Fig. 11 shows the SSFR for star-forming galaxies in the simulations as a function of galaxy stellar mass at the present day and into the future. The colour coding is the same as in Fig. 8. The horizontal dotted lines correspond to the SSFR cut (10⁻² Gyr⁻¹) used to separate star forming from passive galaxies in our Universe. For comparison, observational data from Gilbank et al. (2010) are shown as well. Furlong et al. (2015) show that the SFFR in the reference simulations at the present day is similar to observations in the local universe, with an offset of 0.3 dex. This is possibly consistent with the systematic uncertainties in the calibration of the observation diagnostics. At low masses there is an increase in SSFR with stellar mass; however, this has been found to be a resolution-dependent effect. Hence, we have plotted the results with lighter coloured lines (similar to Schaye et al. 2015). The models without feedback from AGN have higher SSFR for M_* > 10¹⁰ M_⊙, whereas the effect of Λ on the SSFR of galaxies is negligible.

Figure 11.

The SSFR of star-forming galaxies at t = 12.5 Gyr, equivalent to redshift |$\hat{z}=0.1$| for an observer at the present time in a ΛCDM universe (left-hand side) and t = 1.5 × t₀ (right-hand side) for the eagle simulation models. The colour coding is the same as in Fig. 8. Observational data from Gilbank et al. (2010) are shown as symbols. The reference and 'No AGN' ΛCDM models at t = 12.5 Gyr are plotted in the right-hand panel for reference (solid and dashed black lines, respectively). The solid curves show the median relation for star-forming galaxies, defined as those with an SSFR above the limit specified by the horizontal dashed line (10⁻² Gyr⁻¹). The faint shaded regions enclose the 10th–90th percentiles for the ΛCDM and EdS Models. Lines are light coloured when the stellar mass falls below that corresponding to 100 baryonic particles, to indicate that resolution effects will be important. The figure shows that the effect of dark energy on the GSMF is negligible. The models without AGN feedback predict a higher SSFR for massive galaxies (M_* > 10¹⁰ M_⊙). The right-hand panel shows that the overall SSFR drops from t = 12.5 Gyr to t = 1.5 × t₀. For the 'No AGN' models, however, the SSFR increases for massive galaxies (M_* > 10¹¹ M_⊙).

Fig. 11 shows the galaxy population property that has the strongest evolution into the future. We find that over the next 6.8 Gyr the SSFR will drop by ≈0.4 dex.

Interestingly, the models without AGN feedback predict an increase of the SSFR of galaxies with M_* > 10¹¹ M_⊙ in the future. The figure shows that the increase in SFR shown in Fig. 8, and the higher number density of massive galaxies in Fig. 10, originate from massive galaxies that rejuvenate in the future. A plausible explanation for this phenomenon is that, according to simple radiative cooling models, in massive galaxies and clusters, a hot gaseous atmosphere should lose energy by the emission of radiation, and if there is no heating mechanism to compensate the cooling (e.g. AGN feedback), cooling flows should form (Fabian 1994; Peterson & Fabian 2006), triggering star formation. This result will be explored in more detail in a follow-up paper.

6 DISCUSSION AND CONCLUSIONS

In this paper, we explored the dependence of the star formation history of the universe on the existence of a cosmological constant and feedback from accreting BHs. We base our results on the eagle simulation code, which has been shown to compare favourably to observational data, and thus, to provide a good description of the formation of galaxies in our Universe. Feedback from supermassive BHs has been shown to be a key ingredient in achieving this match by suppressing star formation in massive haloes (e.g. Bower et al. 2006; Harrison 2017), while the accelerating expansion rate of the Universe suppresses the accretion rates of haloes at late times (e.g. Jenkins et al. 1998; Huterer et al. 2015). Our study allows us to assess the relative importance of these ingredients.

The universes that we consider are indistinguishable at early times. They share a common epoch of equality and recombination, and have equal amplitudes and spectrum of density fluctuations at early times. We take care to compare the evolution of models with equivalent starting points, and to demonstrate that the simulation code correctly scales the different values of the present-day expansion rate (Hubble parameter). When comparing the universes, it is important that we compare properties at a fixed cosmic time. Since the processes of stellar (and biological) evolution provide a common clock, independent of the large-scale cosmological expansion, these provide an astrophysically relevant comparison.

We have also developed an analytic model derived from a simple relation of star formation to halo mass accretion rate. Despite its simplicity, the model reproduces the overall shape of evolution of the cosmic SFR density. The model and the simulations allow us to explore the effect of a cosmological constant term on the cosmic SFR density.

Our main conclusions are as follows:

• We find that the existence of the cosmological constant has little impact on the star formation history of the Universe. The SFR is suppressed by |${\approx }15\hbox{ per cent}$| at the present time, and we find that the properties of galaxies are almost indistinguishable in the two universes.

• To explore whether this is due to the relatively recent dominance of the dark energy density in our Universe, we continued the simulations 6.8 Gyr into the future. Even after this time, the comoving SFR densities differ only by |${\approx }40\hbox{ per cent}$|⁠. Clearly, the cosmological constant has only a marginal effect on the stellar content of the Universe.

• Using the analytic model, we can recognize that the existence of the peak in the SFR density results from the interaction of the star formation efficiency (set by the ISM physics) which limits the rate at which the galaxy can respond to convert in-falling material into stars, the relative abundance of efficiently star-forming haloes (i.e. of masses ≈10¹² M_⊙), and only at late times, the slowing growth rates of haloes due to the cosmological constant.

• By extrapolating fits to the evolution of the comoving SFR density into the future, we show that, in our Universe, more than |${\approx }88\hbox{ per cent}$| of the stars that will ever be produced, have already been formed by the present cosmic time. In the absence of dark energy, only |${\approx }15\hbox{ per cent}$| more stellar mass would have been formed in the same time. The difference is small, bringing into question whether the ‘coincidence problem’ (the comparable energy densities of matter and dark energy) can be explained by an anthropic argument: the existence of dark energy (at the observed value) has negligible impact on the existence of observers or the ability of humanity to observe the cosmos. In Barnes et al. (2018), we explore this argument in more detail by considering a wider range of Λ values, and determining the likelihood distribution of possible Λ values conditioning the existence of observers.

• In comparison, the existence of BHs has a major impact on the Universe. In the absence of AGN feedback, the comoving SFR density is enhanced by a factor of 2.5 at the present day.

• Even in a universe without BHs or dark energy, we find that the comoving SFR density peaks at 3.5 Gyr (z ≈ 2 according to a present-day observer in our Universe). The decline in star formation is however slower at more recent times.

• For hypothetical universes without feedback from accreting BHs, there is a comeback of SFR, which increases again in the future. This effect originates from massive galaxies (M_* > 10¹¹ M_⊙) that rejuvenate as there is no heating mechanism to compensate the cooling, in turn, triggering star formation.

Footnotes

1

Of course, an observer in a Λ = 0 universe would measure a different angular power spectrum in the cosmic microwave background after 13.8 Gyr, because of the very different expansion history of the Universe at later times.

2

|$\hat{z}=127$| was the reference simulation’s starting redshift.

3

http://www.eaglesim.org

4

http://www.eaglesim.org/database.php

5

The eagle simulation does not include a boost factor the accretion rate of BHs to account for an unresolved clumping factor.

6

BH feedback efficiency left unchanged from Booth & Schaye (2009).

7

The CAMB input parameter file and the linear power spectrum are available at http://eagle.strw.leidenuniv.nl/

8

As we are interested in total mass of stars produced by the universe, equation (31) ignores stellar mass loss.

ACKNOWLEDGEMENTS

We are grateful to all members of the Virgo Consortium and the eagle collaboration, who have contributed to the development of the codes and simulations used here, as well as to the people who helped with the analysis. We thank Alejandro Benitez-Llambay for the visualization software py-sphviewer (Benitez-Llambay 2015).

This work was supported by the Science and Technology Facilities Council (grant number ST/P000541/1); European Research Council (grant number GA 267291 'Cosmiway'), and the Netherlands Organization for Scientific Research (NWO, VICS grant 639.043.409). RAC is a Royal Society University Research Fellow.

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (http://www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. We acknowledge PRACE for awarding us access to the Curie machine based in France at TGCC, CEA, Bruyères-le-Châtel.

Jaime Salcido gratefully acknowledges the financial support from the Mexican Council for Science and Technology (CONACyT), fellow no. 218259.