© The Authors 2024

1. Introduction

Baryonic matter occurs as interstellar gas and stellar material in galaxies. The relationship between interstellar gas, stars, and metals is a central issue for understanding the cosmic baryon cycle in galaxies (Péroux & Howk 2020; Tortora et al. 2022; Mercado et al. 2021). Stars form from the collapse of interstellar gas clouds consisting mainly of hydrogen and helium evolve with time and therefore undergo nuclear fusion reactions and create metals, that is, elements heavier than helium (Carroll & Ostlie 2017; Clayton 1983; Phillips 1999; Lang 2013). It is well established that external pressure on gas clouds reduces the critical cloud mass to facilitate the collapse of the Bonnor-Ebert mass (Ebert 1955; Bonnor 1956), which is considerably smaller than the Jeans mass for the spontaneous collapse of gas clouds. Therefore, it is interesting to investigate both the spontaneous star formation of gas clouds and the triggered star formation due to the interaction of the gas with already formed stars. During the stellar evolution, part of the metal-enriched stellar matter is continuously fed back into the interstellar medium (ISM) by stellar winds and outflows (Dave et al. 2011; Tortora et al. 2022; Sadoun et al. 2016). Moreover, at the termination of stellar fusion reactions the powerful novae and supernovae resulting from the fatal star’s collapse eject further stellar material into the ISM, whereas the central cores of the original stars (depending on the stellar mass) end up as locked-in matter in the form of white dwarfs, neutron stars, and black holes (Portinari et al. 1998; Shapiro & Teukolsky 1983; Lang 1999; Frank et al. 2002; Woosley & Janka 2005; Hansen et al. 2004). While the locked-in matter can no longer feed the gaseous component, the supersonic nova and supernova outflows can trigger efficient star formation in the ambient ISM by providing strong enough external pressure (Low & Klessen 2004; Elmegreen 1997; Hennebelle & Chabrier 2011; Tielens 2005; Naoz & Perna 2014; Cousin et al. 2015; Maoz & Graur 2017).

The purpose of the present study is to describe the linear and nonlinear temporal evolution of the interstellar gas and stellar matter in galaxies by a compartmental model. The nonlinearity stems from the triggered star formation process. Compartmental models are very successful in analyzing and/or predicting the temporal development of infectious diseases and epidemics since their invention about a hundred years ago by Kermack & McKendrick (1927) and the subsequent refinement by Kendall (1956). Here, individuals from a considered population with many N persons are assigned to different compartments. In the simple susceptible-infectious-recovered/removed (SIR) model these compartments are S (susceptible), I (infectious), and R (recovered/removed), respectively. Time-dependent infection and recovery rates then regulate the transitions between the compartments S → I and I → R, respectively (for reviews, see Hethcote 1989, 2000 and Estrada 2020).

In the present manuscript, we restricted the analysis to the case of spontaneous star formation, whereas the additional influence of triggered star formation will be the subject in future papers of this series. The organization of this manuscript is as follows. In Sect. 2, we introduced the general compartmental description of the cosmic baryonic matter by formulating the relevant dynamical equations including spontaneous and triggered star formation. In Sect. 3, the dynamical equations were solved analytically for the case of negligible spontaneous star formation. The analytical solutions were applied to the present-day (redshift z = 0) gaseous fraction of the Universe and the observed cosmic star formation history in Sect. 4. The comparison of the theoretical results to the observations is the subject of Sect. 5 (in the special case of neglected stellar feedback) and Sect. 6 (including stellar feedback). The summary and conclusion in Sect. 7 complete the manuscript.

2. Gas-star-locked-in-matter (GSL) compartmental model for the baryonic matter cycle

We considered the total system of baryonic matter either in stars or in interstellar gas and introduce the compartments G (gas), S (stars), and L (locked-in matter in white dwarfs, neutron stars, and black holes). G(t), S(t), and L(t) denote the relative fractions of baryonic matter in the three compartments, respectively, as a function of time t, so that the sum constraint

$$\begin{array}{c}\hfill G(t)+S(t)+L(t)=1\end{array}$$(1)

holds at all times t after the beginning of the baryonic evolution at time t = t₀. Equation (1) adopts a “closed-box” system, where the total amount of baryonic matter is conserved and contained in at least one of the three fractions. This assumption is justified if we apply Eq. (1) to the whole Universe. However, if one applies Eq. (1) to individual galaxies or clusters of galaxies, one has to account for losses or gains of matter in the particular system. This is caused, for example, by the ejection of gas out of the galaxy by supernovae and black-hole feedback and the accretion of fresh gas from its host environment atmosphere, respectively (see the review by Tumlinson et al. 2017). For an application to such individual systems, the dynamical models considered in Pandya et al. (2023) are more appropriate. Nevertheless, we think that our investigation by using coarse-grained averages of gas fractions from a system of many individual astrophysical systems and enforcing the necessary conservation laws is fruitful and worth pursuing. It can provide new insights into the fundamental physical processes controlling the cosmological baryonic matter cycle when comparing with observations similarly averaged over a large enough number of individual astrophysical objects (see Sect. 4 below).

The temporal evolution of the three fractions G(t), S(t), and L(t) is controlled by the spontaneous star formation rate (SFR) β(t)G(t); the triggered star formation from the interaction of gas and stars with the assumed rate a(t)G(t)S(t); the continuous feedback rate b(t)S(t) of stellar to gaseous matter; and the formation rate c(t)S(t) of white dwarfs, neutron stars, and black holes from stellar evolution. The nonzero feedback rate is essential for the existence of gaseous heavier metals beyond helium and lithium from processed stellar baryons in the Universe. Throughout this paper, the term “stellar feedback” refers to gas returned to the ISM from luminous stars; this is different from the standard usage of this term in the literature, in which it refers to gas ejection out of the galaxy and/or the prevention of its accretion from outside. Similarly, our term “stellar evolution”, referring to the averaged formation rate of locked-in matter from the evolution of luminous stars, is a gross simplification compared to galaxy evolution models based on standard stellar evolution calculations for the production rate of white dwarfs, neutron stars, and black holes for given stellar initial mass functions.

The dynamical equations (Fig. 1) for the three fractions then read

$$\begin{array}{cc}\hfill \frac{\mathrm{d}G}{\mathrm{d}t}& =-a(t)G(t)S(t)-\beta (t)G(t)+b(t)S(t),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}S}{\mathrm{d}t}& =a(t)G(t)S(t)+\beta (t)G(t)-b(t)S(t)-c(t)S(t),\hfill \end{array}$$(3)

Fig. 1. Schematic diagram for nonlinear dynamical exchange (Eqs. (2)–(4)) between gas (G), star (S), and locked-in (L) baryonic matter fractions, regulated by the four potentially time-dependent rates a(t), b(t), c(t), and β(t).

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}L}{\mathrm{d}t}=c(t)S(t),\end{array}$$(4)

respectively. The spontaneous SFR starts in the truly metal-free primordial gas at the recombination era at z = 1100, during which charged electrons and protons became bound to form electrically neutral hydrogen atoms, so that the gas could decouple from the global expansion and begin to cool efficiently and contract (Klessen & Glover 2023). As initial condition, we adopted

$$\begin{array}{c}\hfill G(t_0)=1,S(t_0)=0,L(t_0)=0,\end{array}$$(5)

as the whole baryonic gaseous matter was subject to spontaneous star formation. The fact that most star formation only occurs at the centers of overdensities in the large-scale dark matter distribution called dark matter halos can be accounted for by lower values of the coarse-grained averaged rates of spontaneous (β(t)) and triggered (a(t)) star formation.

The initial time t₀ corresponds to the redshift z = 1100. The first stars are believed to have formed at z ∼ 30, since prior to z ∼ 30 the gas would not have cooled and grown overdense enough to collapse into stars (Klessen & Glover 2023). Our later results, however, indicated a rather minor influence of the exact starting time of the first stars because we compared with observations at redshifts below z = 8 and hence much smaller than 30 or 1100.

The introduction of gaseous and stellar fractions as the basic dynamical variables has the advantage that the results can be applied to a variety of astrophysical systems (Tinsley 1980) including individual galaxies, clusters of galaxies, and the whole Universe. In each particular case, one has to multiply the derived fractions as a function of time with the respective total mass of baryonic matter.

On the other hand, our approach is simplified as we only considered space-averaged rates for star formation, stellar feedback, and stellar evolution, without accounting for details such as the initial stellar mass function. Of course, these averaged rates can be varied for applications to different astrophysical objects. Our approach is different but bears some similarities with the earlier investigated bathtub (Somerville & Dave 2015) or gas-regulator (Peng & Maiolino 2014) models for galaxy evolution considering the competition among the flow of gas into galaxies, conversion to stars by in situ star formation, and ejection out of galaxies by stellar feedback. This includes individual galaxies, clusters of galaxies, and the whole Universe.

We were particularly interested in the formation rate $\mathring{J}(t)$ of new stars as a function of time:

$$\begin{array}{c}\hfill \mathring{J}(t)=\frac{\mathrm{d}J}{\mathrm{d}t}=a(t)G(t)S(t)+\beta (t)G(t).\end{array}$$(6)

Obviously, the first term represents the rate from triggered star formation, whereas the second term is the rate from spontaneous star formation.

As mentioned above, we restricted our analysis here to the case of negligible triggered star formation. In this case, exact analytical solutions of the dynamical equations of the GSL model (1)–(5) could be derived. The derived exact analytical solutions hold for stationary rates as well as for the case of the same time dependence of all rates.

3. Negligible triggered star formation

Here, for a(t) = 0, the GSL model Eqs. (2)–(4) simplify to the following linear equations:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}G}{\mathrm{d}t}& =-\beta (t)G(t)+b(t)S(t),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}S}{\mathrm{d}t}& =\beta (t)G(t)-b(t)S(t)-c(t)S(t),\hfill \end{array}$$(8)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}L}{\mathrm{d}t}=c(t)S(t).\end{array}$$(9)

Introducing the dimensionless time variable

$$\begin{array}{c}\hfill \stackrel{\sim }{t}={\displaystyle {\int }_{t_0}^t}\mathrm{d}\xi \beta (\xi ),\end{array}$$(10)

allowed us to eliminate β(t) from the equations and to thus solve them for arbitrary time-dependent rates β(t). Equations (7)–(9) become

$$\begin{array}{cc}\hfill \frac{\mathrm{d}G}{\mathrm{d}\stackrel{\sim }{t}}& =-G(\stackrel{\sim }{t})+r_1(\stackrel{\sim }{t})S(\stackrel{\sim }{t}),\hfill \end{array}$$(11a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}S}{\mathrm{d}\stackrel{\sim }{t}}& =G(\stackrel{\sim }{t})-[r_1(\stackrel{\sim }{t})+r_2(\stackrel{\sim }{t})]S(\stackrel{\sim }{t}),\hfill \end{array}$$(11b)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}L}{\mathrm{d}\stackrel{\sim }{t}}& =r_2(\stackrel{\sim }{t})S(\stackrel{\sim }{t}),\hfill \end{array}$$(11c)

$$\begin{array}{cc}\hfill 1& =G(\stackrel{\sim }{t})+S(\stackrel{\sim }{t})+L(\stackrel{\sim }{t}),\hfill \end{array}$$(11d)

with the dimensionless ratios

$$\begin{array}{c}\hfill r_1(\stackrel{\sim }{t}(t))=\frac{b(t)}{\beta (t)},r_2(\stackrel{\sim }{t}(t))=\frac{c(t)}{\beta (t)}.\end{array}$$(12)

3.1. Stationary ratios

In the following, we adopted stationary ratios r₁= const. and r₂= const., so our solution allows for arbitrary but given time-dependent spontaneous rates β(t).

The first equation, (11a) readily provides

$$\begin{array}{c}\hfill S(\stackrel{\sim }{t})=\frac{1}{r_1}[\frac{\mathrm{d}G(\stackrel{\sim }{t})}{\mathrm{d}\stackrel{\sim }{t}}+G(\stackrel{\sim }{t})],\end{array}$$(13)

which upon insertion into Eq. (11b) yields

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}G(\stackrel{\sim }{t})}{\mathrm{d}\stackrel{\sim }{t}{}^2}+(1+r_1+r_2)\frac{\mathrm{d}G(\stackrel{\sim }{t})}{\mathrm{d}\stackrel{\sim }{t}}+r_{2}G(\stackrel{\sim }{t})=0.\end{array}$$(14)

Equation (14) is solved by

$$\begin{array}{c}\hfill G(\stackrel{\sim }{t})=e^{-\mu \stackrel{\sim }{t}}[c_{-}{e}^{\nu \stackrel{\sim }{t}}+c_{+}{e}^{-\nu \stackrel{\sim }{t}}],\end{array}$$(15)

with

$$\begin{array}{c}\hfill \mu =\frac{1+r_1+r_2}{2},\nu =\sqrt{{\mu }^2-r_2},\end{array}$$(16)

and with the two integration constants c_±, implying the following for Eq. (13)

$$\begin{array}{c}\hfill S(\stackrel{\sim }{t})=\frac{e^{-\mu \stackrel{\sim }{t}}}{r_1}[c_{-}(1-\mu +\nu )e^{\nu \stackrel{\sim }{t}}+c_{+}(1-\mu -\nu )e^{-\nu \stackrel{\sim }{t}}].\end{array}$$(17)

The initial conditions (5), that is $G(\stackrel{\sim }{t}=0)=1$ and $S(\stackrel{\sim }{t}=0)=0,$ provide

$$\begin{array}{c}\hfill c_{\mp }=\frac{\nu \mp (1-\mu )}{2\nu }\end{array}$$(18)

so,

$$\begin{array}{c}\hfill S(\stackrel{\sim }{t})=\frac{1}{\nu }{e}^{-\mu \stackrel{\sim }{t}}sinh\nu t.\end{array}$$(19)

The dimensionless formation rate from spontaneous star formation only is given by

$$\begin{array}{c}\hfill j(\stackrel{\sim }{t})=\frac{\mathrm{d}J}{\mathrm{d}\stackrel{\sim }{t}}=G(\stackrel{\sim }{t})=e^{-\mu \stackrel{\sim }{t}}[cosh(\nu \stackrel{\sim }{t})-\frac{1-\mu }{\nu }sinh(\nu \stackrel{\sim }{t})]\end{array}$$(20)

and is related to Eq. (6) in this case as $\mathring{J}(t)=\beta (t)j(\stackrel{\sim }{t})$. Consequently,

$$\begin{array}{c}\hfill L(\stackrel{\sim }{t})=1-G(\stackrel{\sim }{t})-S(\stackrel{\sim }{t})=1-e^{-\mu \stackrel{\sim }{t}}[cosh(\nu \stackrel{\sim }{t})+\frac{\mu }{\nu }sinh\left(\nu \stackrel{\sim }{t}\right)].\end{array}$$(21)

For infinitely large times $\stackrel{\sim }{t}\to \infty$ we note that $G_{\infty }=G(\stackrel{\sim }{t}=\infty )=0$ and $S_{\infty }=S(\stackrel{\sim }{t}=\infty )=0$, whereas L_∞ = 1. While $G(\stackrel{\sim }{t})$ and $L(\stackrel{\sim }{t})$ monotonically decrease and increase, respectively, the stellar fraction $S(\stackrel{\sim }{t})$ first increases from zero to the maximum value,

$$\begin{array}{c}\hfill S_{\max }=\frac{{(\mu -\nu )}^{\frac{\mu -\nu }{2\nu }}}{{(\mu +\nu )}^{\frac{\mu +\nu }{2\nu }}}\end{array}$$(22)

at

$$\begin{array}{c}\hfill {\stackrel{\sim }{t}}_{\mathrm{S},\max }=\frac{1}{2\nu }ln\frac{\mu +\nu }{\mu -\nu },\end{array}$$(23)

and the S approaches zero over the course of time: S_∞ = 0.

The new star formation rate (6) for the stationary SFR coefficient β₀ in this case is

$$\begin{array}{cc}\hfill \mathring{J}(t)& ={\beta }_{0}G(\stackrel{\sim }{t}={\beta }_0(t-t_0))\hfill \\ \hfill & ={\beta }_0{e}^{-\mu {\beta }_0(t-t_0)}[cosh\left(\nu {\beta }_0(t-t_0)\right)-\frac{1-\mu }{\nu }sinh(\nu {\beta }_0(t-t_0))].\hfill \end{array}$$(24)

3.2. Negligible stellar feedback

In the special case of negligible stellar feedback b(t) = 0 = r₁, one obtains

$$\begin{array}{cc}\hfill j(\stackrel{\sim }{t},r_1=0)& =G(\stackrel{\sim }{t},r_1=0)\hfill \\ \hfill & =e^{-\frac{1+r_2}{2}\stackrel{\sim }{t}}[cosh\frac{1-r_2}{2}\stackrel{\sim }{t}-sinh\frac{1-r_2}{2}\stackrel{\sim }{t}]=e^{-\stackrel{\sim }{t}},\hfill \end{array}$$(25)

corresponding to $\mathring{J}(t,b=0)={\beta }_0{e}^{-{\beta }_0(t-t_0)}$. Likewise,

$$\begin{array}{c}\hfill S(\stackrel{\sim }{t},r_1=0)=\frac{2}{1-r_2}{e}^{-\frac{1+r_2}{2}\stackrel{\sim }{t}}sinh\left(\frac{1-r_2}{2}\stackrel{\sim }{t}\right)=\frac{e^{-r_2\stackrel{\sim }{t}}-e^{-\stackrel{\sim }{t}}}{1-r_2}\end{array}$$(26)

and

$$\begin{array}{c}\hfill L(\stackrel{\sim }{t},r_1=0)=1-\frac{e^{-r_2\stackrel{\sim }{t}}-r_2{e}^{-\stackrel{\sim }{t}}}{1-r_2}·\end{array}$$(27)

While $j(\stackrel{\sim }{t})=G(\stackrel{\sim }{t})$ and $L(\stackrel{\sim }{t})$ are monotonically decreasing or increasing, respectively, in this case one finds that $S(\stackrel{\sim }{t})$ rises first to the maximum value S_max = r₂^{r₂/(1 − r₂)} at ${\stackrel{\sim }{t}}_{\max }=(ln{r}_2)/(r_2-1)$ and then decreases to the final value of S_∞ = 0.

In the special case of r₂ = 1, Eqs. (25)–(27) simplify to

$$\begin{array}{cc}\hfill j(\stackrel{\sim }{t},r_1=0,r_2=1)& =G(\stackrel{\sim }{t},r_1=0)=e^{-\stackrel{\sim }{t}},\hfill \\ \hfill S(\stackrel{\sim }{t},r_1=0,r_2=1)& =\stackrel{\sim }{t}{e}^{-\stackrel{\sim }{t}},\hfill \\ \hfill L(\stackrel{\sim }{t},r_1=0,r_2=1)& =1-(1+\stackrel{\sim }{t})e^{-\stackrel{\sim }{t}}.\hfill \end{array}$$(28)

Here, $S(\stackrel{\sim }{t})$ first rises to the maximum value of S_max = e⁻¹ at ${\stackrel{\sim }{t}}_{\max }=1$ and then decreases to the final value, S_∞ = 0.

3.3. Negligible stellar evolution

In the special case of negligible stellar evolution, c(t) = 0 = r₂, one obtains

$$\begin{array}{c}\hfill j(\stackrel{\sim }{t},c=0)=G(\stackrel{\sim }{t},c=0)=\frac{1}{1+r_1}[r_1+e^{-(1+r_1)\stackrel{\sim }{t}}],\end{array}$$(29)

corresponding to

$$\begin{array}{c}\hfill \mathring{J}(t,c=0)=\frac{{\beta }_0}{b_0+{\beta }_0}[b_0+{\beta }_0{e}^{-(b_0+{\beta }_0)(t-t_0)}],\end{array}$$(30)

and $L(\stackrel{\sim }{t},c=0)=0$. In Fig. 2, we compared the derived exact analytic solutions with the numerical solutions of Eqs. (7)–(9) for several illustrative examples. For the numerical code, we used a single-step solver based on a second-order modified Rosenbrock formula, implemented by Shampine & Reichelt (1997) as ode23s in Matlab^TM. The excellent agreement proves the validity of our derivations.

Fig. 2. (a) Gas (G), (b) stars (S), and (c) locked-in matter (L) compartments as well as (d) the formation rate $\mathring{J}(t)$ of new stars in the course of time t, at zero star formation from the interaction between gas and stars, a(t) = 0. We show five choices for b0 and c0, while β0 = 1 and η = 10−5. The numerical solution of Eqs. (7)–(9) coincides exactly with the analytical solution provided by Eqs. (19)–(21), and (24).

4. Cosmic star formation history

As an astrophysical application, we considered the present-day (redshift z = 0) gas/star matter ratio and the cosmic star formation history (SFH) of the Universe, which is proportional to the star formation rate (6) as a function of real time. The present-day gas/star matter ratio G(z = 0)/[S(z = 0)+L(z = 0)] = G(z = 0)/[1 − G(z = 0)] in our Milky Way is about 0.1, and a little higher in galaxy clusters. If we take this ratio as characteristic for the whole Universe and allow for a possible contribution of hidden baryons in the intergalactic gas, a value of

$$\begin{array}{c}\hfill G(z=0)\simeq 0.10\pm 0.05\end{array}$$(31)

seems reasonable (Bland-Hawthorn & Gerhard 2016).

With respect to the cosmic SFH, far-ultraviolet and -infrared measurements have indicated (Madau & Dickinson 2014) that the best-fit SFR density per time t is

$$\begin{array}{cc}\hfill \psi (z)& =\frac{0.015{(1+z)}^{2.7}}{1+{\left[\frac{(1+z)}{2.9}\right]}^{5.6}}{\mathrm{M}}_{\odot }{\text{yr}}^{-1}{\text{Mpc}}^{-3}\hfill \\ \hfill & =\frac{3.22·{10}^{-47}{(1+z)}^{2.7}}{1+{\left[\frac{(1+z)}{2.9}\right]}^{5.6}}\text{kg}{\text{s}}^{-1}{\text{m}}^{-3},\hfill \end{array}$$(32)

where the time t(z) has been expressed as a function of redshift z using a flat ΛCDM Friedmann cosmology (Ryden 2017; Weinberg 2008; Carroll et al. 1992; Peebles 1993), and M_⊙ = 1.989 ⋅ 10³⁰ kg in SI-units¹. The data used for the best-fit (32) are given in Table 1 of Madau & Dickinson (2014). The SFR density (32) attains its maximum value ψ_max = 0.133 M_⊙ yr⁻¹ Mpc⁻³ at redshift z = 1.8632, which is practically constant at small redshifts well below unity, and it decreases at large redshifts ψ(z ≫ 2.9)≃1.25 ⋅ 10⁻⁴⁴(1 + z)^−2.9 kg s⁻¹ m⁻³.

Also available in Table 2 and Fig. 11 of Madau & Dickinson (2014) is the integrated stellar density ρ_*(z), which according to Madau & Dickinson (2014) fulfills dρ_*(z)/dt = (1 − R)ψ(z) and is therefore more explicitly given by

$$\begin{array}{c}\hfill {\rho }_{\ast }(z)=(1-R){\displaystyle {\int }_z^{\infty }}\mathrm{d}z\prime \psi (z\prime )\left|\frac{\mathrm{d}t\prime }{\mathrm{d}z\prime }\right|,\end{array}$$(33)

where according to Madau & Dickinson (2014)R is the return factor or the mass fraction of newly born stars; this is put back into the interstellar gas that has been treated as a free fit parameter. As we argue below, within the compartmental model considered here, 1 − R = 1 − L(t(z)) had to be replaced by the time- or redshift-dependent term 1 − L(t) with the fraction L(t) of locked-in stellar matter.

It is worth mentioning earlier works (such as Somerville et al. 2015 and references therein) that have attempted to reproduce the cosmic star formation history and stellar-mass-density evolution jointly with semi-analytic models (see their Fig. 17). The approach followed here is different with respect to the simplicity of the ordinary differential Eqs. (1)–(5), which allowed us to infer analytic constraints (see Sects. 5.1 and 5.2).

4.1. Timescales

We first discuss the relation of our dynamical evolution timescale $t=\stackrel{\sim }{t}/{\beta }_0$ for stationary spontaneous SFRs with the look-back time and redshift z in a flat ΛCDM Friedmann cosmology. In a flat ΛCDM Friedmann cosmology with Ω_m = Ω₀ = 0.3 and H₀ = 70 h₇₀ km s⁻¹ Mpc⁻¹ = 2.268 ⋅ 10⁻¹⁸ h₇₀ s⁻¹, the look-back time T_L is related to the redshift according to Peebles (1993), Hogg (2000), Boylan-Kolchin & Weisz (2021)

$$\begin{array}{cc}\hfill T_L(z)& =H_0^{-1}{\displaystyle {\int }_1^{1+z}}\frac{\mathrm{d}x}{x\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{x}^3}}\hfill \\ \hfill & =\frac{2}{3H_0}{\displaystyle {\int }_1^{{(1+z)}^{3/2}}}\frac{\mathrm{d}u}{u\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{u}^2}}\hfill \\ \hfill & =\frac{2}{3\sqrt{1-{\mathrm{\Omega }}_0}{H}_0}ln\frac{(1+\sqrt{1-{\mathrm{\Omega }}_0}){(1+z)}^{\frac{3}{2}}}{\sqrt{1-{\mathrm{\Omega }}_0}+\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{(1+z)}^3}}·\hfill \\ \hfill \end{array}$$(34)

The look-back time T_L from now (redshift z = 0) back to redshift z refers to the time that would have elapsed on the clock of an observer moving with the Hubble flow between these redshifts. The look-back time (34) corresponds to the age of the universe T₀ given by

$$\begin{array}{cc}\hfill T_0& =H_0^{-1}{\displaystyle {\int }_1^{\infty }}\frac{\mathrm{d}x}{x\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{x}^3}}\hfill \\ \hfill & =\frac{2}{3\sqrt{1-{\mathrm{\Omega }}_0}{H}_0}ln\frac{1+\sqrt{1-{\mathrm{\Omega }}_0}}{\sqrt{{\mathrm{\Omega }}_0}}·\hfill \end{array}$$(35)

The relation between the real time t, the look-back time T_L(z), and the age of the Universe T₀ is given by

$$\begin{array}{cc}\hfill t& =T_0-T_L(z)=H_0^{-1}{\displaystyle {\int }_{1+z}^{\infty }}\frac{\mathrm{d}x}{x\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{x}^3}}\hfill \\ \hfill & =\frac{2}{3\sqrt{1-{\mathrm{\Omega }}_0}{H}_0}\hfill \\ \hfill & \times ln\left[\frac{(1+\sqrt{1-{\mathrm{\Omega }}_0})(\sqrt{1-{\mathrm{\Omega }}_0}+\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{(1+z)}^3})}{\sqrt{{\mathrm{\Omega }}_0}(1+\sqrt{1-{\mathrm{\Omega }}_0}){(1+z)}^{\frac{3}{2}}}\right],\hfill \end{array}$$(36)

implying

$$\begin{array}{c}\hfill t_0=T_0-T_L(z=1100).\end{array}$$(37)

Consequently,

$$\begin{array}{c}\hfill t-t_0=T_L(1100)-T_L(z);\end{array}$$(38)

so, for a constant, spontaneous SFR coefficient β(t) = β₀, we obtained

$$\begin{array}{cc}\hfill \stackrel{\sim }{t}(z)& ={\beta }_0[T_L(1100)-T_L(z)]=\frac{2{\beta }_0}{3\sqrt{1-{\mathrm{\Omega }}_0}{H}_0}\hfill \\ \hfill & \times ln\left[\frac{{(1101)}^{\frac{3}{2}}(\sqrt{1-{\mathrm{\Omega }}_0}+\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{(1+z)}^3})}{{(1+z)}^{\frac{3}{2}}(\sqrt{1-{\mathrm{\Omega }}_0}+\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{(1101)}^3})}\right]\hfill \\ \hfill & =\frac{2{\beta }_0}{3\sqrt{1-{\mathrm{\Omega }}_0}{H}_0}ln\left[{\left(\frac{1101}{1+z}\right)}^{\frac{3}{2}}\frac{1+\sqrt{1+\frac{3}{7}{(1+z)}^3}}{1+\sqrt{1+\frac{3}{7}{(1101)}^3}}\right]\hfill \\ \hfill & =3.513·{10}^{17}\frac{{\beta }_0}{h_{70}}[0.424+ln\frac{1+\sqrt{1+\frac{3}{7}{(1+z)}^3}}{{(1+z)}^{\frac{3}{2}}}]\text{s}.\hfill \end{array}$$(39)

For later use, we note that

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\stackrel{\sim }{t}}{\mathrm{d}z}& =-{\beta }_0\frac{\mathrm{d}{T}_L(z)}{\mathrm{d}z}=-\frac{{\beta }_0}{H_0(1+z)\sqrt{1-{\mathrm{\Omega }}_0+{\mathrm{\Omega }}_0{(1+z)}^3}}\hfill \\ \hfill & =-\frac{1.5{\beta }_0{\tau }_f}{(1+z)\sqrt{{(1+z)}^3+\frac{7}{3}}}\simeq -\frac{3}{2}{\beta }_0{\tau }_f{(1+z)}^{-5/2},\hfill \end{array}$$(40)

with

$$\begin{array}{c}\hfill {\tau }_f=\frac{2}{3H_0\sqrt{{\mathrm{\Omega }}_0}}=5.37·{10}^{17}{h}_{70}^{-1}\text{s}.\end{array}$$(41)

Accordingly, the reduced time (39) for all values of z is well approximated by

$$\begin{array}{c}\hfill \stackrel{\sim }{t}(z)\simeq {(1+z)}^{-3/2}{\beta }_0{\tau }_f.\end{array}$$(42)

In Fig. 3, we show the variation of Eq. (39) as a function of the redshift z. Obviously, large redshifts correspond to early dynamical times $\stackrel{\sim }{t}$, whereas small redshifts correspond to late dynamical times. Also shown in Fig. 3 is the approximation (42).

Fig. 3. Semi-logarithmic $\stackrel{\sim }{t}(z)/{\beta }_0{\tau }_f$ (black solid line) versus z according to Eq. (39) with τf from Eq. (41). Dashed green curve shows the approximate $\stackrel{\sim }{t}(z)/{\beta }_0{\tau }_f={(1+z)}^{-3/2}$; i.e., Eq. (42).

In the following two sections, we calculated the SFR density and the stellar density for the case of spontaneous star formation without and with stellar feedback using our results from Sect. 3. We did not consider triggered star formation here. We adopted a flat Universe with Ω_tot = 1 so that the critical density would not evolve with redshift and is given by

$$\begin{array}{c}\hfill {\rho }_0(z)={\rho }_0=1.36·{10}^{11}{h}_{70}^2{\mathrm{M}}_{\odot }{\text{Mpc}}^{-3}.\end{array}$$(43)

Consequently, with Ω_m = 0.3 the matter density of the Universe is

$$\begin{array}{c}\hfill {\rho }_m=0.3{\rho }_0=4.1·{10}^{10}{h}_{70}^2{\mathrm{M}}_{\odot }{\text{Mpc}}^{-3}.\end{array}$$(44)

The matter density (44) includes contributions from dark matter and from baryonic matter, where the latter is a factor 8 smaller than the dark-matter contribution according to Davis et al. (2014). Hence, the baryonic matter density is given by

$$\begin{array}{c}\hfill {\rho }_b={\rho }_m/8=5.1·{10}^9{h}_{70}^2{\mathrm{M}}_{\odot }{\text{Mpc}}^{-3}=3.45·{10}^{-28}{h}_{70}^2\text{kg}{\text{m}}^{-3}.\end{array}$$(45)

In general, $G(\stackrel{\sim }{t})=j(\stackrel{\sim }{t})$ is the theoretical cosmic SFR density and the integrated stellar density as a function of redshift, ${\psi }_{\mathrm{GLS}}(z)\mathrm{d}z={\rho }_b\mathring{J}(t)\mathrm{d}\stackrel{\sim }{t}$, or, equivalently,

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{GSL}}(z)& ={\beta }_0{\rho }_{b}j(\stackrel{\sim }{t}(z))\left|\frac{\mathrm{d}\stackrel{\sim }{t}}{\mathrm{d}z}\right|\simeq \frac{3{\beta }_0^2{\rho }_b{\tau }_f}{2}\frac{j(\stackrel{\sim }{t}(z))}{{(1+z)}^{5/2}}\hfill \\ \hfill & \simeq A_1\frac{j(\stackrel{\sim }{t}(z))}{{(1+z)}^{5/2}},\hfill \end{array}$$(46)

with $A_1=3{\beta }_0^2{\rho }_b{\tau }_f/2=2.78·{10}^{-10}{\beta }_0^2{h}_{70}\text{kg}{\text{m}}^{-3}\text{s}$.

Using Eq. (46) thus provides for the integrated stellar mass density according to Eq. (33)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{GSL}}^{\ast }(z)& ={\beta }_0^{-1}{\displaystyle {\int }_z^{1100}}\mathrm{d}z\prime [1-L(z\prime )]{\psi }_{\mathrm{GS}}(z\prime )\left|\frac{\mathrm{d}\stackrel{\sim }{t}\prime }{\mathrm{d}z\prime }\right|\hfill \\ \hfill & \simeq \frac{3}{2}{A}_1{\tau }_f{\displaystyle {\int }_z^{1100}}\mathrm{d}z\prime [1-L(z\prime )]\frac{j(\stackrel{\sim }{t}(z\prime ))}{{(1+z\prime )}^5}\hfill \\ \hfill & =B_1{\displaystyle {\int }_z^{1100}}\mathrm{d}z\prime [1-L(z\prime )]\frac{j(\stackrel{\sim }{t}(z\prime ))}{{(1+z\prime )}^5},\hfill \end{array}$$(47)

with B₁ = 2.24 ⋅ 10⁸β₀² kg m⁻³ and where we have replaced 1 − R = 1 − L(z). This factor enters the calculation of ρ_GSL^*(z), but not the calculation of the SFR density (46) because all stars are born as luminous main-sequence stars. The fraction L(z) resulted from the stellar evolution to locked-in stellar matter, so only the fraction 1 − L(z) contributes to the observed integrated density of luminous stars.

4.2. Observational constraints

In the following sections, we compared the predictions of the GSL model without triggered star formation with the observations. We required the GSL model to be in agreement with the following five observational constraints.

First, the observed peak SFR density has large error bars,

$$\begin{array}{cc}\hfill \psi (z_E)& =0.{178}_{-0.044}^{+0.372}{\mathrm{M}}_{\odot }{\text{yr}}^{-1}{\text{Mpc}}^{-3}\hfill \\ \hfill & =(3.{83}_{-0.95}^{+8.00})·{10}^{-46}\text{kg}{\text{s}}^{-1}{\text{m}}^{-3},\hfill \end{array}$$(48)

and occurs in the redshift ranges of 1.62–1.88 (Dahlen et al. 2007) and 1.7–2.5 (Cucciati et al. 2012).

Secondly, we required the following observed peak redshift:

$$\begin{array}{c}\hfill z_E=2.0\pm 1.0.\end{array}$$(49)

As a third constraint, we used the observed integrated stellar mass density at z = 0 (Perez-Gonzalez et al. 2008; Moustakas et al. 2013):

$$\begin{array}{c}\hfill {\rho }^{\ast }(0)=(5.{62}_{-1.35}^{+1.79})·{10}^8\frac{{\mathrm{M}}_{\odot }}{{\text{Mpc}}^3}=(3.{80}_{-0.91}^{+1.21})·{10}^{-29}\text{kg}{\text{m}}^{-3.}\end{array}$$(50)

As a fourth constraint, we required (dψ/dz)_z = 0 > 0 in order to have at least one maximum of ψ(z) at positive z. As a fifth, more stringent constraint, we required ψ_GSL(z) to exhibit exactly one maximum within the redshift range z ∈ [0, 8]. While the first to fourth constraints could be dealt with analytically, the fifth constraint was verified by solving a highly nonlinear equation, Eq. (93), numerically.

5. Results for model without stellar evolution

We start with the spontaneous star formation without stellar evolution (a₀ = c₀ = 0). The omission of stellar evolution is justified as the evolution of most luminous stars occurs on a timescale on the order of the Hubble time (Portinari et al. 1998). We note that for c₀ = 0 the fraction (21) of nonluminous locked-in stellar matter vanishes.

5.1. SFR density

We employed the expression (29) for $j(\stackrel{\sim }{t})$. With $\stackrel{\sim }{t}(z)$ from Eq. (42), the SFR density (46) becomes

$$\begin{array}{c}\hfill {\psi }_{\mathrm{GS}}(z)=\frac{A_1}{1+r_1}\frac{r_1+e^{-{ϵ}_0{(1+z)}^{-3/2}}}{{(1+z)}^{5/2}},\end{array}$$(51)

where we introduced

$$\begin{array}{c}\hfill {ϵ}_0=(1+r_1){\beta }_0{\tau }_f=(b_0+{\beta }_0){\tau }_f=5.37·{10}^{17}(b_0+{\beta }_0)h_{70}^{-1}.\end{array}$$(52)

In terms of

$$\begin{array}{c}\hfill Z={ϵ}_0{(1+z)}^{-3/2}\end{array}$$(53)

the Eq. (51) reads

$$\begin{array}{c}\hfill {\psi }_{\mathrm{GS}}(z)=\frac{A_1}{(1+r_1){ϵ}_0^{5/3}}F(Z),F(Z)=Z^{5/3}(r_1+e^{-Z}).\end{array}$$(54)

The first derivative of the function F(Z) is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}F(Z)}{\mathrm{d}Z}=\frac{5}{3}{Z}^{2/3}[r_1+(1-0.6Z)e^{-Z}],\end{array}$$(55)

which is positive for all values of r₁ > 0.6 e^−5/3 ≈ 0.042, so the function F(Z) has no local extrema in this case.

In order to have a minimum of F(Z), corresponding to an extremum of ψ_GS(z), the ratio

$$\begin{array}{c}\hfill r_1=\frac{b_0}{{\beta }_0}<0.042\end{array}$$(56)

has to be rather small. For these values of r₁, we could already conclude that the gas fraction at infinitely long times,

$$\begin{array}{c}\hfill G(t=\infty )<\frac{r_1}{1+r_1}=\frac{b_0}{b_0+{\beta }_0}<0.040\end{array}$$(57)

is below 4%.

Provided restriction (56) holds, two extrema occur at

$$\begin{array}{c}\hfill Z_E=\frac{5}{3}-W(-\frac{5r_1}{3}{e}^{5/3})=\frac{5}{3}-W(-8.824r_1)\end{array}$$(58)

in terms of the two branches of the Lambert function discussed in Appendix A, where Z_E = ϵ₀(1 + z_E)^−3/2. The maximum of ψ_GS is determined by z_E employing the principal branch W₀, while W₋₁ determines the position of the minimum. We had to require that the minimum does not exist for positive redshift values z; that is

$$\begin{array}{c}\hfill \frac{5}{3}-W_{-1}(-\frac{5r_1{e}^{5/3}}{3})>{ϵ}_0,\end{array}$$(59)

corresponding to

$$\begin{array}{c}\hfill r_1<(\frac{3{ϵ}_0}{5}-1)e^{-{ϵ}_0},\end{array}$$(60)

which, since r₁ > 0, requires ϵ₀ > 5/3, or with Eq. (52) that

$$\begin{array}{c}\hfill {\beta }_0>3.10·{10}^{-18}\frac{h_{70}}{1+r_1}\text{Hz},\end{array}$$(61)

as first constraint on the value of the spontaneous SFR coefficient β₀. The maximum then occurs at the peak redshift:

$$\begin{array}{c}\hfill z_E={[\frac{5}{3{ϵ}_0}-\frac{1}{{ϵ}_0}{W}_0(-\frac{5r_1{e}^{5/3}}{3})]}^{-2/3}-1.\end{array}$$(62)

For small r₁ ≪ 1, the Taylor expansion yielded

$$\begin{array}{c}\hfill z_E\simeq \frac{(3-2r_1{e}^{5/3}){ϵ}_0^{2/3}}{{75}^{1/3}}-1\simeq {(0.6{ϵ}_0)}^{2/3}-1.\end{array}$$(63)

The peak SFR density (51) evaluated to

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{GS}}(z_E)& =\frac{3A_1}{5(1+r_1){ϵ}_0^{5/3}}{Z}_E^{8/3}{e}^{-Z_E}\hfill \\ \hfill & =-\frac{A_1{r}_1}{(1+r_1){ϵ}_0^{5/3}}\frac{{[\frac{5}{3}-W_0(-\frac{5r_1{e}^{5/3}}{3})]}^{8/3}}{W_0(-\frac{5r_1{e}^{5/3}}{3})}\hfill \\ \hfill & =-\frac{A_1{r}_1{ϵ}_0}{(1+r_1)W_0(-\frac{5r_1{e}^{5/3}}{3}){(1+z_E)}^4}\hfill \\ \hfill & =-\frac{5.37·{10}^{17}{A}_1{r}_1{\beta }_0{h}_{70}^{-1}}{W_0(-\frac{5r_1{e}^{5/3}}{3}){(1+z_E)}^4}\simeq 1.69·{10}^7\frac{{\beta }_0^3}{{(1+z_E)}^4}\hfill \\ \hfill & =3.46·{10}^{-40}\frac{{\beta }_0^{1/3}{h}_{70}^{8/3}}{{(1+r_1)}^{8/3}}\text{kg}{\text{s}}^{-1}{\text{m}}^{-3},\hfill \end{array}$$(64)

where we used the expansion W₀(|x|≪1)≃x and Eq. (63). Equating (64) and (48) yielded

$$\begin{array}{c}\hfill {\beta }_0=(1.{36}_{-0.78}^{+38.61})·{10}^{-18}{\left(\frac{1+r_1}{h_{70}}\right)}^8\text{Hz}.\end{array}$$(65)

This β₀-value² is partially consistent with the lower limit (61), especially for large values of the observed peak SFR density (48). Consequently, according to Eqs. (52) and (63) we obtained

$$\begin{array}{c}\hfill {ϵ}_0=(0.{73}_{-0.42}^{+20.73}){\left(\frac{1+r_1}{h_{70}}\right)}^9,z_E<4.49{\left(\frac{1+r_1}{h_{70}}\right)}^6;\end{array}$$(66)

this large upper limit on z_E, however, results from adopting large values of the observed peak SFR density (49) within its large error bars. The upper constraint on z_E is consistent with the observational requirement (49).

Consequently, for large values of the observed peak SFR as well as the spontaneous formation rate (65), and values of r₁ below the small value (56), the GS model with spontaneous star formation indeed provided a rough but far from optimal fit to the observed SFR density (see Fig. 4).

Fig. 4. GS model with zero triggered star formation (a0 = c0 = 0) with r1 = 0.001 and ${\stackrel{\sim }{\beta }}_0=2$ for reasons described in the text. (a) ψGS(z) and (b) ρ*(z) obtained numerically (black solid line) and analytically (green dashed line; Eqs. (54) and (68), respectively). Symbols with error bars are data for ψ(z) and ρ*(z) tabulated by Madau & Dickinson (2014).

5.2. Stellar density

Here, we calculated the integrated stellar density for the spontaneous star formation model. With the help of

$$\begin{array}{c}\hfill {\displaystyle {\int }^z}\mathrm{d}z\frac{e^{-\alpha {(1+z)}^{-3/2}}}{{(1+z)}^5}=\frac{2}{3{\alpha }^{8/3}}\mathrm{\Gamma }[\frac{8}{3},\frac{\alpha }{{(1+z)}^{3/2}}],\end{array}$$(67)

the integrated stellar mass density (47) with L(z) = 0 becomes

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{GS}}^{\ast }(z)& =\frac{B_1}{1+r_1}{\displaystyle {\int }_z^{1100}}\mathrm{d}z\prime \frac{r_1+e^{-{ϵ}_0{(1+z\prime )}^{-3/2}}}{{(1+z\prime )}^5}\hfill \\ \hfill & =\frac{B_1}{1+r_1}(\frac{r_1}{4}(\frac{1}{{(1+z)}^4}-{1100}^{-4})\hfill \\ \hfill & +\frac{2}{3{ϵ}_0^{8/3}}[\mathrm{\Gamma }(\frac{8}{3},{(1100)}^{3/2}{ϵ}_0)-\mathrm{\Gamma }(\frac{8}{3},\frac{{ϵ}_0}{{(1+z)}^{3/2}})])\hfill \\ \hfill & \simeq \frac{2.24·{10}^8{\beta }_0^2}{4(1+r_1)}[\frac{r_1}{{(1+z)}^4}+\frac{8}{3}\frac{\gamma (\frac{8}{3},\frac{{ϵ}_0}{{(1+z)}^{3/2}})}{{ϵ}_0^{8/3}}]\text{kg}{\text{m}}^{-3},\hfill \\ \hfill \end{array}$$(68)

where γ is the lower incomplete gamma function. The stellar mass density thus monotonically decreases from the value

$$\begin{array}{c}\hfill {\rho }_{\mathrm{GS}}^{\ast }(0)=\frac{5.60·{10}^7{\beta }_0^2}{1+r_1}(r_1+\frac{8}{3}\frac{\gamma (\frac{8}{3},{ϵ}_0)}{{ϵ}_0^{8/3}})\text{kg}{\text{m}}^{-3}\end{array}$$(69)

in a power-law fashion at large redshifts z ≫ z_E,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{GSL}}^{\ast }(z)\simeq \frac{B_1}{4{(1+z)}^4},\end{array}$$(70)

since γ(s, ϵ ≪ 1)≃s⁻¹ϵ^s to the lowest order in ϵ. Equating the result (69) with the observed value (50) provides

$$\begin{array}{c}\hfill \frac{{\stackrel{\sim }{\beta }}_0^2}{1+r_1}(r_1+\frac{8}{3}\frac{\gamma (\frac{8}{3},{ϵ}_0)}{{ϵ}_0^{8/3}})=(6.{79}_{-1.70}^{+2.16})·{10}^{-3},\end{array}$$(71)

where we scaled

$$\begin{array}{c}\hfill {\beta }_0={10}^{-17}{\stackrel{\sim }{\beta }}_0\end{array}$$(72)

in units of 10⁻¹⁷ Hz. The function $\gamma (8/3,5.37{\stackrel{\sim }{\beta }}_0)$ is shown in Fig. 5 and varies from 0 to the finite value Γ(8/3)≈1.504 at large arguments. For the nominal values r₁ = 0 and h₇₀ = 1, we obtained the following from Eq. (71)

$$\begin{array}{c}\hfill \frac{{\stackrel{\sim }{\beta }}_0^{2/3}}{\gamma (8/3,5.37{\stackrel{\sim }{\beta }}_0)}\simeq \frac{{\stackrel{\sim }{\beta }}_0^{2/3}}{1.504}=4.{44}_{-1.07}^{+1.48},\end{array}$$(73)

Fig. 5. Lower incomplete gamma function, $\gamma (8/3,5.37{\stackrel{\sim }{\beta }}_0)$, involved in the analytic expressions for ρGS*(z) and ρGSL*(z), versus ${\stackrel{\sim }{\beta }}_0$. Particular values are 0.0487, 0.21, 0.88, 1.40, and 1.50 for ${\stackrel{\sim }{\beta }}_0=0.1$, 0.2, 0.5, 1.0, and 2, respectively.

which leads to

$$\begin{array}{c}\hfill {\stackrel{\sim }{\beta }}_0=17.{26}_{-5.85}^{+9.31}.\end{array}$$(74)

The two determinations (65) and (74) of the rate of spontaneous star formation, resulting from reproducing the peak value of the SFR density and the present-day integrated stellar mass density within their partially large uncertainties, exclude each other by at a factor of at least 2.8. For example, it is not possible to explain these two constraints with the GS model with one single ${\stackrel{\sim }{\beta }}_0$ value.

Adopting ${\stackrel{\sim }{\beta }}_0=2.0$ as a best compromise value, we compared the SFR density and the integrated stellar densities for r₁ = 0.001 in Fig. 4. This choice implied z_E = 2.46 as peak redshift. Only because of the large error bars in the peak SFR density are the theoretical results of the simplified GS model marginally consistent with the observations. Figure 4 also illustrates the basic dilemma of this simplified model very well. Because of the dependencies ${\psi }_{\mathrm{GS}}(z_E)\propto {\stackrel{\sim }{\beta }}_0^{1/3}$ in Eq. (64) and ${\rho }_{\mathrm{GS}}^{\ast }(0)\propto {\stackrel{\sim }{\beta }}_0^{-2/3}$ in Eq. (69) on the only free parameter, ${\stackrel{\sim }{\beta }}_0$, it is obvious that by decreasing ${\stackrel{\sim }{\beta }}_0$ from its compromise value 2 would decrease ψ_GS(z) to smaller values, but it would increase ρ_GS^*(z) at the same time. Alternatively, increasing ${\stackrel{\sim }{\beta }}_0$ from its compromise value 2 would decrease ρ_GS^*(z) to smaller values, but at the same time it would make ψ_GS(z) even larger. With the compromise value ${\stackrel{\sim }{\beta }}_0=2.0$ and the ratio r₁ = 0.04, the present-day gas fraction is

$$\begin{array}{c}\hfill G(r_2=0,z=0)=3.85·{10}^{-2},\end{array}$$(75)

which is basically identical to the gas fraction at infinite time G_∞, and in marginal agreement with the observation (31).

5.3. Special case: Negligible stellar feedback (r₁ = 0) and stellar evolution (r₂ = 0)

The results of the last subsections included the GS model without stellar feedback (r₁ = 0) and stellar evolution (r₂ = 0) as a special case. As in this case some of the theoretical results are severely simplified and thus more transparent, we include a detailed analysis of this case here. The SFR density (51) simplified to

$$\begin{array}{c}\hfill {\psi }_{\mathrm{GS}}(z)=A_1\frac{e^{-ϵ{(1+z)}^{-3/2}}}{{(1+z)}^{5/2}}\end{array}$$(76)

with

$$\begin{array}{c}\hfill ϵ={ϵ}_0(r_1=0)=5.37{\stackrel{\sim }{\beta }}_0{h}_{70}^{-1}.\end{array}$$(77)

The density (76) attains its maximum value,

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{GS}}(z_E)& ={\left(\frac{5}{3ϵ}\right)}^{5/3}{A}_1{e}^{-\frac{5}{3}}\hfill \\ \hfill & =7.47·{10}^{-46}{\stackrel{\sim }{\beta }}_0^{1/3}{h}_{70}^{8/3}\text{kg}{\text{s}}^{-1}{\text{m}}^{-3}\hfill \end{array}$$(78)

at the peak redshift

$$\begin{array}{c}\hfill z_E={(0.6ϵ)}^{2/3}-1=2.18{\left(\frac{{\stackrel{\sim }{\beta }}_0}{h_{70}}\right)}^{2/3}-1.\end{array}$$(79)

The two observational constraints (49) and (48) for z_E and the peak SFR density then provided the constraints

$$\begin{array}{c}\hfill {\stackrel{\sim }{\beta }}_0=(1.{61}_{-0.73}^{+0.87})h_{70}\end{array}$$(80)

and

$$\begin{array}{c}\hfill {\stackrel{\sim }{\beta }}_0=(0.{14}_{-0.08}^{+3.84})h_{70}^{-8},\end{array}$$(81)

respectively, which are consistent which each other due to the large uncertainty in the peak SFR density.

The integrated stellar mass density (68) simplifies to

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{GS}}^{\ast }(z)& =B_1{\displaystyle {\int }_z^{1100}}\mathrm{d}z\prime \frac{e^{ϵ{(1+z\prime )}^{-3/2}}}{{(1+z\prime )}^5}\simeq \frac{2B_1}{3{ϵ}^{8/3}}\gamma (\frac{8}{3},\frac{ϵ}{{(1+z)}^{3/2}})\hfill \\ \hfill & \simeq \frac{1.49·{10}^{-26}{\stackrel{\sim }{\beta }}_0^2}{{ϵ}^{8/3}}\gamma (\frac{8}{3},\frac{ϵ}{{(1+z)}^{3/2}})\text{kg}{\text{m}}^{-3}\hfill \\ \hfill & =\frac{1.69·{10}^{-28}{h}_{70}^{8/3}}{{\stackrel{\sim }{\beta }}_0^{2/3}}\gamma (\frac{8}{3},\frac{5.37{\stackrel{\sim }{\beta }}_0}{h_{70}{(1+z)}^{3/2}})\text{kg}{\text{m}}^{-3}.\hfill \end{array}$$(82)