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

Preprint

Article

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

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

Reinhard Schlickeiser

Martin Kröger^*

Reinhard Schlickeiser

Martin Kröger^*

This version is not peer-reviewed

Submitted:

14 August 2024

Posted:

15 August 2024

You are already at the latest version

Alerts

Abstract

Context: 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. Aims: The competition of spontaneous star formation, stellar feedback and stellar evolution is theoretically investigated to understand the baryonic matter cycle including luminous baryonic matter in main-sequence stars and weakly luminous matter in white dwarfs, neutron stars and black holes (referred to as locked-in matter). Of particular interest is the understanding of the cosmic star formation history and the present-day gas fraction with compartmental models. Methods: For stationary rates of spontaneous star formation, continuous stellar feedback and stellar evolution exact analytical solutions of the time evolution of the fractions of gaseous, luminous stellar and locked-in stellar matter are derived. The accuracy of the analytical solutions is proven by the favorite comparison with the exact numerical solutions of the dynamical equations. Results: 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. Without stellar evolution the observations cannot be reproduced very well. 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. The gas, luminous star and locked-in stellar matter fractions indicate that the vast majority of the baryons in the present-day universe resides in the form of locked-in stellar matter in white dwarfs, neutron stars and black holes.

Keywords:

Subject: Physical Sciences - Astronomy and Astrophysics

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 (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

G (t) + S (t) + L (t) = 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 (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 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{matrix} \frac{d G}{d t} & = & - a (t) G (t) S (t) - β (t) G (t) + b (t) S (t), \end{matrix}

(2)

\begin{matrix} \frac{d S}{d t} & = & a (t) G (t) S (t) + β (t) G (t) - b (t) S (t) - c (t) S (t), \end{matrix}

(3)

and

\frac{d L}{d t} = c (t) S (t),

(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 (t_{0}) = 1, S (t_{0}) = 0, L (t_{0}) = 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

\overset{˚}{J} (t)

of new stars as a function of time,

\overset{˚}{J} (t) = a (t) G (t) S (t) + β (t) G (t) .

(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 (t) = 0

, the GSL-model Equations (2)–(4) simplify to the linear equations

\begin{matrix} \frac{d G}{d t} & = & - β (t) G (t) + b (t) S (t), \end{matrix}

(7)

\begin{matrix} \frac{d S}{d t} & = & β (t) G (t) - b (t) S (t) - c (t) S (t), \end{matrix}

(8)

and

\begin{matrix} \frac{d L}{d t} & = & c (t) S (t) . \end{matrix}

(9)

Introducing the dimensionless time variable

\tilde{t} = \int_{t_{0}}^{t} d ξ β (ξ),

(10)

Equations (7)–(9) become

\begin{matrix} \frac{d G}{d \tilde{t}} & = & - G (\tilde{t}) + r_{1} (\tilde{t}) S (\tilde{t}), \end{matrix}

(11a)

\begin{matrix} \frac{d S}{d \tilde{t}} & = & G (\tilde{t}) - [r_{1} (\tilde{t}) + r_{2} (\tilde{t})] S (\tilde{t}), \end{matrix}

(11b)

\begin{matrix} \frac{d L}{d \tilde{t}} & = & r_{2} (\tilde{t}) S (\tilde{t}), \end{matrix}

(11c)

\begin{matrix} 1 & = & G (\tilde{t}) + S (\tilde{t}) + L (\tilde{t}), \end{matrix}

(11d)

with the dimensionless ratios

r_{1} (\tilde{t} (t)) = \frac{b (t)}{β (t)}, r_{2} (\tilde{t} (t)) = \frac{c (t)}{β (t)} .

(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

β (t)

The first Equation (11a) readily provides

S (\tilde{t}) = \frac{1}{r_{1}} [\frac{d G (\tilde{t})}{d \tilde{t}} + G (\tilde{t})],

(13)

which upon insertion in Equation (11b) yields

\frac{d^{2} G (\tilde{t})}{d {\tilde{t}}^{2}} + (1 + r_{1} + r_{2}) \frac{d G (\tilde{t})}{d \tilde{t}} + r_{2} G (\tilde{t}) = 0 .

(14)

Equation (14) is solved by

G (\tilde{t}) = e^{- μ \tilde{t}} [c_{-} e^{ν \tilde{t}} + c_{+} e^{- ν \tilde{t}}],

(15)

with

μ = \frac{1 + r_{1} + r_{2}}{2}, ν = \sqrt{μ^{2} - r_{2}}

(16)

and with the two integration constants

c_{\pm}

, implying for Equation (13)

S (\tilde{t}) = \frac{e^{- μ \tilde{t}}}{r_{1}} [c_{-} (1 - μ + ν) e^{ν \tilde{t}} + c_{+} (1 - μ - ν) e^{- ν \tilde{t}}] .

(17)

The initial conditions (5), i.e.

G (z = 0) = 1

and

S (z = 0) = 0

provide

c_{\mp} = \frac{ν \mp (1 - μ)}{2 ν}

(18)

so that

\begin{matrix} j (\tilde{t}) & = & G (\tilde{t}) = e^{- μ \tilde{t}} [cosh (ν \tilde{t}) - \frac{1 - μ}{ν} sinh (ν \tilde{t})] \end{matrix}

(19)

and

S (\tilde{t}) = \frac{1}{ν} e^{- μ \tilde{t}} sinh ν t

(20)

Consequently,

\begin{matrix} L (\tilde{t}) & = & 1 - G (\tilde{t}) - S (\tilde{t}) = 1 - e^{- μ \tilde{t}} [cosh (ν \tilde{t}) + \frac{μ}{ν} sinh (ν \tilde{t})] \end{matrix}

(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 (\tilde{t})

and

L (\tilde{t})

monotonically decrease and increase, respectively, the stellar fraction

S (\tilde{t})

first increases from zero to the maximum value

S_{\max} = \frac{{(μ - ν)}^{\frac{μ - ν}{2 ν}}}{{(μ + ν)}^{\frac{μ + ν}{2 ν}}}

(22)

{\tilde{t}}_{S, \max} = \frac{1}{2 ν} ln \frac{μ + ν}{μ - ν},

(23)

and the S approaches zero in the course of time,

S_{\infty} = 0

The new star formation rate (6) for stationary SFR coefficient

β_{0}

in this case is

\begin{matrix} \overset{˚}{J} (t) & = & β_{0} G (\tilde{t} = β_{0} (t - t_{0})) \\ = & β_{0} e^{- μ β_{0} (t - t_{0})} [cosh (ν β_{0} (t - t_{0})) - \frac{1 - μ}{ν} sinh (ν β_{0} (t - t_{0}))] . \end{matrix}

(24)

3.2. Negligible Stellar Feedback

In the special case of negligible stellar feedback

b (t) = 0 = r_{1}

one obtains

\begin{matrix} j (\tilde{t}, r_{1} = 0) & = & G (\tilde{t}, r_{1} = 0) \\ = & e^{- \frac{1 + r_{2}}{2} \tilde{t}} [cosh \frac{1 - r_{2}}{2} \tilde{t} - sinh \frac{1 - r_{2}}{2} \tilde{t}] = e^{- \tilde{t}}, \end{matrix}

(25)

corresponding to

\overset{˚}{J} (t, b = 0) = β_{0} e^{- β_{0} (t - t_{0})}

. Likewise

S (\tilde{t}, r_{1} = 0) = \frac{2}{1 - r_{2}} e^{- \frac{1 + r_{2}}{2} \tilde{t}} sinh (\frac{1 - r_{2}}{2} \tilde{t}) = \frac{e^{- r_{2} \tilde{t}} - e^{- \tilde{t}}}{1 - r_{2}},

(26)

and

L (\tilde{t}, r_{1} = 0) = 1 - \frac{e^{- r_{2} \tilde{t}} - r_{2} e^{- \tilde{t}}}{1 - r_{2}} .

(27)

While

j (\tilde{t}) = G (\tilde{t})

and

L (\tilde{t})

in this case are monotonically decreasing or increasing, respectively, one finds that

S (\tilde{t})

rises first to the maximum value

S_{\max} = r_{2}^{r_{2} / (1 - r_{2})}

{\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{matrix} j (\tilde{t}, r_{1} = 0, r_{2} = 1) & = & G (\tilde{t}, r_{1} = 0) = e^{- \tilde{t}}, \\ S (\tilde{t}, r_{1} = 0, r_{2} = 1) & = & \tilde{t} e^{- \tilde{t}}, \\ L (\tilde{t}, r_{1} = 0, r_{2} = 1) & = & 1 - (1 + \tilde{t}) e^{- \tilde{t}} . \end{matrix}

(28)

Here

S (\tilde{t})

rises first to the maximum value

S_{\max} = e^{- 1}

{\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 (t) = 0 = r_{2}

one obtains

j (\tilde{t}, c = 0) = G (\tilde{t}, c = 0) = \frac{1}{1 + r_{1}} [r_{1} + e^{- (1 + r_{1}) \tilde{t}}],

(29)

corresponding to

\overset{˚}{J} (t, c = 0) = \frac{β_{0}}{b_{0} + β_{0}} [b_{0} + β_{0} e^{- (b_{0} + β_{0}) (t - t_{0})}],

(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

{Matlab}^{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) / [S (z = 0) + L (z = 0)] = G (z = 0) / [1 - G (z = 0)]

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) ≃ 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{matrix} ψ (z) & = & \frac{0.015 {(1 + z)}^{2.7}}{1 + {[\frac{(1 + z)}{2.9}]}^{5.6}} M_{⊙} {yr}^{- 1} {Mpc}^{- 3} \\ = & \frac{3.22 \cdot 10^{- 47} {(1 + z)}^{2.7}}{1 + {[\frac{(1 + z)}{2.9}]}^{5.6}} kg s^{- 1} m^{- 3}, \end{matrix}

(32)

where the time

t (z)

has been expressed as a function of redshift z using a flat

Λ

CDM Friedmann cosmology [4,32,38,46], and where

M_{⊙} = 1.989 \cdot 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

ψ_{\max} = 0.133

M_{⊙}

{yr}^{- 1}

{Mpc}^{- 3}

at redshift

z = 1.8632

, is practically constant at small redshifts well below unity, and it decreases at large redshifts

ψ (z ≫ 2.9) ≃ 1.25 \cdot 10^{- 44} {(1 + z)}^{- 2.9}

{kg}^{- 1}

m^{- 3}

Also available in Table 2 and Figure 11 of [27] is the integrated stellar density

ρ_{*} (z)

which according to [27] fulfills

d ρ_{*} (z) / d t = (1 - R) ψ (z)

and is therefore more explicitly given by

\begin{matrix} ρ_{*} (z) & = & (1 - R) \int_{z}^{\infty} d z^{'} ψ (z^{'}) |\frac{d t^{'}}{d z^{'}}|, \end{matrix}

(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 (t (z))

has to be replaced by the time- or redshift dependent term

1 - L (t)

with the fraction

L (t)

of locked-in stellar matter.

We first discuss the relation of our dynamical evolution time scale

t = \tilde{t} / β_{0}

for stationary spontaneous SFRs with the lookback time and redshift z in a flat

Λ

CDM Friedmann cosmology. In a flat

Λ

CDM Friedmann cosmology with

Ω_{m} = Ω_{0} = 0.3

and

H_{0} = 70 h_{70}

s^{- 1}

{Mpc}^{- 1} = 2.268 \cdot 10^{- 18} h_{70}

s^{- 1}

the lookback time

T_{L}

is related to the redshift according to [2,19,32]

\begin{matrix} T_{L} (z) & = & H_{0}^{- 1} \int_{1}^{1 + z} \frac{d x}{x \sqrt{1 - Ω_{0} + Ω_{0} x^{3}}} \\ = & \frac{2}{3 H_{0}} \int_{1}^{{(1 + z)}^{3 / 2}} \frac{d u}{u \sqrt{1 - Ω_{0} + Ω_{0} u^{2}}} \\ = & \frac{2}{3 \sqrt{1 - Ω_{0}} H_{0}} ln \frac{(1 + \sqrt{1 - Ω_{0}} {(1 + z)}^{\frac{3}{2}}}{\sqrt{1 - Ω_{0}} + \sqrt{1 - Ω_{0} + Ω_{0} {(1 + z)}^{3}}} . \end{matrix}

(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{matrix} T_{0} & = & H_{0}^{- 1} \int_{1}^{\infty} \frac{d x}{x \sqrt{1 - Ω_{0} + Ω_{0} x^{3}}} \\ = & \frac{2}{3 \sqrt{1 - Ω_{0}} H_{0}} ln \frac{1 + \sqrt{1 - Ω_{0}}}{\sqrt{Ω_{0}}} . \end{matrix}

(35)

The relation between the real time t, the lookback time

T_{L} (z)

and the age of the universe

T_{0}

is given by

\begin{matrix} t & = & T_{0} - T_{L} (z) = H_{0}^{- 1} \int_{1 + z}^{\infty} \frac{d x}{x \sqrt{1 - Ω_{0} + Ω_{0} x^{3}}} \\ = & \frac{2}{3 \sqrt{1 - Ω_{0}} H_{0}} \times \\ ln [\frac{(1 + \sqrt{1 - Ω_{0}} (\sqrt{1 - Ω_{0}} + \sqrt{1 - Ω_{0} + Ω_{0} {(1 + z)}^{3}}}{\sqrt{Ω_{0}} (1 + \sqrt{1 - Ω_{0}} {(1 + z)}^{\frac{3}{2}}}], \end{matrix}

(36)

implying

t_{0} = T_{0} - T_{L} (z = 1100) .

(37)

Consequently,

t - t_{0} = T_{L} (1100) - T_{L} (z),

(38)

so that for a constant spontaneous SFR coefficient

β (t) = β_{0}

one obtains

\begin{matrix} \tilde{t} (z) & = & β_{0} [T_{L} (1100) - T_{L} (z)] = \frac{2 β_{0}}{3 \sqrt{1 - Ω_{0}} H_{0}} \times \\ ln [\frac{{(1101)}^{\frac{3}{2}} (\sqrt{1 - Ω_{0}} + \sqrt{1 - Ω_{0} + Ω_{0} {(1 + z)}^{3}}}{{(1 + z)}^{\frac{3}{2}} (\sqrt{1 - Ω_{0}} + \sqrt{1 - Ω_{0} + Ω_{0} {(1101)}^{3}}}] \\ = & \frac{2 β_{0}}{3 \sqrt{1 - Ω_{0}} H_{0}} ln [{(\frac{1101}{1 + z})}^{\frac{3}{2}} \frac{1 + \sqrt{1 + \frac{3}{7} {(1 + z)}^{3}}}{1 + \sqrt{1 + \frac{3}{7} {(1101)}^{3}}}] \\ = & 3.513 \cdot 10^{17} \frac{β_{0}}{h_{70}} [0.424 + ln \frac{1 + \sqrt{1 + \frac{3}{7} {(1 + z)}^{3}}}{{(1 + z)}^{\frac{3}{2}}}] s \end{matrix}

(39)

For later use we note that

\begin{matrix} \frac{d \tilde{t}}{d z} & = & - β_{0} \frac{d T_{L} (z)}{d z} = - \frac{β_{0}}{H_{0} (1 + z) \sqrt{1 - Ω_{0} + Ω_{0} {(1 + z)}^{3}}} \\ = & - \frac{1.5 β_{0} τ_{f}}{(1 + z) \sqrt{{(1 + z)}^{3} + \frac{7}{3}}} ≃ - \frac{3}{2} β_{0} τ_{f} {(1 + z)}^{- 5 / 2}, \end{matrix}

(40)

with

τ_{f} = \frac{2}{3 H_{0} \sqrt{Ω_{0}}} = 5.37 \cdot 10^{17} h_{70}^{- 1} s .

(41)

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

\tilde{t} (z) ≃ {(1 + z)}^{- 3 / 2} β_{0} τ_{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

Ω_{tot} = 1

so that the critical density does not evolve with redshift and is given by

ρ_{0} (z) = ρ_{0} = 1.36 \cdot 10^{11} h_{70}^{2} M_{⊙} {Mpc}^{- 3} .

(43)

Consequently, with

Ω_{m} = 0.3

the matter density of the universe then is

ρ_{m} = 0.3 ρ_{0} = 4.1 \cdot 10^{10} h_{70}^{2} M_{⊙} {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

ρ_{b} = ρ_{m} / 8 = 5.1 \cdot 10^{9} h_{70}^{2} M_{⊙} {Mpc}^{- 3} = 3.45 \cdot 10^{- 28} h_{70}^{2} kg m^{- 3} .

(45)

For general

G (\tilde{t}) = j (\tilde{t})

the theoretical cosmic SFR density and the integrated stellar density as a function of redshift

ψ_{GLS} (z) d z = ρ_{b} \overset{˚}{J} (t) d \tilde{t}

, or equivalently,

\begin{matrix} ψ_{GSL} (z) & = & β_{0} ρ_{b} j (\tilde{t} (z)) |\frac{d \tilde{t}}{d z}| ≃ \frac{3 β_{0}^{2} ρ_{b} τ_{f}}{2} \frac{j (\tilde{t} (z))}{{(1 + z)}^{5 / 2}} \\ ≃ & A_{1} \frac{j (\tilde{t} (z))}{{(1 + z)}^{5 / 2}}, \end{matrix}

(46)

with

A_{1} = 3 β_{0}^{2} ρ_{b} τ_{f} / 2 = 2.78 \cdot 10^{- 10} β_{0}^{2} h_{70} kg m^{- 3} s

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

\begin{matrix} ρ_{GSL}^{*} (z) & = & β_{0}^{- 1} \int_{z}^{1100} d z^{'} [1 - L (z^{'})] ψ_{GS} (z^{'}) |\frac{d {\tilde{t}}^{'}}{d z^{'}}| \\ ≃ & \frac{3}{2} A_{1} τ_{f} \int_{z}^{1100} d z^{'} [1 - L (z^{'})] \frac{j (\tilde{t} (z^{'}))}{{(1 + z^{'})}^{5}} \\ = & B_{1} \int_{z}^{1100} d z^{'} [1 - L (z^{'})] \frac{j (\tilde{t} (z^{'}))}{{(1 + z^{'})}^{5}}, \end{matrix}

(47)

with

B_{1} = 2.24 \cdot 10^{8} β_{0}^{2}

m^{- 3}

, 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)

results from the stellar evolution to locked-in stellar matter, so that only the fraction

1 - L (z)

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{matrix} ψ (z_{E}) & = & 0 . 178_{- 0.044}^{+ 0.372} M_{⊙} {yr}^{- 1} {Mpc}^{- 3} \\ = & (3 . 83_{- 0.95}^{+ 8.00} \cdot 10^{- 46} kg s^{- 1} m^{- 3}, \end{matrix}

(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]

ρ^{*} (0) = (5 . 62_{- 1.35}^{+ 1.79} \cdot 10^{8} \frac{M_{⊙}}{{Mpc}^{3}} = (3 . 80_{- 0.91}^{+ 1.21} \cdot 10^{- 29} kg m^{- 3}

(50)

As fourth constraint we demand that

{(d ψ / d z)}_{z = 0} > 0

in order to have at least one maximum of

ψ (z)

at positive z. As fifth more stringent constraint we require that

ψ_{GSL} (z)

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 (\tilde{t})

. With

\tilde{t} (z)

from Equation (42), the SFR density (46) becomes

\begin{matrix} ψ_{GS} (z) & = & \frac{A_{1}}{1 + r_{1}} \frac{r_{1} + e^{- ϵ_{0} {(1 + z)}^{- 3 / 2}}}{{(1 + z)}^{5 / 2}}, \end{matrix}

(51)

where we introduce

ϵ_{0} = (1 + r_{1}) β_{0} τ_{f} = (b_{0} + β_{0}) τ_{f} = 5.37 \cdot 10^{17} (b_{0} + β_{0}) h_{70}^{- 1} .

(52)

In terms of

Z = ϵ_{0} {(1 + z)}^{- 3 / 2}

(53)

Equation (51) reads

ψ_{GS} (z) = \frac{A_{1}}{(1 + r_{1}) ϵ_{0}^{5 / 3}} F (Z), F (Z) = Z^{5 / 3} (r_{1} + e^{- Z}) .

(54)

The first derivative of the function

F (Z)

is given by

\frac{d F (Z)}{d Z} = \frac{5}{3} Z^{2 / 3} [r_{1} + (1 - 0.6 Z) e^{- Z}],

(55)

which is positive for all values of

r_{1} > 0.6 e^{- 5 / 3} \approx 0.042

, so that 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

r_{1} = \frac{b_{0}}{β_{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) < \frac{r_{1}}{1 + r_{1}} = \frac{b_{0}}{b_{0} + β_{0}} < 0.040

(57)

is below 4 percent.

Provided restriction (56) holds, two extrema occur at

Z_{E} = \frac{5}{3} - W (- \frac{5 r_{1}}{3} e^{5 / 3}) = \frac{5}{3} - W (- 8.824 r_{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

ψ_{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.

\frac{5}{3} - W_{- 1} (- \frac{5 r_{1} e^{5 / 3}}{3}) > ϵ_{0},

(59)

corresponding to

r_{1} < (\frac{3 ϵ_{0}}{5} - 1) e^{- ϵ_{0}},

(60)

which, since

r_{1} > 0

, requires

ϵ_{0} > 5 / 3

or with Equation (52) that

β_{0} > 3.10 \cdot 10^{- 18} \frac{h_{70}}{1 + r_{1}} Hz,

(61)

as first constraint on the value of the spontaneous SFR coefficient

β_{0}

. The maximum then occurs at the peak redshift

z_{E} = {[\frac{5}{3 ϵ_{0}} - \frac{1}{ϵ_{0}} W_{0} (- \frac{5 r_{1} e^{5 / 3}}{3})]}^{- 2 / 3} - 1 .

(62)

For small

r_{1} ≪ 1

, Taylor expansion yields

z_{E} ≃ \frac{(3 - 2 r_{1} e^{5 / 3} ϵ_{0}^{2 / 3}}{75^{1 / 3}} - 1 ≃ {(0.6 ϵ_{0})}^{2 / 3} - 1

(63)

The peak SFR density evaluates to

\begin{matrix} ψ_{GS} (z_{E}) & = & \frac{3 A_{1}}{5 (1 + r_{1}) ϵ_{0}^{5 / 3}} Z_{E}^{8 / 3} e^{- Z_{E}} \\ = & - \frac{A_{1} r_{1}}{(1 + r_{1}) ϵ_{0}^{5 / 3}} \frac{{[\frac{5}{3} - W_{0} (- \frac{5 r_{1} e^{5 / 3}}{3})]}^{8 / 3}}{W_{0} (- \frac{5 r_{1} e^{5 / 3}}{3})} \\ = & - \frac{A_{1} r_{1} ϵ_{0}}{(1 + r_{1}) W_{0} (- \frac{5 r_{1} e^{5 / 3}}{3}) {(1 + z_{E})}^{4}} \\ = & - \frac{5.37 \cdot 10^{17} A_{1} r_{1} β_{0} h_{70}^{- 1}}{W_{0} (- \frac{5 r_{1} e^{5 / 3}}{3}) {(1 + z_{E})}^{4}} ≃ 1.69 \cdot 10^{7} \frac{β_{0}^{3}}{{(1 + z_{E})}^{4}} \\ = & 3.46 \cdot 10^{- 40} \frac{β_{0}^{1 / 3} h_{70}^{8 / 3}}{{(1 + r_{1})}^{8 / 3}} kg s^{- 1} m^{- 3}, \end{matrix}

(64)

where we used the expansion

W_{0} (| x | ≪ 1) ≃ x

and Equation (63). Equating (64) and (48) yields

β_{0} = (1 . 36_{- 0.78}^{+ 38.61} \cdot 10^{- 18} {(\frac{1 + r_{1}}{h_{70}})}^{8} Hz .

(65)

This

β_{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} {(\frac{1 + r_{1}}{h_{70}})}^{9}, z_{E} < 4.49 {(\frac{1 + r_{1}}{h_{70}})}^{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} d z \frac{e^{- α {(1 + z)}^{- 3 / 2}}}{{(1 + z)}^{5}} = \frac{2}{3 α^{8 / 3}} Γ [\frac{8}{3}, \frac{α}{{(1 + z)}^{3 / 2}}]

(67)

the integrated stellar mass density (47) with

L (z) = 0

becomes

\begin{matrix} ρ_{GS}^{*} (z) & = & \frac{B_{1}}{1 + r_{1}} \int_{z}^{1100} d z^{'} \frac{r_{1} + e^{- ϵ_{0} {(1 + z^{'})}^{- 3 / 2}}}{{(1 + z^{'})}^{5}} \\ = & \frac{B_{1}}{1 + r_{1}} (\frac{r_{1}}{4} (\frac{1}{{(1 + z)}^{4}} - 1100^{- 4} \\ + & \frac{2}{3 ϵ_{0}^{8 / 3}} [Γ (\frac{8}{3}, {(1100)}^{3 / 2} ϵ_{0}) - Γ (\frac{8}{3}, \frac{ϵ_{0}}{{(1 + z)}^{3 / 2}}]) \\ ≃ & \frac{2.24 \cdot 10^{8} β_{0}^{2}}{4 (1 + r_{1})} [\frac{r_{1}}{{(1 + z)}^{4}} + \frac{8}{3} \frac{γ (\frac{8}{3}, \frac{ϵ_{0}}{{(1 + z)}^{3 / 2}}}{ϵ_{0}^{8 / 3}}] kg m^{- 3}, \end{matrix}

(68)

where

γ

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

ρ_{GS}^{*} (0) = \frac{5.60 \cdot 10^{7} β_{0}^{2}}{1 + r_{1}} (r_{1} + \frac{8}{3} \frac{γ (\frac{8}{3}, ϵ_{0})}{ϵ_{0}^{8 / 3}} kg m^{- 3}

(69)

in a power-law like fashion at large redshifts

z ≫ z_{E}

ρ_{GSL}^{*} (z) ≃ \frac{B_{1}}{4 {(1 + z)}^{4}},

(70)

since

γ (s, ϵ ≪ 1) ≃ s^{- 1} ϵ^{s}

to lowest order in

ϵ

. Equating the result (69) with the observed value (50) provides

\frac{{\tilde{β}}_{0}^{2}}{1 + r_{1}} (r_{1} + \frac{8}{3} \frac{γ (\frac{8}{3}, ϵ_{0})}{ϵ_{0}^{8 / 3}} = (6 . 79_{- 1.70}^{+ 2.16}) \cdot 10^{- 3},

(71)

where we scale

β_{0} = 10^{- 17} {\tilde{β}}_{0}

(72)

in units of

10^{- 17}

Hz. The function

γ (8 / 3, 5.37 {\tilde{β}}_{0})

is shown in Figure 4 and varies from 0 to the finite value

Γ (8 / 3) \approx 1.504

at large arguments. For the nominal values

r_{1} = 0

and

h_{70} = 1

one obtains from Equation (71)

\frac{{\tilde{β}}_{0}^{2 / 3}}{γ (8 / 3, 5.37 {\tilde{β}}_{0})} ≃ \frac{{\tilde{β}}_{0}^{2 / 3}}{1.504} = 4 . 44_{- 1.07}^{+ 1.48},

(73)

which leads to

{\tilde{β}}_{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{β}}_{0}

-value these two constraints with the GS-model.

Adopting as a best compromise value

{\tilde{β}}_{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

ψ_{GS} (z_{E}) \propto {\tilde{β}}_{0}^{1 / 3}

in Equation (64) and

ρ_{GS}^{*} (0) \propto {\tilde{β}}_{0}^{- 2 / 3}

in Equation (69) on the only free parameter

{\tilde{β}}_{0}

, it is obvious that by decreasing

{\tilde{β}}_{0}

from its compromise value 2 would decrease

ψ_{GS} (z)

to smaller values, but would increase at the same time

ρ_{GS}^{*} (z)

. Alternatively, increasing

{\tilde{β}}_{0}

from its compromise value 2 would decrease

ρ_{GS}^{*} (z)

to smaller values, but at the same time would make

ψ_{GS} (z)

even larger.

With this compromise value

{\tilde{β}}_{0} = 2.0

and the ratio

r_{1} = 0.04

the present-day gas fraction is

G (r_{2} = 0, z = 0) = 3.85 \cdot 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

ψ_{GS} (z) = A_{1} \frac{e^{- ϵ {(1 + z)}^{- 3 / 2}}}{{(1 + z)}^{5 / 2}}

(76)

with

ϵ = ϵ_{0} (r_{1} = 0) = 5.37 {\tilde{β}}_{0} h_{70}^{- 1} .

(77)

The density (76) attains its maximum value

\begin{matrix} ψ_{GS} (z_{E}) & = & {(\frac{5}{3 ϵ})}^{5 / 3} A_{1} e^{- \frac{5}{3}} \\ = & 7.47 \cdot 10^{- 46} {\tilde{β}}_{0}^{1 / 3} h_{70}^{8 / 3} kg s^{- 1} m^{- 3} \end{matrix}

(78)

at the peak redshift

z_{E} = {(0.6 ϵ)}^{2 / 3} - 1 = 2.18 {(\frac{{\tilde{β}}_{0}}{h_{70}})}^{2 / 3} - 1

(79)

The two observational constrains (49) and (48) for

z_{E}

and the peak SFR density then provide the constraints

{\tilde{β}}_{0} = (1 . 61_{- 0.73}^{+ 0.87}) h_{70}

(80)

and

{\tilde{β}}_{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{matrix} ρ_{GS}^{*} (z) & = & B_{1} \int_{z}^{1100} d z^{'} \frac{e^{ϵ {(1 + z^{'})}^{- 3 / 2}}}{{(1 + z^{'})}^{5}} ≃ \frac{2 B_{1}}{3 ϵ^{8 / 3}} γ (\frac{8}{3}, \frac{ϵ}{{(1 + z)}^{3 / 2}}) \\ ≃ & \frac{1.49 \cdot 10^{- 26} {\tilde{β}}_{0}^{2}}{ϵ^{8 / 3}} γ (\frac{8}{3}, \frac{ϵ}{{(1 + z)}^{3 / 2}}) kg m^{- 3} \\ = & \frac{1.69 \cdot 10^{- 28} h_{70}^{8 / 3}}{{\tilde{β}}_{0}^{2 / 3}} γ (\frac{8}{3}, \frac{5.37 {\tilde{β}}_{0}}{h_{70} {(1 + z)}^{3 / 2}}) kg m^{- 3} . \end{matrix}

(82)

The stellar mass density thus monotonically decreases from the value

\begin{matrix} ρ_{GS}^{*} (0) & = & \frac{1.69 \cdot 10^{- 28} h_{70}^{8 / 3}}{{\tilde{β}}_{0}^{2 / 3}} γ (\frac{8}{3}, 5.37 \frac{{\tilde{β}}_{0}}{h_{70}}) kg m^{- 3} \\ = & \frac{2.54 \cdot 10^{- 28} h_{70}^{8 / 3}}{{\tilde{β}}_{0}^{2 / 3}} kg m^{- 3} \end{matrix}

(83)

in a power-law like fashion at large redshifts

z ≫ z_{E}

ρ_{GSL}^{*} (z) ≃ \frac{B_{1}}{4 {(1 + z)}^{4}},

(84)

since

γ (8 / 3, x) = (3 / 8) x^{8 / 3} + O (x^{10 / 3}

. The estimate (74) for

{\tilde{β}}_{0}

from the result (83) with the observed value (50) for the nominal values and

h_{70} = 1

remains unchanged, i.e.

{\tilde{β}}_{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{β}}_{0}

-value these three constraints with the simplified GS-model. Again a compromise value of

{\tilde{β}}_{0} = 2.0

is suggested. With this compromise value

{\tilde{β}}_{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 \cdot 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} \neq 0

) besides spontaneous star formation (

β_{0} \neq 0

) and stellar evolution

c_{0} \neq 0

by using (19), instead of (29), for

j (\tilde{t})

6.1. Theoretical Results

With

\tilde{t} (z)

from Equation (42) and Equation (19) for

j (\tilde{t})

the SFR density (46) becomes

\begin{matrix} ψ_{GSL} (z) & = & \frac{A_{1}}{{(1 + z)}^{5 / 2}} e^{- μ \tilde{t}} [cosh (ν \tilde{t}) - \frac{1 - μ}{ν} sinh (ν \tilde{t})] \\ = & \frac{A_{1}}{ϵ^{5 / 3}} F (\tilde{t}), \end{matrix}

(87)

with

μ = (1 + r_{1} + r_{2}) / 2

ν = \sqrt{μ^{2} - r_{2}}

ϵ = β_{0} τ_{f} = 5.37 {\tilde{β}}_{0} h_{70}^{- 1}

, as before in Equation (77), and the function3

F (\tilde{t}) = {\tilde{t}}^{5 / 3} e^{- μ \tilde{t}} [cosh (ν \tilde{t}) + \frac{μ - 1}{ν} sinh (ν \tilde{t})] .

(88)

Likewise with Equation (21) one obtains

1 - L (z) = e^{- μ \tilde{t}} [cosh ν \tilde{t} + \frac{μ}{ν} sinh ν \tilde{t}]

(89)

with

\tilde{t} (z) = ϵ {(1 + z)}^{- 3 / 2}

. Consequently, the integrated stellar density (47) is given by

\begin{matrix} ρ_{GSL}^{*} (z) & = & B_{1} \int_{z}^{1100} d z^{'} \frac{e^{- 2 μ \tilde{t} (z^{'})}}{{(1 + z^{'})}^{5}} [cosh ν \tilde{t} (z^{'}) \\ + \frac{μ}{ν} sinh ν \tilde{t} (z^{'})] [cosh ν \tilde{t} (z^{'}) + \frac{μ - 1}{ν} sinh ν \tilde{t} (z^{'})] \\ ≃ & \frac{B_{1}}{6 ν^{2} ϵ^{8 / 3}} \int_{0}^{\frac{ϵ}{{(1 + z)}^{3 / 2}}} d \tilde{t} {\tilde{t}}^{5 / 3} e^{- 2 μ \tilde{t}} [(ν + μ) e^{ν \tilde{t}} + (ν - μ) e^{- ν \tilde{t}}] \\ \times & [(ν + μ - 1) e^{ν \tilde{t}} + (ν + 1 - μ) e^{- ν \tilde{t}}] \\ = & \frac{B_{1}}{6 ν^{2} ϵ^{8 / 3}} \int_{0}^{\frac{ϵ}{{(1 + z)}^{3 / 2}}} d \tilde{t} {\tilde{t}}^{5 / 3} e^{- 2 μ \tilde{t}} [2 (ν^{2} + μ (1 - μ)) \\ + + (ν + μ) (ν + μ - 1) e^{2 ν \tilde{t}} + (ν - μ) (ν + 1 - μ) e^{- 2 ν \tilde{t}}] \\ = & \frac{B_{1}}{6 ν^{2} {(2 ϵ)}^{8 / 3}} [\frac{2 (ν^{2} + μ (1 - μ))}{μ^{8 / 3}} γ (\frac{8}{3}, \frac{2 μ ϵ}{{(1 + z)}^{3 / 2}} \\ + + \frac{(ν + μ) (ν + μ - 1)}{{(μ - ν)}^{8 / 3}} γ (\frac{8}{3}, \frac{2 (μ - ν) ϵ}{{(1 + z)}^{3 / 2}}) \\ + \frac{(ν - μ) (ν + 1 - μ)}{{(μ + ν)}^{8 / 3}} γ (\frac{8}{3}, \frac{2 (μ + ν) ϵ}{{(1 + z)}^{3 / 2}}] . \end{matrix}

(90)

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

G (z) = \frac{1}{2 ν} [(ν + μ - 1) e^{\frac{(ν - μ) ϵ}{{(1 + z)}^{3 / 2}}} + (ν + 1 - μ) e^{- \frac{(ν + μ) ϵ}{{(1 + z)}^{3 / 2}}}] .

(91)

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

\begin{matrix} \frac{d F (\tilde{t})}{d \tilde{t}} & = & \frac{{\tilde{t}}^{2 / 3} e^{- μ \tilde{t}} cosh ν \tilde{t}}{ν} {\frac{5}{3} [ν + (μ - 1) tanh (ν \tilde{t})] \\ + \tilde{t} [(μ - r_{2}) tanh (ν \tilde{t}) - ν]}, \end{matrix}

(92)

so that a maximum occurs at

{\tilde{t}}_{E}

given by the solution of the transcendental equation

[ν - (μ - r_{2}) tanh (ν {\tilde{t}}_{E})] {\tilde{t}}_{E} = \frac{5}{3} [ν + (μ - 1) tanh (ν {\tilde{t}}_{E})] .

(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

ν = μ = (1 + r_{1}) / 2

and with

(1 + r_{1}) {\tilde{t}}_{E} = Z

, with Z defined in Equation (53), Equation (93) simplifies to

0.6 Z [1 - tanh \frac{Z}{2}] = 1 + r_{1} + (r_{1} - 1) tanh \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{matrix} (5 . 715_{- 1.369}^{+ 1.820} \frac{ν {\tilde{β}}_{0}^{2 / 3}}{h_{70}^{8 / 3}}) & = & \frac{2 (ν^{2} + μ (1 - μ))}{μ^{8 / 3}} γ (\frac{8}{3}, \frac{10.74 μ {\tilde{β}}_{0}}{h_{70}}) \\ + \frac{(ν + μ) (ν + μ - 1)}{{(μ - ν)}^{8 / 3}} γ (\frac{8}{3}, \frac{10.74 (μ - ν) {\tilde{β}}_{0}}{h_{70}}) \\ + \frac{(ν - μ) (ν + 1 - μ)}{{(μ + ν)}^{8 / 3}} γ (\frac{8}{3}, \frac{10.74 (μ + ν) {\tilde{β}}_{0}}{h_{70}}) . \end{matrix}

(95)

Adopting values of

{\tilde{β}}_{0} / h_{70} > 1

and assuming

| μ - ν | > 0.5

, all incomplete

γ

-functions are well approximated by

1.504

, so that the constraint (95) simplifies to

\begin{matrix} (4 . 421_{- 0.910}^{+ 1.230} ν \frac{{\tilde{β}}_{0}^{2 / 3}}{h_{70}^{8 / 3}} & = & \frac{2 (ν^{2} + μ (1 - μ))}{μ^{8 / 3}} + \frac{(ν + μ) (ν + μ - 1)}{{(μ - ν)}^{8 / 3}} \\ + \frac{(ν - μ) (ν + 1 - μ)}{{(μ + ν)}^{8 / 3}} . \end{matrix}

(96)

Secondly, the observed peak redshift (49) yields for

{\tilde{t}}_{E} = \frac{ϵ}{{(1 + z_{E})}^{3 / 2}} = (1 . 03_{- 0.36}^{+ 0.87}) {\tilde{β}}_{0} h_{70}^{- 1} .

(97)

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

(3 . 83_{- 0.95}^{+ 8.00}) \cdot 10^{- 46} = 1.69 \cdot 10^{- 45} {\tilde{β}}_{0}^{1 / 3} h_{70}^{8 / 3} F ({\tilde{t}}_{E})

(98)

\begin{matrix} 0 . 227_{- 0.056}^{+ 0.473} = {\tilde{β}}_{0}^{1 / 3} h_{70}^{8 / 3} F ({\tilde{t}}_{E}) = [{(1 . 03_{- 0.36}^{+ 0.87}]}^{5 / 3} {\tilde{β}}_{0}^{2} h_{70} \\ \times exp [- (1 . 03_{- 0.36}^{+ 0.87}) μ {\tilde{β}}_{0} h_{70}^{- 1}] (cosh ((1 . 03_{- 0.36}^{+ 0.87}) ν {\tilde{β}}_{0} h_{70}^{- 1}) \\ + \frac{μ - 1}{ν} sinh ((1 . 03_{- 0.36}^{+ 0.87}) ν {\tilde{β}}_{0} h_{70}^{- 1}) \end{matrix}

(99)

with

F ({\tilde{t}}_{E}) = {\tilde{t}}_{E}^{5 / 3} e^{- μ {\tilde{t}}_{E}} [cosh (ν {\tilde{t}}_{E}) + \frac{μ - 1}{ν} sinh (ν {\tilde{t}}_{E})]

(100)

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

As fourth constraint we demand that

{(d ψ / d z)}_{z = 0} > 0

in order to have a maximum of

ψ (z)

. With

\frac{d ψ_{GSL} (z)}{d z} = - \frac{2 A_{1}}{3 ϵ {\tilde{t}}^{5 / 3}} \frac{d F (\tilde{t})}{d \tilde{t}}

(101)

and Equation (92) one has to demand that

ν + (μ - 1) tanh (ν ϵ) < \frac{3 ϵ}{5} [ν - (μ - r_{2}) tanh (ν ϵ)],

(102)

leading, after insertion of Equation (77), to

\begin{matrix} tanh (\frac{5.37 ν {\tilde{β}}_{0}}{h_{70}}) & < & \frac{ν (3.22 \frac{{\tilde{β}}_{0}}{h_{70}} - 1)}{μ - 1 + 3.22 \frac{{\tilde{β}}_{0}}{h_{70}} (μ - r_{2})} . \end{matrix}

(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{β}}_{0} \in [0.5, 2.7]

the compatibility of the five constraints on

ψ_{GSL} (z_{E})

z_{E}

ρ_{GSL}^{*} (0)

, and

{d ψ (z) / d z |}_{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

ρ_{GSL}^{*} (0)

. The green regions on the rightest panel (labeled summary) for the values of

{\tilde{β}}_{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{β}}_{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{β}}_{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{β}}_{0} = 0.9, 1.0, 1.4, 2.0

) within the blue (three constraints fulfilled for

{\tilde{β}}_{0} = 0.5, 0.7

) areas of the corresponding panels in Figure 5. Although the agreement is not perfect for any

({\tilde{β}}_{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{β}}_{0}

values we also determine the ranges of the present-day gas fraction from Equation (91)

G (0) = \frac{1}{2 ν} [(ν + μ - 1) e^{(ν - μ) ϵ} + (ν + 1 - μ) e^{- (ν + μ) ϵ}] .

(104)

Using all admissible (

r_{1}, r_{2}

) pairs we obtain

G (0) \in [0.007, 0.014]

for

{\tilde{β}}_{0} = 0.9

G (0) \in [0.0047, 0.0205]

for

{\tilde{β}}_{0} = 1

G (0) \in [2.2 \times 10^{- 5}, 0.0077]

for

{\tilde{β}}_{0} = 2

, and

G (0) \in [1.48 \times 10^{- 6}, 0.0041]

for

{\tilde{β}}_{0} = 2.5

. All these values are within precision below about 2 percent.

Figure 5. For each

{\tilde{β}}_{0}

value there is one line. Each plot versus

r_{1}

and

r_{2}

(double-logarithmic). From left to right: Constraint 1:

ψ_{GSL} (z_{E}) = 3 . 83_{- 0.95}^{+ 8} \cdot 10^{- 46}

s^{- 1}

m^{- 3}

according to Equation (48), constraint 2:

z_{E} = 2_{- 1}^{+ 1}

(49), constraint 3:

ρ_{GSL}^{*} (0) = 3 . 80_{- 0.91}^{+ 1.21} \cdot 10^{- 29}

m^{- 3}

(50), constraint 4: the sign of

ψ_{GSL}^{'} (0)

, constraint 5: the existence of a single extremum in

ψ_{GSL} (z)

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

ψ_{GSL} (z = 0) > 0

. The summary displays the amount of simultaneously fulfilled criteria. For

{\tilde{β}}_{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{β}}_{0}

value there is one line. Each plot versus

r_{1}

and

r_{2}

(double-logarithmic). From left to right: Constraint 1:

ψ_{GSL} (z_{E}) = 3 . 83_{- 0.95}^{+ 8} \cdot 10^{- 46}

s^{- 1}

m^{- 3}

according to Equation (48), constraint 2:

z_{E} = 2_{- 1}^{+ 1}

(49), constraint 3:

ρ_{GSL}^{*} (0) = 3 . 80_{- 0.91}^{+ 1.21} \cdot 10^{- 29}

m^{- 3}

(50), constraint 4: the sign of

ψ_{GSL}^{'} (0)

, constraint 5: the existence of a single extremum in

ψ_{GSL} (z)

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

ψ_{GSL} (z = 0) > 0

. The summary displays the amount of simultaneously fulfilled criteria. For

{\tilde{β}}_{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

ψ (z)

(symbols) and

ψ_{GSL} (z)

according to Equation (87) (family of gray lines) as well as (Right column) reported

ρ_{*} (z)

(symbols) and

ρ_{GSL}^{*} (z)

according to Equation (90) (family of gray lines) for

h_{70} = 1

and various

{\tilde{β}}_{0}

(increasing from top to bottom). The families are created using 1000 randomly chosen, but eligible (

r_{1}, r_{2}

) pairs for given

{\tilde{β}}_{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{β}}_{0} = 0.5

and

{\tilde{β}}_{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

ψ (z)

(symbols) and

ψ_{GSL} (z)

according to Equation (87) (family of gray lines) as well as (Right column) reported

ρ_{*} (z)

(symbols) and

ρ_{GSL}^{*} (z)

according to Equation (90) (family of gray lines) for

h_{70} = 1

and various

{\tilde{β}}_{0}

(increasing from top to bottom). The families are created using 1000 randomly chosen, but eligible (

r_{1}, r_{2}

) pairs for given

{\tilde{β}}_{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{β}}_{0} = 0.5

and

{\tilde{β}}_{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

ψ_{GLS} (z_{E}) \propto h_{70}^{8 / 3}

ρ_{GLS}^{*} (0) \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

s^{- 1}

{Mpc}^{- 1}

or even of

H_{0} = 60

s^{- 1}

{Mpc}^{- 1}

indeed provides better agreements of our theoretical results in Figure with the observations as in this case

ψ_{GLS} (z_{E})

and

ρ_{GLS} (0)

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

s^{- 1}

{Mpc}^{- 1}

(corresponding to

h_{70} = 0.857

) for

{\tilde{β}}_{0} = 0.8

0.9

, and

1.0

(from top to bottom).

Figure 7. Same as Figure 6 for

H_{0} = 60

km s

s^{- 1}

{Mpc}^{- 1}

(corresponding to

h_{70} = 0.857

) for

{\tilde{β}}_{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{matrix} F (\tilde{t}) & = & {\tilde{t}}^{5 / 3} e^{- \tilde{t}} \\ 1 - L (\tilde{t}) & = & \frac{e^{- r_{2} \tilde{t}} - r_{2} e^{- \tilde{t}}}{1 - r_{2}} \end{matrix}

(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 (\tilde{t})

has a single maximum at

{\tilde{t}}_{E} = 5 / 3

, so that the SFR density

ψ_{GLS} (z)

attains its maximum

ψ_{GLS} (z_{E}) = {(\frac{5}{3 ϵ})}^{5 / 3} A_{1} e^{- \frac{5}{3}} = 7.47 \cdot 10^{- 46} {\tilde{β}}_{0}^{1 / 3} h_{70}^{8 / 3} kg s^{- 1} m^{- 3}

(106)

at the peak redshift

1 + z_{E} = {(0.6 ϵ)}^{2 / 3} = 2.18 {\tilde{β}}_{0}^{2 / 3} h_{70}^{- 2 / 3} .

(107)

Consequently,

\frac{ψ_{GLS} (z)}{ψ_{GLS} (z_{E})} = {(\frac{1 + z_{E}}{1 + z})}^{5 / 2} e^{- ϵ [{(1 + z)}^{- 3 / 2} - {(1 + z_{E})}^{- 3 / 2}]} .

(108)

Likewise with

μ = (1 + r_{2}) / 2

and

ν = | 1 - r_{2} | / 2

in this case Equation (90) simplifies to

ρ_{GSL}^{*} (z, r_{1} = 0, r_{2}) = \frac{2 B_{1}}{3 (1 - r_{2}) {(2 ϵ)}^{8 / 3}} \times [{(\frac{2}{1 + r_{2}})}^{8 / 3} γ (\frac{8}{3}, \frac{(1 + r_{2}) ϵ}{{(1 + z)}^{3 / 2}}) - r_{2} γ (\frac{8}{3}, \frac{2 ϵ}{{(1 + z)}^{3 / 2}})] .

(109)

and in the particular case

r_{2} = 1

ρ_{GSL}^{*} (z, r_{1} = 0, r_{2} = 1) ≃ \frac{2 B_{1}}{3 ϵ^{8 / 3}} \int_{0}^{\frac{ϵ}{{(1 + z)}^{3 / 2}}} d \tilde{t} {\tilde{t}}^{5 / 3} (1 + \tilde{t}) e^{- 2 \tilde{t}} = \frac{B_{1}}{3 {(2 ϵ)}^{8 / 3}} [2 γ (\frac{8}{3}, \frac{2 ϵ}{{(1 + z)}^{3 / 2}} + γ (\frac{11}{3}, \frac{2 ϵ}{{(1 + z)}^{3 / 2}})],

(110)

while Equations (109) and (110) imply for

r_{2} \neq 1

\begin{matrix} ρ_{GSL}^{*} (z = 0, r_{1} = 0) = \frac{2.660 \cdot 10^{- 29} h_{70}^{8 / 3}}{(1 - r_{2}) {\tilde{β}}_{0}^{2 / 3}} \times \\ [{(\frac{2}{1 + r_{2}})}^{8 / 3} γ (\frac{8}{3}, 5.37 (1 + r_{2})) - r_{2} γ (\frac{8}{3}, 10.74)] \\ ≃ & \frac{4.00 \cdot 10^{- 29} [{(\frac{2}{1 + r_{2}})}^{8 / 3} - r_{2}] h_{70}^{8 / 3}}{(1 - r_{2}) {\tilde{β}}_{0}^{2 / 3}} kg m^{- 3}, \end{matrix}

(111)

whereas for

r_{2} = 1

ρ_{GSL}^{*} (0, r_{1} = 0, r_{2} = 1)) ≃ 9.31 \cdot 10^{- 29} h_{70}^{8 / 3} {\tilde{β}}_{0}^{- 2 / 3} kg m^{- 3},

(112)

where we approximated both incomplete gamma functions by

Γ (8 / 3) = 1.504

and

Γ (11 / 3) = (8 / 3) Γ (8 / 3)

. Note that with the limit

lim_{r_{2} \to 1} \frac{[{(\frac{2}{1 + r_{2}})}^{8 / 3} - r_{2}]}{1 - r_{2}} = lim_{x \to 0} [\frac{{(\frac{2}{2 - x})}^{8 / 3} + x - 1}{x}] = 7 / 3

(113)

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

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

{\tilde{β}}_{0} = (1 . 61_{- 0.73}^{+ 0.88}) h_{70},

(114)

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

(115)

and

{\tilde{β}}_{0} = (1 . 080_{- 0.367}^{+ 0.548}) {[\frac{{(\frac{2}{1 + r_{2}})}^{8 / 3} - r_{2}}{1 - r_{2}}]}^{3 / 2} h_{70}^{4}

(116)

Whereas the first two constraints (114) and (115) restrict the value of

{\tilde{β}}_{0}

, we illustrate in Figure 8 the allowed range of values values of the parameters

{\tilde{β}}_{0}

and

r_{2}

in the case

r_{1} = 0

provided by constraint (116) from the observed

ρ^{*} (0)

For

r_{1} = 0

, in the limit

r_{2} \to \infty

\begin{matrix} ρ_{GSL} (z) & = & \frac{2 B_{1}}{3 ϵ^{8 / 3}} \int_{0}^{\frac{ϵ}{{(1 + z)}^{3 / 2}}} d \tilde{t} {\tilde{t}}^{5 / 3} e^{- 2 \tilde{t}}, \\ = & \frac{2 B_{1}}{3 {(2 ϵ)}^{8 / 3}} γ (\frac{8}{3}, \frac{2 ϵ}{{(1 + z)}^{3 / 2}}), \end{matrix}

(117)

so that

ρ_{GSL}^{*} (0) = 4.0 \cdot 10^{- 29} h_{70}^{8 / 3} {\tilde{β}}_{0}^{- 2 / 3}

m^{- 3}

. This latter expression agrees with the values for

ρ_{GSL} (0)

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{matrix} G (z) & = & e^{\frac{- 5.37 {\tilde{β}}_{0} h_{70}^{- 1}}{{(1 + z)}^{3 / 2}}}, \\ S (z) & = & \frac{e^{\frac{- 5.37 {\tilde{β}}_{0} h_{70}^{- 1}}{{(1 + z)}^{3 / 2}}} - e^{\frac{- 5.37 r_{2} {\tilde{β}}_{0} h_{70}^{- 1}}{{(1 + z)}^{3 / 2}}}}{r_{2} - 1}, \\ L (z) & = & 1 - \frac{r_{2} e^{\frac{- 5.37 {\tilde{β}}_{0} h_{70}^{- 1}}{{(1 + z)}^{3 / 2}}} - e^{\frac{- 5.37 r_{2} {\tilde{β}}_{0} h_{70}^{- 1}}{{(1 + z)}^{3 / 2}}}}{r_{2} - 1} . \end{matrix}

(118)

Adopting for

h_{70} = 1

the value

{\tilde{β}}_{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{β}}_{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 \cdot 10^{- 3}

S (z = 0) = 8.19 \cdot 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 (t)

as well as in main-sequence stars with the fraction

S (t)

. The third compartment with fraction

L (t)

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 (

β (t) G (t)

), of stellar feedback (

b (t) S (t)

) of stellar to gaseous matter, and of the formation (

c (t) S (t)

) 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 ξ β (ξ)

for arbitrarily but given time-dependent spontaneous SFR coefficient

β (t)

, as well as the ratios

r_{1} = b (t) / β (t)

and

r_{2} = c (t) / β (t)

, 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

Λ

CDM Friedmann cosmology the relationship between the reduced time variable

\tilde{t} (z)

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

β_{0}

. The best compromise value

β_{0} = 2.0 \cdot 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.

Acknowledgments

R.S. gratefully acknowledges the institutional support by the Astrophysics Group headed by Prof. Dr. Wolfgang Duschl and Prof. Dr. Sebastian Wolf at the Institut für Theoretische Physik und Astrophysik of the Christian-Albrechts-Universität in Kiel, Germany.

Appendix A. Lambert Functions

In this manuscript we encounter transcendental equations of the form

y = A ln (y) + B,

(A1)

with arbitrary values of A and C. Substituting

y = e^{- x}

then yields for Equation (A1)

e^{- x} = - A (x - \frac{B}{A}),

(A2)

which can be solved (see Appendix G in [23] as

x = \frac{B}{A} + W (- \frac{e^{- B / A}}{A})

(A3)

in terms of Lambert functions. Consequently, the solution of Equation (A1) is given by

y = e^{- x} = - A W (- \frac{e^{- B / A}}{A}),

(A4)

where we used Lambert’s equation

z = W (z) e^{W (z)}

(A5)

defining the Lambert function

W (z)

. The Lambert Equation (A5) can only be solved for real-valued z if

z \geq - 1 / e

. One obtains the principal branch

W_{0} (z)

for arguments

z \geq 0

and the two branches

W_{0} (z)

and

W_{- 1} (z)

- 1 / e \leq z \leq 0

. For these negative real arguments the principal branch has values

W_{0} (z) \in [- 1, 0]

and the lower branch

W_{- 1} (z) < - 1

has values smaller than

- 1

References

Bonnor, W. B. Monthly Notices R. Astron. Soc 1956, 16, 351. [Google Scholar]
Boylan-Kolchin, M.; Weisz, D.R. Monthly Notices R. Astron. Soc. 2021, 505, 2764. [Google Scholar] [CrossRef]
Carroll, B.W.; Ostlie, D.A. An Introduction to Modern Astrophysics, 2nd ed.; Cambridge University Press: Cambridge, 2017. [Google Scholar]
Carroll, S.M.; Press, W.H.; Turner, E.L. Ann. Rev. Astron. Astrophys. 1992, 30, 499. [CrossRef]
Clayton, D. D. Principles of Stellar Evolution and Nucleosynthesis; University of Chicago Press: Chicago, 1983. [Google Scholar]
Cousin, M.; Lagache, G.; Bethermin, M.; Guiderdoni, B. Astron. Astrophys. 2015, 575.
Cucciati, O., Tresse, L., Ilbert, O., et al. Astron. Astrophys 2012, 539, A31. [CrossRef]
Dahlen, T., Mobasher, B., Dickinson, M., et al. Astrophys. J. 2007.
Dave, R.; Finlator, K.; Oppenheimer, B.D. Monthly Notices R. Astron. Soc. 2011, 416, 1354. [Google Scholar]
Davis, A.J.; Khochfar, S.; Dalla, V.C. Monthly Notices R. Astron. Soc. 2014, 443, 985. [Google Scholar] [CrossRef]
Ebert, R. Z. Astrophysik 1955, 37, 215. [Google Scholar]
Elmegreen, B.G. Astrophys. J. 1997, 477, 196. [Google Scholar]
Estrada, E. Phys. Rep. 2020, 869, 1. [CrossRef] [PubMed]
Frank, J.; King, A.; Raine, D. Accretion Power in Astrophysics, 3rd ed.; Cambridge University Press: Cambridge, 2002. [Google Scholar]
Hansen, C.J.; Kawaler, S.D.; Trimble, V. Stellar Interiors: Physical Principles, Structure, and Evolution, 2nd ed.; S: (New York, 2004. [Google Scholar]
Hennebelle, P.; Chabrier, G. Astrophys. J. Lett. 2011, 743, L29. [Google Scholar]
Hethcote, H.W. Biomathematics. In Applied Mathematical Ecology; Levin, S.A., Hallam, T.G., Gross, L.J., Eds.; Springer: Berlin, 1989; Vol. 18, pp. 119–144. [Google Scholar]
Hethcote, H.W. SIAM. Rev 2000, 42, 599. [Google Scholar]
Hogg, D.W. 2000, arXiv:astro-ph/9905116.
Kendall, D.G. Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability; University of California Press: Berkeley, CA, USA, 1956; p. 149. [Google Scholar]
Kermack, W.O.; McKendrick, A.G. Proc. R. Soc. A 1927, 115, 700.
Klessen, R.S.; Glover, S.C. Annu. Rev. Astron. Astrophys. 2023, 61, 65. [CrossRef]
Kröger, M.; Schlickeiser, R. J. Phys. A 2020, 53, 505601. [CrossRef]
Lang, K.R. Astrophysical Formulae: Space, Time, Matter, and Cosmology, 3rd ed.; Springer: Berlin, 1999. [Google Scholar]
Lang, K.R. The Life and Death of Stars; Cambridge University Press: Cambridge, 2013. [Google Scholar]
Low, M. M. M.; Klessen, R.S. Control of Star Formation by Supersonic Turbulence. Rev. Mod. Phys. 2004, 76, 125–194. [Google Scholar]
Madau, P.; Dickinson, M. Annu. Rev. Astron. Astrophys 2014, 52, 415. [CrossRef]
Maoz, D.; Graur, O. Astrophys. J. 2017, 848, 25. [CrossRef]
Mercado, F. J., Bullock, J. S., Boylan-Kolchin, M., et al. Monthly Notices R. Astron. Soc. 2021, 501, 5121. [Google Scholar]
Moustakas, J., Coil, A. L., Aird, J., et al. Astrophys. J. 2013. [Google Scholar]
Naoz, S.; Perna, R. Astrophys. J. 2014. [Google Scholar]
Peebles, P. J. E. Principles of Physical Cosmology; Princeton University Press: Princeton, NJ, United States, 1993. [Google Scholar]
Peng, Y.; Maiolino, R. Monthly Notices R. Astron. Soc. 2014, 443, 3643. [Google Scholar] [CrossRef]
Perez-Gonzalez, P. G., Rieke, G. H., Villar, V., et al. Astrophys. J. 2008, 675, 234. [Google Scholar]
Péroux, C.; Howk, J.C. Annu. Rev. Astron. Astrophys. 2020, 58, 363. [CrossRef]
Phillips, A. The Physics of Stars; John Wiley & Sons: Chichester, 1999. [Google Scholar]
Portinari, L.; Chiosi, C.; Bressan, A. Astron. Astrophys 1998, 334, 505. [Google Scholar]
Ryden, B. Introduction to Cosmology, 2nd ed.; C: (Cambridge, 2017. [Google Scholar]
Sadoun, R.; Shlosman, I.; Choi, J.-H.; Romano-Diaz, E. Astrophys. J. 2016, 829. [Google Scholar]
Shampine, L.F.; Reichelt, M. SIAM. J. Sci. Comput. 1997, 18, 1. [Google Scholar]
Shapiro, S.L.; Teukolsky, S.A. Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects; Wiley-Interscience: New York, 1983. [Google Scholar]
Somerville, R.S.; Dave, R. Annu. Rev. Astron. Astrophys. 2015, 53, 51. [Google Scholar] [CrossRef]
Tielens, A.G.G.M. The Physics and Chemistry of the Interstellar Medium; Cambridge University Press: Cambridge, 2005. [Google Scholar]
Tinsley, B. M. Fundament. Cosmic Phys. 1980, 5, 287. [Google Scholar]
Tortora, C.; Hunt, L.K.; Ginolfi, M. Astron. Astrophys 2022, 657. [Google Scholar]
Weinberg, S. Cosmology; Oxford University Press: Oxford, 2008. [Google Scholar]
Woosley, S. E.; Janka, H.-T. Nat. Phys. 2005, 1, 147. [CrossRef]

1	Note that $M_{⊙} {yr}^{- 1} {Mpc}^{- 3} = 2.15 \cdot 10^{- 45} kg s^{- 1} m^{- 3}$ .
2	We emphasize the extremely sensitive dependence on the value of the Hubble constant $h_{70}^{- 8}$ and also on the ratio $r_{1}$ as with the upper limit (45) for $r_{1}$ the factor ${(1 + r_{1})}^{8} < 1.39$ .
3	To evaluate $F (\tilde{t})$ for large arguments $ν \tilde{t} ≫ 1$ (implying $μ \tilde{t} ≫ 1$ ) numerically, one can use $F (\tilde{t}) \approx (μ + ν - 1) {\tilde{t}}^{5 / 3} e^{- (μ - ν) \tilde{t}} / 2 ν$ .

Figure 1. (a) Gas (G), (b) stars (S) and (c) locked-in matter (L) compartments as well as (d) the formation rate

\overset{˚}{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

. Shown are five choices for

b_{0}

and

c_{0}

, while

β_{0} = 1

and

η = 10^{- 5}

. The numerical solution of Equations (7)–(9) exactly coincides with the analytical solution provided by Equations (19), (20), (21), and (24).

Figure 1. (a) Gas (G), (b) stars (S) and (c) locked-in matter (L) compartments as well as (d) the formation rate

\overset{˚}{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

. Shown are five choices for

b_{0}

and

c_{0}

, while

β_{0} = 1

and

η = 10^{- 5}

. The numerical solution of Equations (7)–(9) exactly coincides with the analytical solution provided by Equations (19), (20), (21), and (24).

Figure 2. Semi-logarithmic

\tilde{t} (z) / β_{0} τ_{f}

(black solid line) versus z according to Equation (39) with

τ_{f}

from Equation (41), Dashed green curve shows the approximant

\tilde{t} (z) / β_{0} τ_{f} = {(1 + z)}^{- 3 / 2}

, i.e., Equation (42).

Figure 2. Semi-logarithmic

\tilde{t} (z) / β_{0} τ_{f}

(black solid line) versus z according to Equation (39) with

τ_{f}

from Equation (41), Dashed green curve shows the approximant

\tilde{t} (z) / β_{0} τ_{f} = {(1 + z)}^{- 3 / 2}

, i.e., Equation (42).

Figure 3. GS-model with zero triggered star formation (

a_{0} = c_{0} = 0

) with

r_{1} = 0.001

and

{\tilde{β}}_{0} = 2

for reasons described in the text. (a)

ψ_{GS} (z)

and (b)

ρ_{*} (z)

obtained numerically (black solid line) and analytically (green dashed line, Equations (54) and (68), respectively). Symbols with error bars are data for

ψ (z)

and

ρ_{*} (z)

tabulated by Madau [27].

Figure 3. GS-model with zero triggered star formation (

a_{0} = c_{0} = 0

) with

r_{1} = 0.001

and

{\tilde{β}}_{0} = 2

for reasons described in the text. (a)

ψ_{GS} (z)

and (b)

ρ_{*} (z)

obtained numerically (black solid line) and analytically (green dashed line, Equations (54) and (68), respectively). Symbols with error bars are data for

ψ (z)

and

ρ_{*} (z)

tabulated by Madau [27].

Figure 4. The lower incomplete gamma function

γ (8 / 3, 5.37 {\tilde{β}}_{0})

, involved in the analytic expressions for

ρ_{GS}^{*} (z)

and

ρ_{GSL}^{*} (z)

, versus

{\tilde{β}}_{0}

. Particular values are

0.0487

0.21

0.88

1.40

, and

1.50

for

{\tilde{β}}_{0} = 0.1

0.2

0.5

1.0

, and 2, respectively.

Figure 4. The lower incomplete gamma function

γ (8 / 3, 5.37 {\tilde{β}}_{0})

, involved in the analytic expressions for

ρ_{GS}^{*} (z)

and

ρ_{GSL}^{*} (z)

, versus

{\tilde{β}}_{0}

. Particular values are

0.0487

0.21

0.88

1.40

, and

1.50

for

{\tilde{β}}_{0} = 0.1

0.2

0.5

1.0

, and 2, respectively.

Figure 8.

{log}_{10} [ρ_{GSL}^{*} (0) {kg}^{- 1} m^{3}]

for

r_{1} = 0

versus

{\tilde{β}}_{0}

and

{log}_{10} (r_{2})

. The reported values for

ρ_{*} (0)

are located between the two red contour lines. The

ρ_{GSL}^{*} (0)

had been evaluated using the last line of Equation (90). For

r_{2} = 1

, the expression (110) could be used. The values of

ρ_{GSL}^{*} (0)

at the uppermost value

r_{2} = 10^{3}

basically agree with the analytic expression

ρ_{GSL}^{*} (0) = 4.0 \cdot 10^{- 29} h_{7}^{8 / 3} {\tilde{β}}_{0}^{- 2 / 3}

m^{- 3}

, that follows from Equation (117).

Figure 8.

{log}_{10} [ρ_{GSL}^{*} (0) {kg}^{- 1} m^{3}]

for

r_{1} = 0

versus

{\tilde{β}}_{0}

and

{log}_{10} (r_{2})

. The reported values for

ρ_{*} (0)

are located between the two red contour lines. The

ρ_{GSL}^{*} (0)

had been evaluated using the last line of Equation (90). For

r_{2} = 1

, the expression (110) could be used. The values of

ρ_{GSL}^{*} (0)

at the uppermost value

r_{2} = 10^{3}

basically agree with the analytic expression

ρ_{GSL}^{*} (0) = 4.0 \cdot 10^{- 29} h_{7}^{8 / 3} {\tilde{β}}_{0}^{- 2 / 3}

m^{- 3}

, that follows from Equation (117).

Figure 9.

G (z)

S (z)

, and

L (z)

for

r_{1} = 0

r_{2} = 6.7

, and

{\tilde{β}}_{0} = 1.2

, according to Equations (118).

Figure 9.

G (z)

S (z)

, and

L (z)

for

r_{1} = 0

r_{2} = 6.7

, and

{\tilde{β}}_{0} = 1.2

, according to Equations (118).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. GSL-Compartmental Model for the Baryonic Matter Cycle

3. Negligible Triggered Star Formation

3.1. Stationary Ratios

3.2. Negligible Stellar Feedback

3.3. Negligible Stellar Evolution

4. Cosmic Star Formation History

4.1. Observational Constraints

5. Results for Model without Stellar Evolution

5.1. SFR Density

5.2. Stellar Density

5.3. Special Case: Negligible Stellar Feedback ( r 1 = 0 ) and Stellar Evolution ( r 2 = 0 )

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

6.1. Theoretical Results

6.2. Constraints on Parameters

6.2.1. Results for standard Hubble constant value h 70 = 1

6.2.2. Smaller Values of the Hubble Constant

6.3. Special Case of Negligible Stellar Feedback

7. Summary and Conclusions

Acknowledgments

Appendix A. Lambert Functions

References

MDPI Initiatives

Important Links

Subscribe

5.3. Special Case: Negligible Stellar Feedback ( $r_{1} = 0$ ) and Stellar Evolution ( $r_{2} = 0$ )

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