Compartmental Description of the Cosmological Baryonic Matter Cycle. I. Competition of Spontaneous Star Formation, Stellar Feedback and Stellar Evolution

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 [29,35,45]. Stars form from the collapse of interstellar gas clouds consisting mainly of hydrogen and helium, evolve with time undergoing nuclear fusion reactions thereby creating metals i.e. elements heavier than helium [3,5,25,36]. It is well established that external pressure on gas clouds reduces the critical cloud mass for collapse to the Bonnor-Ebert mass [1,11] being 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 as well as 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 [9,39,45]. Moreover, at the termination of stellar fusion reactions the powerful novae and supernovae resulting from the fatal star’s collapse eject further stellar material to 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 [14,15,24,37,41,47]. While the locked-in matter no longer can fed the gaseous component, the supersonic nova and supernova outflows can trigger efficient star formation in the ambient ISM by providing strong enough external pressure [6,12,16,26,28,31,43].

It is the purpose of the present study 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 analysing and/or predicting the temporal development of infection diseases and epidemics since their invention about a hundred years ago by [21] and the subsequent refinement by [20]. Here persons from a considered population with many N persons are assigned to different compartments. In the simple 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\to I$ and $I\to R$, respectively (for reviews see [13,17,18].

In the present manuscript we restrict the analysis to the case of spontaneous star formation only, 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 Section 2 we introduce the general compartmental description of the cosmic baryonic matter by formulating the relevant dynamical equations including spontaneous as well as triggered star formation. In Section 3. the dynamical equations are solved analytically for the case of negligible spontaneous star formation. The analytical solutions are applied to the present-day (redshift $z=0$) gaseous fraction of the universe and the observed cosmic star formation history (Section 4. The comparison of the theoretical results with the observations are the subjects of Section 5 (in the special case of neglected stellar feedback) and Section 6 (including stellar feedback). The summary and conclusion (Section 7 completes the manuscript.

2. GSL-Compartmental Model for the Baryonic Matter Cycle

We consider 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\left(t\right)$, $S\left(t\right)$ and $L\left(t\right)$ denote the relative fractions of baryonic matter in the three compartments, respectively, as a function of time t, so that the sum constraint

$$G\left(t\right)+S\left(t\right)+L\left(t\right)=1$$

(1)

holds at all times t after the begin of the baryonic evolution at time $t=t_0$.

The temporal evolution of the three fractions $G\left(t\right)$, $S\left(t\right)$ and $L\left(t\right)$ is controlled by the spontaneous star formation rate (SFR) $\beta \left(t\right)G\left(t\right)$, the triggered star formation from the interaction of gas and stars with the assumed rate $a\left(t\right)G\left(t\right)S\left(t\right)$, the continuous feedback rate $b\left(t\right)S\left(t\right)$ of stellar to gaseous matter, and the formation rate $c\left(t\right)S\left(t\right)$ of white dwarfs, neutron stars and black holes from stellar evolution. The non-zero feedback rate is essential for the existence of gaseous heavier metals beyond helium and lithium from processed stellar baryons in the universe.

The dynamical equations for the three fractions then read

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dG}{dt}}& =& -a\left(t\right)G\left(t\right)S\left(t\right)-\beta \left(t\right)G\left(t\right)+b\left(t\right)S\left(t\right),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& a\left(t\right)G\left(t\right)S\left(t\right)+\beta \left(t\right)G\left(t\right)-b\left(t\right)S\left(t\right)-c\left(t\right)S\left(t\right),\hfill \end{array}$$

(3)

and

$${\displaystyle \frac{dL}{dt}}=c\left(t\right)S\left(t\right),$$

(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 [22]. As initial condition we adopt

$$G\left(t_0\right)=1,S\left(t_0\right)=0,L\left(t_0\right)=0,$$

(5)

as the whole baryonic gaseous matter was subject to spontaneous star formation. The initial time $t_0$ corresponds to the redshift $z=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 [44] 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 consider 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 [42] or gas regulator [33] models for galaxy evolution considering the competition between the flow of gas into galaxies, conversion into stars by in-situ star formation and ejection out of galaxies by stellar feedback. including individual galaxies, clusters of galaxies and the whole universe.

We are particularly interested in the formation rate $\mathring{J}\left(t\right)$ of new stars as a function of time,

$$\mathring{J}\left(t\right)=a\left(t\right)G\left(t\right)S\left(t\right)+\beta \left(t\right)G\left(t\right).$$

(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 restrict 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) can 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\left(t\right)=0$, the GSL-model Equations (2)–(4) simplify to the linear equations

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dG}{dt}}& =& -\beta \left(t\right)G\left(t\right)+b\left(t\right)S\left(t\right),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \beta \left(t\right)G\left(t\right)-b\left(t\right)S\left(t\right)-c\left(t\right)S\left(t\right),\hfill \end{array}$$

(8)

and

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dL}{dt}}& =& c\left(t\right)S\left(t\right).\hfill \end{array}$$

(9)

Introducing the dimensionless time variable

$$\tilde{t}={\int }_{t_0}^{t}d\xi \beta \left(\xi \right),$$

(10)

Equations (7)–(9) become

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dG}{d\tilde{t}}}& =& -G\left(\tilde{t}\right)+r_1\left(\tilde{t}\right)S\left(\tilde{t}\right),\hfill \end{array}$$

(11a)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{d\tilde{t}}}& =& G\left(\tilde{t}\right)-[r_1\left(\tilde{t}\right)+r_2\left(\tilde{t}\right)]S\left(\tilde{t}\right),\hfill \end{array}$$

(11b)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dL}{d\tilde{t}}}& =& r_2\left(\tilde{t}\right)S\left(\tilde{t}\right),\hfill \end{array}$$

(11c)

$$\begin{array}{ccc}\hfill 1& =& G\left(\tilde{t}\right)+S\left(\tilde{t}\right)+L\left(\tilde{t}\right),\hfill \end{array}$$

(11d)

with the dimensionless ratios

$$r_1\left(\tilde{t}\left(t\right)\right)={\displaystyle \frac{b\left(t\right)}{\beta \left(t\right)}},r_2\left(\tilde{t}\left(t\right)\right)={\displaystyle \frac{c\left(t\right)}{\beta \left(t\right)}}.$$

(12)

3.1. Stationary Ratios

In the following we adopt stationary ratios $r_1=$ const. and $r_2=$ const., so that our solution allows for arbitrary but given time dependent spontaneous rates $\beta \left(t\right)$.

The first Equation (11a) readily provides

$$S\left(\tilde{t}\right)={\displaystyle \frac{1}{r_1}}\left[{\displaystyle \frac{dG\left(\tilde{t}\right)}{d\tilde{t}}}+G\left(\tilde{t}\right)\right],$$

(13)

which upon insertion in Equation (11b) yields

$${\displaystyle \frac{d^{2}G\left(\tilde{t}\right)}{d{\tilde{t}}^2}}+(1+r_1+r_2){\displaystyle \frac{dG\left(\tilde{t}\right)}{d\tilde{t}}}+r_{2}G\left(\tilde{t}\right)=0.$$

(14)

Equation (14) is solved by

$$G\left(\tilde{t}\right)=e^{-\mu \tilde{t}}\left[c_{-}{e}^{\nu \tilde{t}}+c_{+}{e}^{-\nu \tilde{t}}\right],$$

(15)

with

$$\mu ={\displaystyle \frac{1+r_1+r_2}{2}},\nu =\sqrt{{\mu }^2-r_2}$$

(16)

and with the two integration constants $c_{\pm }$, implying for Equation (13)

$$S\left(\tilde{t}\right)={\displaystyle \frac{e^{-\mu \tilde{t}}}{r_1}}\left[c_{-}\left(1-\mu +\nu \right)e^{\nu \tilde{t}}+c_{+}\left(1-\mu -\nu \right)e^{-\nu \tilde{t}}\right].$$

(17)

The initial conditions (5), i.e. $G(z=0)=1$ and $S(z=0)=0$ provide

$$c_{\mp }={\displaystyle \frac{\nu \mp (1-\mu )}{2\nu }}$$

(18)

so that

$$\begin{array}{ccc}\hfill j\left(\tilde{t}\right)& =& G\left(\tilde{t}\right)=e^{-\mu \tilde{t}}[cosh\left(\nu \tilde{t}\right)-{\displaystyle \frac{1-\mu }{\nu }}sinh\left(\nu \tilde{t}\right)]\hfill \end{array}$$

(19)

and

$$S\left(\tilde{t}\right)={\displaystyle \frac{1}{\nu }}{e}^{-\mu \tilde{t}}sinh\nu t$$

(20)

Consequently,

$$\begin{array}{ccc}\hfill L\left(\tilde{t}\right)& =& 1-G\left(\tilde{t}\right)-S\left(\tilde{t}\right)=1-e^{-\mu \tilde{t}}\left[cosh\left(\nu \tilde{t}\right)+{\displaystyle \frac{\mu }{\nu }}sinh\left(\nu \tilde{t}\right)\right]\hfill \end{array}$$

(21)

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

$$S_{\max }={\displaystyle \frac{{(\mu -\nu )}^{\frac{\mu -\nu }{2\nu }}}{{(\mu +\nu )}^{\frac{\mu +\nu }{2\nu }}}}$$

(22)

at

$${\tilde{t}}_{\mathrm{S},\max }={\displaystyle \frac{1}{2\nu }}ln{\displaystyle \frac{\mu +\nu }{\mu -\nu }},$$

(23)

and the S approaches zero in the course of time, $S_{\infty }=0$.

The new star formation rate (6) for stationary SFR coefficient ${\beta }_0$ in this case is

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

(24)

3.2. Negligible Stellar Feedback

In the special case of negligible stellar feedback $b\left(t\right)=0=r_1$ one obtains

$$\begin{array}{ccc}\hfill j(\tilde{t},r_1=0)& =& G(\tilde{t},r_1=0)\hfill \\ & =\hfill & e^{-{\displaystyle \frac{1+r_2}{2}}\tilde{t}}\left[cosh{\displaystyle \frac{1-r_2}{2}}\tilde{t}-sinh{\displaystyle \frac{1-r_2}{2}}\tilde{t}\right]=e^{-\tilde{t}},\hfill \end{array}$$

(25)

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

$$S(\tilde{t},r_1=0)={\displaystyle \frac{2}{1-r_2}}{e}^{-{\displaystyle \frac{1+r_2}{2}}\tilde{t}}sinh\left({\displaystyle \frac{1-r_2}{2}}\tilde{t}\right)={\displaystyle \frac{e^{-r_2\tilde{t}}-e^{-\tilde{t}}}{1-r_2}},$$

(26)

and

$$L(\tilde{t},r_1=0)=1-{\displaystyle \frac{e^{-r_2\tilde{t}}-r_2{e}^{-\tilde{t}}}{1-r_2}}.$$

(27)

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

In the special case of $r_2=1$ Equations (25)–(27) simplify to

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

(28)

Here $S\left(\tilde{t}\right)$ rises first to the maximum value $S_{\max }=e^{-1}$ at ${\tilde{t}}_{\max }=1$ and then decreases to the final value $S_{\infty }=0$.

3.3. Negligible Stellar Evolution

In the special case of negligible stellar evolution $c\left(t\right)=0=r_2$ one obtains

$$j(\tilde{t},c=0)=G(\tilde{t},c=0)={\displaystyle \frac{1}{1+r_1}}\left[r_1+e^{-(1+r_1)\tilde{t}}\right],$$

(29)

corresponding to

$$\mathring{J}(t,c=0)={\displaystyle \frac{{\beta }_0}{b_0+{\beta }_0}}\left[b_0+{\beta }_0{e}^{-(b_0+{\beta }_0)(t-t_0)}\right],$$

(30)

and $L(\tilde{t},c=0)=0$. In Figure 1 we compare the derived exact analytic solutions with the numerical solutions of Equations (7)–(9) for several illustrative examples. For the numerical code we use a single-step solver based on a modified Rosenbeck formula of order 2, implemented by [40] as ode23s in ${\mathrm{Matlab}}^{\mathrm{TM}}$. The excellent agreement proofs the validity of our derivations.

4. Cosmic Star Formation History

As astrophysical application we consider 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)/\left[S\right(z=0)+L(z=0\left)\right]=G(z=0)/[1-G(z=0\left)\right]$ in our Milky Way is about 0.1 and a bit 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

$$G(z=0)\simeq 0.10\pm 0.05$$

(31)

seems to be reasonable.

With respect to the cosmic SFH far ultraviolet and infrared measurements have indicated [27] that the best fit SFR density per time t is

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

(32)

where the time $t\left(z\right)$ has been expressed as a function of redshift z using a flat $\Lambda$CDM Friedmann cosmology [4,32,38,46], and where $M_{\odot }=1.989·{10}^{30}$ kg in SI-units1. The data used for the best fit (32) are given in Table 1 of [27]. The SFR density (32) attains its maximum value ${\psi }_{\max }=0.133$ ${\mathrm{M}}_{\odot }$ ${\mathrm{yr}}^{-1}$ ${\mathrm{Mpc}}^{-3}$ at redshift $z=1.8632$, is practically constant at small redshifts well below unity, and it decreases at large redshifts $\psi (z\gg 2.9)\simeq 1.25·{10}^{-44}{(1+z)}^{-2.9}$ ${\mathrm{kg}}^{-1}$ ${\mathrm{m}}^{-3}$.

Also available in Table 2 and Figure 11 of [27] is the integrated stellar density ${\rho }_{*}\left(z\right)$ which according to [27] fulfills $d{\rho }_{*}\left(z\right)/dt=(1-R)\psi \left(z\right)$ and is therefore more explicitly given by

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

(33)

where according to [27] R is the return factor or the mass fraction of newly born stars that is put back into the interstellar gas which has been treated as free fit parameter. As we will argue below within the compartmental model considered here $1-R=1-L\left(t\right(z\left)\right)$ has to be replaced by the time- or redshift dependent term $1-L\left(t\right)$ with the fraction $L\left(t\right)$ of locked-in stellar matter.

We first discuss the relation of our dynamical evolution time scale $t=\tilde{t}/{\beta }_0$ for stationary spontaneous SFRs with the lookback time and redshift z in a flat $\Lambda$CDM Friedmann cosmology. In a flat $\Lambda$CDM Friedmann cosmology with ${\Omega }_m={\Omega }_0=0.3$ and $H_0=70h_{70}$ km ${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}=2.268·{10}^{-18}{h}_{70}$ ${\mathrm{s}}^{-1}$ the lookback time $T_L$ is related to the redshift according to [2,19,32]

$$\begin{array}{ccc}\hfill T_L\left(z\right)& =& H_0^{-1}{\int }_1^{1+z}{\displaystyle \frac{dx}{x\sqrt{1-{\Omega }_0+{\Omega }_0{x}^3}}}\hfill \\ & =& {\displaystyle \frac{2}{3H_0}}{\int }_1^{{(1+z)}^{3/2}}{\displaystyle \frac{du}{u\sqrt{1-{\Omega }_0+{\Omega }_0{u}^2}}}\hfill \\ & =& {\displaystyle \frac{2}{3\sqrt{1-{\Omega }_0}{H}_0}}ln{\displaystyle \frac{(1+\sqrt{1-{\Omega }_0}{(1+z)}^{\frac{3}{2}}}{\sqrt{1-{\Omega }_0}+\sqrt{1-{\Omega }_0+{\Omega }_0{(1+z)}^3}}}.\hfill \end{array}$$

(34)

The lookback 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 lookback time (34) corresponds to the age of the universe $T_0$ given by

$$\begin{array}{ccc}\hfill T_0& =& H_0^{-1}{\int }_1^{\infty }{\displaystyle \frac{dx}{x\sqrt{1-{\Omega }_0+{\Omega }_0{x}^3}}}\hfill \\ & =& {\displaystyle \frac{2}{3\sqrt{1-{\Omega }_0}{H}_0}}ln{\displaystyle \frac{1+\sqrt{1-{\Omega }_0}}{\sqrt{{\Omega }_0}}}.\hfill \end{array}$$

(35)

The relation between the real time t, the lookback time $T_L\left(z\right)$ and the age of the universe $T_0$ is given by

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

(36)

implying

$$t_0=T_0-T_L(z=1100).$$

(37)

Consequently,

$$t-t_0=T_L\left(1100\right)-T_L\left(z\right),$$

(38)

so that for a constant spontaneous SFR coefficient $\beta \left(t\right)={\beta }_0$ one obtains

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

(39)

For later use we note that

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d\tilde{t}}{dz}}& =& -{\beta }_0{\displaystyle \frac{d{T}_L\left(z\right)}{dz}}=-{\displaystyle \frac{{\beta }_0}{H_0(1+z)\sqrt{1-{\Omega }_0+{\Omega }_0{(1+z)}^3}}}\hfill \\ & =& -{\displaystyle \frac{1.5{\beta }_0{\tau }_f}{(1+z)\sqrt{{(1+z)}^3+\frac{7}{3}}}}\simeq -{\displaystyle \frac{3}{2}}{\beta }_0{\tau }_f{(1+z)}^{-5/2},\hfill \end{array}$$

(40)

with

$${\tau }_f={\displaystyle \frac{2}{3H_0\sqrt{{\Omega }_0}}}=5.37·{10}^{17}{h}_{70}^{-1}\mathrm{s}.$$

(41)

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

$$\tilde{t}\left(z\right)\simeq {(1+z)}^{-3/2}{\beta }_0{\tau }_f.$$

(42)

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

In the following two sections we calculate the SFR density and the stellar density for the case of spontaneous star formation without and with stellar feedback using our results from the earlier Section 3. We do not consider triggered star formation here. We adopt a flat universe with ${\Omega }_{\mathrm{tot}}=1$ so that the critical density does not evolve with redshift and is given by

$${\rho }_0\left(z\right)={\rho }_0=1.36·{10}^{11}{h}_{70}^2{M}_{\odot }{\mathrm{Mpc}}^{-3}.$$

(43)

Consequently, with ${\Omega }_m=0.3$ the matter density of the universe then is

$${\rho }_m=0.3{\rho }_0=4.1·{10}^{10}{h}_{70}^2{M}_{\odot }{\mathrm{Mpc}}^{-3}.$$

(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 [10]. Hence the baryonic matter density is given by

$${\rho }_b={\rho }_m/8=5.1·{10}^9{h}_{70}^2{M}_{\odot }{\mathrm{Mpc}}^{-3}=3.45·{10}^{-28}{h}_{70}^2\mathrm{kg}{\mathrm{m}}^{-3}.$$

(45)

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

$$\begin{array}{ccc}\hfill {\psi }_{\mathrm{GSL}}\left(z\right)& =& {\beta }_0{\rho }_{b}j\left(\tilde{t}\left(z\right)\right)\left|{\displaystyle \frac{d\tilde{t}}{dz}}\right|\simeq {\displaystyle \frac{3{\beta }_0^2{\rho }_b{\tau }_f}{2}}{\displaystyle \frac{j\left(\tilde{t}\left(z\right)\right)}{{(1+z)}^{5/2}}}\hfill \\ & \simeq & A_1{\displaystyle \frac{j\left(\tilde{t}\left(z\right)\right)}{{(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}\mathrm{kg}{\mathrm{m}}^{-3}\mathrm{s}$.

Using Equation (46) then provides for the integrated stellar mass density according to Equation (33)

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

(47)

with $B_1=2.24·{10}^8{\beta }_0^2$ kg ${\mathrm{m}}^{-3}$, where we have replaced $1-R=1-L\left(z\right)$. This factor enters the calculation of ${\rho }_{\mathrm{GSL}}^{*}\left(z\right)$ but not the calculation of the SFR density (46) because all stars are born as luminous main sequence stars. The fraction $L\left(z\right)$ results from the stellar evolution to locked-in stellar matter, so that only the fraction $1-L\left(z\right)$ contributes to the observed integrated density of luminous stars.

4.1. Observational Constraints

In the following Sections. we compare the predictions the GSL-model without triggered star formation with the observations. We require for the GSL-model to be in accord with the following five observational constraints:

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

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

(48)

and occurs in the redshift range $1.62-1.88$[8] and $1.7-2.5$ [7].

Secondly, we require for the observed peak redshift

$$z_E=2.0\pm 1.0.$$

(49)

As third constraint we use the observed integrated stellar mass density at $z=0$ [30,34]

$${\rho }^{*}\left(0\right)=(5.{62}_{-1.35}^{+1.79}·{10}^8{\displaystyle \frac{M_{\odot }}{{\mathrm{Mpc}}^3}}=(3.{80}_{-0.91}^{+1.21}·{10}^{-29}\mathrm{kg}{\mathrm{m}}^{-3}$$

(50)

As fourth constraint we demand that ${(d\psi /dz)}_{z=0}>0$ in order to have at least one maximum of $\psi \left(z\right)$ at positive z. As fifth more stringent constraint we require that ${\psi }_{\mathrm{GSL}}\left(z\right)$ exhibits exactly one maximum within the redshift range $z\in [0,8]$. While constraints 1 to 4 can be dealt with analytically, the fifth constraint is verified by solving a highly nonlinear equation, the later Equation (93), numerically.

5. Results for Model without Stellar Evolution

We start with the spontaneous star formation without stellar evolution ($a_0=c_0=0$). The omission of stellar evolution is justified as the evolution of most luminous stars occurs on a time scale of the order of the Hubble time [37]. Note that for $c_0=0$ the fraction (21) of non-luminous locked-in stellar matter vanishes.

5.1. SFR Density

We employ the expression (29) for $j\left(\tilde{t}\right)$. With $\tilde{t}\left(z\right)$ from Equation (42), the SFR density (46) becomes

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

(51)

where we introduce

$${ϵ}_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}.$$

(52)

In terms of

$$Z={ϵ}_0{(1+z)}^{-3/2}$$

(53)

Equation (51) reads

$${\psi }_{\mathrm{GS}}\left(z\right)={\displaystyle \frac{A_1}{(1+r_1){ϵ}_0^{5/3}}}F\left(Z\right),F\left(Z\right)=Z^{5/3}(r_1+e^{-Z}).$$

(54)

The first derivative of the function $F\left(Z\right)$ is given by

$${\displaystyle \frac{dF\left(Z\right)}{dZ}}={\displaystyle \frac{5}{3}}{Z}^{2/3}[r_1+(1-0.6Z)e^{-Z}],$$

(55)

which is positive for all values of $r_1>0.6e^{-5/3}\approx 0.042$, so that the function $F\left(Z\right)$ has no local extrema in this case.

In order to have a minimum of $F\left(Z\right)$, corresponding to an extremum of ${\psi }_{\mathrm{GS}}\left(z\right)$, the ratio

$$r_1={\displaystyle \frac{b_0}{{\beta }_0}}<0.042$$

(56)

has to be rather small. For these values of $r_1$ we can conclude already that the gas fraction at infinitely large times

$$G(t=\infty )<{\displaystyle \frac{r_1}{1+r_1}}={\displaystyle \frac{b_0}{b_0+{\beta }_0}}<0.040$$

(57)

is below 4 percent.

Provided restriction (56) holds, two extrema occur at

$$Z_E={\displaystyle \frac{5}{3}}-W\left(-{\displaystyle \frac{5r_1}{3}}{e}^{5/3}\right)={\displaystyle \frac{5}{3}}-W(-8.824r_1)$$

(58)

in terms of the two branches of the Lambert function discussed in Appendix A, where $Z_E={ϵ}_0{(1+z_E)}^{-3/2}$. The maximum of ${\psi }_{\mathrm{GS}}$ is determined by $z_E$ employing the principal branch $W_0$, while $W_{-1}$ determines the position of the minimum. We need to require that the minimum does not exist for positive redshift values z, i.e.

$${\displaystyle \frac{5}{3}}-W_{-1}\left(-{\displaystyle \frac{5r_1{e}^{5/3}}{3}}\right)>{ϵ}_0,$$

(59)

corresponding to

$$r_1<\left({\displaystyle \frac{3{ϵ}_0}{5}}-1\right)e^{-{ϵ}_0},$$

(60)

which, since $r_1>0$, requires ${ϵ}_0>5/3$ or with Equation (52) that

$${\beta }_0>3.10·{10}^{-18}{\displaystyle \frac{h_{70}}{1+r_1}}\mathrm{Hz},$$

(61)

as first constraint on the value of the spontaneous SFR coefficient ${\beta }_0$. The maximum then occurs at the peak redshift

$$z_E={\left[{\displaystyle \frac{5}{3{ϵ}_0}}-{\displaystyle \frac{1}{{ϵ}_0}}{W}_0\left(-{\displaystyle \frac{5r_1{e}^{5/3}}{3}}\right)\right]}^{-2/3}-1.$$

(62)

For small $r_1\ll 1$, Taylor expansion yields

$$z_E\simeq {\displaystyle \frac{(3-2r_1{e}^{5/3}{ϵ}_0^{2/3}}{{75}^{1/3}}}-1\simeq {\left(0.6{ϵ}_0\right)}^{2/3}-1$$

(63)

The peak SFR density evaluates to

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

(64)

where we used the expansion $W_0\left(\right|x|\ll 1)\simeq x$ and Equation (63). Equating (64) and (48) yields

$${\beta }_0=(1.{36}_{-0.78}^{+38.61}·{10}^{-18}{\left({\displaystyle \frac{1+r_1}{h_{70}}}\right)}^8\mathrm{Hz}.$$

(65)

This ${\beta }_0$-value2 is partially consistent with the lower limit (61) especially for large values of the observed peak SFR density (48). It implies according to Equations (52) and (63)

$${ϵ}_0=(0.{73}_{-0.42}^{+20.73}{\left({\displaystyle \frac{1+r_1}{h_{70}}}\right)}^9,z_E<4.49{\left({\displaystyle \frac{1+r_1}{h_{70}}}\right)}^6,$$

(66)

however, only if one adopts 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_1$ below the small value (56), the GS-model with spontaneous star formation indeed provides a rough but far from optimal fit to the observed SFR density (see Figure 3).

5.2. Stellar Density

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

$${\int }^{z}dz{\displaystyle \frac{e^{-\alpha {(1+z)}^{-3/2}}}{{(1+z)}^5}}={\displaystyle \frac{2}{3{\alpha }^{8/3}}}\Gamma \left[{\displaystyle \frac{8}{3}},{\displaystyle \frac{\alpha }{{(1+z)}^{3/2}}}\right]$$

(67)

the integrated stellar mass density (47) with $L\left(z\right)=0$ becomes

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

(68)

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

$${\rho }_{\mathrm{GS}}^{*}\left(0\right)={\displaystyle \frac{5.60·{10}^7{\beta }_0^2}{1+r_1}}(r_1+{\displaystyle \frac{8}{3}}{\displaystyle \frac{\gamma (\frac{8}{3},{ϵ}_0)}{{ϵ}_0^{8/3}}}\mathrm{kg}{\mathrm{m}}^{-3}$$

(69)

in a power-law like fashion at large redshifts $z\gg z_E$,

$${\rho }_{\mathrm{GSL}}^{*}\left(z\right)\simeq {\displaystyle \frac{B_1}{4{(1+z)}^4}},$$

(70)

since $\gamma (s,ϵ\ll 1)\simeq s^{-1}{ϵ}^s$ to lowest order in $ϵ$. Equating the result (69) with the observed value (50) provides

$${\displaystyle \frac{{\tilde{\beta }}_0^2}{1+r_1}}(r_1+{\displaystyle \frac{8}{3}}{\displaystyle \frac{\gamma (\frac{8}{3},{ϵ}_0)}{{ϵ}_0^{8/3}}}=(6.{79}_{-1.70}^{+2.16})·{10}^{-3},$$

(71)

where we scale

$${\beta }_0={10}^{-17}{\tilde{\beta }}_0$$

(72)

in units of ${10}^{-17}$ Hz. The function $\gamma (8/3,5.37{\tilde{\beta }}_0)$ is shown in Figure 4 and varies from 0 to the finite value $\Gamma (8/3)\approx 1.504$ at large arguments. For the nominal values $r_1=0$ and $h_{70}=1$ one obtains from Equation (71)

$${\displaystyle \frac{{\tilde{\beta }}_0^{2/3}}{\gamma (8/3,5.37{\tilde{\beta }}_0)}}\simeq {\displaystyle \frac{{\tilde{\beta }}_0^{2/3}}{1.504}}=4.{44}_{-1.07}^{+1.48},$$

(73)

which leads to

$${\tilde{\beta }}_0=17.{26}_{-5.85}^{+9.31}.$$

(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 least a factor 2.8. E.g. it is not possible to explain with one single ${\tilde{\beta }}_0$-value these two constraints with the GS-model.

Adopting as a best compromise value ${\tilde{\beta }}_0=2.0$ we compare in Figure 3 the SFR density and the integrated stellar densities for $r_1=0.001$. This choice implies as peak redshift $z_E=2.46$. Only because of the large error bars in the peak SFR density the theoretical results of the simplified GS-model are marginally consistent with the observations. This Figure 3 also illustrates the basic dilemma of this simplified model very well. Because of the dependencies ${\psi }_{\mathrm{GS}}\left(z_E\right)\propto {\tilde{\beta }}_0^{1/3}$ in Equation (64) and ${\rho }_{\mathrm{GS}}^{*}\left(0\right)\propto {\tilde{\beta }}_0^{-2/3}$ in Equation (69) on the only free parameter ${\tilde{\beta }}_0$, it is obvious that by decreasing ${\tilde{\beta }}_0$ from its compromise value 2 would decrease ${\psi }_{\mathrm{GS}}\left(z\right)$ to smaller values, but would increase at the same time ${\rho }_{\mathrm{GS}}^{*}\left(z\right)$. Alternatively, increasing ${\tilde{\beta }}_0$ from its compromise value 2 would decrease ${\rho }_{\mathrm{GS}}^{*}\left(z\right)$ to smaller values, but at the same time would make ${\psi }_{\mathrm{GS}}\left(z\right)$ even larger.

With this compromise value ${\tilde{\beta }}_0=2.0$ and the ratio $r_1=0.04$ the present-day gas fraction is

$$G(r_2=0,z=0)=3.85·{10}^{-2},$$

(75)

basically identical to the gas fraction at infinite time $G_{\infty }$, and in marginal agreement with the observation (31).

5.3. Special Case: Negligible Stellar Feedback ($r_1=0$) and Stellar Evolution ($r_2=0$)

The results of the last subsections include as special case the GS-model without stellar feedback ($r_1=0$) and stellar evolution ($r_2=0$). 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) simplifies to

$${\psi }_{\mathrm{GS}}\left(z\right)=A_1{\displaystyle \frac{e^{-ϵ{(1+z)}^{-3/2}}}{{(1+z)}^{5/2}}}$$

(76)

with

$$ϵ={ϵ}_0(r_1=0)=5.37{\tilde{\beta }}_0{h}_{70}^{-1}.$$

(77)

The density (76) attains its maximum value

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

(78)

at the peak redshift

$$z_E={\left(0.6ϵ\right)}^{2/3}-1=2.18{\left({\displaystyle \frac{{\tilde{\beta }}_0}{h_{70}}}\right)}^{2/3}-1$$

(79)

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

$${\tilde{\beta }}_0=(1.{61}_{-0.73}^{+0.87})h_{70}$$

(80)

and

$${\tilde{\beta }}_0=(0.{14}_{-0.08}^{+3.84})h_{70}^{-8},$$

(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}{ccc}\hfill {\rho }_{\mathrm{GS}}^{*}\left(z\right)& =& B_1{\int }_z^{1100}d{z}^{\prime }{\displaystyle \frac{e^{ϵ{(1+z^{\prime })}^{-3/2}}}{{(1+z^{\prime })}^5}}\simeq {\displaystyle \frac{2B_1}{3{ϵ}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}})\hfill \\ & \simeq & {\displaystyle \frac{1.49·{10}^{-26}{\tilde{\beta }}_0^2}{{ϵ}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}})\mathrm{kg}{\mathrm{m}}^{-3}\hfill \\ & =& {\displaystyle \frac{1.69·{10}^{-28}{h}_{70}^{8/3}}{{\tilde{\beta }}_0^{2/3}}}\gamma \left({\displaystyle \frac{8}{3}},{\displaystyle \frac{5.37{\tilde{\beta }}_0}{h_{70}{(1+z)}^{3/2}}}\right)\mathrm{kg}{\mathrm{m}}^{-3}.\hfill \end{array}$$

(82)

The stellar mass density thus monotonically decreases from the value

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{GS}}^{*}\left(0\right)& =& {\displaystyle \frac{1.69·{10}^{-28}{h}_{70}^{8/3}}{{\tilde{\beta }}_0^{2/3}}}\gamma \left({\displaystyle \frac{8}{3}},5.37{\displaystyle \frac{{\tilde{\beta }}_0}{h_{70}}}\right)\mathrm{kg}{\mathrm{m}}^{-3}\hfill \\ & =& {\displaystyle \frac{2.54·{10}^{-28}{h}_{70}^{8/3}}{{\tilde{\beta }}_0^{2/3}}}\mathrm{kg}{\mathrm{m}}^{-3}\hfill \end{array}$$

(83)

in a power-law like fashion at large redshifts $z\gg z_E$,

$${\rho }_{\mathrm{GSL}}^{*}\left(z\right)\simeq {\displaystyle \frac{B_1}{4{(1+z)}^4}},$$

(84)

since $\gamma (8/3,x)=(3/8)x^{8/3}+\mathcal{O}(x^{10/3}$. The estimate (74) for ${\tilde{\beta }}_0$ from the result (83) with the observed value (50) for the nominal values and $h_{70}=1$ remains unchanged, i.e.

$${\tilde{\beta }}_0=17.{26}_{-5.85}^{+9.31}.$$

(85)

The two determinations by Equations (80)–(81) and Equation (85) of the rate of spontaneous star formation, resulting from reproducing the peak redshift and peak value of the SFR density and the present-day integrated stellar mass density within their partially large uncertainties, exclude each other. E.g. it is not possible to explain with one single ${\tilde{\beta }}_0$-value these three constraints with the simplified GS-model. Again a compromise value of ${\tilde{\beta }}_0=2.0$ is suggested. With this compromise value ${\tilde{\beta }}_0=2.0$ and the present-day gas fraction in the case of neglected stellar feedback is

$$G(r_1=0,r_2=0,z=0)=e^{-10.74}=2.17·{10}^{-5},$$

(86)

another strong disagreement with the observational constraint (31) of this model.

6. SFR and Stellar Density for Spontaneous Star Formation with Stellar Feedback

Here we extend the analysis of the previous section by including additionally stellar feedback ($b_0\ne 0$) besides spontaneous star formation (${\beta }_0\ne 0$) and stellar evolution $c_0\ne 0$ by using (19), instead of (29), for $j\left(\tilde{t}\right)$.

6.1. Theoretical Results

With $\tilde{t}\left(z\right)$ from Equation (42) and Equation (19) for $j\left(\tilde{t}\right)$ the SFR density (46) becomes

$$\begin{array}{ccc}\hfill {\psi }_{\mathrm{GSL}}\left(z\right)& =& {\displaystyle \frac{A_1}{{(1+z)}^{5/2}}}{e}^{-\mu \tilde{t}}[cosh\left(\nu \tilde{t}\right)-{\displaystyle \frac{1-\mu }{\nu }}sinh\left(\nu \tilde{t}\right)]\hfill \\ & =& {\displaystyle \frac{A_1}{{ϵ}^{5/3}}}F\left(\tilde{t}\right),\hfill \end{array}$$

(87)

with $\mu =(1+r_1+r_2)/2$, $\nu =\sqrt{{\mu }^2-r_2}$, $ϵ={\beta }_0{\tau }_f=5.37{\tilde{\beta }}_0{h}_{70}^{-1}$, as before in Equation (77), and the function3

$$F\left(\tilde{t}\right)={\tilde{t}}^{5/3}{e}^{-\mu \tilde{t}}[cosh\left(\nu \tilde{t}\right)+{\displaystyle \frac{\mu -1}{\nu }}sinh\left(\nu \tilde{t}\right)].$$

(88)

Likewise with Equation (21) one obtains

$$1-L\left(z\right)=e^{-\mu \tilde{t}}\left[cosh\nu \tilde{t}+{\displaystyle \frac{\mu }{\nu }}sinh\nu \tilde{t}\right]$$

(89)

with $\tilde{t}\left(z\right)=ϵ{(1+z)}^{-3/2}$. Consequently, the integrated stellar density (47) is given by

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{GSL}}^{*}\left(z\right)& =& B_1{\int }_z^{1100}d{z}^{\prime }{\displaystyle \frac{e^{-2\mu \tilde{t}\left(z^{\prime }\right)}}{{(1+z^{\prime })}^5}}[cosh\nu \tilde{t}\left(z^{\prime }\right)\hfill \\ & & +{\displaystyle \frac{\mu }{\nu }}sinh\nu \tilde{t}\left(z^{\prime }\right)][cosh\nu \tilde{t}\left(z^{\prime }\right)+{\displaystyle \frac{\mu -1}{\nu }}sinh\nu \tilde{t}\left(z^{\prime }\right)]\hfill \\ & \simeq & {\displaystyle \frac{B_1}{6{\nu }^2{ϵ}^{8/3}}}{\int }_0^{{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}}}d\tilde{t}{\tilde{t}}^{5/3}{e}^{-2\mu \tilde{t}}[(\nu +\mu )e^{\nu \tilde{t}}+(\nu -\mu )e^{-\nu \tilde{t}}]\hfill \\ & \times & [(\nu +\mu -1)e^{\nu \tilde{t}}+(\nu +1-\mu )e^{-\nu \tilde{t}}]\hfill \\ & =& {\displaystyle \frac{B_1}{6{\nu }^2{ϵ}^{8/3}}}{\int }_0^{{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}}}d\tilde{t}{\tilde{t}}^{5/3}{e}^{-2\mu \tilde{t}}[2({\nu }^2+\mu (1-\mu ))\hfill \\ & & ++(\nu +\mu )(\nu +\mu -1)e^{2\nu \tilde{t}}+(\nu -\mu )(\nu +1-\mu )e^{-2\nu \tilde{t}}]\hfill \\ & =& {\displaystyle \frac{B_1}{6{\nu }^2{\left(2ϵ\right)}^{8/3}}}[{\displaystyle \frac{2({\nu }^2+\mu (1-\mu ))}{{\mu }^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{2\mu ϵ}{{(1+z)}^{3/2}}}\hfill \\ & & ++{\displaystyle \frac{(\nu +\mu )(\nu +\mu -1)}{{(\mu -\nu )}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{2(\mu -\nu )ϵ}{{(1+z)}^{3/2}}})\hfill \\ & & +{\displaystyle \frac{(\nu -\mu )(\nu +1-\mu )}{{(\mu +\nu )}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{2(\mu +\nu )ϵ}{{(1+z)}^{3/2}}}].\hfill \end{array}$$

(90)

Moreover, the gas fraction (19) as a function of redshift reads

$$G\left(z\right)={\displaystyle \frac{1}{2\nu }}[(\nu +\mu -1)e^{{\displaystyle \frac{(\nu -\mu )ϵ}{{(1+z)}^{3/2}}}}+(\nu +1-\mu )e^{-{\displaystyle \frac{(\nu +\mu )ϵ}{{(1+z)}^{3/2}}}}].$$

(91)

The first derivative of the function (88) is given by

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dF\left(\tilde{t}\right)}{d\tilde{t}}}& =& {\displaystyle \frac{{\tilde{t}}^{2/3}{e}^{-\mu \tilde{t}}cosh\nu \tilde{t}}{\nu }}\{{\displaystyle \frac{5}{3}}\left[\nu +(\mu -1)tanh\left(\nu \tilde{t}\right)\right]\hfill \\ & & +\tilde{t}\left[(\mu -r_2)tanh\left(\nu \tilde{t}\right)-\nu \right]\},\hfill \end{array}$$

(92)

so that a maximum occurs at ${\tilde{t}}_E$ given by the solution of the transcendental equation

$$\left[\nu -(\mu -r_2)tanh\left(\nu {\tilde{t}}_E\right)\right]{\tilde{t}}_E={\displaystyle \frac{5}{3}}\left[\nu +(\mu -1)tanh\left(\nu {\tilde{t}}_E\right)\right].$$

(93)

As an aside we note that for $r_1=r_2=0$ the Equation (93) correctly reduces to ${\tilde{t}}_E(r_1=0,r_2=0)=5/3$ in agreement with Section 5.3. Likewise for $r_2=0$ one obtains $\nu =\mu =(1+r_1)/2$ and with $(1+r_1){\tilde{t}}_E=Z$, with Z defined in Equation (53), Equation (93) simplifies to

$$0.6Z\left[1-tanh{\displaystyle \frac{Z}{2}}\right]=1+r_1+(r_1-1)tanh{\displaystyle \frac{Z}{2}},$$

(94)

reproducing with $tanh(Z/2)=(1-e^{-Z}/(1+e^{-Z}$ correctly the determining Equation (54) for the extrema in the case of negligible stellar feedback.

6.2. Constraints on Parameters

First, the observed present-day integrated stellar mass density (50) provides from Equation (90)

$$\begin{array}{ccc}\hfill (5.{715}_{-1.369}^{+1.820}{\displaystyle \frac{\nu {\tilde{\beta }}_0^{2/3}}{h_{70}^{8/3}}})& =& {\displaystyle \frac{2({\nu }^2+\mu (1-\mu ))}{{\mu }^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{10.74\mu {\tilde{\beta }}_0}{h_{70}}})\hfill \\ & & +{\displaystyle \frac{(\nu +\mu )(\nu +\mu -1)}{{(\mu -\nu )}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{10.74(\mu -\nu ){\tilde{\beta }}_0}{h_{70}}})\hfill \\ & & +{\displaystyle \frac{(\nu -\mu )(\nu +1-\mu )}{{(\mu +\nu )}^{8/3}}}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{10.74(\mu +\nu ){\tilde{\beta }}_0}{h_{70}}}).\hfill \end{array}$$

(95)

Adopting values of ${\tilde{\beta }}_0/h_{70}>1$ and assuming $|\mu -\nu |>0.5$, all incomplete $\gamma$-functions are well approximated by $1.504$, so that the constraint (95) simplifies to

$$\begin{array}{ccc}\hfill (4.{421}_{-0.910}^{+1.230}\nu {\displaystyle \frac{{\tilde{\beta }}_0^{2/3}}{h_{70}^{8/3}}}& =& {\displaystyle \frac{2({\nu }^2+\mu (1-\mu ))}{{\mu }^{8/3}}}+{\displaystyle \frac{(\nu +\mu )(\nu +\mu -1)}{{(\mu -\nu )}^{8/3}}}\hfill \\ & & +{\displaystyle \frac{(\nu -\mu )(\nu +1-\mu )}{{(\mu +\nu )}^{8/3}}}.\hfill \end{array}$$

(96)

Secondly, the observed peak redshift (49) yields for

$${\tilde{t}}_E={\displaystyle \frac{ϵ}{{(1+z_E)}^{3/2}}}=(1.{03}_{-0.36}^{+0.87}){\tilde{\beta }}_0{h}_{70}^{-1}.$$

(97)

As third constraint we infer from Equations (48) and (87)

$$(3.{83}_{-0.95}^{+8.00})·{10}^{-46}=1.69·{10}^{-45}{\tilde{\beta }}_0^{1/3}{h}_{70}^{8/3}F\left({\tilde{t}}_E\right)$$

(98)

or

$$\begin{array}{ccc}& & 0.{227}_{-0.056}^{+0.473}={\tilde{\beta }}_0^{1/3}{h}_{70}^{8/3}F\left({\tilde{t}}_E\right)=[{(1.{03}_{-0.36}^{+0.87}]}^{5/3}{\tilde{\beta }}_0^2{h}_{70}\hfill \\ & & \times exp[-(1.{03}_{-0.36}^{+0.87})\mu {\tilde{\beta }}_0{h}_{70}^{-1}](cosh\left((1.{03}_{-0.36}^{+0.87})\nu {\tilde{\beta }}_0{h}_{70}^{-1}\right)\hfill \\ & & +{\displaystyle \frac{\mu -1}{\nu }}sinh((1.{03}_{-0.36}^{+0.87})\nu {\tilde{\beta }}_0{h}_{70}^{-1})\hfill \end{array}$$

(99)

with

$$F\left({\tilde{t}}_E\right)={\tilde{t}}_E^{5/3}{e}^{-\mu {\tilde{t}}_E}[cosh\left(\nu {\tilde{t}}_E\right)+{\displaystyle \frac{\mu -1}{\nu }}sinh\left(\nu {\tilde{t}}_E\right)]$$

(100)

according to Equation (88), and where we inserted Equation (97).

As fourth constraint we demand that ${(d\psi /dz)}_{z=0}>0$ in order to have a maximum of $\psi \left(z\right)$. With

$${\displaystyle \frac{d{\psi }_{\mathrm{GSL}}\left(z\right)}{dz}}=-{\displaystyle \frac{2A_1}{3ϵ{\tilde{t}}^{5/3}}}{\displaystyle \frac{dF\left(\tilde{t}\right)}{d\tilde{t}}}$$

(101)

and Equation (92) one has to demand that

$$\nu +(\mu -1)tanh\left(\nu ϵ\right)<{\displaystyle \frac{3ϵ}{5}}[\nu -(\mu -r_2)tanh\left(\nu ϵ\right)],$$

(102)

leading, after insertion of Equation (77), to

$$\begin{array}{ccc}\hfill tanh\left({\displaystyle \frac{5.37\nu {\tilde{\beta }}_0}{h_{70}}}\right)& <& {\displaystyle \frac{\nu \left(3.22\frac{{\tilde{\beta }}_0}{h_{70}}-1\right)}{\mu -1+3.22\frac{{\tilde{\beta }}_0}{h_{70}}(\mu -r_2)}}.\hfill \end{array}$$

(103)

6.2.1. Results for standard Hubble constant value $h_{70}=1$

In Figure 5 we calculate for the wide range of values of ${\tilde{\beta }}_0\in [0.5,2.7]$ the compatibility of the five constraints on ${\psi }_{\mathrm{GSL}}\left(z_E\right)$, $z_E$, ${\rho }_{\mathrm{GSL}}^{*}\left(0\right)$, and ${d\psi \left(z\right)/dz|}_{z=0}$ in the $r_1$–$r_2$-parameter plane over 9 decades of values each. Obviously, the tightest constraints on $r_1$ and $r_2$ are provided by the constraints on $z_E$ and ${\rho }_{\mathrm{GSL}}^{*}\left(0\right)$. The green regions on the rightest panel (labeled summary) for the values of ${\tilde{\beta }}_0\in [0.9,2.7]$ indicates the range of parameters $r_1$ and $r_2$ where all five constraints are fulfilled. parameter-range Although there is no set of fit parameters $({\tilde{\beta }}_0,r_1,r_2)$, where all four parameter constraints are fulfilled. This is certainly a significant improvement of the GSL-model with spontaneous star formation over the simplified models considered in Section 5.

To illustrate this improvement we compare in Figure for six different values of ${\tilde{\beta }}_0$ the redshift dependencies of the SFR density and the integrated stellar density, for all pairs of ($r_1,r_2$) values residing either within the green (all five constraints fulfilled for ${\tilde{\beta }}_0=0.9,1.0,1.4,2.0$) within the blue (three constraints fulfilled for ${\tilde{\beta }}_0=0.5,0.7$) areas of the corresponding panels in Figure 5. Although the agreement is not perfect for any $({\tilde{\beta }}_0,r_1,r_2)$ triplet, it is certainly an improvement as compared to the model ignoring stellar feedback in Section 5. Even with three free parameters it is still not possible to reproduce the observations nearly perfectly with the GSL-model for spontaneous star formation only. It defines a challenge for future work to demonstrate that the inclusion of the triggered star formation process can remedy the situation significantly.

For each of these four ${\tilde{\beta }}_0$ values we also determine the ranges of the present-day gas fraction from Equation (91)

$$G\left(0\right)={\displaystyle \frac{1}{2\nu }}[(\nu +\mu -1)e^{(\nu -\mu )ϵ}+(\nu +1-\mu )e^{-(\nu +\mu )ϵ}].$$

(104)

Using all admissible ($r_1,r_2$) pairs we obtain $G\left(0\right)\in [0.007,0.014]$ for ${\tilde{\beta }}_0=0.9$, $G\left(0\right)\in [0.0047,0.0205]$ for ${\tilde{\beta }}_0=1$, $G\left(0\right)\in [2.2\times {10}^{-5},0.0077]$ for ${\tilde{\beta }}_0=2$, and $G\left(0\right)\in [1.48\times {10}^{-6},0.0041]$ for ${\tilde{\beta }}_0=2.5$. All these values are within precision below about 2 percent.

Figure 5. For each ${\tilde{\beta }}_0$ value there is one line. Each plot versus $r_1$ and $r_2$ (double-logarithmic). From left to right: Constraint 1: ${\psi }_{\mathrm{GSL}}\left(z_E\right)=3.{83}_{-0.95}^{+8}·{10}^{-46}$ kg ${\mathrm{s}}^{-1}$ ${\mathrm{m}}^{-3}$ according to Equation (48), constraint 2: $z_E=2_{-1}^{+1}$ (49), constraint 3: ${\rho }_{\mathrm{GSL}}^{*}\left(0\right)=3.{80}_{-0.91}^{+1.21}·{10}^{-29}$ kg ${\mathrm{m}}^{-3}$ (50), constraint 4: the sign of ${\psi }_{\mathrm{GSL}}^{\prime }\left(0\right)$, constraint 5: the existence of a single extremum in ${\psi }_{\mathrm{GSL}}\left(z\right)$ for positive z, and the summary in the last column. The gray-shaded region highlights the parameter regime for which a quantity is within the allowed range, or ${\psi }_{\mathrm{GSL}}(z=0)>0$. The summary displays the amount of simultaneously fulfilled criteria. For ${\tilde{\beta }}_0\in [0.9,2.7]$ all five criteria are fulfilled over ranges of $r_1$ and $r_2$ values.

Figure 5. For each ${\tilde{\beta }}_0$ value there is one line. Each plot versus $r_1$ and $r_2$ (double-logarithmic). From left to right: Constraint 1: ${\psi }_{\mathrm{GSL}}\left(z_E\right)=3.{83}_{-0.95}^{+8}·{10}^{-46}$ kg ${\mathrm{s}}^{-1}$ ${\mathrm{m}}^{-3}$ according to Equation (48), constraint 2: $z_E=2_{-1}^{+1}$ (49), constraint 3: ${\rho }_{\mathrm{GSL}}^{*}\left(0\right)=3.{80}_{-0.91}^{+1.21}·{10}^{-29}$ kg ${\mathrm{m}}^{-3}$ (50), constraint 4: the sign of ${\psi }_{\mathrm{GSL}}^{\prime }\left(0\right)$, constraint 5: the existence of a single extremum in ${\psi }_{\mathrm{GSL}}\left(z\right)$ for positive z, and the summary in the last column. The gray-shaded region highlights the parameter regime for which a quantity is within the allowed range, or ${\psi }_{\mathrm{GSL}}(z=0)>0$. The summary displays the amount of simultaneously fulfilled criteria. For ${\tilde{\beta }}_0\in [0.9,2.7]$ all five criteria are fulfilled over ranges of $r_1$ and $r_2$ values.

Figure 6. (Left column) Experimentally reported $\psi \left(z\right)$ (symbols) and ${\psi }_{\mathrm{GSL}}\left(z\right)$ according to Equation (87) (family of gray lines) as well as (Right column) reported ${\rho }_{*}\left(z\right)$ (symbols) and ${\rho }_{\mathrm{GSL}}^{*}\left(z\right)$ according to Equation (90) (family of gray lines) for $h_{70}=1$ and various ${\tilde{\beta }}_0$ (increasing from top to bottom). The families are created using 1000 randomly chosen, but eligible ($r_1,r_2$) pairs for given ${\tilde{\beta }}_0$ (specified in the panels). The eligible pairs reside in the green regions of the corresponding panel in Fig. Figure 5. Notice that for ${\tilde{\beta }}_0=0.5$ and ${\tilde{\beta }}_0=0.7$ no green regions exist: i.e. not all constraints can be fulfilled. For those two values, the eligible pairs are from the blue regions.

Figure 6. (Left column) Experimentally reported $\psi \left(z\right)$ (symbols) and ${\psi }_{\mathrm{GSL}}\left(z\right)$ according to Equation (87) (family of gray lines) as well as (Right column) reported ${\rho }_{*}\left(z\right)$ (symbols) and ${\rho }_{\mathrm{GSL}}^{*}\left(z\right)$ according to Equation (90) (family of gray lines) for $h_{70}=1$ and various ${\tilde{\beta }}_0$ (increasing from top to bottom). The families are created using 1000 randomly chosen, but eligible ($r_1,r_2$) pairs for given ${\tilde{\beta }}_0$ (specified in the panels). The eligible pairs reside in the green regions of the corresponding panel in Fig. Figure 5. Notice that for ${\tilde{\beta }}_0=0.5$ and ${\tilde{\beta }}_0=0.7$ no green regions exist: i.e. not all constraints can be fulfilled. For those two values, the eligible pairs are from the blue regions.

6.2.2. Smaller Values of the Hubble Constant

We have noticed before the strong dependency of the constraints on the adopted value of the Hubble constant with ${\psi }_{\mathrm{GLS}}\left(z_E\right)\propto h_{70}^{8/3}$, ${\rho }_{\mathrm{GLS}}^{*}\left(0\right)\propto h_{70}^{8/3}$ and $1+z_E\propto h_{70}^{-2/3}$. Lowering the values of the Hubble constant to values of $H_0=63$ km ${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}$ or even of $H_0=60$ km ${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}$ indeed provides better agreements of our theoretical results in Figure with the observations as in this case ${\psi }_{\mathrm{GLS}}\left(z_E\right)$ and ${\rho }_{\mathrm{GLS}}\left(0\right)$ are reduced by the factors 0.755 and 0.663, respectively, whereas the peak redshift increases by the factors 1.07 and 1.11, respectively. In Figure 7 we show the comparison for $h_{70}=0.857$.

Figure 7. Same as Figure 6 for $H_0=60$ km s${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}$ (corresponding to $h_{70}=0.857$) for ${\tilde{\beta }}_0=0.8$, $0.9$, and $1.0$ (from top to bottom).

Figure 7. Same as Figure 6 for $H_0=60$ km s${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}$ (corresponding to $h_{70}=0.857$) for ${\tilde{\beta }}_0=0.8$, $0.9$, and $1.0$ (from top to bottom).

6.3. Special Case of Negligible Stellar Feedback

The general investigation of the case with stellar evolution ($r_2>0$) has indicated that the special case of negligible stellar evolution ($r_1=0$) also meets all constraints. For completeness we consider this special case here in more detail. We readily find from Equations (88) and (89) that

$$\begin{array}{ccc}\hfill F\left(\tilde{t}\right)& =& {\tilde{t}}^{5/3}{e}^{-\tilde{t}}\hfill \\ \hfill 1-L\left(\tilde{t}\right)& =& {\displaystyle \frac{e^{-r_2\tilde{t}}-r_2{e}^{-\tilde{t}}}{1-r_2}}\hfill \end{array}$$

(105)

which also hold for $r_2=1$, where $L(r_1=0,r_2=1,\tilde{t})=(1+\tilde{t})e^{-\tilde{t}}$. The function $F\left(\tilde{t}\right)$ has a single maximum at ${\tilde{t}}_E=5/3$, so that the SFR density ${\psi }_{\mathrm{GLS}}\left(z\right)$ attains its maximum

$${\psi }_{\mathrm{GLS}}\left(z_E\right)={\left({\displaystyle \frac{5}{3ϵ}}\right)}^{5/3}{A}_1{e}^{-{\displaystyle \frac{5}{3}}}=7.47·{10}^{-46}{\tilde{\beta }}_0^{1/3}{h}_{70}^{8/3}\mathrm{kg}{\mathrm{s}}^{-1}{\mathrm{m}}^{-3}$$

(106)

at the peak redshift

$$1+z_E={\left(0.6ϵ\right)}^{2/3}=2.18{\tilde{\beta }}_0^{2/3}{h}_{70}^{-2/3}.$$

(107)

Consequently,

$${\displaystyle \frac{{\psi }_{\mathrm{GLS}}\left(z\right)}{{\psi }_{\mathrm{GLS}}\left(z_E\right)}}={\left({\displaystyle \frac{1+z_E}{1+z}}\right)}^{5/2}{e}^{-ϵ[{(1+z)}^{-3/2}-{(1+z_E)}^{-3/2}]}.$$

(108)

Likewise with $\mu =(1+r_2)/2$ and $\nu =|1-r_2|/2$ in this case Equation (90) simplifies to

$$\begin{gathered} {\rho }_{\mathrm{GSL}}^{*}(z,r_1=0,r_2)={\displaystyle \frac{2B_1}{3(1-r_2){\left(2ϵ\right)}^{8/3}}} \\ \times \left[{\left({\displaystyle \frac{2}{1+r_2}}\right)}^{8/3}\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{(1+r_2)ϵ}{{(1+z)}^{3/2}}})-r_2\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{2ϵ}{{(1+z)}^{3/2}}})\right]. \end{gathered}$$

(109)

and in the particular case $r_2=1$

$$\begin{gathered} {\rho }_{\mathrm{GSL}}^{*}(z,r_1=0,r_2=1)\simeq {\displaystyle \frac{2B_1}{3{ϵ}^{8/3}}}{\int }_0^{{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}}}d\tilde{t}{\tilde{t}}^{5/3}(1+\tilde{t})e^{-2\tilde{t}} \\ ={\displaystyle \frac{B_1}{3{\left(2ϵ\right)}^{8/3}}}\left[2\gamma ({\displaystyle \frac{8}{3}},{\displaystyle \frac{2ϵ}{{(1+z)}^{3/2}}}+\gamma ({\displaystyle \frac{11}{3}},{\displaystyle \frac{2ϵ}{{(1+z)}^{3/2}}})\right], \end{gathered}$$

(110)

while Equations (109) and (110) imply for $r_2\ne 1$

$$\begin{array}{ccc}& & {\rho }_{\mathrm{GSL}}^{*}(z=0,r_1=0)={\displaystyle \frac{2.660·{10}^{-29}{h}_{70}^{8/3}}{(1-r_2){\tilde{\beta }}_0^{2/3}}}\times \hfill \\ & & \left[{\left({\displaystyle \frac{2}{1+r_2}}\right)}^{8/3}\gamma ({\displaystyle \frac{8}{3}},5.37(1+r_2))-r_2\gamma ({\displaystyle \frac{8}{3}},10.74)\right]\hfill \\ & \simeq & {\displaystyle \frac{4.00·{10}^{-29}[{\left(\frac{2}{1+r_2}\right)}^{8/3}-r_2]h_{70}^{8/3}}{(1-r_2){\tilde{\beta }}_0^{2/3}}}\mathrm{kg}{\mathrm{m}}^{-3},\hfill \end{array}$$

(111)

whereas for $r_2=1$

$${\rho }_{\mathrm{GSL}}^{*}(0,r_1=0,r_2=1))\simeq 9.31·{10}^{-29}{h}_{70}^{8/3}{\tilde{\beta }}_0^{-2/3}\mathrm{kg}{\mathrm{m}}^{-3},$$

(112)

where we approximated both incomplete gamma functions by $\Gamma (8/3)=1.504$ and $\Gamma (11/3)=(8/3)\Gamma (8/3)$. Note that with the limit

$$\underset{r_2\to 1}{lim}{\displaystyle \frac{[{\left(\frac{2}{1+r_2}\right)}^{8/3}-r_2]}{1-r_2}}=\underset{x\to 0}{lim}\left[{\displaystyle \frac{{\left(\frac{2}{2-x}\right)}^{8/3}+x-1}{x}}\right]=7/3$$

(113)

one can infer directly Equation (112) from Equation (111).

The three constraints (48)–(50) then provide

$${\tilde{\beta }}_0=(1.{61}_{-0.73}^{+0.88})h_{70},$$

(114)

$${\tilde{\beta }}_0=(0.{135}_{-0.78}^{+3.836})h_{70}^{-8},$$

(115)

and

$${\tilde{\beta }}_0=(1.{080}_{-0.367}^{+0.548}){\left[{\displaystyle \frac{{\left(\frac{2}{1+r_2}\right)}^{8/3}-r_2}{1-r_2}}\right]}^{3/2}{h}_{70}^4$$

(116)

Whereas the first two constraints (114) and (115) restrict the value of ${\tilde{\beta }}_0$, we illustrate in Figure 8 the allowed range of values values of the parameters ${\tilde{\beta }}_0$ and $r_2$ in the case $r_1=0$ provided by constraint (116) from the observed ${\rho }^{*}\left(0\right)$.

For $r_1=0$, in the limit $r_2\to \infty$

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{GSL}}\left(z\right)& =& {\displaystyle \frac{2B_1}{3{ϵ}^{8/3}}}{\int }_0^{{\displaystyle \frac{ϵ}{{(1+z)}^{3/2}}}}d\tilde{t}{\tilde{t}}^{5/3}{e}^{-2\tilde{t}},\hfill \\ & =& {\displaystyle \frac{2B_1}{3{\left(2ϵ\right)}^{8/3}}}\gamma \left({\displaystyle \frac{8}{3}},{\displaystyle \frac{2ϵ}{{(1+z)}^{3/2}}}\right),\hfill \end{array}$$

(117)

so that ${\rho }_{\mathrm{GSL}}^{*}\left(0\right)=4.0·{10}^{-29}{h}_{70}^{8/3}{\tilde{\beta }}_0^{-2/3}$ kg ${\mathrm{m}}^{-3}$. This latter expression agrees with the values for ${\rho }_{\mathrm{GSL}}\left(0\right)$ shown in Figure 8 at $r_2={10}^3$.

Using Equations (25)–(27) and Equation (42) one obtains in the case of negligible stellar feedback for the redshift dependencies of the three fractions

$$\begin{array}{ccc}\hfill G\left(z\right)& =& e^{{\displaystyle \frac{-5.37{\tilde{\beta }}_0{h}_{70}^{-1}}{{(1+z)}^{3/2}}}},\hfill \\ \hfill S\left(z\right)& =& {\displaystyle \frac{e^{\frac{-5.37{\tilde{\beta }}_0{h}_{70}^{-1}}{{(1+z)}^{3/2}}}-e^{\frac{-5.37r_2{\tilde{\beta }}_0{h}_{70}^{-1}}{{(1+z)}^{3/2}}}}{r_2-1}},\hfill \\ \hfill L\left(z\right)& =& 1-{\displaystyle \frac{r_2{e}^{\frac{-5.37{\tilde{\beta }}_0{h}_{70}^{-1}}{{(1+z)}^{3/2}}}-e^{\frac{-5.37r_2{\tilde{\beta }}_0{h}_{70}^{-1}}{{(1+z)}^{3/2}}}}{r_2-1}}.\hfill \end{array}$$

(118)

Adopting for $h_{70}=1$ the value ${\tilde{\beta }}_0=1.2$ as best choice from constraints (114) and (115) then provides $r_2=6.7$ from the third constraint (116). In Figure 9 we show the redshift dependencies of the three fractions (118) for the parameter set $({\tilde{\beta }}_0=1.2,r_1=0,r_2=6.7)$. For the present-day values at $z=0$ we then find $G(z=0)=4.67·{10}^{-3}$, $S(z=0)=8.19·{10}^{-4}$. The vast majority (more than 99.4 percent) of the baryons in the the present universe reside in form of locked-in stellar matter in white dwarfs, neutron stars and black holes.

7. Summary and Conclusions

The compartmental description, well-known from the description of infection diseases and epidemics, is applied to describe the temporal evolution of the baryonic matter in interstellar gas and stars. The introduction of gaseous and stellar fractions of the total baryonic matter as the basic dynamical variables is advantageous because it allows to apply the description to a variety of astrophysical systems.

In this first paper of a series the competition of spontaneous star formation, stellar feedback and stellar evolution is theoretically investigated in order to understand the baryonic matter cycle whereas the inclusion of also the triggered star formation process will be the subject of the second paper of this series. Luminous baryonic matter occurs as interstellar and intergalactic gas with the fraction $G\left(t\right)$ as well as in main-sequence stars with the fraction $S\left(t\right)$. The third compartment with fraction $L\left(t\right)$ denotes weakly luminous matter in white dwarfs, neutron stars and black holes (referred to as locked-in matter) which have no significant stellar feedback to the gaseous matter compartment. The temporal evolution of the three fractions are controlled by the respective rates of spontaneous star formation ($\beta \left(t\right)G\left(t\right)$), of stellar feedback ($b\left(t\right)S\left(t\right)$) of stellar to gaseous matter, and of the formation ($c\left(t\right)S\left(t\right)$) of white dwarfs, neutron stars and black holes from stellar evolution.

By introducing the dimensionless reduced time variable $\tilde{t}={\int }_{t_0}^{t}d\xi \beta \left(\xi \right)$ for arbitrarily but given time-dependent spontaneous SFR coefficient $\beta \left(t\right)$, as well as the ratios $r_1=b\left(t\right)/\beta \left(t\right)$ and $r_2=c\left(t\right)/\beta \left(t\right)$, the derived exact solutions of the dynamical equations hold for stationary rates as well as for the case of the same time dependence of all rates. The accuracy of the analytical solutions is proven by the favorite comparison with the exact numerical solutions of the dynamical equations.

Of particular interest is the understanding of the cosmic star formation history and the present-day gas fraction with compartmental models. For a flat $\Lambda$CDM Friedmann cosmology the relationship between the reduced time variable $\tilde{t}\left(z\right)$ and the cosmological redshift z is used to calculate the redshift dependence of the cosmological star formation rate, the integrated stellar density and the present-day gas fraction from the derived gaseous fraction determining the formation rate of news stars.

The comparison with the observed cosmological star formation rate and the integrated stellar density indicates that the simplified GS-model ignoring stellar evolution cannot explain the observations reasonably well with one single stationary spontaneous SFR coefficient ${\beta }_0$. The best compromise value ${\beta }_0=2.0·{10}^{-17}$ Hz implies a present-day gas fraction in the universe of about $3.85$ percent.

However, the situation is considerably improved in the full GSL-model including stellar evolution, although the agreement with observations is not perfect. Here the observed cosmological star formation rate and the integrated stellar density as a function of redshift are reasonably well explained by the compartmental model without triggered star formation by the competition of spontaneous star formation and stellar evolution whereas the influence of stellar feedback is less important. The action of stellar evolution provides a significant redshift dependent reduction factor when calculating the integrated stellar density from the star formation rate. Then the fits to the observation allow us conclusions on the relative importance of spontaneous star formation, stellar evolution and feedback in the early universe after the recombination era until today. Lowering the value of the Hubble constant to smaller than the nominal one ($h_{70}=1$) improves the agreement of the GSL-Model with the observations. The present-day gas, luminous star and locked-in stellar matter fractions indicate that the vast majority of the baryons (more than 99.4 percent) resides in the form of locked-in stellar matter in white dwarfs, neutron stars and black holes. These objects are the astrophysical sites for any hidden baryons.

The non-perfect agreement of the GSL-model with only spontaneous star formation process defines a challenge for future work to demonstrate that the inclusion of the triggered star formation process can remedy the situation significantly. In conclusion, this work has demonstrated that compartmental models of the type introduced here lead to new and original insights on the cosmological baryonic matter cycle in the universe.