Dynamical Reason for a Cyclic Universe

Abstract

By analyzing the energy-momentum tensor and equations of state of ideal gas, scalar, spinor and vector potential in detail, we find that the total mass density of all matter is always positive, and the initial total pressure is negative. Under these conditions, by qualitatively analyzing the global behavior of the dynamical equation of cosmological model, we get the following results: (i) $K=1$, namely, the global spatial structure of the universe should be a three-dimensional sphere $S^3$; (ii) $0\le \Lambda <{10}^{-24}{\mathrm{ly}}^{-2}$, the cosmological constant should be zero or an infinitesimal; (iii) $a\left(t\right)>0$, the initial singularity of the universe is unreachable, and the evolution of the universe should be cyclic in time. Since the matter components considered are quite complete and the proof is very elementary and strict, these conclusions are quite reliable in logic and compatible with all observational data. Obviously, these conclusions will be very helpful to correct some popular misconceptions and bring great convenience to further research other problems in cosmology such as the properties of dark matter and dark energy. In addition, the macroscopic Lagrangian of fluid model is derived.

1. Introduction

In cosmology, we have two important constants to be determined. They are cosmic curvature K and cosmological constant $\Lambda$. Some characteristic parameters of the universe, such as the age $\mathcal{T}$, Hubble constant $H_0$, total mass density ${\Omega }_{\mathrm{tot}}$, etc., have been measured to high accuracy [1,2,3,4,5,6]. To determine the cosmic curvature K, the usual approach is to transform the Friedmann equation into an algebraic equation ${\Omega }_K\equiv K{\dot{a}}^{-2}={\Omega }_{\mathrm{tot}}-1$. In theory, $K=0,\pm 1$ can be judged by contrasting observational data ${\Omega }_{\mathrm{tot}}>1,=1$ or $<1$. From the observations we have ${\Omega }_K=-0.0020\pm 0.0047$, which is very close to the case of flat space. Considering measurement errors, it is hard to determine what type the space is. In fact, no matter what the case of spatial curvature is, for a young universe, it is easy to calculate that we always have ${\Omega }_{\mathrm{tot}}\approx 1$. So this criterion is rather ambiguous.

Cosmological constant $\Lambda$ has a dramatic history. Since Einstein introduced $\Lambda$ in 1917 to get a static and closed universe, the debate over whether $\Lambda$ is zero or not has been repeated many times [7,8,9]. Now dark matter and dark energy are attracting the attention of scientists around the world, becoming the hottest topic, challenging traditional standard models of particle and cosmology. The usual description of dark matter and dark energy uses the equation of state $P=w\rho$ and $w=w\left(a\right)$ or $w=w\left(z\right)$, and many specific models have been obtained by fitting the observed data [7,8,9,10,11,12,13,14]. However, the problem is far from being solved [15].

In [12], by introducing the potential function $V\left(a\right)$, the Friedmann equations with some known dark energy models are converted into Hamiltonian dynamics, and the evolving trajectory is analyzed to explain the accelerating expansion of the universe. The literatures [16,17,18,19,20,21,22] provide some similar discussion on specific gravitational sources. The nonlinear scalar fields were discussed in reference [16,17,18] to obtain the cyclic universe model. In [19], a set of precise cyclic solutions with ordinary dust and radiation are obtained, and in [20], the exact solutions of ghost and electromagnetic fields are derived. The quantized nonlinear spinor fields and trajectories are calculated and a cyclic solution $a\left(t\right)$ is solved in [21,22]. In [23,24], the authors use $({\Omega }_K,a)$ phase plane to discuss the dynamical behavior of the universe, and conclude that a cyclical universe is reasonable from a dynamical systems perspective, and requires in addition to standard cosmological assumptions, only two conditions: (i) the spatial sections must have positive spatial curvature $K=+1$, and (ii) the late time effective cosmological “constant” must decay fast enough as a function of the scale factor. Both of these conditions are consistent with all current observations to date. In 2008, M. Novello and S. E. Perez Bergliaffa reviewed the general features of nonsingular universes and cyclic universes, discussed the mechanisms behind the bounce, and analyzed examples of solutions that implement these mechanisms [25].

In [26], after analyzing the Planck Legacy 2018 release, the authors found enhanced gravitational lens amplitude in the power spectrum of the cosmic microwave background (CMB), which is different from the data of the standard $\Lambda$CDM model. The cosmic space is more like a giant balloon that is closed and curved. The main task of the Planck satellite is to detect the tiny fluctuation of the CMB temperature. The study of the fluctuation of the CMB temperature is the key to uncover the relevant cosmological models and parameters. This information defines the expansion, composition and origin of the cosmic large-scale structure. The authors use ‘closed universe’ to explain this anomalous effect. The spectra are now more inclined to a positive curvature greater than 99% confidence level. The positive curvature can explain the anomalous amplitude of the gravitational lens.

This paper is a development of [21,27]. In [21], the calculations showed the term $\Lambda a^4$ is a higher order infinitesimal as $a\to 0$, so its influence on the behavior of the early universe is negligible. Some nonlinear effects of the matter fields should be important for the early behavior. The calculations show that the nonlinear potential of spinors leads to a higher order term $a^{-6}$, which is equivalent to negative pressure and becomes dominant factor for the behavior of the early universe. In [27], we showed that the global behavior of Friedmann equation is mainly determined by two hypotheses of positive mass density and initial negative pressure. In order to verify the reasonability of the two hypotheses, some contents were added continually, and we obtain the present paper.

Under two ordinary hypotheses for functions of state, namely the total mass density of all matter is always positive, and the initial total pressure is negative, by qualitatively analyzing the global behavior of the dynamical equation of cosmological model, we get the following results: the global spatial structure of the universe should be a closed 3-d sphere $S^3$, the cosmological constant should be zero or a tiny number and the initial singularity of the universe is unreachable, and the evolution of universe should be cyclic in time. Obviously, these conclusions will be very helpful to correct some popular misconceptions and bring great convenience in further researching other problems in cosmology such as the properties of dark matter and dark energy.

2. Energy-Momentum Tensor of Matter

In order to understand the mystery and dynamically behavior of the universe, Energy-momentum tensor(EMT) and equation of state(EOS) of all kinds of matter are crucial factors. In this paper, we establish EMT and EOS of some elementary components according to credible theories in detail. Cosmology contains a variety of contents, so it is necessary to clarify the conventions and notations frequently used in the paper at first. The line element of the space-time is given by [28]

$$\begin{array}{c}\hfill d\mathbf{x}={\gamma }_{\mu }d{x}^{\mu }={\gamma }_a\delta X^a,\end{array}$$

(1)

where ${\gamma }^a$ and ${\gamma }^{\mu }$ are tetrad expressed by Dirac matrices

$$\begin{array}{c}\hfill {\gamma }^a=\left(\begin{array}{cc}0& {\tilde{\sigma }}^a\\ {\sigma }^a& 0\end{array}\right),{\gamma }^{\mu }=f_a^{\mu }{\gamma }^a,{\gamma }_{\mu }=f_{\mu }^a{\gamma }_a,\end{array}$$

(2)

which satisfies the $C{\ell }_{1,3}$ Clifford algebra

$$\begin{array}{c}\hfill {\gamma }_a{\gamma }_b+{\gamma }_b{\gamma }_a=2{\eta }_{ab},{\eta }_{ab}=\mathrm{diag}(1,-1,-1,-1),{\gamma }_{\mu }{\gamma }_{\nu }+{\gamma }_{\nu }{\gamma }_{\mu }=2g_{\mu \nu }.\end{array}$$

(3)

The Pauli matrices are expressed by

$${\sigma }^a=\left\{\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\right\},$$

(4)

$$\begin{array}{ccc}\hfill {\tilde{\sigma }}^a& =& ({\sigma }^0,-{\sigma }^1,-{\sigma }^2,-{\sigma }^3),{\sigma }^{\mu }=f_a^{\mu }{\sigma }^a,{\tilde{\sigma }}_{\mu }=f_{\mu }^a{\tilde{\sigma }}_a,\hfill \end{array}$$

(5)

where $a,\mu \in \{0,1,2,3\}$. We take $c=1$ as unit of speed. For index, we use the Greek characters for curvilinear coordinates, Latin characters for local Minkowski coordinates and $\{i,j,k,l\}$ for spatial indices. The Einstein’s summation is used.

Similarly to the case of metric $g_{\mu \nu }$, the definition of Ricci tensor can also differ by a negative sign. We take the definition as follows

$$\begin{array}{c}\hfill R_{\mu \nu }\equiv {\partial }_{\alpha }{\Gamma }_{\mu \nu }^{\alpha }-{\partial }_{\mu }{\Gamma }_{\nu \alpha }^{\alpha }-{\Gamma }_{\mu \beta }^{\alpha }{\Gamma }_{\nu \alpha }^{\beta }+{\Gamma }_{\mu \nu }^{\alpha }{\Gamma }_{\alpha \beta }^{\beta },R=g^{\mu \nu }{R}_{\mu \nu }.\end{array}$$

Since the Lagrangian of a compound physical system is a superposable real scalar, in cosmology the Lagrangian of main kinds of gravitational sources is generally given by

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{2\kappa }(R-2\Lambda )+{\mathcal{L}}_{\mathrm{m}},{\mathcal{L}}_{\mathrm{m}}={\mathcal{L}}_{\varphi }+{\mathcal{L}}_p+{\mathcal{L}}_{\psi }+{\mathcal{L}}_A+{\mathcal{L}}_{\Phi }+\cdots ,\end{array}$$

(6)

in which $\kappa =8\pi G$, $\Lambda$ is the cosmological constant, ${\mathcal{L}}_{\mathrm{m}}$ the total Lagrangian of all matter.

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\varphi }& =& \frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V\left(\varphi \right),V\sim \frac{m^2}{2k}{\varphi }^{2k}\left(\right|\varphi |\to \infty ),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathcal{L}}_p& =& -\sum _n{m}_n\sqrt{1-v_n^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n),v_n^2=g^{00}{\overline{g}}_{kl}\frac{d{X}_n^k}{dt}\frac{d{X}_n^l}{dt},\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\psi }& =& \sum _n\left(\Re \langle {\psi }_n^{+}{\alpha }^{\mu }{\widehat{p}}_{n|\mu }{\psi }_n\rangle -{\psi }_n^{+}{\Omega }_{\mu }{\widehat{S}}^{\mu }{\psi }_n-m_n{\psi }_n^{+}{\gamma }^0{\psi }_n+N_n\right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathcal{L}}_A& =& -\frac{1}{2}{\nabla }_{\mu }{A}_{\nu }{\nabla }^{\mu }{A}^{\nu },{\mathcal{L}}_{\Phi }=\frac{1}{2}({\nabla }_{\mu }{\Phi }_{\nu }{\nabla }^{\mu }{\Phi }^{\nu }-b^2{\Phi }_{\mu }{\Phi }^{\mu }),\hfill \end{array}$$

(10)

where $\varphi$ is the global slow-roll scalar field. ${\mathcal{L}}_p$ is the Lagrangian of ideal dust, whose statistical average is perfect fluid model. $v_n$ is the speed of n-th particle in usual sense. ${\mathcal{L}}_{\psi }$ is the Lagrangian for spinors with nonlinear potential $N_n$, electromagnetic potential $A^{\mu }$ and strong short distance interaction ${\Phi }^{\mu }$, in which the momentum operator

$$\begin{array}{c}\hfill {\widehat{p}}_k^{\mu }=i(\hslash {\partial }^{\mu }+{{\rm Y}}^{\mu })-e_k{A}^{\mu }-s_k{\Phi }^{\mu },e_k=\{0,\pm e\},s_k=\{0,s\},\end{array}$$

(11)

where ${\alpha }^{\mu }$ is current operator and ${\widehat{S}}^{\mu }$ spin operator. They are defined respectively by

$$\begin{array}{c}\hfill {\alpha }^{\mu }=\mathrm{diag}({\sigma }^{\mu },{\tilde{\sigma }}^{\mu }),{\widehat{S}}^{\mu }=\frac{1}{2}\mathrm{diag}({\sigma }^{\mu },-{\tilde{\sigma }}^{\mu }).\end{array}$$

(12)

${{\rm Y}}_{\mu }\in {\Lambda }^1$ is Keller connection and ${\Omega }^{\alpha }$ is Gu-Nester potential. They are calculated by [28,29,30,31].

$$\begin{array}{c}\hfill {{\rm Y}}_{\mu }=\frac{1}{2}{f}_a^{\nu }({\partial }_{\mu }{f}_{\nu }^a-{\partial }_{\nu }{f}_{\mu }^a),{\Omega }^{\alpha }=\frac{1}{2}{f}_d^{\alpha }{f}_a^{\mu }{f}_b^{\nu }{\partial }_{\mu }{f}_{\nu }^e{ϵ}^{abcd}{\eta }_{ce}.\end{array}$$

(13)

If the metric of the space-time can be orthogonalized, we have ${\Omega }_{\mu }\equiv 0$, so the spin-gravity coupling term ${\widehat{S}}_{\mu }{\Omega }^{\mu }=0$ in cosmology. $\Re \langle \rangle$ means taking real part of a complex number. In this paper we take the following simplest nonlinear potential as example to show its dynamical effects in cosmology

$$\begin{array}{c}\hfill N_n=\frac{1}{2}\tilde{w}{\check{\gamma }}_n^2,\check{\gamma }={\psi }_n^{+}{\gamma }^0{\psi }_n.\end{array}$$

(14)

The Hamiltonian formalism and classical mechanics can be clearly described only in the Gu’s natural coordinate system [32]

$$\begin{array}{c}\hfill d{\mathbf{x}}^2=g_{00}d{t}^2-{\overline{g}}_{kl}d{x}^{k}d{x}^l,\sqrt{g}=\sqrt{g_{00}}\sqrt{\overline{g}},\overline{g}=det\left({\overline{g}}_{kl}\right),\end{array}$$

(15)

where $d\tau =\sqrt{g_{00}}dt$ defines the Newton’s realistic cosmic time, which is different from the proper time of a particle $d{s}_k=\sqrt{1-v_k^2}d\tau$. The Dirac-$\delta$ is defined as

$$\begin{array}{c}\hfill {\int }_{{\mathbb{R}}^3}f\left(\overrightarrow{x}\right){\delta }^3(\overrightarrow{x}-\overrightarrow{X})\sqrt{\overline{g}}{d}^{3}x=f\left(\overrightarrow{X}\right),\forall f\in C_c^{\infty }\left({\mathbb{R}}^3\right),\end{array}$$

which means ${\delta }^3(\overrightarrow{x}-\overrightarrow{X})\propto {\left(\sqrt{\overline{g}}\right)}^{-1}$. Only scalar, spinor and vector fields can construct a proper Lagrangian [22], so (7)–(10) become the representation of all kinds of matter. In cosmology, to study them clearly is enough for theoretical analysis.

Variation of the Lagrangian (6) with respect to $g_{\mu \nu }$, we get Einstein’s field equation

$$\begin{array}{c}\hfill G^{\mu \nu }+\Lambda g^{\mu \nu }+\kappa T^{\mu \nu }=0,G^{\mu \nu }\equiv R^{\mu \nu }-\frac{1}{2}{g}^{\mu \nu }R=-\frac{\delta \left(R\sqrt{g}\right)}{\sqrt{g}\delta g_{\mu \nu }}.\end{array}$$

(16)

where $\frac{\delta }{\delta g_{\mu \nu }}$ is the Euler derivative, and $T^{\mu \nu }$ is EMT of all matter defined by

$$\begin{array}{ccc}\hfill T^{\mu \nu }& =& -2\frac{\delta \left({\mathcal{L}}_{\mathrm{m}}\sqrt{g}\right)}{\sqrt{g}\delta g_{\mu \nu }}=-2\left(\frac{\partial {\mathcal{L}}_{\mathrm{m}}}{\partial g_{\mu \nu }}-({\partial }_{\alpha }+{\Gamma }_{\alpha \gamma }^{\gamma })\frac{\partial {\mathcal{L}}_{\mathrm{m}}}{\partial \left({\partial }_{\alpha }{g}_{\mu \nu }\right)}\right)-g^{\mu \nu }{\mathcal{L}}_{\mathrm{m}}.\hfill \end{array}$$

(17)

By calculation, we have [22,31,33]

$$\begin{array}{ccc}\hfill T^{\mu \nu }& =& T_{\varphi }^{\mu \nu }+T_p^{\mu \nu }+T_{\psi }^{\mu \nu }+T_A^{\mu \nu }+T_{\Phi }^{\mu \nu }+\cdots ,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill T_{\varphi }^{\mu \nu }& =& {\partial }^{\mu }\varphi {\partial }^{\nu }\varphi -g^{\mu \nu }[\frac{1}{2}{g}^{\alpha \beta }{\partial }_{\alpha }\varphi {\partial }_{\beta }\varphi -V\left(\varphi \right)],\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill T_p^{\mu \nu }& =& \sum _n{m}_n{u}_n^{\mu }{u}_n^{\nu }\sqrt{1-v_n^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n),u_n^{\mu }=\frac{d}{ds}{X}_n^{\mu },\hfill \end{array}$$

(20)

and

$$\begin{array}{ccc}\hfill T_{\psi }^{\mu \nu }& =& \sum _n\left(\frac{1}{2}\Re \langle {\psi }_n^{+}({\alpha }^{\mu }{\widehat{p}}_n^{\nu }+{\alpha }^{\nu }{\widehat{p}}_n^{\mu }-2S_{ab}^{\mu \nu }{f}_c^{\alpha }{\eta }^{bc}{\alpha }^a{\widehat{p}}_{n|\alpha }){\psi }_n\rangle +g^{\mu \nu }{N}_n\right),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill T_A^{\mu \nu }& =& -{\nabla }^{\mu }{A}^{\alpha }{\nabla }^{\nu }{A}_{\alpha }+\frac{1}{2}{g}^{\mu \nu }{\nabla }_{\beta }{A}_{\alpha }{\nabla }^{\beta }{A}^{\alpha },\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill T_{\Phi }^{\mu \nu }& =& {\nabla }^{\mu }{\Phi }^{\alpha }{\nabla }^{\nu }{\Phi }_{\alpha }-\frac{1}{2}{g}^{\mu \nu }{\nabla }_{\beta }{\Phi }_{\alpha }{\nabla }^{\beta }{\Phi }^{\alpha }+\frac{1}{2}{b}^2({\Phi }^{\mu }{\Phi }^{\nu }+g^{\mu \nu }{\Phi }_{\alpha }{\Phi }^{\alpha }).\hfill \end{array}$$

(23)

In (21) $S_{ab}^{\mu \nu }=f_a^{\{\mu }{f}_b^{\nu \}}\mathrm{sign}(a-b)=-S_{ba}^{\mu \nu }$ for $LU$ decomposition of metric.

For classical approximation of (21), we have [30,31]

$$\begin{array}{c}\hfill {\psi }^{+}{\alpha }^{\nu }\psi \to u^{\nu }\sqrt{1-v^2}{\delta }^3(\overrightarrow{x}-\overrightarrow{X}),{\widehat{p}}^{\mu }\psi \to m{u}^{\mu }\psi ,N\to w\sqrt{1-v^2}{\delta }^3(\overrightarrow{x}-\overrightarrow{X}).\end{array}$$

(24)

Substituting (24) into (21) and noticing $S_{ab}^{\mu \nu }=-S_{ba}^{\mu \nu }$, we have

$$\begin{array}{c}\hfill \Re \langle {\psi }^{+}{S}_{ab}^{\mu \nu }{f}_c^{\alpha }{\eta }^{cb}{\alpha }^a{\widehat{p}}_{\alpha }\psi \rangle \to m{S}_{ab}^{\mu \nu }{u}^a{u}^b\sqrt{1-v^2}{\delta }^3(\overrightarrow{x}-\overrightarrow{X})=0.\end{array}$$

(25)

So, we get the classical approximation for EMT of dark spinor with self-interaction

$$\begin{array}{c}\hfill T_{\psi }^{\mu \nu }\to \sum _n\left(m_n{u}_n^{\mu }{u}_n^{\nu }+w_n{g}^{\mu \nu }\right)\sqrt{1-v_n^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n),\end{array}$$

(26)

$w_n\propto {\left(\sqrt{\overline{g}}\right)}^{-1}>0$ acts like negative pressure and takes the place of cosmological constant $\Lambda$ in Equation (16). If $w_n=0$, the spinor moves along geodesic like a mass-point. Some researchers studied varying cosmological ‘constant’ [34,35,36,37,38,39,40,41], but a directly varying $\Lambda \left(t\right)$ violates covariance of dynamics. The following calculation discloses that the main physical origin of a covariant $\Lambda \left(a\right)$ may be just the nonlinear potential ${\sum }_n{N}_n\sim a^{-6}$.

For energy-momentum tensor, we have the following useful theorem, which means the energy-momentum of any independent system is conserved respectively.

Theorem 1.

Assume matter consists of two subsystems I and II, namely ${\mathcal{L}}_{\mathrm{m}}={\mathcal{L}}_I\left(\varphi \right)+{\mathcal{L}}_{II}\left(\psi \right)$, then we have

$$\begin{array}{c}\hfill T^{\mu \nu }=T_I^{\mu \nu }+T_{II}^{\mu \nu }.\end{array}$$

(27)

If the subsystems I and II have not interaction with each other, namely,

$$\begin{array}{c}\hfill \frac{\delta }{\delta \psi }{\mathcal{L}}_I\left(\varphi \right)=\frac{\delta }{\delta \varphi }{\mathcal{L}}_{II}\left(\psi \right)=0,\end{array}$$

(28)

then the two subsystems have independent energy-momentum conservation laws respectively,

$$\begin{array}{c}\hfill T_{I;\nu }^{\mu \nu }=0,T_{II;\nu }^{\mu \nu }=0.\end{array}$$

(29)

Proof. 

By the definition of EMT (17), the variation $\delta g:{\mathcal{L}}_{\mathrm{m}}\mapsto T^{\mu \nu }$ is a linear mapping, so (27) holds. By (28), the variables $\varphi$ and $\psi$ have decoupling dynamic equations. Since the dynamics of variables is sufficient condition of energy-momentum conservation law, we can derive $T_{I;\nu }^{\mu \nu }=0$ from dynamic equation of $\varphi$, and $T_{II;\nu }^{\mu \nu }=0$ from dynamic equation of $\psi$ independently, so (29) holds. The proof is completed. □

By Theorem 1, we find in (7) the slow-roll scalar $\varphi$ have not interaction with other components of matter, so we have $T_{\varphi ;\nu }^{\mu \nu }=0$. In (8), each particle of ideal gas has not interaction with other components, so each particle satisfies energy-momentum conservation law. For k-th particle, by (20) we have

$$\begin{array}{ccc}\hfill 0& =& {\int }_{{\mathbb{R}}^3}{T}_{k;\nu }^{\mu \nu }\sqrt{g}{d}^{3}x={\int }_{{\mathbb{R}}^3}\left({\partial }_{\nu }\left(T_k^{\mu \nu }\sqrt{g}\right)+{\Gamma }_{\alpha \beta }^{\mu }{T}_k^{\alpha \beta }\sqrt{g}\right)d^{3}x\hfill \\ & \to & \frac{d}{dt}\left(\left(m_k{u}_k^{\mu }{u}_k^0\right)\sqrt{1-v_k^2}\sqrt{g_{00}}\right)+{\Gamma }_{\alpha \beta }^{\mu }\left(m_k{u}_k^{\alpha }{u}_k^{\beta }\right)\sqrt{1-v_k^2}\sqrt{g_{00}}.\hfill \end{array}$$

(30)

By $\frac{ds}{dt}=\sqrt{1-v_k^2}\sqrt{g_{00}}={\left(u_k^0\right)}^{-1}$, we get geodesic equation for k-th particle

$$\begin{array}{c}\hfill \frac{d}{ds}{u}_k^{\mu }+{\Gamma }_{\alpha \beta }^{\mu }{u}_k^{\alpha }{u}_k^{\beta }=0.\end{array}$$

(31)

Therefore, for ideal mass point, energy momentum conservation law is equivalent to dynamics. For perfect fluid model, $T_{;\nu }^{\mu \nu }=0$ also gives dynamical equations, but an additional equation of state $P=P\left(\rho \right)$ is needed and some information of the system is lost in the simplification of fluid.

If $e_k\ne 0$ or $s_k\ne 0$, spinor ${\psi }_k$ interacts with $A^{\mu }$ or ${\Phi }^{\mu }$, so we should take it as one system. For an isolated particle, the classical approximation of static $A^{\mu }$ or ${\Phi }^{\mu }$ of the particle can be treated as an additional mass $\delta m_k$ due to linearity of $A^{\mu }$ and ${\Phi }^{\mu }$. We calculate $\delta m_k$ in the next section. The propagating $A^{\mu }$ is photon, which can be treated as massless particles. Thus except the global scalar $\varphi$, the classical approximation of EMT for other fields is usually depicted by

$$\begin{array}{c}\hfill T^{\mu \nu }\to (\rho +P)U^{\mu }{U}^{\nu }+(W-P)g^{\mu \nu }.\end{array}$$

(32)

The statistical expectation of nonlinear potentials is function of state $W\left(\rho \right)$ which acts like negative pressure. By (14) we find $W_{\psi }\sim {\left({\delta }^3\left(\overrightarrow{x}\right)\right)}^2$ is scale dependent, so we have $W_{\psi }\propto \rho {\left(\sqrt{\overline{g}}\right)}^{-1}\propto {\left(\sqrt{\overline{g}}\right)}^{-2}$.

The following calculation discloses that the function of state $W\left(\rho \right)$ may be partly the physical origin of cosmological ‘constant’ $\Lambda$. By the energy-momentum conservation law $T_{;\nu }^{\mu \nu }=0$ of (32), we have equation of motion for dark spinor gas

$$\begin{array}{c}\hfill (\rho +P)U^{\nu }{U}_{;\nu }^{\mu }=(g^{\mu \nu }-U^{\mu }{U}^{\nu }){\partial }_{\nu }(P-W).\end{array}$$

(33)

The Equation (33) is non-geodesic, which is quite different from the geodesic equation $U^{\nu }{U}_{;\nu }^{\mu }=0$ for ordinary particles, thus the dark halo in a galaxy is automatically separated from ordinary matter during galaxy formation [31,42].

3. Equation of State in Cosmology

3.1. Basic Relations and Assumptions for EMT

In cosmology, we have Friedmann-Lemaitre-Robertson-Walker(FLRW) metric

$$d{s}^2=d{\tau }^2-a{\left(\tau \right)}^2\left(d{r}^2+\mathcal{S}{\left(r\right)}^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2)\right),$$

(34)

in which $\mathcal{S}\left(r\right)=sinr,r,sinhr$ correspond to $K=1,0,-1$ respectively. By symmetry we get

$$\begin{array}{ccc}\hfill G_{\tau }^{\tau }& =& -\frac{3}{a^2}(a^{\prime 2}+K),R=\frac{6}{a^2}(a{a}^{″}+a^{\prime 2}+K),\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill G_r^r& =& G_{\theta }^{\theta }=G_{\phi }^{\phi }=-\frac{1}{a^2}(2a{a}^{″}+a^{\prime 2}+K),a^{\prime }=\frac{d}{d\tau }a\left(\tau \right).\hfill \end{array}$$

(36)

For all gravitational sources ${\mathcal{L}}_{\mathrm{m}}$ we have EMT in traditional form

$$\begin{array}{c}\hfill T_{\nu }^{\mu }=\mathrm{diag}({\rho }_{\mathrm{m}},-P_{\mathrm{m}},-P_{\mathrm{m}},-P_{\mathrm{m}}).\end{array}$$

(37)

The dynamics of cosmology (16) and $T_{;\nu }^{\mu \nu }=0$ are overdetermined due to high symmetry of FLRW metric. It is easy to check, among all dynamical equations, only the following equations are independent equations, but the spatial components of the Einstein’s field equations can be derived from these equations. So we have dynamics for FLRW universe as

$$\begin{array}{ccc}\hfill G_{\tau }^{\tau }+\Lambda +\kappa T_{\tau }^{\tau }=0& \iff & a^{\prime 2}=\frac{1}{3}\Lambda a^2-K+\frac{\kappa }{3}{\rho }_{\mathrm{m}}{a}^2.\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill T_{i;\nu }^{\mu \nu }=0& \iff & d\left({\rho }_i{a}^3\right)=-P_{i}d{a}^3,(\forall i),\hfill \end{array}$$

(39)

where index i means independent subsystem. All other equations can be derived from the above equations. Together with initial values and equations of state $P_i=P_i\left({\rho }_i\right)$, the above equations are closed and the solution is uniquely determined.

The analysis in the next section shows that the asymptotic property of P as $a\to +0$ is crucial for the fate of the universe, so we give detailed discussion on this problem. By (39) we have energy conservation law for all matter

$$\begin{array}{c}\hfill d\left({\rho }_{\mathrm{m}}{a}^3\right)=-P_{\mathrm{m}}d{a}^3,{\rho }_{\mathrm{m}}=\sum _i{\rho }_i,P_{\mathrm{m}}=\sum _i{P}_i.\end{array}$$

(40)

By $G_r^r+\Lambda +\kappa T_r^r=0$ or derivative of (38) together with (40), we get second order equation for $a\left(\tau \right)$ as

$$\begin{array}{c}\hfill a^{″}=\frac{1}{6}a[2\Lambda -\kappa ({\rho }_{\mathrm{m}}+3P_{\mathrm{m}})].\end{array}$$

(41)

If $a^{″}>0$ as $a\to +0$, then we have

$$P_{\mathrm{m}}<\frac{1}{3\kappa }(2\Lambda -\kappa {\rho }_{\mathrm{m}}),(a\to +0).$$

(42)

Since ${\rho }_{\mathrm{m}}\to \infty$ as $a\to +0$, by (42) we certainly have $P_{\mathrm{m}}<0$ as $a\to +0$. Therefore, the existence of inflation in the early universe implies the initial negative pressure [11,12,13,14].

According to the above analysis, we have the following conclusions for total EMT (37) in general:

$1^{\circ }$ The total mass-energy density is always positive, namely

$${\rho }_{\mathrm{m}}>0,(\forall a>+0).$$

(43)

$2^{\circ }$ The total pressure $P_{\mathrm{m}}<0$ when the universe is small enough, namely

$$P_{\mathrm{m}}<0,(a\to +0).$$

(44)

These two conditions are the fundamental hypotheses of the following dynamical analysis, which are decisive for the evolution of the universe. Next, we give more discussions for the reasonability of the assumptions according to EOS of various materials.

3.2. Equation of State of Ideal gases

At first, we examine the ideal gases. To solve the geodesic of the particles, we prove a useful theorem.

Theorem 2.

Assume the line element of a manifold has the following form

$$\begin{array}{c}\hfill d{s}^2=a{\left(t\right)}^{2}d{t}^2+g_{\mu \nu }\left(t\right)d{x}^{\mu }d{x}^{\nu },\mu ,\nu \in \{1,2,\cdots ,n\},\end{array}$$

(45)

then we have the first integral for the geodesic in the manifold,

$$\begin{array}{c}\hfill \frac{d}{ds}t=\frac{1}{a}\sqrt{1-g^{\mu \nu }{b}_{\mu }{b}_{\nu }},\frac{d}{ds}{x}^{\mu }=g^{\mu \nu }{b}_{\nu },\end{array}$$

(46)

where $b_{\mu }$ are integral constants.

Proof. 

The nonzero Christoffel symbols are ${\Gamma }_{tt}^t$, ${\Gamma }_{\mu \nu }^t$ and

$$\begin{array}{c}\hfill {\Gamma }_{t\nu }^{\mu }={\Gamma }_{\nu t}^{\mu }=\frac{1}{2}{g}^{\mu \alpha }{g}_{\alpha \nu }^{\prime },(\mu ,\nu \ne 0),\end{array}$$

(47)

where $g_{\alpha \nu }^{\prime }=\frac{d}{dt}{g}_{\alpha \nu }$, $g^{\mu \alpha }{g}_{\alpha \nu }={\delta }_{\nu }^{\mu }$. Since $\frac{d}{ds}t$ can be directly derived from (45), it is unnecessary to calculate ${\Gamma }_{tt}^t$ and ${\Gamma }_{\mu \nu }^t$. Denoting ${\dot{x}}^{\mu }=\frac{d}{ds}{x}^{\mu }$, the geodesic equation of $x^{\mu }$ is given by

$$\begin{array}{c}\hfill {\ddot{x}}^{\mu }+2{\Gamma }_{0\nu }^{\mu }\dot{t}{\dot{x}}^{\nu }={\ddot{x}}^{\mu }+g^{\mu \alpha }{g}_{\alpha \nu }^{\prime }\dot{t}{\dot{x}}^{\nu }=0.\end{array}$$

(48)

Multiplying $g_{\mu \beta }{\dot{t}}^{-1}$ on both sides we get

$$\begin{array}{c}\hfill g_{\mu \beta }\frac{d}{dt}{\dot{x}}^{\mu }+\frac{d}{dt}\left(g_{\beta \nu }\right){\dot{x}}^{\nu }=\frac{d}{dt}\left(g_{\mu \beta }{\dot{x}}^{\mu }\right)=0.\end{array}$$

(49)

The solution is given by

$$\begin{array}{c}\hfill g_{\mu \nu }{\dot{x}}^{\nu }=b_{\mu },{\dot{x}}^{\mu }=g^{\mu \nu }{b}_{\nu },\end{array}$$

(50)

where $b_{\nu }$ is integral constant. By (45) and (50) we get

$$\begin{array}{c}\hfill {\dot{t}}^2=a^{-2}(1-g_{\mu \nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu })=a^{-2}(1-g^{\mu \nu }{b}_{\mu }{b}_{\nu }),\end{array}$$

(51)

so (46) holds. The proof is completed. □

In the FLRW space-time (34), if a particle moves along r coordinate, then the first integral of its geodesic in this subspace-time $(\tau ,r)$ can be solved as

$$\begin{array}{c}\hfill \frac{dr}{ds}=\frac{C}{a^2},\frac{d\tau }{ds}=\frac{1}{a}\sqrt{a^2+C^2},\end{array}$$

(52)

where C is a constant determined by initial data. By (52) we get the drifting speed of n-th particle in usual sense

$$\begin{array}{c}\hfill v_n\equiv \frac{adr}{d\tau }=\frac{b_n}{\sqrt{a^2+b_n^2}},\sqrt{1-v_n^2}=\frac{a}{\sqrt{a^2+b_n^2}}.\end{array}$$

(53)

So the momentum of the particle satisfies

$$\begin{array}{c}\hfill p\left(\tau \right)a\left(\tau \right)=p\left({\tau }_0\right)a\left({\tau }_0\right)=C,p=\frac{m_n{v}_n}{\sqrt{1-v_n^2}},\end{array}$$

(54)

where $m_n$ is the proper mass of the particle. By the symmetry of the space-time, (54) holds for all particles. For the massless photons, it is easy to check that the wavelength $\lambda \left(\tau \right)\propto a\left(\tau \right)$, so their momentum p also satisfy (54); that is to say, (54) is also suitable for photon gas.

Now we establish the relation between the state functions $({\rho }_p,P_p)$ of the ideal gases and scale factor a. By (54) we have $p_n^2=C_n{a}^{-2}$, where $C_n$ are constants determined by initial data at $\tau ={\tau }_0$. Then on the one hand, for all particles we have the average square momentum as

$$\begin{array}{c}\hfill {\overline{p}}^2=\frac{C_0}{a^2},\end{array}$$

(55)

where $C_0$ is a constant only determined by initial data at ${\tau }_0$.

On the other hand, the relativistic relation between momentum p and the kinetic energy $\mathcal{K}$ is given by

$$\begin{array}{c}\hfill p^2=\mathcal{K}(\mathcal{K}+2m),\mathcal{K}=\frac{m}{\sqrt{1-v^2}}-m,\end{array}$$

(56)

so ${\overline{p}}^2$ can be also calculated by statistics. Assuming the distribution of kinetic energy $\mathcal{K}$ of the ideal gases is given by $d\mathcal{P}=\mathcal{F}\left(\mathcal{K}\right)d\mathcal{K}$. Then it satisfies the equation of moments

$$\begin{array}{c}\hfill {\int }_0^{\infty }\mathcal{F}\left(\mathcal{K}\right)d\mathcal{K}=1,{\int }_0^{\infty }\mathcal{K}d\mathcal{P}=\frac{3}{2}kT,{\int }_0^{\infty }{\mathcal{K}}^{2}d\mathcal{P}=\frac{3}{2\sigma }{\left(kT\right)}^2,\end{array}$$

(57)

where the second formula is the definition of temperature T of the gases, $k=1.38\times {10}^{-23}$ J/K is the Boltzmann constant, $\sigma \sim \frac{2}{5}$ is a constant reflecting the distribution function of particles. In the case of Maxwell distribution [31], we have $\sigma =\frac{2}{5}$.

By (56) and the moments (57) we have

$$\begin{array}{c}\hfill {\overline{p}}^2=\frac{1}{N}\sum _{n=1}^N{\int }_0^{\infty }{p}_n^2\mathcal{F}\left({\mathcal{K}}_n\right)d{\mathcal{K}}_n=\frac{3}{2\sigma }kT(kT+2\sigma \overline{m}),\end{array}$$

(58)

where $\overline{m}=\frac{1}{N}{\sum }_n{m}_n$ is average proper mass of all particles. Comparing (58) with (55), we get the relation $T\left(a\right)$,

$$\begin{array}{c}\hfill kT=\frac{\sigma \overline{m}}{a}(\sqrt{a^2+b^2}-a),a=\frac{\sigma b\overline{m}}{\sqrt{kT(kT+2\sigma \overline{m})}},\end{array}$$

(59)

where $b>0$ is a constant determined by initial data. We take a as independent variable in statistical calculation.

In microscopic view, the energy momentum tensor of ideal gas is given by (20). The average energy-momentum tensor takes the form (37). By (20), in average sense we have [31]

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\rho }_p\equiv T_{0,p}^0=\sum _n{m}_n{g}_{00}{u}_n^0{u}_n^0\sqrt{1-v_n^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n)=\sum _n\frac{m_n}{\sqrt{1-v_n^2}}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n),\hfill \\ P_p\equiv -\frac{1}{3}{T}_{k,p}^k=-\frac{1}{3}\sum _n{m}_n{g}_{kl}{u}_n^k{u}_n^l\sqrt{1-v_n^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n)=\frac{1}{3}\sum _n\frac{m_n{v}_n^2}{\sqrt{1-v_n^2}}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n).\hfill \end{array}\right.\end{array}$$

(60)

By (20) and (59) we have

$$\begin{array}{ccc}\hfill {\rho }_p& =& \frac{1}{V}\int \sum _{X_n\in V}\frac{m_n}{\sqrt{1-v_n^2}}d\mathcal{P}=\frac{1}{a^3\Omega }\sum _{X_n\in \Omega }\int (m_n+{\mathcal{K}}_n)d\mathcal{P}\hfill \\ & =& \frac{1}{a^3\Omega }\sum _{X_n\in \Omega }(\overline{m}+\frac{3}{2}kT)=\frac{\varrho }{a^3}\left(1+\frac{3\sigma }{2a}(\sqrt{a^2+b^2}-a)\right),\hfill \end{array}$$

(61)

where $V=a^3\Omega$, and $\varrho =\frac{1}{\Omega }{\sum }_{X_n\in \Omega }{m}_n$ is the angular density of proper mass, which is a constant. Substituting (61) into energy conservation law (39) we get

$$\begin{array}{c}\hfill P_p=\frac{\sigma \varrho b^2}{2a^4\sqrt{a^2+b^2}}=\frac{NkT}{V}\left(1-\frac{kT}{2(\sigma \overline{m}+kT)}\right).\end{array}$$

(62)

While $a\to +0$ we find ${\rho }_p\sim P_p\to C{a}^{-4}$. Ideal gases cannot provide negative pressure.

In the above derivation, the FLRW metric (34) is used as piston-cylinder system to drive the ideal gas, and we can prove that the elastic collisions of particles have no influence on these results [31]. So the results are actually valid in general cases. By (59) we have relation

$$\begin{array}{c}\hfill a=\frac{\sigma b}{\sqrt{J(J+2\sigma )}},J\equiv \frac{kT}{\overline{m}{c}^2}=\frac{\sigma }{a}(\sqrt{a^2+b^2}-a),\end{array}$$

(63)

where J is dimensionless temperature. The above results conclude the following theorem.

Theorem 3.

For relativistic ideal gas, we have the EOS as

$$\begin{array}{c}\hfill P_p=\frac{\overline{\rho }kT}{\overline{m}{c}^2}\left(1-\frac{kT}{2(\sigma \overline{m}{c}^2+kT)}\right),\overline{m}\equiv \frac{1}{N}\sum _{n=1}^N{m}_n,\overline{\rho }\equiv \frac{N\overline{m}}{V}\propto \frac{1}{a^3}.\end{array}$$

(64)

For adiabatic process the functions of state satisfy parameter equation

$$\begin{array}{c}\hfill \overline{\rho }={\rho }_0{\left[J(J+2\sigma )\right]}^{\frac{3}{2}},{\rho }_p=\overline{\rho }(1+\frac{3}{2}J),P_p=\overline{\rho }\frac{J(J+2\sigma )}{2(J+\sigma )},\end{array}$$

(65)

where J acts as independent variable, ${\rho }_0=\varrho {\left(\sigma b\right)}^{-3}$ is a constant of density dimension.

The above derivation is compatible with relativity and includes the driving effect of gravity. In the case of low temperature, Theorem 3 gives the EOS for the adiabatic monatomic gas

$$\begin{array}{c}\hfill P\dot{=}\rho J\dot{=}{P}_0{\rho }^{\frac{5}{3}},(J\ll 1,\mathrm{or}kT\ll \overline{m}),\end{array}$$

(66)

which is identical to the empirical law of thermodynamics. Letting $J\to \infty$ or $\overline{m}\to 0$, we get $P\to \frac{1}{3}\rho$ and the Stefan-Boltzmann’s law $\rho \propto T^4$. Thus the above results are automatically suitable for photons, and the Stefan-Boltzmann’s law is also valid for the ultra-relativistic particles. In general relativity, all processes occur automatically, and ${\rho }_0$ is independent of any practical process. An EOS omitting the driving effect of gravity is invalid in general relativity.

3.3. Asymptotic Behavior of Scalar Field in the Early Universe

For scalar field $\varphi$, by (19) we have

$$\begin{array}{c}\hfill {\mathcal{L}}_{\varphi }=\frac{1}{2}{\varphi }^{\prime 2}-V,{\rho }_{\varphi }=\frac{1}{2}{\varphi }^{\prime 2}+V,P_{\varphi }=\frac{1}{2}{\varphi }^{\prime 2}-V.\end{array}$$

(67)

In the case $a\to 0$, we have $\left|\varphi \right|\to \infty$, so we only need to consider the highest order term in $V\left(\varphi \right)$. For potential $V=\frac{1}{2k}{m}^2{\varphi }^{2k},(k\ge 1)$, by variation with respect to $\varphi$ we get dynamical equation

$$\begin{array}{c}\hfill {\varphi }^{″}+\frac{3a^{\prime }}{a}{\varphi }^{\prime }+m^2{\varphi }^{2k-1}=0.\end{array}$$

(68)

(68) can be also derived by substituting (67) into energy conservation law (39). If the universe has an initial singularity, we set $a\left(0\right)=0$. For low temperature particles, we have ${\rho }_p{a}^3\sim \mathrm{const}$. By Friedmann Equation (38), we find $a\ge a_0\sqrt[3]{{\tau }^2}$ while $\tau \to +0$. So we assume in general

$$\begin{array}{c}\hfill a\to a_0{\tau }^j,\varphi \to \frac{C_0}{{\tau }^i},(j\ge \frac{2}{3},i>0).\end{array}$$

(69)

About the range of j, Barrow made some analysis in [43]. In the mode of accelerating expansion $a^{″}>0$, we should have $j>1$. Substituting (69) into (68) we get solution in the case of $k>1$,

$$\begin{array}{ccc}\hfill i& =& \frac{1}{k-1},m^2=\frac{3j(k-1)-k}{C_0^{2(k-1)}{(k-1)}^2},j>\frac{k}{3(k-1)};\hfill \end{array}$$

(70)

$$\begin{array}{ccc}\hfill {\rho }_{\varphi }& \to & \frac{3j{C}_0^2}{2k(k-1)}{\tau }^{-\frac{2k}{k-1}},P_{\varphi }\to -\frac{3j(k-1)-2k}{2k{(k-1)}^2}{C}_0^2{\tau }^{-\frac{2k}{k-1}}.\hfill \end{array}$$

(71)

Substituting $a\to a_0{\tau }^j$ into (71), we get

$$\begin{array}{c}\hfill P_{\varphi }\to \left(\frac{2k}{3(k-1)}-j\right)\frac{C^2}{a^n},(2<n<6,a\to +0).\end{array}$$

(72)

By (72), we find the nonlinear scalar $\varphi$ provides initial negative pressure if $j>\frac{2k}{3(k-1)}$. We have $j>\frac{4}{3}$ if $k=2$, $j>1$ if $k=3$ and $j>\frac{2}{3}$ if $k\to \infty$.

In the case of $k=1$, (68) is a linear equation of $\varphi$, and its solution is an oscillating function of $\tau$. For asymptotic solution of (68) as $\tau \to +0$, $|{\varphi }^{\prime }|$ is a higher order infinity than $\left|\varphi \right|$, so we certainly have ${\varphi }^{\prime 2}>m^2{\varphi }^2$ and then $P>0$. The linear scalar $\varphi$ cannot provides initial negative pressure. However, all asymptotic solutions are similar to the following one with $j=\frac{4}{3}$,

$$\begin{array}{c}\hfill \varphi \to \frac{A}{{\tau }^3}[sin(m\tau +\chi )-m\tau cos(m\tau +\chi )],\end{array}$$

(73)

where $A,0<\chi <\frac{\pi }{2}$ are constants. The corresponding pressure is given by

$$\begin{array}{c}\hfill P\to \frac{9A^2}{{\tau }^8}\left[{sin}^2(\chi +m\tau )-m\tau sin\left[2(\chi +m\tau )\right]+\mathrm{O}\left({\tau }^2\right)\right],\end{array}$$

(74)

which is also an oscillating function with $P\to +\infty (\tau \to +0)$. Of course, P becomes negative when $\tau$ is large enough. Substituting the above results into (67) we get

$$\begin{array}{c}\hfill {\rho }_{\varphi }\dot{=}{P}_{\varphi }\to \frac{C_1^2}{{\tau }^{6j}}=\frac{C_2^2}{a^6}>0,(a\to +0).\end{array}$$

(75)

3.4. Equation of State of Spinor Gas

The above calculation shows that perfect fluid and scalar model cannot provide negative pressure. Now we examine the dark spinor with nonlinear potential. In this case, the classical approximation of spinor EMT is given by (26). The spinor moves approximately along geodesic, i.e., for nonlinear spinors the equivalence principle holds approximately, so the functions of state $({\rho }_{\psi },P_{\psi })$ should approximately equal $({\rho }_p,P_p)$ given by (61) and (62). We only need to derive the function $W_{\psi }$. By (26), the definition of $W_{\psi }$ in microscopic view reads

$$\begin{array}{c}\hfill W_{\psi }\equiv \frac{1}{V}{\int }_V\sum _n{w}_n{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_n)\sqrt{1-v_n^2}dV=\frac{1}{V}\sum _{X_n\in V}{w}_n\sqrt{1-v_n^2},\end{array}$$

(76)

The statistical expectation of $\sqrt{1-v_n^2}$ cannot be directly calculated according to moments (57). By (53), we get

$$\begin{array}{c}\hfill \sqrt{1-v_n^2}=\frac{b_n}{\sqrt{a^2+b_n^2}}=\frac{d}{da}\frac{a}{\sqrt{1-v_n^2}}=\frac{1}{m_n}\frac{d}{da}\left[a(m_n+{\mathcal{K}}_n)\right].\end{array}$$

(77)

Substituting it into (76) and using (59), we get expectation value

$$\begin{array}{ccc}\hfill W_{\psi }& =& \frac{1}{V}\sum _{X_n\in V}\int w_n\sqrt{1-v_n^2}d\mathcal{P}=\frac{1}{V}\sum _{X_n\in V}\frac{w_n}{m_n}\frac{d}{da}\int a(m_n+{\mathcal{K}}_n)d\mathcal{P}\hfill \\ & =& \frac{1}{V}\sum _{X_n\in V}\frac{w_n}{m_n}\frac{d}{da}\left(a(m_n+\frac{3}{2}kT)\right)\hfill \\ & =& \frac{1}{V}\sum _{X_n\in V}\left(w_n-\frac{3w_n}{2m_n}\sigma \overline{m}(1-\frac{a}{\sqrt{a^2+b^2}})\right)\hfill \end{array}$$

(78)

$$\begin{array}{ccc}& =& \overline{\rho }\left(\frac{\overline{w}}{\overline{m}}-\frac{3\mu \sigma kT}{2(\sigma \overline{m}+kT)}\right),\hfill \end{array}$$

(79)

where the average parameters are defined by

$$\begin{array}{c}\hfill \overline{w}=\frac{1}{N}\sum _{n=1}^N{w}_n\to \frac{C_1}{a^3},\mu =\frac{1}{N}\sum _{n=1}^N\frac{w_n}{m_n}\to \frac{C_2}{a^3},(a\to +0).\end{array}$$

(80)

By (63) and (65) we obtain

$$\begin{array}{c}\hfill W_{\psi }=\overline{\rho }\left(\frac{\overline{w}}{\overline{m}}-\frac{3\mu \sigma J}{2(J+\sigma )}\right)=\frac{\varrho }{a^3}\left(\frac{\overline{w}}{\overline{m}}-\frac{3\mu \sigma (\sqrt{a^2+b^2}-a)}{2\sqrt{a^2+b^2}}\right)\to \frac{C}{a^6}.\end{array}$$

(81)

For nonlinear spinors, by (32) we have

$$\begin{array}{c}\hfill T_{\tau |\psi }^{\tau }={\rho }_{\psi }+W_{\psi },T_{r|\psi }^r=T_{\theta |\psi }^{\theta }=T_{\phi |\psi }^{\phi }=W_{\psi }-P_{\psi }.\end{array}$$

(82)

So its EOS in cosmology is given by [30,31]

$$\begin{array}{c}\hfill w_{\psi }\left(a\right)\equiv -\frac{T_{r|\psi }^r}{T_{\tau |\psi }^{\tau }}=\frac{P_{\psi }-W_{\psi }}{{\rho }_{\psi }+W_{\psi }}\to -1,(a\to +0).\end{array}$$

(83)

The above results suggest that dark nonlinear spinors may be the real dark energy and dark matter, which determines the large scale structure of the universe. $W_{\psi }$ acts like dark matter in energy density $T_{\tau |\psi }^{\tau }={\rho }_{\psi }+W_{\psi }$, but acts like dark energy in ‘pressure’ $\frac{1}{3}{T}_{k|\psi }^k=W_{\psi }-P_{\psi }$. We confirm this conclusion in the next section.

Now we consider the classical approximation of EMT (22) of electromagnetic field. We have two types of energy-momentum for $A^{\mu }$. The propagating photon can be directly treated as a particle with $P=\frac{1}{3}\rho$. The stationary electromagnetic field of ${\psi }_k$ satisfies

$$\begin{array}{c}\hfill {\partial }_{\mu }{\partial }^{\mu }{A}_k^{\nu }=e_k{\check{\alpha }}_k^{\nu },{\check{\alpha }}_k^{\nu }\equiv {\psi }_k^{+}{\alpha }^{\nu }{\psi }_k.\end{array}$$

(84)

By the principle of superposition, we have $A^{\mu }={\sum }_k{A}_k^{\mu }$. (84) in the comoving coordinate system with natural boundary condition $(A^{\mu }\to 0,r\to \infty )$ is static electromagnetic field, which can be solved by means of Green function. Then the general solution of (84) can be derived from static solution by local Lorentz transformation. We calculate the classical approximation of EMT (22) concretely.

The stationary field $A^{\mu }$ has only a very tiny geometrical structure like a Dirac-$\delta$, so we can omit the effect of curved space-time in comoving frame and calculate it in tangent Minkowski space-time. For simplicity, we take one particle as example to calculate EMT, and omit the particle index temporarily. By local Lorentz transformation we have [31,44]

$$\begin{array}{ccc}\hfill T^{ab}\left(\mathbf{x}\right)& =& \sqrt{1-v^2}{L}_c^a{L}_d^b{\overline{T}}^{cd}\left(\mathbf{X}\right),d{x}^a=L_b^{a}d{X}^b.\hfill \end{array}$$

(85)

where $\mathbf{X}={\gamma }_a{X}^a$ is the comoving coordinate of the particle, $L_b^a$ is Lorentz transformation, and ${\overline{T}}^{cd}\left(\mathbf{X}\right)$ is the static EMT of the spinor $\psi$.

In the comoving coordinate system, we have classical approximation

$$\begin{array}{c}\hfill {\overline{T}}^{ab}\left(\mathbf{X}\right)\to {\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^3}\left(-{\partial }^a{A}_c{\partial }^b{A}^c+\frac{1}{2}{\eta }^{ab}{\partial }_c{A}_d{\partial }^c{A}^d\right)d^{3}X.\end{array}$$

(86)

For static charge, we have ${\partial }_0{A}^a\equiv 0$, and then ${\overline{T}}^{0k}\equiv 0$,

$$\begin{array}{ccc}\hfill {\overline{T}}^{00}\left(\mathbf{X}\right)& \to & -\frac{1}{2}{\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^3}(\nabla A_d)·(\nabla A^d)d^{3}X\hfill \\ & =& \frac{1}{2}{\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^3}{A}_a\Delta A^a{d}^{3}X=-\frac{1}{2}{\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^3}e{A}_a{\check{\alpha }}^a{d}^{3}X\hfill \\ & =& -\frac{e^2}{8\pi }{\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^6}\frac{{\check{\alpha }}_a\left(\overrightarrow{\xi }\right){\check{\alpha }}^a\left(\overrightarrow{X}\right)}{|\overrightarrow{\xi }-\overrightarrow{X}|}{d}^3\xi d^{3}X=-m_e{\delta }^3\left(\overrightarrow{X}\right),\hfill \end{array}$$

(87)

in which $m_e$ is the proper electromagnetic mass of the particle defined by

$$\begin{array}{c}\hfill m_e=\frac{e^2}{8\pi }{\int }_{{\mathbb{R}}^6}\frac{{\check{\alpha }}_a\left(\overrightarrow{\xi }\right){\check{\alpha }}^a\left(\overrightarrow{X}\right)}{|\overrightarrow{\xi }-\overrightarrow{X}|}{d}^3\xi d^{3}X.\end{array}$$

(88)

Different from the nonlinear potential, $m_e$ is a constant independent of scale.

By symmetry of ${\check{\alpha }}^a\left(\overrightarrow{X}\right)$, we have ${\overline{T}}^{kl}=0,\left(k\ne l\right)$, so ${\overline{T}}^{ab}$ is diagonal.

$$\begin{array}{c}\hfill {\overline{T}}_a^a\left(\mathbf{X}\right)\to {\delta }^3\left(\overrightarrow{X}\right){\int }_{{\mathbb{R}}^3}\left(-{\partial }_a{A}_b{\partial }^a{A}^b+2{\partial }_c{A}_d{\partial }^c{A}^d\right)d^{3}X=-2{\overline{T}}_0^0=2m_e{\delta }^3\left(\overrightarrow{X}\right).\end{array}$$

(89)

Thus, in comoving coordinate system we get average classical approximation of EMT as

$$\begin{array}{c}\hfill {\overline{T}}_b^a\left(\mathbf{X}\right)\to m_e{\delta }^3\left(\overrightarrow{X}\right)\mathrm{diag}\left(-1,1,1,1\right),\end{array}$$

(90)

which provided negative energy and negative pressure. Converting (91) into general form we get

$$\begin{array}{c}\hfill T_A^{\mu \nu }\left(\mathbf{x}\right)\to \sum _k{m}_{ek}(-2u_k^{\mu }{u}_k^{\nu }+g^{\mu \nu })\sqrt{1-v_k^2}{\delta }^3(\overrightarrow{x}-{\overrightarrow{X}}_k).\end{array}$$

(91)

From ${\widehat{p}}_{\psi }^{\mu }$ defined in (11), we find ${\psi }_k$ and $A_k^{\mu }$ are not independent systems, and some EMT components of electromagnetic field are included in $T_{\psi }^{\mu \nu }$ given by (21). Therefore, the EMT (91) of stationary $A^{\mu }$ can be merged into $T_{\psi }^{\mu \nu }$, and we don’t need to calculate it separately. The situation of field ${\Phi }^{\mu }$ is the same of $A^{\mu }$, we don’t discuss it any more.

From the discussion of this section, we get some important knowledge about the total mass density ${\rho }_{\mathrm{m}}=T_0^0$ and the total pressure $P_{\mathrm{m}}=-\frac{1}{3}{T}_k^k$. Except for the global scalar field $\varphi$, the classical approximation of EMT of all other matter has a standard formalism (32), and the EOS is algebraic equation. By local structure in the universe such as galaxies, the uniform scalar $\varphi$ is obviously much weaker than other kind components of matter. While $a\to +0$, $P_{\mathrm{m}}$ is controlled by function $W\sim -\mathrm{O}\left(a^{-6}\right)$ [30,31], and we have $P_{\mathrm{m}}<0$ as $a\to +0$. In cosmology, although we call $P_{\mathrm{m}}$ pressure, but it is actually a variable including the potentials of all fields, so the negative pressure is reasonable [11,21,22].

3.5. About the Macroscopic Lagrangian of Fluid

By (8) and (60) we get relation for ideal gases in comoving coordinate system

$$\begin{array}{c}\hfill {\mathcal{L}}_p=3P_p-{\rho }_p,\end{array}$$

(92)

where ${\rho }_p$ is energy density acting as dynamical variable, EOS is $P_p=P\left({\rho }_p\right)$. Variation of (92) with respect to $\rho$ results in $P=\frac{1}{3}\rho$, which shows too many degrees of freedom of the system have been lost in the transformation from (8) into (92). The usual treatment requires some auxiliary conditions such as the first law of thermodynamics. However, this is unreasonable, because the Lagrangian is the most primitive constraint on a system, and can only determine the energy conservation law of the system by Lagrangian, but not the contrary.

Next we discuss this problem. For perfect fluid, its microscopic Lagrangian has the following form,

$$\begin{array}{c}\hfill {\mathcal{L}}_f=-\sum _n{m}_n\sqrt{1-v_n^2}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)-\mathcal{P},\end{array}$$

(93)

where $\mathcal{P}$ represents various scalar potentials related to specific interactions such as W in (32). These potentials may be complex but often easy to handle. Now we mainly explore how to reduce (93) to a reasonable macroscopic formalism.

Assuming the fluid has a mean velocity field $\overrightarrow{V}$ relative to the Gu’s natural frame (15), so we have a local Lorentz transformation in the co-moving frame ${\overline{x}}^{\mu }$ as

$$\begin{array}{c}\hfill \sqrt{1-v_n^2}=\sqrt{1-{\overline{v}}_n^2}\frac{d\overline{t}}{d\tau }=\sqrt{1-{\overline{v}}_n^2}\frac{\sqrt{1-{\overrightarrow{V}}^2}}{1-|\overrightarrow{V}|{\overline{v}}_{n}cos\theta },\end{array}$$

in which ${\overrightarrow{V}}^2={\overline{g}}_{kl}{V}^k{V}^l$ and $\overrightarrow{V}·{\overrightarrow{v}}_n=|\overrightarrow{V}|{\overline{v}}_{n}cos\theta$. For mass density (60), we have

$$\begin{array}{ccc}\hfill {\rho }_p& =& \sum _n\frac{m_n}{\sqrt{1-v_n^2}}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)=\frac{1}{\sqrt{1-{\overrightarrow{V}}^2}}\sum _n\frac{m_n(1-|\overrightarrow{V}|{\overline{v}}_{n}cos\theta )}{\sqrt{1-{\overline{v}}_n^2}}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\hfill \\ & =& \frac{1}{\sqrt{1-{\overrightarrow{V}}^2}}\sum _n{m}_n\left(\frac{1}{4\pi }{\int }_{\Omega }\frac{(1-|\overrightarrow{V}|{\overline{v}}_{n}cos\theta )sin\theta d\theta d\phi }{\sqrt{1-{\overline{v}}_n^2}}\right)\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\hfill \\ & =& \frac{1}{\sqrt{1-{\overrightarrow{V}}^2}}\sum _n\frac{m_n}{\sqrt{1-{\overline{v}}_n^2}}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)=\frac{\rho }{\sqrt{1-{\overrightarrow{V}}^2}},\hfill \end{array}$$

where $\rho$ is the co-moving energy density. Obviously, this formula is the right mass energy density of a star moving with velocity $\overrightarrow{V}$ relative to the natural coordinate system.

By (8), in average sense we have

$$\begin{array}{ccc}\hfill -{\mathcal{L}}_p& =& \sum _n{m}_n\sqrt{1-v_n^2}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)=\sqrt{1-{\overrightarrow{V}}^2}\sum _n{m}_n\frac{\sqrt{1-{\overline{v}}_n^2}}{1-|\overrightarrow{V}|{\overline{v}}_{n}cos\theta }\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\hfill \\ & =& \sqrt{1-{\overrightarrow{V}}^2}\sum _n{m}_n\left(\frac{1}{4\pi }{\int }_{\Omega }\frac{\sqrt{1-{\overline{v}}_n^2}sin\theta d\theta d\phi }{1-|\overrightarrow{V}|{\overline{v}}_{n}cos\theta }\right)\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\hfill \\ & =& \sqrt{1-{\overrightarrow{V}}^2}\sum _n{m}_n\frac{\sqrt{1-{\overline{v}}_n^2}}{|\overrightarrow{V}|{\overline{v}}_n}ln\sqrt{\frac{1+|\overrightarrow{V}|{\overline{v}}_n}{1-|\overrightarrow{V}|{\overline{v}}_n}}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n).\hfill \end{array}$$

(94)

We find $\overline{\rho }\sqrt{1-{\overrightarrow{V}}^2}<-{\mathcal{L}}_p<\rho \sqrt{1-{\overrightarrow{V}}^2}$. For the static ideal dust without other interaction, ${\mathcal{L}}_p$ should have the structure of (92) as $\overrightarrow{V}=0$. In order to get the right equation of motion, we assume

$$\begin{array}{c}\hfill {\mathcal{L}}_f=-\rho \sqrt{1-{\overrightarrow{V}}^2}+\tilde{P}-\mathcal{P},\end{array}$$

(95)

where $\tilde{P}$ is a function of P to be determined, and the particle index ‘${}_p$’ is omitted for simplicity. $\mathcal{P}$ still represents other potentials. From (95) we find $\tilde{P}$ is a scalar, but $\rho$ is not. As a matter of fact, energy and energy density indeed are not scalars. (95) seems to be a reasonable model for realistic fluid.

For Lagrangian (95) of fluid, we have action as

$$\begin{array}{ccc}\hfill S_f& =& \int \left(-\rho \sqrt{1-V^2}+\tilde{P}-\mathcal{P}\right)\sqrt{g}{d}^{4}x\hfill \\ & =& -\int \rho \sqrt{g_{\mu \nu }d{X}^{\mu }d{X}^{\nu }}\sqrt{\overline{g}}{d}^{3}x+\int (\tilde{P}-\mathcal{P})\sqrt{g}{d}^{4}x,\hfill \end{array}$$

(96)

where $d{X}^0=dt$. If $\mathcal{P}$ does not manifestly depend on $g_{\mu \nu }$, variation of (96) with respect to $g_{\mu \nu }$ we get

$$\begin{array}{ccc}\hfill \delta S_f& =& -\frac{1}{2}\int \frac{\rho d{X}^{\mu }d{X}^{\nu }}{\sqrt{g_{\mu \nu }d{X}^{\mu }d{X}^{\nu }}}\delta g_{\mu \nu }\sqrt{\overline{g}}{d}^{3}x-\frac{1}{2}\int \rho \sqrt{g_{\mu \nu }d{X}^{\mu }d{X}^{\nu }}{g}^{kl}\delta g_{kl}\sqrt{\overline{g}}{d}^{3}x\hfill \\ & & +\frac{1}{2}\int (\tilde{P}-\mathcal{P})g^{\mu \nu }\delta g_{\mu \nu }\sqrt{g}{d}^{4}x\hfill \\ & =& -\frac{1}{2}\int \rho \sqrt{1-V^2}{U}^{\mu }{U}^{\nu }\delta g_{\mu \nu }\sqrt{g}{d}^{4}x-\frac{1}{2}\int \rho \sqrt{1-V^2}{g}^{kl}\delta g_{kl}\sqrt{g}{d}^{4}x\hfill \\ & & +\frac{1}{2}\int (\tilde{P}-\mathcal{P})g^{\mu \nu }\delta g_{\mu \nu }\sqrt{g}{d}^{4}x,\hfill \end{array}$$

(97)

in which

$$\begin{array}{c}\hfill U^{\mu }=\frac{d{X}^{\mu }}{d\tau }\frac{1}{\sqrt{1-V^2}}=\frac{d{X}^{\mu }}{\sqrt{g_{00}}dt}\frac{1}{\sqrt{1-V^2}},U_{\alpha }{U}^{\alpha }=g_{\mu \nu }{U}^{\mu }{U}^{\nu }=1.\end{array}$$

By (97) we obtain the EMT of the fluid as

$$\begin{array}{c}\hfill T_f^{\mu \nu }=-2\frac{\delta \left({\mathcal{L}}_f\sqrt{g}\right)}{\sqrt{g}\delta g_{\mu \nu }}=\left\{\begin{array}{c}\rho \sqrt{1-V^2}{U}^0{U}^0+(\mathcal{P}-\tilde{P})g^{00},\hfill \\ \rho \sqrt{1-V^2}(U^k{U}^l+g^{kl})+(\mathcal{P}-\tilde{P})g^{kl}.\hfill \end{array}\right.\end{array}$$

(98)

The Lagrangian and EMT of fluid can be derived only in natural coordinate system (15). Comparing (98) with (37) we obtain

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{m}}& =& T_0^0=\frac{\rho }{\sqrt{1-V^2}}+\mathcal{P}-\tilde{P},\hfill \end{array}$$

(99)

$$\begin{array}{ccc}\hfill P_{\mathrm{m}}& =& -\frac{1}{3}{T}_k^k=\frac{\rho V^2}{3\sqrt{1-V^2}}-\rho \sqrt{1-V^2}-\mathcal{P}+\tilde{P}.\hfill \end{array}$$

(100)

For static ideal dust we have $\overrightarrow{V}=0,\mathcal{P}=0$. According to (100) we should set $\tilde{P}=P$, where P is the normal pressure defined by (60). Therefore we finally obtain the Lagrangian of fluid as

$$\begin{array}{c}\hfill {\mathcal{L}}_f=-\rho \sqrt{1-V^2}+P-\mathcal{P}.\end{array}$$

(101)

(100) includes another source of negative pressure $-\rho \sqrt{1-V^2}$, which is a spatial effect. The symmetry between time and space is broken by the factor $\sqrt{1-{\overrightarrow{V}}^2}$.

For the microscopic Lagrangian (8) of ideal dust, we have action

$$\begin{array}{ccc}\hfill S_p& =& -\int \sum _n{m}_n\sqrt{g_{\mu \nu }d{X}_n^{\mu }d{X}_n^{\nu }}\left(\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\sqrt{\overline{g}}\right)d^{3}x\hfill \\ & =& -\int \sum _n{m}_n\sqrt{g_{\mu \nu }d{X}_n^{\mu }d{X}_n^{\nu }}\widehat{\delta }(\overrightarrow{x}-{\overrightarrow{X}}_n)d^{3}x,\hfill \end{array}$$

(102)

where ${\int }_{R^3}f\left(x\right)\widehat{\delta }(\overrightarrow{x}-{\overrightarrow{X}}_n)d^{3}x=f\left({\overrightarrow{X}}_n\right)$, so $\widehat{\delta }\left(\overrightarrow{x}\right)=\delta \left(\overrightarrow{x}\right)\sqrt{\overline{g}}$ is a function independent of $g_{\mu \nu }$. Variation of (102) with respect to $g_{\mu \nu }$ we get

$$\begin{array}{ccc}\hfill \delta S_p& =& -\frac{1}{2}\int \sum _n\frac{m_{n}d{X}^{\mu }d{X}^{\nu }}{\sqrt{g_{\mu \nu }d{X}^{\mu }d{X}^{\nu }}}\delta g_{\mu \nu }\widehat{\delta }(\overrightarrow{x}-{\overrightarrow{X}}_n)d^{3}x\hfill \\ & =& -\frac{1}{2}\int \sum _n{m}_n{u}_n^{\mu }{u}_n^{\nu }\sqrt{1-v_n^2}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n)\delta g_{\mu \nu }\sqrt{g}{d}^{4}x.\hfill \end{array}$$

Thus we have EMT of the dust as

$$\begin{array}{c}\hfill T_p^{\mu \nu }=-2\frac{\delta {\mathcal{L}}_p}{\delta g_{\mu \nu }}=\sum _n{m}_n{u}_n^{\mu }{u}_n^{\nu }\sqrt{1-v_n^2}\delta (\overrightarrow{x}-{\overrightarrow{X}}_n).\end{array}$$

(103)

For (103) we find that the symmetry between time and space holds due to $\widehat{\delta }(\overrightarrow{x}-{\overrightarrow{X}}_n)$.

4. Dynamical Constraints on $\mathit{K}$ and $\Lambda$

Since the Friedmann equation is a dynamical equation, it is hard to determine its constants by static analysis. In [12,18,23,24], under respective assumptions the authors have provided some dynamical analysis for the behavior of the universe. Here we also qualitatively analyze the dynamical properties of the Friedmann equation according to the above results. Under the assumptions of positive mass density (43) and negative initial pressure (44), we find that the universe cannot reach the initial singularity, as well as the parameters $K=1$ and $\Lambda \simeq 0$. That is to say, the spatial structure of the universe is closed 3 dimensional sphere $S^3$, and the cosmological constant is likely to be zero. Besides, the universe should be cyclic in time. Obviously, these conclusions will help to correct some popular misconceptions and bring great convenience to further study the properties of other problems in cosmology such as dark matter and dark energy.

From the previous analysis, we find that the scale factor $a\left(\tau \right)\sim \mathrm{O}\left({\tau }^j\right)$ is nonanalytic at origin, which increases difficulty for analysis, so we adopt the conformal FLRW metric,

$$d{s}^2=a{\left(t\right)}^2\left(d{t}^2-d{r}^2-\mathcal{S}{\left(r\right)}^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2)\right).$$

(104)

Then the Friedmann Equation (38) becomes

$$a^{\prime 2}=-K{a}^2+\frac{1}{3}\Lambda a^4+\frac{1}{3}\kappa {\rho }_{\mathrm{m}}{a}^4.$$

(105)

where $a^{\prime }=\frac{d}{dt}a\left(t\right)$. (105) can be rewritten as

$$a^{\prime 2}=F\left(a\right),F\left(a\right)\equiv 2\overline{R}a-K{a}^2+\frac{1}{3}\Lambda a^4+X\left(a\right),$$

(106)

where $\overline{R}$ corresponds to the total conformal density of proper mass,

$$\overline{R}\equiv \underset{a\to \infty }{lim}\frac{1}{6}\kappa {\rho }_{\mathrm{m}}{a}^3\approx \frac{1}{6}\kappa \varrho$$

(107)

which is a constant. (107) means that the total conformal energy of the universe is bounded. By (106) we find $\overline{R}$ has length dimension, it is the average scale of the universe [21].

Comparing (106) with (105), we get relation

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{m}}& =& \frac{3}{\kappa a^4}\left(F\left(a\right)+K{a}^2-\frac{1}{3}\Lambda a^4\right)\hfill \end{array}$$

(108)

$$\begin{array}{ccc}& =& \frac{3}{\kappa a^4}\left[2\overline{R}a+X\left(a\right)\right].\hfill \end{array}$$

(109)

$X\left(a\right)$ corresponds to the rest and unknown parts of ${\rho }_{\mathrm{m}}$, which satisfies $X\left(a\right)\sim -\mathrm{O}\left(a^{-2}\right)$ as $a\to +0$ and $X\left(a\right)/a\to 0$ as $a\to \infty$. In physics, $X\left(a\right)$ is mainly determined by nonlinear potential ${\sum }_k{N}_k$. The specific form of $X\left(a\right)$ is not important for qualitative analysis, only its asymptotic properties as $a\to +0$ have influence on the following discussion. The property of the solution of (106) can be clearly discussed by means of phase trajectories.

Substituting (109) into energy conservation law (39), we get the total pressure as

$$\begin{array}{c}\hfill P_{\mathrm{m}}=-\frac{1}{\kappa a^4}\left[X^{\prime }\left(a\right)a-X\left(a\right)\right].\end{array}$$

(110)

We find $P_{\mathrm{m}}$ is irrelative to $\overline{R}$. Since the derivatives of pressure and potential correspond to ordinary forces which should be finite, so $P_{\mathrm{m}}$ should be at least continuous. Then by (110) we have at least $X\left(a\right)\in C^1$. Consequently, by the definition of $F\left(a\right)$ in (106), we also have $F\left(a\right)\in C^1$.

The following discussion is based on Friedamnn Equation (106) as well as two assumptions (43) and (44), namely, the positive total energy density ${\rho }_{\mathrm{m}}>0$ and negative initial pressure $P_{\mathrm{m}}\left(a\right)<0,(a\to +0)$. Clearly the two assumptions are compatible with observational facts [24].

Theorem 4.

For function $F\left(a\right)\in C^1$ defined in (106), if ${lim}_{a\to 0}X\left(a\right)\ne 0$ and condition (44) holds, we have

$$F\left(a\right)<0,(a\to +0).$$

(111)

Proof. 

In the case of $0<\left|X\right(0\left)\right|<\infty$, by (110) and condition (44), we have

$$\begin{array}{c}\hfill P_{\mathrm{m}}\to \frac{X\left(0\right)}{\kappa a^4}<0,(a\to +0).\end{array}$$

(112)

Consequently, by the definition of $F\left(a\right)$ in (106), we get

$$F\left(0\right)=X\left(0\right)<0.$$

(113)

In the case of $X\to X_0{a}^{-n},(a\to +0,n>0)$, again by (110) and condition (44), we have

$$\begin{array}{c}\hfill P_{\mathrm{m}}\to \frac{(n+1)X_0}{\kappa a^{4+n}},(a\to +0).\end{array}$$

(114)

by $P_{\mathrm{m}}<0$ we find $X_0<0$. According to the definition of $F\left(a\right)$ in (106), we get

$$F\to \frac{X_0}{a^n}<0,(a\to +0).$$

(115)

Then we prove (111) holds in all cases. □

The above theorem implies the following important conclusion:

$a>0$, the evolution of the universe cannot reach the initial singularity.

Now we check $a>0$. For the solution of Friedmann Equation (106), we have $F\left(a\right)=a^{\prime 2}\ge 0$. By (111) and the continuity of $F\left(a\right)$, the equation $F\left(a\right)=0$ certainly have a positive root $0<a_0\ll \overline{R}$.

If $F\left(a\right)=0$ has only this positive root $a_0$, by $F\left(a\right)\in C^1$, then $F\left(a\right)$ can be expressed as

$$F\left(a\right)=(a-a_0)A\left(a\right),(A>0,\forall a\ge a_0).$$

(116)

If $F\left(a\right)=0$ has a series of different positive roots $0<a_0<a_1<a_2<\cdots$, $F\left(a\right)$ for the practical universe should be simply connected and then can be expressed as

$$F\left(a\right)=(a-a_0)(a_1-a)B\left(a\right),(B>0,a_0\le a\le a_1).$$

(117)

Since Friedmann equation is an equation in average sense, the multiple roots are meaningless in physics. The connected phase trajectories $a\sim a^{\prime }$ of dynamical Equation (106) with (116) or (117) are displayed in Figure 1, in which we have set the average scale $\overline{R}=1$. (117) corresponds to the cyclic cosmological model, and (116) to the bouncing one. We set the time origin $t=0$ at the turning point $a\left(0\right)=a_0$. Form Figure 1 we learn clearly the initial singularity is absent, i.e., a cannot reaches 0 point.

Substituting (116) or (117) into (108) and letting $a=a_0$, in both cases we have

$$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}\left(a_0\right)=\frac{3}{\kappa a_0^2}\left(K-\frac{1}{3}\Lambda a_0^2\right)>0.\end{array}$$

(118)

Since $\Lambda \ge 0$ in cosmology, by (117) we certainly have $K=1$ due to ${\rho }_{\mathrm{m}}>0$. Then we get another conclusion:

$K=1$, the space of the universe is a closed 3 dimensional sphere$S^3$.

In the cyclic closed case (117), we have an estimation of upper bound for the cosmological constant $\Lambda$. Substituting (117) into (108) and letting $a=a_1$, by (43) we have

$$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}\left(a_1\right)=\frac{3}{\kappa a_1^2}\left(1-\frac{1}{3}\Lambda a_1^2\right)>0.\end{array}$$

(119)

Since $\overline{R}a$ forms the main part of mass-energy density at present time, which can be estimated by observational data [21], and $\left|X\left(a\right)\right|\ll \overline{R}a$ as $a\to \infty$ can be omitted. For $\Lambda \ge 0$ we have $a_1\ge 2\overline{R}$ and estimation

$$\begin{array}{c}\hfill 0\le \Lambda <\frac{3}{a_1^2}\le \frac{3}{4{\overline{R}}^2}\sim {10}^{-24}{\mathrm{ly}}^{-2}.\end{array}$$

(120)

So for the cosmological constant $\Lambda$ in a cyclic and closed universe, we get the third conclusion:

$0\le \Lambda <{10}^{-24}{\mathrm{ly}}^{-2}$, the cosmological constant is an infinitesimal.

This estimation is less than the present observational data. This difference can be explained by the potentials W in energy-momentum tensor, which is a fast decaying $\Lambda \left(a\right)$ in Friedmann equation. Therefore setting constant $\Lambda =0$ is a good choice in cosmology.

For a bouncing cosmological model, while $a\to \infty$, the behavior of Friedmann Equation (106) is controlled by dominant term $\frac{1}{3}\Lambda a^4$, and the fast decaying term $X\left(a\right)/a\to 0$ can be omitted. In this case, to be clear, (106) can be replaced by the following dimensionless Hamiltonian-like equation

$$\begin{array}{c}\hfill \frac{1}{2}{q}^{\prime 2}+V\left(q\right)=H,V\left(q\right)\equiv \frac{1}{6}\lambda q^4-\frac{1}{54}q(q+6){(q-3)}^2,\end{array}$$

(121)

in which

$$\begin{array}{c}\hfill q\left(t\right)=\frac{a\left(t\right)}{\overline{R}},\lambda =\frac{1}{9}-\Lambda {\overline{R}}^2.\end{array}$$

(122)

For solution of Friedmann equation we have $H\equiv 0$, which is conserved.

In (121), $q\left(t\right)$ is equivalent to the coordinate of a unit mass, $V\left(q\right)$ potential, and H energy. The potential function $V\left(q\right)$ and phase trajectories $(q^{\prime },q)$ is displayed by Figure 2 and Figure 3. The solution is closed if $\lambda >0\iff \Lambda <\frac{1}{9}{\overline{R}}^{-2}\sim {10}^{-25}{\mathrm{ly}}^{-2}$, but bouncing if $\lambda \le 0\iff \Lambda \ge \frac{1}{9}{\overline{R}}^{-2}$. As we can see from Figure 2 and Figure 3, $\lambda$ or $\Lambda$ has only influence on the behavior of a fully developed universe but has little influence on a small universe. On the contrary, the function of $X\left(a\right)$ can prevent $a\left(t\right)$ from reaching the initial singularity but has little influence on a fully developed universe. For closed universe, we have the second root $2\overline{R}\le a_1<3\overline{R}$. In contrast $\Lambda <\frac{1}{9}{\overline{R}}^{-2}$ with (120), we find the estimation of $\Lambda$ by criterion ${\rho }_{\mathrm{m}}>0$ is larger. In bouncing case [25], we have ${\left(\frac{d}{d\tau }a\right)}^2\to C^2{a}^2\iff a\to C_{0}cosh\left(C\tau \right)$ as $\tau \to \pm \infty$. In this case, the solution of Friedmann equation globally exists and the asymptotic behavior is controlled by term $\frac{1}{3}\Lambda a^4$. For the big rip model, we need an even higher order term in Friedmann equation, which is unreasonable in physics.

The bouncing model with closed space is inconsistent with the isotropy and homogeneity of the present universe, because the universe should be heavily anisotropy and inhomogeneity before the turning point $t<0$ due to the lack of initial causality among remote parts, and some information should be kept today. However, this problem is absent in a cyclic and closed universe, because the anisotropy and inhomogeneity in a period will be polished by the subsequent evolution due to the entropy increasing principle. In a closed system, the dissipative effect will work when the inhomogeneity of the matter distribution is greater than the controlling degree of internal interaction of matter, it reduces the inhomogeneity to match the scale of the interactive force. From the above analysis we find $\Lambda a^4$ is purely a trouble term without any practical purpose.

5. Towards a Realistic Cosmological Model

In this section, we discuss the cosmological model with nonlinear potentials; that is, we examine the influence of the term $X\left(a\right)\sim \mathrm{O}\left(a^{-2}\right)$ on the Friedmann Equations (105)–(109). Based on the analysis of Section 2 and Section 3, we find that in (106), X has an influence on the behavior of the universe at the early age, and $\Lambda$ only has an influence on the very old universe. $\Lambda$ has not any clear purpose in physics, so we set $\Lambda =0$. In the conformal coordinate system (104), using (106) and $K=1$, the dynamics for a cyclic and closed universe is generally given by (105) and (117) as follows

$$\begin{array}{ccc}\hfill a^{\prime 2}& =& -a^2+\frac{\kappa }{3}{\rho }_{\mathrm{m}}{a}^4=-a^2+2\overline{R}a+c_0+\frac{c_1}{a}+\frac{c_2}{a^2}\hfill \\ & =& (a-\alpha )(\omega -a)\left(\frac{{(a-\delta )}^2+{\epsilon }^2}{a^2}\right),\hfill \end{array}$$

(123)

where $\alpha >0$ is the minimum scale and $\omega$ is the maximum scale, $\delta ,\epsilon$ are two parameters that reflect the influence of nonlinear potentials on the universe. and the parameters $\left|\delta \right|\sim \epsilon \sim \alpha \ll \omega$ all have length dimension. In fact, (117) only required ${(a-\delta )}^2\pm {\epsilon }^2>0$ in the domain $\alpha \le a\le \omega$. However, $\pm {\epsilon }^2$ is a high order small term and has no influence on the following analysis and calculation, and taking the term $-{\epsilon }^2$ adds a lot of cases for judgment, so only the case ${(a-\delta )}^2+{\epsilon }^2>0$ is needed for the following discussion.

For the convenience of discussion, we define a dimensionless function

$$\begin{array}{c}\hfill Q\equiv \frac{1}{a^2}\sqrt{(a-\alpha )(\omega -a)\left[{(a-\delta )}^2+{\epsilon }^2\right]}.\end{array}$$

For $\alpha \ll a\ll \omega$, we have the approximation $0<\alpha /a\ll 1$ and $0<a/\omega \ll 1$ in the calculation. The age of the universe is calculated by

$$\begin{array}{c}\hfill \mathcal{T}={\int }_{\alpha }^A\frac{a}{a^{\prime }}da={\int }_{\alpha }^A\frac{da}{Q\left(a\right)},\end{array}$$

where $\alpha \ll A\ll \omega$ is the present scale of the universe. The solution is the elliptic integrals. Denote the half period of the universe by

$$\begin{array}{c}\hfill {\mathcal{T}}_m={\int }_{\alpha }^{\omega }\frac{da}{Q\left(a\right)}.\end{array}$$

By the approximation $\alpha \ll A\ll \omega$, we have

$$\begin{array}{ccc}\hfill \mathcal{T}& =& {\mathcal{T}}_m-{\int }_A^{\omega }\frac{a^{2}da}{\sqrt{(a-\alpha )(\omega -a)\left[{(a-\delta )}^2+{\epsilon }^2\right]}},\hfill \\ & \dot{=}& [{\mathcal{T}}_m-\frac{\pi }{2}(\alpha +2\delta +\omega )]+(\alpha +2\delta )B+\frac{1}{3}(2\omega -\alpha -2\delta )B^3,\hfill \end{array}$$

where B is defined by

$$\begin{array}{c}\hfill B=\sqrt{\frac{A}{\omega -A}}\ll 1,\omega =(1+B^{-2})A.\end{array}$$

(124)

So we have the approximation

$$\begin{array}{ccc}\hfill {\mathcal{T}}_m& \dot{=}& \frac{\pi }{2}(\alpha +2\delta +\omega ),\hfill \end{array}$$

(125)

$$\begin{array}{ccc}\hfill \mathcal{T}& \dot{=}& (\alpha +2\delta )B+\frac{1}{3}(2\omega -\alpha -2\delta )B^3.\hfill \end{array}$$

(126)

In the above calculation, we used $a\ge A$, $x=\alpha /a\ll 1$ and

$$\begin{array}{ccc}\hfill Q& =& \sqrt{{\alpha }^3}{\left(\sqrt{(\omega -a)x(1-x)\left[{\epsilon }^2{x}^2+{(\alpha -\delta x)}^2\right]}\right)}^{-1}\hfill \\ & \to & \frac{1}{8\sqrt{{\alpha }^{3}x(\omega -a)}}\left(8{\alpha }^2+4\alpha (\alpha +2\delta )x+\mathrm{O}\left(x^2\right)\right)\hfill \\ & =& \frac{1}{2a\sqrt{a(\omega -a)}}\left((\alpha +2\delta )a+2a^2+\mathrm{O}\left(x^2\right)\right).\hfill \end{array}$$

then the integral can be calculated.

The Hubble parameter can be calculated as

$$\begin{array}{ccc}\hfill H& =& \frac{\dot{a}}{a}=\frac{a^{\prime }\left(t\right)}{a^2}=\frac{Q\left(a\right)}{a},\hfill \end{array}$$

(127)

$$\begin{array}{ccc}\hfill H_0& \dot{=}& \frac{1}{AB}\left(1-\frac{1}{2A}(\alpha +2\delta )-\frac{1}{8A^2}({\alpha }^2-4\alpha \delta -4{\epsilon }^2)\right).\hfill \end{array}$$

(128)

Denoting the critical density

$$\begin{array}{c}\hfill {\rho }_c=\frac{3}{8\pi G}{H}_0^2\sim 8\times {10}^{-30}\left({\mathrm{g}/\mathrm{cm}}^3\right),\end{array}$$

we have the dimensionless mass density

$$\begin{array}{ccc}\hfill \Omega & =& \frac{{\rho }_{\mathrm{m}}}{{\rho }_c}=1+{\left(\frac{a}{a^{\prime }}\right)}^2=1+Q^{-2}\hfill \\ & =& 1+B^2\left(1+\frac{1}{A}(\alpha +2\delta )+\frac{1}{A^2}({\alpha }^2+2\alpha \delta +3{\delta }^2-{\epsilon }^2)\right).\hfill \end{array}$$

(129)

In the second-order approximation, we find that $\alpha +2\delta$ acts as one parameter and the influence of $\epsilon$ disappears. This overcomes the shortage of useful observational data. Omitting the higher order terms in (128) and (129), and making the transformation

$$\begin{array}{c}\hfill C=2-\frac{\alpha +2\delta }{A},(0<C<2),\alpha +2\delta =(2-C)A,\end{array}$$

(130)

we obtain 3 equations for 3 parameters $A,B,C$.

$$\begin{array}{ccc}\hfill \mathcal{T}=\frac{1}{3}AB\left((8-3C)+C{B}^2\right),& & H_0=\frac{C}{2AB},\hfill \end{array}$$

(131)

$$\begin{array}{ccc}\hfill \Omega =1+(3-C)B^2.& & \hfill \end{array}$$

(132)

Theoretically, we can solve for $A,B,C$ from the above equations using observational data $\mathcal{T},H_0,\Omega$. By (131), we obtain

$$\begin{array}{c}\hfill (3-B^2)C^2-8C+6\mathcal{T}{H}_0=0,\end{array}$$

(133)

where $0<B\ll 1$. Substituting (133) into (132), we obtain the equation for C,

$$\begin{array}{c}\hfill \mathcal{T}{H}_0=\frac{C}{3-C}(C\delta \Omega +3C^2-17C+24),\delta \Omega =\Omega -1.\end{array}$$

(134)

The curves of $(C,\mathcal{T}{H}_0)$ are displayed in Figure 4. By the definition (130), we find that C is an increasing function of A or $\mathcal{T}$, so only the increasing part of the curves are meaningful in physics.

$\mathcal{T}{H}_0$ has a maximum. Omitting $\mathrm{O}\left(\delta {\Omega }^3\right)$, we have

$$\begin{array}{c}\hfill \mathcal{T}{H}_m=\frac{8}{9}+\frac{8}{45}\delta \Omega +\frac{392}{5625}\delta {\Omega }^2,C_m=\frac{4}{3}+\frac{28}{75}\delta \Omega +\frac{756}{3125}\delta {\Omega }^2.\end{array}$$

For simplicity, we take B as the independent parameter instead of $\delta \Omega$. From (133), we obtain the solution

$$\begin{array}{c}\hfill C=\frac{1}{3-B^2}\left(4-\sqrt{16-6\mathcal{T}{H}_0(3-B^2)}\right).\end{array}$$

(135)

Substituting it into (134) and noting $B\ll 1$, we obtain

$$\begin{array}{c}\hfill B^2=\frac{3\delta \Omega }{5+\sqrt{2(8-9\mathcal{T}{H}_0)}}.\end{array}$$

(136)

Substituting it into (135), we obtain C determined by observational data, then by (131), we obtain A. Other parameters can be calculated by

$$\begin{array}{c}\hfill \omega =A(1+B^{-2}),\alpha +2\delta =A(2-C),{\mathcal{T}}_m=\frac{\pi }{2}A(3-C+B^{-2}).\end{array}$$

(137)

From the observational data set [3,4,5,6], we have the following typical range of parameter values [4],

$$\begin{array}{c}\hfill \mathcal{T}=10\sim 14\mathrm{Gyr},h=0.50\sim 0.75,\delta \Omega ={10}^{-4}\sim {10}^{-2},\end{array}$$

where h is the dimensionless Hubble parameter

$$\begin{array}{c}\hfill H_0=100h\mathrm{km}·{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}=0.1023h{\mathrm{Gyr}}^{-1}.\end{array}$$

Based on the equations, we calculate the typical values for all parameters, some results are displayed in

Table 1. Correlation Among Cosmic Parameters.

$\Omega -1$	Equation (129)	0.00010	0.0010	0.0025	0.0050	0.0050	0.0075	0.0003
$\mathcal{T}$, Gyr	Cosmic Age	12.0000	12.000	13.000	13.000	11.000	11.000	12.000
h	Hubble’s	0.60000	0.6000	0.6500	0.6500	0.7000	0.7500	0.7000
$H_0$, 1/Gyr	Hubble’s	0.06138	0.0614	0.0665	0.0665	0.0716	0.0767	0.0716
$\mathcal{T}·H_0$	Equation (131)	0.73656	0.7366	0.8644	0.8644	0.7877	0.8440	0.8593
B	Equation (136)	0.00671	0.0212	0.0364	0.0515	0.0486	0.0618	0.0114
C	Equation (135)	0.78137	0.7813	1.1110	1.1097	0.8828	1.0314	1.0900
A, Gly	Present Scale	948.07	299.78	229.56	162.14	126.82	108.83	665.23
$\alpha +2\delta$, Gly	Equation (137)	1155.35	365.34	204.09	144.35	141.68	105.42	605.33
$\omega$, Gly	Largest Scale	21,035,047	665,393	173,573	61,381	53,810	28,643	5,082,600
${\mathcal{T}}_m$, Gyr	Half Period	33,043,667	1,045,773	272,969	96,644	84,747	45,157	7,984,698

.

By calculation, we have

$$\begin{array}{c}\hfill a\ddot{a}\dot{=}\frac{3-2C}{2B^2}>0,\ddot{a}=\frac{d^{2}a}{d{\tau }^2}.\end{array}$$

This means that the expansion of the universe is still accelerating. Using (123), we obtain the relation ${\rho }_{\mathrm{m}}\left(a\right)$ as

$$\begin{array}{ccc}\hfill {\rho }_{\mathrm{m}}& =& \frac{3}{\kappa a^4}\left(a^2+(a-\alpha )(\omega -a)\frac{{(a-\delta )}^2+{\epsilon }^2}{a^2}\right)\hfill \end{array}$$

(138)

$$\begin{array}{ccc}& \to & \frac{3\omega }{\kappa a^3}\left(1+(\alpha +2\delta )(\frac{1}{\omega }-\frac{1}{a})\right).\hfill \end{array}$$

(139)

Substituting (138) into the definition of X (109), we have the relation $X\left(a\right)$ as

$$\begin{array}{ccc}\hfill X& =& \omega \left(-(\alpha +2\delta )+(2\alpha \delta +{\delta }^2+{\epsilon }^2)\frac{1}{a}-\alpha ({\delta }^2+{\epsilon }^2)\frac{1}{a^2}\right)\hfill \\ & & -(2\alpha \delta +{\delta }^2+{\epsilon }^2)\frac{1}{a}+\alpha ({\delta }^2+{\epsilon }^2)\frac{1}{a^2}.\hfill \end{array}$$

(140)

Substituting it into (110), we obtain the relation $P_{\mathrm{m}}\left(a\right)$ as

$$\begin{array}{c}\hfill P_{\mathrm{m}}=\frac{\omega }{\kappa a^4}\left(-(\alpha +2\delta )+(2\alpha \delta +{\delta }^2+{\epsilon }^2)(\frac{2}{a}-\frac{1}{\omega })\right)+\mathrm{O}\left(\frac{1}{a^6}\right).\end{array}$$

(141)

From (138) and (141), we obtain the equation of state $w\left(a\right)$ as

$$\begin{array}{ccc}\hfill w& =& \frac{P_{\mathrm{m}}}{{\rho }_{\mathrm{m}}}\to -(\alpha +2\delta )\frac{1}{3a}+\mathrm{O}\left(\frac{1}{a^2}\right)\hfill \\ & \dot{=}& -\frac{(2-C)A}{3a}=-\frac{1}{3}(2-C)(1+z)<0,\hfill \end{array}$$

(142)

where z is the cosmological redshift.

From Table 1, we find that the structural parameters are sensitive to mass density $\delta \Omega =\Omega -1$: the larger the mass, the smaller the universe, and the rate is very large. The difficulty to detect the little value of $\delta \Omega$ is a good news for humanity, as it reflects that the universe is still very young, and we can expect a better future for humanity. Relatively, the parameters are not very sensitive to the age $\mathcal{T}$ and the Hubble parameter $H_0$. Since measuring $\delta \Omega$ is very difficult and ambiguous, we need a more efficient method to measure cosmic data. Obviously, the above model is hardly contradictory against observational data, this is because the above calculating results are all compatible with observational facts and the effective observational data are still not enough to determine all the coefficients in (123).

The cyclic and closed model can provide more natural explanations for cosmological problems. From (123), we learn that there is a cradle of rebirth $t=0$ rather than an initial singularity. At this time, the drifting speed of the particles or the cosmic temperature takes the maximum, and the volume of the Universe takes the minimum. This situation will result in violent collisions among galaxies and the particles. Such collisions will recover the vitality of the dead stars and combined particles. Besides, the nonlinear potential will become dominant at this time, and the property of the particles will change. Just like the rebirth of the Phoenix, the Universe gets her rebirth at the turning point.

6. Discussion and Conclusions

Since the above derivations are all elementary and reliable, and all concepts have clearly physical meanings, so the conclusions should be quite credible. However, the above conclusions contradict the singularity theorems, this is because some preconditions of these theorems are invalid in physical world [25]. For example, the existence of negative pressure or potential is ignored in the energy condition. Besides, we have only unique realistic simultaneous Cauchy surface in the space-time [30,32], but the derivation of Raychaudhury equation unconsciously assumes and uses the future properties of the space-time. So this equation cannot be generally used for dynamical analysis. A realistic singularity in Nature is actually contradictory and incomprehensible, so the absence of singularity in Nature is a basic principle in physics.

To sum up, by studying the properties of the EMT and EOS of various physical fields, especially the EMT and EOS of nonlinear spinors, and qualitatively analyzing the dynamical behavior of the general Friedmann equation and logical relations between parameters, we get some definite constraints on $(K,\Lambda )$. We find that only the cyclic and closed cosmological model with a tiny or vanishing $\Lambda$ is natural and reasonable in physics. The other cases include nonphysical effects or logical contradictions. Such constraints will be helpful for the research of some other issues in cosmology. In some sense, we restored Heraclitus’ ancient faith: “The world, an entity out of everything, was created by neither gods nor men, but was, is and will be eternally living fire, regularly becoming ignited and regularly becoming extinguished”.