Cold Collapse and Bounce of a FLRW Cloud

1. Introduction

In 1970, Hawking and Penrose proved that singularities cannot be avoided in GR [1]. Due to the success of these singularity theorems, the solutions and ideas of non-singular or bouncing universes have been overlooked in the last 30 years of the 20th century [2]. A bounce of collapsing matter at supranuclear densities naturally avoids the formation of a singularity. To allow for such a bounce in this class of solutions, Einstein gravity would have to violate the strong energy conditions (SEC). For a perfect fluid with energy density $\rho$ and pressure P, the SEC corresponds to $\rho +3P\ge 0$. This SEC violation only started to be accepted in the arising models of inflation [3,4]. The discovery of positive acceleration of the universe [5,6,7,8,9,10,11] proved the existence of SEC violation, and the models of bouncing cosmologies received new interest [2]. These models aim to solve the problems encountered by the standard model of cosmology, among which are: i) the singularity problem: the existence of a singularity means that the spacetime is geodesically past incomplete; ii) the flatness problem: the critical value of matter-energy density of the universe needs to be extremely fine-tuned, any small deviation from this value will lead to a non-flat universe; and iii) the horizon problem: we see a uniform CMB sky where superhorizon scales appear (these are scales larger than the Hubble horizon at the time when the CMB photons were emmitted). How is this possible if the structures were not in casual contact in the past? Alternatively, new theories of modified gravity and the concept of quantum gravity [12,13] were developed. However, the latter is incomplete and cannot be easily testable (see [14] for a review).

One realization of the bouncing universe is the Black Hole (BH) Universe or BHU model presented in [15,16]. In this model, some of the problems faced by the standard model are avoided without the need to introduce new theories of modified gravity or Dark Energy (DE). In the BHU model, our observable universe is in the interior of its corresponding BH event horizon, and the $\mathsf{\Lambda }$-term in the standard cosmology arises as a boundary to the evolution of the system inside the BH event horizon [17,18]. The initial condition of such a universe could be a tenuous, cold, but large Friedmann-Lemaitre-Robertson-Walker (FLRW) cloud containing the total mass of our observable universe. This cloud collapses under its own gravity into a BH. For such a large mass, the density at the time of the formation of the BH is very small, and the pressureless collapse inside the BH event horizon proceeds as if there was no BH formed. Note that we use the expression cold to describe a fluid at zero temperature, resulting in zero thermal pressure. We use both expressions equally. Therefore, if not stated otherwise, pressureless refers to zero thermal pressure.

The gravitational collapse occurs homologously on a shorter time scale than thermal interaction of particles could create any significant thermal pressure. Thus, it is reasonable to assume that the cloud fluid stays cold during the entire collapse process. This means that if nothing else occurs, the density of the cloud continues to increase during the collapse until a singularity is formed.

However, there are effects that can slow down, if not even bring the cold collapse to a halt. Whenever one of the fermionic constituents of the cloud reaches its quantum ground state, the Pauli Exclusion Principle creates a degeneracy pressure independent of temperature. Assuming a cloud of neutral hydrogen, first the electrons would become degenerate and create a pressure that could stabilize configurations up to the Chandrasekhar mass limit ($\sim 1.4$M_⊙). Typical maximum densities allowed in this way are lower than $\rho \lesssim {10}^7$g cm⁻³ At higher densities, reached for larger initial masses, the electrons recombine with protons, and the collapse continues. Similarly, at nuclear saturation densities, the collapse could be stopped due to the quantum pressure of the ground state of nuclear matter. The maximum mass allowed for this to happen sensitively depends on the equation of state (EoS) of nuclear matter [19,20], but it should be between $2-3$ M_⊙.

Similar scenarios are actually observed in nature at the end of the life of massive stars. The fate of these objects depends on the final mass of their collapsing cores. Stars can lose a significant amount of their initial mass during their different burning stages through stellar winds or at the end of their lifetime by ejecting the outer layers during the collapse of the central cores. For stars with zero-age main-sequence masses of up to about $\sim 8$M_⊙, the final cores reach at most the Chandrasekhar mass ($\sim 1.4$M_⊙). Therefore, these stars end their lives as White Dwarfs. The cores of more massive stars ($\sim 8$M_⊙$\lesssim M\lesssim 40$M_⊙) eventually become unstable and collapse once the Chandrasekhar mass is reached. This collapse is stopped once the collapsing core of the progenitor of $\sim 1.4$M_⊙ reaches nuclear saturation density, and matter cannot be compressed further. The still-infalling outer layers bounce back from the proto-neutron star (PNS) surface and create an outward-moving shock wave. However, this shock wave stalls because of energy losses by neutrinos and dissociation of the infalling nuclei. The details of how this shock revives and the bounce finally leads to an exploding supernova are a fascinating topic of current astrophysical research (see [21,22,23,24,25,26] for reviews). What is of interest to the current work is the existence of the bounce of the infalling matter. Note that the densities in neutron stars and in the nucleus of an atom are comparable, despite the former having a mass ${10}^{57}$ times larger than the latter. Having this in mind, the step from stellar size objects we study here with tens to thousand solar masses to the mass of the universe ($\sim {10}^{22}$ M_⊙, Gaztañaga [27]) appears less exaggerated.

However, we have to note that there are essential differences between the cold collapse we are interested in and a supernova explosion. Supernova are `hot’ events. This not only means that there are photons present that mediate thermal pressure but also that there are plenty of nuclear reactions occurring simultaneously. Depending on the conditions of temperature, pressure, and density, different reactions create and destroy nuclei, releasing copious amounts of neutrinos that may interact with matter in different locations. The outward-moving shock increases entropy and, hence, the temperature of the infalling matter, which leads to the dissociation of the nuclei into nucleons. The composition of matter depends sensitively on all the conditions and on the dynamics of the entire process. Fortunately, most of these complications disappear once we make the reasonable assumption of a cold collapse of neutral hydrogen.

The scenarios discussed above describe what happens for density distributions whose total mass is sufficiently low, and the bounce occurs before the respective distribution is confined within its gravitational radius. The larger the total mass, the lower the density at which a BH forms, and consequently, for more massive density distributions, (supra)nuclear densities only appear inside the corresponding horizon. This brings us back to our BHU model, where we investigate the collapse of a very massive cloud, for which the cold collapse occurs within the black hole itself. Given the enormous mass, one might assume that the collapse would not stop at nuclear saturation density, potentially leading to even higher densities.

In the BHU model (Gaztañaga [27]), the working hypothesis is that some kind of quantum degeneracy pressure or a more energetic ground state of matter could be reached at some supranuclear density. For masses as large as our universe, the threshold density could be much larger than the nuclear saturation density. In analogy to a supernova, the matter then may not be compressed beyond this threshold, the collapse halts abruptly, and the infalling matter bounces back. Consequently, the propagation of matter changes from collapse to expansion, potentially providing an alternative explanation for cosmic inflation and the related dark energy. A mathematical singularity, which is inevitable in non-bouncing GR models, can be avoided with such a bounce. After the Big Bang, the superhorizon perturbations from the collapsing phase reenter the horizon as the expansion proceeds and seed structure formation (Gaztañaga [28]), solving the horizon problem. As we see in the following sections, the curvature value in the GR framework (k) goes to a finite value during the bounce, but the spacetime flattens out as the expansion proceeds and density reduces. Therefore, the flatness problem is solved in a manner similar to that of inflation.

As a first step towards a model for our universe, we first need to understand the mechanism of the cold collapse and bouncing. We begin our study by considering masses in the range of $5-1000$M_⊙, where numerical results of the formation of BH are available. Traditionally, simulations done in GR usually i) stop when a trapped surface or an apparent horizon forms [for example [29], ii) `excise the central part of the evolution inside the Apparent Horizon (AH) [30,31,32,33,34], or iii) replace the BH interior with analytic solutions, for example, `punctures [35,36,37]. What happens inside the black hole is generally not of interest, since no information can come out of the black hole and does not influence the exterior dynamics. By ignoring the interior evolution, the singularity arising from the BH formation is effectively excluded, and the exterior evolution beyond the BH formation can be followed.

In contrast, here we are interested in what happens when the matter continues to collapse inside the BH. [15] showed how a spherically symmetric FLRW cloud collapses to form a BHU (called the Black Hole Universe). The collapse continues inside the BH. This is the starting point of our study, and we investigate the dynamics of the GR solution for the collapse of such a uniform, pressure-less spherical mass distribution with total mass M. For this configuration, the relativistic solution is mathematically equivalent to the corresponding Newtonian solution for a system with net zero total energy. We further assert that this kind of mathematical equivalence exists also for cases with nonzero pressure. To investigate the collapses and bounces of uniform spherical clouds of different masses quantitatively, we simulate their evolution numerically in Section 3. For simplicity, we assume a polytropic EoS. In particular, we investigate the initial conditions and EoS parameters required for the bounce and the properties of the latter. We conclude in Section 5.

2. Collapse of a Uniform Spherical Mass Distribution

To find the analytical equations that describe the dynamics of the collapse of a uniform spherical mass distribution with total mass M, we closely follow the approach presented in [27]. The spherically symmetric metric $g_{\mu \nu }$ in polar coordinates with a proper time $d{x}^{\mu }=(d\tau ,d\chi ,d\theta ,d\varphi )$ can be expressed by the line element [38,39,40]:

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-d{\tau }^2+e^{\lambda (\tau ,\chi )}d{\chi }^2+r^2(\tau ,\chi )d{\mathsf{\Omega }}^2.$$

(1)

Here, we use units of speed of light $c=1$, and the time and radial coordinates $(\tau ,\chi )$ are comoving with the matter (we use t to refer to rest frame or Newtonian coordinates). The corresponding metric is sometimes called the Lemaitre-Tolman [38,41] or the Lemaitre-Tolman-Bondi (LTB) metric. This metric is localized around a central reference point in space, which we have set as the origin ($\overrightarrow{r}=0$) for simplicity.

For an observer moving with a perfect fluid, the energy-momentum tensor is diagonal: $T_{\mu }^{\nu }=diag[-\rho ,P,P,P]$, where $\rho =\rho (\tau ,\chi )$ is the relativistic energy density and $P=P(\tau ,\chi )$ is the pressure. The off-diagonal terms of the field equation vanish, for example, $8\pi G{T}_0^1=G_0^1=0$. These equations can be reformulated into:

$$\left({\partial }_{\tau }\lambda \right)\left({\partial }_{\chi }r\right)=2{\partial }_{\tau }\left({\partial }_{\chi }r\right).$$

(2)

This equation can be solved with the ansatz: $e^{\lambda }=C{\left({\partial }_{\chi }r\right)}^2$, where ${\partial }_{\tau }C=0$. The particular case $C=1$ corresponds to a flat geometry:

$$d{s}^2=-d{\tau }^2+{\left({\partial }_{\chi }r\right)}^{2}d{\chi }^2+r^{2}d{\mathsf{\Omega }}^2.$$

(3)

Consequently, there is only one function required to be solved: $r=r(\tau ,\chi )$ and r corresponds to the radial proper distance. For a perfect fluid in the comoving frame:

$$T_{\mu }^{\nu }=diag[-\rho ,P,P,P]=\left(\begin{array}{cccc}-\rho & 0& 0& 0\\ 0& P& 0& 0\\ 0& 0& P& 0\\ 0& 0& 0& P\end{array}\right).$$

(4)

The field equations for r are then:

$$\begin{array}{ccc}\hfill H^2& \equiv & r_H^{-2}\equiv {\left({\displaystyle \frac{\dot{r}}{r}}\right)}^2={\displaystyle \frac{2GM(\tau ,\chi )}{r^3}},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill M(\tau ,\chi )& \equiv & {\int }_0^{\chi }\rho 4\pi r^2\left({\partial }_{\chi }r\right)d\chi ,\hfill \end{array}$$

(6)

where M is the mass inside $\chi$ at time $\tau$ and the dot here is $\dot{r}\equiv {\partial }_{\tau }r$. See Appendix A for a discussion on the Newtonian and relativistic masses. The remaining field equations correspond to energy conservation, which we discuss later in Equation (13). When $\rho =\rho \left(\tau \right)$ is uniform, we have $M=4\pi r^3\rho /3$ and:

$$r=a\left(\tau \right)\chi ,$$

(7)

where $a\left(\tau \right)$ denotes the scale factor of the collapsing cloud. Thus, the metric in Equation (3) reproduces the flat (global) FLRW metric:

$$d{s}^2=-d{\tau }^2+a^2\left(d{\chi }^2+{\chi }^{2}d{\mathsf{\Omega }}^2\right).$$

(8)

Equation (5) can then be expressed as $3H^2\left(\tau \right)=8\pi G\rho \left(\tau \right)$, and corresponds to a global solution. Since $\rho \left(\tau \right)$ does not depend on $\chi$, any point can be chosen to be the center $\overrightarrow{r}=0$. In general, non-flat geometries or more sophisticated global topologies can be reproduced with a more general metric, that is, $C\ne 1$. When considering spherical perturbations within a flat FLRW background, it is analytically appropriate to model the perturbation using a closed FLRW topology (i.e., with positive curvature, $k=+1$) in general relativity (GR). This approach accounts for the compensation of the surrounding background around the spherical perturbation (see Section 2.3). In particular, the positive local curvature in the FLRW metric reflects the curvature induced by the local overdensity. This helps maintain the isotropy and homogeneity of the perturbed solution, as it mirrors the curvature effects arising from the spherical symmetry of the perturbation.

2.1. Uniform Cloud

The Equations (5) and (6) can also be used to solve nonhomogeneous configurations. The simplest example is an expanding (or collapsing) uniform spherical mass distribution inside a comoving coordinate ${\chi }_{*}$, embedded in vacuum:

$$\rho (\tau ,\chi )=\left\{\begin{array}{cc}\rho \left(\tau \right)\hfill & \mathrm{for}\chi \le {\chi }_{*}\hfill \\ 0\hfill & \mathrm{for}\chi >{\chi }_{*}\hfill \end{array}\right..$$

(9)

Consequently, we have

$$M(\tau ,\chi )=\left\{\begin{array}{cc}{\displaystyle \frac{4}{3}}\pi r^3\rho \left(\tau \right)\hfill & \mathrm{for}\chi \le {\chi }_{*}\hfill \\ M_T\hfill & \mathrm{for}\chi >{\chi }_{*}\hfill \end{array}\right..$$

(10)

Here, we introduced

$$M_T\equiv M(\chi ={\chi }_{*})={\displaystyle \frac{4}{3}}\pi R^3\rho \left(\tau \right),$$

(11)

where $R\equiv a{\chi }_{*}$ is the boundary or surface of the cloud. Inside $\chi <{\chi }_{*}$, this solution is identical to the infinite FLRW solution $r=a\chi$ with the same $\rho$. Any off-center comoving observer inside ${\chi }_{*}$ still sees the same homogeneous and isotropic metric (Gaztañaga [15]). The Hubble-Lemaitre expansion law: $\dot{r}=H\left(\tau \right)r$ of Equation (5) with $3H^2\left(\tau \right)=8\pi G\rho \left(\tau \right)$ is also the same in both cases within $\chi <{\chi }_{*}$. The only difference is that for a mass distribution like Equation (9), the FLRW solution only holds inside; outside ${\chi }_{*}$, there is empty space. Consequently, the total mass is finite $M_T<\infty$.

The equivalence is given by Birkhoff’s Theorem (see Johansen and Ravndal [42], Faraoni and Atieh [43]), which states that a sphere cut out from an infinite uniform distribution has the same spherical symmetry as the total distribution. Thus, the FLRW metric is both a solution to a global homogeneous (i.e. infinite $M_T$) uniform background and also to the inside of a local (finite mass) uniform sphere centred around one particular point. The local solution is called the FLRW cloud (Gaztañaga [15]).

2.2. Newtonian Equivalence

The total energy of a fluid with zero total energy $E=0$, uniform density and zero pressure in a Newtonian framework can be expressed as:

$$\begin{array}{c}\hfill E=0={\displaystyle \frac{1}{2}}{\dot{\tilde{r}}}^2-{\displaystyle \frac{GM(t,\tilde{r})}{\tilde{r}}}\\ \hfill \Rightarrow {\dot{\tilde{r}}}^2\left(t\right)={\displaystyle \frac{8\pi G}{3}}{\rho }_{newt}\left(\tau \right){\tilde{r}}^2\left(\tau \right),\end{array}$$

(12)

where $\tilde{r}$ is here the radius of the spherically symmetric mass distribution and ${\rho }_{newt}$ is the usual Newtonian mass density. This Equation (12) is equivalent to Equation (5) when replacing the surface of the relativistic cloud at $\tilde{r}\to r=R$. Solving Equation (5) for $R=R\left(\tau \right)$ requires two more equations. One is provided by the EoS $P=P\left(\rho \right)$, and the other is the continuity equation for $\rho$. For a perfect fluid, the latter can be derived from GR energy conservation ${\nabla }_{\mu }{T}^{\mu \nu }=0$. The corresponding time component ($\mu =0$) for a uniform distribution gives:

$${\displaystyle \frac{d\rho }{d\tau }}={\displaystyle \frac{\partial \rho }{\partial \tau }}+3H(\rho +P)=0.$$

(13)

The radial equation ($\mu =1$) is:

$${\displaystyle \frac{\partial p(\tau ,\chi )}{\partial \chi }}=0,$$

(14)

which expresses the uniformity of the pressure in the comoving frame. We recall that we have chosen a comoving proper time $\tau$ so that $g_{00}=-1$. The angular equations $\mu =2,3$ are readily satisfied for a spherically symmetric fluid.

For vanishing pressure $P=0$, the energy conservation Equations (13) and (14) also apply to a uniform Newtonian cloud with M given by Equation (10) and replacing $r\to \tilde{r}$ and $\rho \to {\rho }_{newt}$. Thus, when written in the appropriate coordinates, the equations of motion describing the fully relativistic GR collapse of the flat FLRW cold cloud (defined by Equations (8) and (9)) are equivalent to the equations of motion of the corresponding Newtonian problem with $E=0$. This correspondence is well-known in the context of GR. Moreover, for the FLRW metric, this property does not depend on the particular choice of an EoS (see [43] for a comprehensive description of the correspondence of GR and Newtonian solutions).

As is well known, speeds greater than c are allowed in Newtonian mechanics and also appear in our simulations (see Section 3). For the description in the GR framework, velocities larger than c pose no problem because they can be interpreted in a novel way. Since we are considering a background solution, the faster than light velocities are related coordinate velocities $\dot{r}$, which do not carry physical information. The Newtonian equivalence for $P=0$ is evident, and the interpretation is straightforward in GR. How to interpret different EoS with $P=P\left(\rho \right)\ne 0$ within the Newtonian treatment in the GR solution requires more work. In any case, the description of the GR solution with the use of a metric does not allow for information to travel faster than c.

2.3. Spherical Collapse

Consider a spherical perturbation of radius R embedded in a background with mean density $\overline{\rho }$. The perturbed region has initiated their gravitational collapse and thus has a higher density $\rho =\overline{\rho }(1+\delta )$, which is initially only slightly higher than $\overline{\rho }$. We assume that the background with mean density $\overline{\rho }$ follows the evolution of a flat universe dominated by matter (i.e., EdS universe). To guarantee mass conservation, the increase in density in $r<R$ has to be compensated for by an area of low or vanishing density. This low density zone (or void) is generated when the high density perturbation has initiated its collapse from its initial radius $r=r_{\mathrm{out}}$ to its current value $r=R$ and left behind an empty space between $R<r<r_{\mathrm{out}}$, see Figure 1 for a graphical representation. The background could be expanding or collapsing at a lower rate. In any case, mass conservation implies:

$$M={\displaystyle \frac{4\pi }{3}}{\left(r_{\mathrm{out}}\right)}^3\overline{\rho }={\displaystyle \frac{4\pi }{3}}{R}^3\overline{\rho }(1+\delta ).$$

(15)

Due to Gauss’s theorem (or the corollary of Birchhoff’s theorem) in GR, see Johansen and Ravndal [42]), both spheres in Figure 1 follow the FLRW* metric (where * denotes the flat FLRW metric with the constraint $r=a\chi <R$), but each has a different expansion/collapse rate $H^2$. The outside has $k=0$, and the inside has $k>0$ (k is known as the curvature: $k=1/{\chi }_k^2$, where ${\chi }_k$ is the comoving radius of the curvature of the metric of Equation (20) below). The exact dynamics resulting from R in GR also agrees with Newton’s law:

$$\ddot{R}=-{\displaystyle \frac{GM}{R^2}}=-{\displaystyle \frac{4\pi }{3}}G\rho R=-{\displaystyle \frac{1}{2}}{\overline{H}}^2(1+\delta )R.$$

(16)

In the last equality, we have used the Friedmann Equation (5), where ${\overline{H}}^2={\overline{H}}_0^2{a}^{-3}$ refers to the background expansion rate: $\dot{r}=\overline{H}r$ and $\ddot{r}=-{\overline{H}}^{2}r/2$, where $r\ge r_{out}$. Combining with the second derivative of Equation (15):

$$\ddot{\delta }+2\overline{H}\dot{\delta }-{\displaystyle \frac{3}{2}}{\overline{H}}^2\delta ={\displaystyle \frac{4}{3}}{\displaystyle \frac{{\dot{\delta }}^2}{1+\delta }}+{\displaystyle \frac{3}{2}}{\overline{H}}^2{\delta }^2.$$

(17)

This equation is exactly the same in Newtonian Cosmology (Newtonian dynamics in expanding coordinates) and in the GR dynamics, as given by the Raychaudhuri equation for an irrotational and shear-free fluid (see Padmanabhan [44]):

$${\nabla }_{\mu }{q}^{\mu }={\displaystyle \frac{d\mathsf{\Theta }}{d\tau }}+{\displaystyle \frac{1}{3}}{\mathsf{\Theta }}^2=-R_{\mu \nu }{u}^{\mu }{u}^{\nu }=-4\pi G(\rho +3P)+\mathsf{\Lambda },$$

(18)

where $q^{\mu }$ is the acceleration four-vector, $\mathsf{\Lambda }$ is the cosmological constant, $\mathsf{\Theta }={\nabla }_{\nu }{u}^{\nu }$ and $u^{\nu }$ is the 4-velocity. The Raychaudhuri equation is the GR equivalent of the Newtonian Poisson equation. It is a purely geometric equation describing the evolution of the dilatation coefficient $\mathsf{\Theta }$ of a bundle of nearby geodesics in proper time $\tau$. Together with the continuity equation (mass conservation), it results in Equation (17). In other words, the Newtonian spherical collapse describes the actual GR dynamics when the tidal effects are neglected.

2.4. Dynamics Close to the Ground State

We assume here that the cold collapse to a classical singularity $\rho \to \infty$ will be stopped by the ground state of a matter degeneracy pressure, similar to what is expected from the Pauli Exclusion Principle in Quantum Mechanics. A ground state is characterized by a system dominated by its potential energy. In Appendix B, we show how for a generic system described by a single degree of freedom, which we represent by an effective scalar field $\psi$ in curved spacetime, such a ground state corresponds to Equation (A16). In this case, we see from the Equation (A14)-(A15) that $\rho =-P$ and the EoS of the fluid $P=P\left(\rho \right)=\omega \rho$ is characterized by $\omega =-1$.1 Combining Equation (13) with Equation (5), the acceleration in the FLRW metric is given by:

$${\displaystyle \frac{\ddot{a}}{a}}=-{\displaystyle \frac{4\pi G}{3}}(\rho +3P),$$

(19)

which should be compared to Equation (18). We are interested in the collapse of cold, non-interacting matter characterized by initially vanishing pressure $P=0$. Once the collapse reaches a density at which matter cannot be compressed further, the pressure rises and causes the bounce. Just before the bounce, $\dot{\psi }<<V$, and from Equations (A14)-(A15) we obtain $P=V\left(\psi \right)=-\rho$. This is also what happens during inflation or DE domination.

When P changes from zero to $-\rho$, $\ddot{a}/a$ in Equation (19) changes from negative (accelerated collapse) to positive, describing a bounce and the subsequent expansion. Note that Equation (19) is independent of the curvature k.

2.5. A Closed FLRW Perturbation

In the spherical collapse solution presented in Figure 1, we have assumed a closed local geometry for the perturbation with $k=1/{\chi }_k^2$. This is the standard approach in perturbation theory (see, [44,45] for example). Such perturbation can be modelled as a closed FLRW metric:

$$d{s}^2=-d{\tau }^2+a^2\left({\displaystyle \frac{d{\chi }^2}{1-k{\chi }^2}}+{\chi }^{2}d{\mathsf{\Omega }}^2\right).$$

(20)

Note, that this metric is only valid for $\chi <{\chi }_k$, i.e. inside the closed perturbation. Outside, we assume that the metric is flat, as in the spherical collapse configuration of Figure 1. In proper coordinates, the curvature radius is $a{\chi }_k$ so that in the limit of large a, we map to the flat case (infinite curvature radius). The Friedman Equation (5) for such GR perturbed spacetime is then given by:

$$H^2+{\displaystyle \frac{k}{a^2}}={\displaystyle \frac{8\pi G}{3}}\rho ,$$

(21)

where the value of k refers to the curvature scale at the time ${\tau }_0=0$. At this time, we choose $a=1$: i.e. $k=1/{\chi }_k^2$, where ${\chi }_k={\chi }_{*}$ is the comoving curvature radius in Equation (9). Since H is related to the collapse/expansion velocity $\dot{r}$, and the latter should vanish right at the time of bounce, we also need $H=0$ at this time. Therefore, at the bounce we can define ${\rho }_B$ as function of $a=a_B$:

$${\rho }_B={\displaystyle \frac{3k}{8\pi G{a}_B^2}}.$$

(22)

As explained in the subsection above, the conditions for the bounce are: change to positive acceleration, which requires $\omega <-1/3$ and is fulfilled by a ground state with $\omega =-1$ and zero velocity ($H=0$) at the bounce itself. When dominated by the ground state, $\omega =-1$, the expansion becomes exponential, but it turns into a regular expansion as soon as the density becomes smaller than the ground state ${\rho }_B$. In this work, we will study the Newtonian collapses of pressureless and polytropic fluids. In both the Newtonian and GR cases, we start from the same asymptotically free initial conditions. In the Newtonian framework, this corresponds to $E=0$, while in GR, it corresponds to $k=0$. The perturbation grows through gravitational instability, and we begin the numerical study at a point where the Newtonian solution already has some initial velocity, corresponding to a finite overdensity. At this stage, the GR equivalent problem has developed $k>0$ due to the finite curvature radius associated with a finite perturbation mass. Note that [15] focus on the junction conditions in Equation (9) for a flat $k=0$ case (which is the right approximation for our expanding universe or a flat expanding perturbation). The corresponding case with curvature during collapse was solved in §12.5.1 of [44]. We start from an initial perturbation with constant density up to $\chi ={\chi }_{*}$ in Equation (9), which means that ${\chi }_{*}={\chi }_k$ in the spherical collapse model. The junction conditions in this case correspond to:

$$r_s=2G{M}_T={(H_0/c)}^2{\chi }_k^3\Rightarrow k={[r_s/c/H_0^2]}^{-2/3}.$$

(23)

After the bounce, the curvature term is diluted by the exponential expansion and the expansion follows the flat case described above.

2.6. White Holes (WH)

When the bounce happens inside the Black Hole (BH), in both the Newtonian and GR solutions, the spherical perturbation expands and crosses back the gravitational radius $r_s$ from the inside out. This corresponds to a WH (White Hole) solution. Although this behavior is natural in the Newtonian framework, [46] has argued that such a solution cannot exist in a GR scenario bouncing inside a BH. The standard argument supporting WH solutions is based on the time symmetry of GR equations: Since the metric is quadratic in time, one can reverse the arrow of time and convert a BH solution into a WH solution. Although this is theoretically true, we cannot arbitrarily change our choice of the time arrow as the solution evolves.

For an external observer outside the black hole (BH), the arrow of time points towards collapse. Once the bounce occurs inside the BH, the solution cannot transform into a WH because, according to the Schwarzschild metric, once the collapsing time arrow is fixed, nothing can escape the Schwarzschild radius $r_s$. The resolution to this conundrum, as detailed in [46], involves imposing a boundary condition on the GR action, given that the internal dynamics are confined within $r_s$. See Appendix C for more details. This modifies Einstein’s field equations of the expanding solutions from Equation (5) to:

$$H^2={\displaystyle \frac{2GM(\tau ,\chi )}{r^3}}+{\displaystyle \frac{1}{r_s^2}}.$$

(24)

This boundary term is only needed for the bouncing cloud that falls inside its gravitational radius $r_s$. In our GR equivalent solution, we have ignored this term because it provides a very small contribution close to the bounce, but it is straightforward to include it if we want more accurate modelling after the bounce. which prevents the expansion from crossing $r_s$. This term corresponds to a reinterpretation of the measured $\mathsf{\Lambda }$ in cosmic expansion as a boundary term in the BHU solution: $r_s=\sqrt{3/\mathsf{\Lambda }}$. According to this approach, a WH solution cannot emerge out of a BH. After the bounce, the expansion inside the BH is initially exponential, dominated by the ground state’s negative pressure. It then transitions to being dominated by the more regular expanding density (which, in a more realistic scenario, will also contain radiation density) and eventually becomes dominated by the $\mathsf{\Lambda }$ boundary term, which is effectively determined by the initial mass $M_T$ of the collapsing perturbation: $r_s=2G{M}_T=\sqrt{3/\mathsf{\Lambda }}$. In the comoving frame, the expansion accelerates, but in the rest (de Sitter) frame, the expansion becomes asymptotically static within $r_s$ as given by de Sitter metric (see Gaztañaga [47] for further details).

3. Numerical Simulations

As we have seen in the previous section, the solution for the collapse of an FLRW cloud is equivalent to solving the Newtonian problem outlined by Equation (12), utilizing Equation (13) for vanishing pressure P. In this section, we make use of this equivalence and investigate the collapse numerically with the CASTRO2 code.

In Section 3.1, we first concentrate on the pressureless collapse, for which we again highlight that the Newtonian formulation is not an approximation of the relativistic problem. Rather, both formulations are mathematically equivalent to each other. In GR, r is interpreted as a surface coordinate $r=a\left(\tau \right)\chi$ rather than a physical distance as in the Newtonian framework. In the Newtonian solution, velocities $\dot{r}>1$ larger than the speed of light are allowed, and they indicate the formation of a Black Hole. In the GR formulation, the coordinate distance r is not related to the propagation of information. In contrast, the world lines of actual physical events or objects are described by Equation (8) or the corresponding metric, and their propagation follows causality by definition. In this context, we note that, generally, no BH is formed in Newtonian simulations. However, we can estimate when a BH would form in GR as the time $t_{\mathrm{BH}}$ when the matter of a given mass is located within its Schwarzschild radius $r_s\equiv 2GM/c^2$ in the Newtonian solution.

Once we have reproduced the pressureless case, we investigated the collapse of a polytropic fluid. Here, the pressure is negligible initially, and the analogy between Newtonian and GR solutions holds approximately. However, at some point, we postulate a strong interaction of quantum particles based on the Pauli principle, which leads to a significant increase of the pressure and finally to a bounce of the infalling matter. We then tentatively discuss the approximate Newtonian-GR equivalence for the second case in Section 4.1

3.1. Pressureless Collapse

In the Black Hole Universe (BHU) interpretation, the observable Universe is expanding inside its own gravitational radius. This suggests that our observable Universe could originate from the collapse and subsequent bounce of a very large cloud. Given that the total mass $M_T$ is extremely large (of the order of $M_T\simeq 5\times {10}^{22}{M}_{\odot }$), the density of the cloud at black hole formation ${\rho }_{BH}$ would be correspondingly very low:

$$\begin{array}{cc}\hfill {\rho }_{BH}& ={\displaystyle \frac{M}{V}}={\displaystyle \frac{3M}{4\pi r_S^3}}={\displaystyle \frac{3r_S^{-2}}{8\pi G}}\simeq 3.9\times {10}^{-48}{\left[{\displaystyle \frac{M}{M_T}}\right]}^2{\displaystyle \frac{M_{\odot }}{{\mathrm{km}}^3}}\hfill \end{array}$$

(25)

For our Universe, this density is on the order of 1 atom per m³ and motivates the idea that such an initial cloud existed within a very large manifold with comparable or even lower density. The original density must have been so low that the intermolecular separation was much larger than any potential source or interaction term, which implies that it was also extremely cold. In such a manifold (which was called the Apollonian Universe by [16,28]), small perturbations could give rise to a new type of primordial stars and black holes in hierarchical structures of different sizes. This perspective leads us to consider that the starting point of the Universe originated from such dilute, cold ($T=0$ K) clumps of different sizes made of non-interacting, neutral hydrogen atoms.

There is no specific reason to assume that this primordial cloud was composed of neutral hydrogen atoms, and other initial compositions are certainly possible. However, we select this scenario for concreteness, as the results do not appear to depend significantly on the exact initial composition as long as the collapse remains cold. Although absolute zero temperature is not physically realizable, assuming a very low initial density makes this a reasonable approximation. Additionally, we assume that the cloud is spherically symmetric, homogeneous, non-rotating, and non-dissipative.

During the homologous collapse of this clump, the atoms feel each other only through their gravitational interaction, and if initially $T=0$K, no temperature will be created. As the collapse proceeds, the dust cloud reaches densities at which the wave functions of electrons in nearby atoms overlap. At this point, electrons become more and more degenerate as a result of the Pauli exclusion principle. Consequently, a pressure component is built up and ionizes the atoms. This phenomenon is completely temperature-independent and is called pressure ionization [49]. This electron degeneracy pressure is responsible for the stability of White Dwarfs and allows hydrostatic equilibrium in stellar cores until they reach masses greater than $\approx 1.44$M_⊙, the famously-known Chandrasekhar mass limit [50]. For more massive, spherically symmetric mass distributions, the electron degeneracy pressure is insufficient to prevent gravitational collapse. In our simulations of collapsing clumps starting at masses of at least $5$M_⊙, which should ensure that electron degeneracy pressure does not affect the collapse significantly. With masses significantly above the Chandrasekhar limit, our configurations are too massive to be influenced. Additionally, we start the collapse at very low densities ($E=0$), and once the densities reach values where electron degeneracy becomes important, the atoms have already acquired a significant infall velocity. Therefore, the collapse cannot be stopped by the electron degeneracy pressure.

Since it is not physically realizable, the particular case of a pressureless fluid is generally not included in hydrodynamical codes and causes numerical problems. However, practically a very small pressure does not affect the evolution of the collapse of a dust cloud. We follow Monchmeyer and Muller [51] and use a gamma-law EoS with an initial pressure

$$P_0<<{\displaystyle \frac{4\pi G}{\gamma }}{\rho }_0^2{r}_0^2,$$

(26)

where $\gamma$ is the relation of specific heats of the gas, and ${\rho }_0$ and $r_0$ are the initial, uniform density and radius of the cloud, respectively. In order to save numerical resources, we do not start the simulations right at the onset of the collapse. We rather approximate the initial phase of the collapse, where pressure can be neglected, with the free fall solution for $P=0$. The velocity at any point of the evolution is given by $u_0=c\sqrt{r_s/r_0}$. The collapse time of the analytical solution is

$${\tau }_c={\displaystyle \frac{2\sqrt{(r_0^3/r_s)}}{3c}}$$

(27)

and will serve as the reference time for our numerical simulations.

We set up our initial data inside a simulation box of radius $r_b$ and place a uniform-density (${\rho }_0$) sphere of radius $r_0<r_b=300\times {10}^6$cm in the center. The exterior of the sphere is filled with a very low-density medium (${\rho }_{\mathrm{amb}}$). In between the two regions, the transition profile is defined as [48]:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{newt}}={\rho }_0-{\displaystyle \frac{{\rho }_0-{\rho }_{\mathrm{amb}}}{2}}\left[1+tanh\left({\displaystyle \frac{r-r_0}{h}}\right)\right],\end{array}$$

(28)

where h is the smoothing length, which we fix to be $h=4\times {10}^4$ cm in our simulations.

We simulate the cold or pressureless collapse of an $8$M_⊙, spherically symmetric, and homogeneous dust cloud starting at $r=100r_s$. We obtain the corresponding values of initial density and collapse time at that radius from the theoretical solution for a pressureless collapse. The corresponding values are given in the first row of values of Table 1. We set the initial pressure to $P_0={10}^{10}$dyne cm⁻²$<<{\displaystyle \frac{4\pi G}{\gamma }}{\rho }_0^2{r}_0^2=2.26\times {10}^{27}$ dyne cm⁻². We stop the simulation at the estimated collapse time $t_c=0.0532$ s. The simulation uses a CFL number of $0.5$, and the radial resolution is set to 10240 grid points. At the outer boundary, we use outflow boundary conditions, while at the center, symmetry conditions are applied. The gravitational potential is calculated using monopole gravity ${\varphi }_G\left(R\right)=-GM/R$, which represents the spherically symmetric component of the gravitational field.

Results

In the upper panel of Figure 2, we plot the density as a function of the radius r. The initial configuration (blue curve, $t=0.0t_c$) is a solid sphere with $r=100r_s$. The sphere has a sharp transition to the low-density exterior. During evolution, the density increases and the radius decreases as expected. The density within the sphere remains roughly constant until the mass is concentrated within $r=r_s$. At this point, the resolution in the center limits the accuracy of our simulation. In addition, unavoidable numerical diffusion leads to a smearing of the initially sharp edge between the high-density sphere and the ambient matter. However, the edge is still qualitatively well defined and visible at all times. The density jump at the edge decreases from 12 orders of magnitude in the beginning to 3 orders at $t=t_c$. Note that the exterior also evolves with time.

In the last time step $t=t_c$, the densities reached supranuclear densities ${\rho }_{\mathrm{newt}}\approx {10}^7{\rho }_0=2.8\times {10}^{15}g/c{m}^3$. As the simulation progresses, the total mass of the system is concentrated in smaller and smaller radii (central panel of Figure 2). Once the radius of the sphere is comparable to $r_s$, the limited numerical resolution causes a loss of some of the mass. We could increase the resolution or use adaptive methods. However, our interest here is the time up to the collapse, which is reproduced well. The formation of a black hole is indicated by a free-fall velocity $\left|u\right|>c$. As we can see in the bottom panel of Figure 2, this is reached around $t=0.99964t_c$ (green curve).

We can compare our simulations with the theoretical results predicted by GR. For this purpose, we plot the radius and the corresponding density of the collapsing spherical cloud as functions of time in Figure 3. From the Friedmann Equation (12), we get the radius $r\left(\tau \right)$:

$$\begin{array}{c}\hfill r=r_0{\left(1-{\displaystyle \frac{\tau }{{\tau }_c}}\right)}^{2/3}={\left({\displaystyle \frac{3}{2}}({\tau }_c-\tau )\right)}^{2/3}.\end{array}$$

(29)

In our simulations, the radius is defined as the location of the edge of the mass density profile, and the density of the sphere is approximated as the average of all densities inside the sphere enclosed by this radius. The comparison of the simulation (magenta dots) with the analytical prediction (blue line) for the radius (density) is shown in the top (bottom) panel of Figure 3. We obtain a good agreement up to short before $t=t_c$. In the theoretical GR solution, r tends to zero, and the density ${\rho }_R$ increases towards infinity. We anticipate that in the presence of non-zero pressure, the GR solution will not be exact anymore, and the collapse time will be longer. This is already reproduced in our numerical results, which include a small but unavoidable pressure. The results of the simulations lag behind the theoretical curve after $t>0.999t_c$.

We have seen that the numerical resolution in our simulations affects the total mass inside r, and the edge of the sphere becomes less sharp with time. In order to take this into account, when considering the formation of BHs in the GR analogy, we take a more conservative approach. Instead of checking whether the entire mass is within its Schwarzschild radius, we check for each radius whether the mass $M_{\mathrm{encl}}\left(r\right)=4\pi {\int }_0^r\rho \left(\tilde{r}\right){\left(\tilde{r}\right)}^{2}d\tilde{r}$ enclosed by a given radius r has reached its Schwarzschild radius $r_s\left(M_{\mathrm{encl}}\right)=r_s\left(r\right)$ or whether it is still more extended $r>r_s\left(r\right)$. In Figure 4, we plot the Schwarzschild radius of the mass distribution with $8$M_⊙ at each radius at different times, both rescaled by the Schwarzschild radius of the entire mass $r_{s0}$. Initially (blue line), the radii are much larger than the corresponding $r_s$. With time, the collapse proceeds, and the radii approach the Schwarzschild radii. Once the coloured curves, representing the evolution, cross the black dashed line where $r=r_s$ a BH would form in the GR case. This has happened for the last two steps shown, and the system’s total mass is located inside its Schwarzschild radius. As we will see later in our simulations, some of the models bounce back before $r=r_s$ is reached. We find that a threshold value of $95\%$ of the mass is a good indicator to see whether the collapse of the mass distribution has reached $r=r_s$ or not.

3.2. Collapse with a Polytropic EoS

Even for a zero-temperature gas, the negligible pressure approximation applied in the previous section breaks down at some point. First, electrons become degenerate and create temperature-independent pressure. By performing simulations with a zero-temperature-white-dwarf EoS, we found that the collapse for sufficiently high masses ($M\ge 5$M_⊙), is basically indistinguishable from the zero pressure case. However, already at nuclear saturation densities, nucleon degeneracy can halt the collapse and initiate a bounce for masses up to several ${10}^{\prime }s$ of M_⊙. This scenario is actually observed at the end of the life of massive stars as supernova explosions. Analogously, we assume that for more massive objects, a yet unknown pressure, based on the Pauli principle for quantum objects at the subnuclear scale, prevents further collapse and leads to a bounce. To simulate this effect, we investigate the behaviour of collapsing clouds with a polytropic EoS. For a perfect fluid, the latter is given by:

$$P=K{\rho }^{\gamma },$$

(30)

where K is the polytropic constant and $n=1/(\gamma -1)$ is the polytropic index. This kind of EoS can be a good approximation for nuclear degenerate matter with $\gamma =2\dots 3$ (see [52]). Though not taking all subtleties of the nucleonic interaction properly into account, the polytropic EOS has the advantage of being simple. Since we also do not know what happens at much higher densities than nuclear saturation density, we can encode all our ignorance on what may happen at supra-nuclear densities into the parameters of a polytropic EoS. We do not aim to explain the realistic bounce of a collapsing cloud of matter under these conditions. However, the polytropic EoS allows us to study qualitatively under which conditions a bounce may happen and what the properties of this bounce are.

For the $\gamma$ values in the range $2\dots 3$, we choose the K parameter of the polytropic EoS such that the pressure resembles a strong repulsion once the nuclear densities are reached. The combinations of $\gamma$ and K are given in the 3rd and 4th columns of Table 1. The EoS are constructed so that they have the same pressure at $\rho ={10}^{14}$g cm⁻³, which is slightly lower than the nuclear saturation density ${\rho }_{ns}$. This choice ensures that the evolution until reaching densities of the order of ${\rho }_{ns}$ is not affected significantly by the pressure and that all simulations are comparable until these densities. The larger $\gamma$, the stiffer the EoS and the stronger the pressure increases with density. So, for larger $\gamma$, we would expect a bounce to occur at lower densities. Similarly, bounces at higher (lower) pressures for a given density can be achieved by choosing a lower (higher) value of K.

For the study of the effect of the different polytropic EoS, we consider three different masses $M=5$M_⊙, 20M_⊙ and 1000M_⊙. The other initial conditions of the mass distributions are chosen similar to the ones used in subSection 3.1. We place a barotropic fluid with given $\gamma$ and K of uniform density ${\rho }_0$ in a sphere of radius $r_0$. This sphere is placed inside a simulation box of radius $r_b>1.1\times r_0$ filled with a very low density, ${\rho }_{\mathrm{amb}}$. The temperature is $T=0K$, and the initial pressure is obtained from the polytropic EoS. The initial quantities are given in Table 1.

To compare our results to the previously obtained ones, we use the time of the collapse (Equation (27)) of a pressureless cloud as a reference. However, now we do not stop the simulations at ${\tau }_c$ and continue the simulation beyond it. For the polytropic EoS, the simulations inevitably create a bounce. The pressure increases with density to arbitrarily high values, which does not allow for further compactification. Consequently, at some point, the pressure obtained during the collapse of any mass configuration will reach such high values that it can balance the gravitation and lead to a bounce. Once the bounce occurs, we may encounter a number of numerical problems, like strong oscillations of the central density or the bounced mass starting to leave our predefined outer grid boundary $r_b$. In both cases, we stop the analysis, even if the simulation continues. The time-stepping is controlled by the CFL condition with a CFL number of $0.5$, and the spatial resolution is set to 10240 grid points. As before, symmetry conditions are used at the centre, and outflow boundary conditions are set at the outer boundary of our grid.

Note that collapses of polytropic mass distributions have already been studied by [53]. They find that for $5/3>\gamma >4/3$, the system can either oscillate and be in stable equilibrium or lead to a bounce. However, the configuration with $\gamma =5/3$ invariably leads to a bounce in the Newtonian case. However, the initial conditions studied by [53] are different from the ones we use here. In our case, the pressure of the whole cloud is determined by the polytropic EoS, whereas Bondi uses a polytropic EoS only for the center of the cloud. Additionally, the equations used in [53] work only outside the gravitational radius. However, our results are in line with the ones found by Bondi, supporting the close correspondence between Newtonian and GR results.

Results

We performed a total of nine simulations consisting of models with three different masses and each with three different polytropic indices (see Table 1). In Figure 5, we plot the central densities as a function of time. The black line is the density prediction of the zero-pressure case. As we can see, the curves of all nine configurations follow the black line until after we reach densities ${\rho }_{\mathrm{newt}}>{\rho }_{\mathrm{ns}}$. Initially, the pressure from the polytropic EoS is still very low, and the evolution is as if there was no pressure. Note that in the figure for 1000 M_⊙, we start the simulation for radii $r<r_s$, such that in GR, the simulations would have started within a BH from the beginning. In all cases, the central densities rise to values higher than ${\rho }_{ns}$, and then, sooner or later, all models bounce back. This bounce can be recognized as a drop in density. Note that an apparent higher density in our polytropic simulations, compared to the pressureless case, is simply caused by a limited number of the time snapshots we used for drawing the coloured curves, while the black line is analytic.

Similarly to what we have discussed before for numerical viscosity, the pressure introduces a complication when determining the radius of the initially sharp edge of the mass distribution. With increasing pressure, the outermost matter falls in slower than in the pressureless case, and the edge smears out. Therefore, when checking whether a potential BH may form in a GR solution for these configurations, we define $r_{95}$ as the radius enclosing $95\%$ of the total mass and plot its ratio with the Schwarzschild radius of this mass $r_{95}/r_{95s0}$ as a function of time in Figure 6. When $r_{95}/r_{95s0}=1$, we expect that a BH forms in GR. In all models, the radius starts to expand after an initial phase with decreasing $r_{95}$. We call the time when this happens bounce time $t_B$. For almost all models, the bounce occurs within the Schwarzschild radius. Only the model with $M=5$M_⊙, and $\gamma =3$ bounces before reaching $r_s$. The bounce happens so early for this model because the pressure grows the fastest for $\gamma$ = 3, and, in the case of a $M=5$M_⊙ sphere, it can counteract the gravitational pull before total mass falls into its gravitational radius $r_S$.

The parameters obtained for $t_B$, density, and bounce pressure for all models are given in Table 2. The densities (and pressures) increase with mass and decrease with increasing $\gamma$. Consequently, the configuration with $M=1000$M_⊙ and $\gamma =2$ reaches the highest densities and pressures. We expect this trend to hold for larger masses. When extrapolating to a configuration with a mass as large as that of the entire universe, we expect to reach densities many orders of magnitude higher than the nuclear saturation density at bounce. This is far beyond any reachable energies we can study experimentally today. Further work in high-energy physics is needed to develop an EoS to describe the matter at those densities.

In addition, higher pressures and densities are reached for a given mass for a lower $\gamma$ value or a higher polytropic index. This is intuitive since a higher $\gamma$ means a higher stiffness. The cloud cannot be deformed as much as in softer EOS, with lower $\gamma$. In addition, the bounce time for a given mass strongly depends on the polytropic index. It takes longer for the bounce to occur for a lower polytropic index, since the pressure increase is slower.

As we have seen in this section, in the Newtonian problem, there is no problem for the matter to expand past its Schwarzschild radius after it had collapsed to within $r<r_s$ (see Figure 6).

4. Discussion

4.1. Newtonian vs GR Equivalence

As discussed in subSection 2.2, the equations of motion describing a GR collapse of a flat FLRW cloud have mathematically equivalent Newtonian equations for a system with zero total energy. In the absence of pressure, the GR solution can be satisfactorily reproduced by Newtonian simulations, as we showed in Section 3.1. When the pressure does not vanish, the situation is more complex. The interpretation of the Newtonian solution in terms of the GR equivalent is not straightforward and requires further work.

In this work, we focus on the dynamics of spherical collapse in Section 2.3, where tidal effects are absent. As argued from Equation (17), the Newtonian spherical collapse accurately captures the full GR dynamics for a perfect fluid with an arbitrary EoS as long as pressure and density are homogeneous. In GR, our aim is to simulate a flat ($k=0$) background, which corresponds to $E=0$ in the Newtonian solution. However, an important distinction arises once we introduce a perturbation. A spherical GR perturbation on a flat background is characterized by a positive local curvature $k=1/{\left(a{\chi }_k\right)}^2$ in the Equation (20). In the corresponding perturbed Newtonian solution, $k\ne 0$ would introduce a constant energy term $E=R_0/{\chi }_k$ (the relationship between the initial radius $R_0$ and the initial curvature radius ${\chi }_k$ at $a=1$). This would pose a problem if we want to compare the Newtonian to a flat background GR solution because, in Newtonian physics, energy conservation is a fundamental principle. In contrast, in GR, there is no direct counterpart. To ensure that the Newtonian simulation aligns with the flat-background GR space asymptotically, we retain the $E=0$ initial condition. This distinction becomes vital when deriving the equivalent GR EoS from a given Newtonian EoS, as the curvature and energy terms are handled differently between the two frameworks.

Given the equivalence between the GR and Newtonian equations, under certain approximations, we can estimate what to expect in GR from the results of the equivalent Newtonian simulation. This kind of mapping allows us to leverage the simplicity of Newtonian simulations while gaining insight into the corresponding GR behaviour, providing a useful approximation for more complex relativistic scenarios.

Let us now estimate the qualitative behavior expected for the GR solution from our Newtonian simulations with a non-zero-pressure EoS. As a first step, we use Newtonian simulations with a polytropic EoS to estimate the scale factor $a\left(t\right)$. Using Equation (7) and the initial condition $a\left(0\right)=1$, we have $a\left(t\right)=R\left(t\right)/R_0$. With $a\left(t\right)$, we can estimate the Hubble rate using Equation (5): $H=\dot{a}\left(t\right)/a\left(t\right)$. Once H is determined, we can calculate the relativistic density $\rho \left(t\right)$ using the Friedmann Equation (21). Finally, the pressure $P\left(t\right)$ can be obtained from the continuity Equation (13). With P and $\rho$ given, the relativistic EoS is determined.

We discuss our findings as exemplary for the model $M=20$M_⊙ and $\gamma =2.5$. Before discussing the results, we note that finding the exact radius $R\left(t\right)$ where the collapsing shell terminates at a given time is tricky once $P\ne 0$. The reason is that once there is pressure, the density does not remain constant within the cloud, and the matter at the outer edge of the shell collapses somewhat slower than in the centre. This is depicted in Figure 7, where we show the density as a function of radius for a few time steps. There is still a region of almost constant density in the very center of the cloud. However, this constant density region no longer accounts for all the matter of the cloud. With subsequent collapse, the mass having constant density decreases further and, at the time of bounce, the density decreases monotonically from the inside out (cyan curve in Figure 7). Actually, this moment of maximum central density is only the onset of the bounce, and most of the matter is still in-falling. The process takes time, and a hump in the density moves outwards and marks the region where bounced matter encounters the still in-falling material. Note that since the density is not constant, the assumption of an FLRW cloud no longer holds, and our results and interpretations are only approximations.

Without a well-defined edge of the density profile, we estimate the equivalent radius of our collapsing non-zero-pressure cloud by using four-thirds of the mass-weighted radius

$$R_{\mathrm{equiv}}={\displaystyle \frac{4}{3}}{\displaystyle \frac{{\int }_0^{r_b}\left(4\pi \rho r\right)r^{2}dr}{M}}.$$

(31)

With this $R_{\mathrm{equiv}}$, we can calculate the scale factor and the Hubble rate. For our 20M_⊙ model Equation (23) gives $k=6.29\times {10}^6$ s⁻². This k corresponds to a curvature radius of ${\chi }_k\simeq 1.19\times {10}^7$ cm, which is lower than the Hubble Horizon $c/H_0$ at the initial condition: $H_0=H(t=0)\simeq 1.69\times {10}^7$ cm. This means that, initially, the cloud is still inside its Hubble horizon, and everything is causally connected.

To estimate the behaviour of the bounce in GR, we plot different representative densities in Figure 8. The dotted blue line represents the exact FLRW solution of the pressureless collapse, and the solid blue line represents the corresponding central density of the Newtonian simulation with the polytropic EoS. As expected, the pressure has a vanishing influence for $t\ll t_c$. However, once approaching $t\sim 0.96t_c$, ${\rho }_{\mathrm{newt},}\mathrm{central}$ decreases due to the bounce of the matter at the centre. As we have seen in Figure 7, the density is not constant, and to approximate the GR solution better, we define a mean density:

$${\overline{\rho }}_{newt}={\displaystyle \frac{3}{4}}{\displaystyle \frac{M_{\mathrm{encl}}}{R_{\mathrm{equiv}}^3}}.$$

(32)

For early times, this ${\overline{\rho }}_{newt}$ is identical to the central density, but the closer the simulations come to the bounce, the more it deviates. It reaches its maximum around $t\sim t_c$, indicating that the bulk of the matter has bounced by that time.

The approximated relativistic energy density has a different behavior (solid black line in the figure). It shows two maxima just before and after $t=t_c$. To understand this result, we show the two contributions to $\rho ={\rho }_1+{\rho }_2$ explicitly. Here, ${\rho }_1=3{\dot{a}}^2{\left(8\pi G{a}^2\right)}^{-1}$ is the kinetic contribution or the term related to the Hubble rate (dashed black line) and ${\rho }_2=3k{\left(8\pi G{a}^2\right)}^{-1}$ (dot-dashed, black line) is related to the curvature. In contrast to the other curves, ${\rho }_1\left(t\right)$ has two peaks. This is caused by the Hubble rate, which follows the dynamics of the collapse. Initially, the cloud collapses slowly, so H is small and negative. The collapse accelerates until the bounce of the very central matter sets in, causing H to slowly decrease to a (negative) minimum. At this point, some of the matter changes direction and, therefore, the magnitude of the Hubble rate (and ${\rho }_1$) decreases. Around $t\sim t_c$, H and ${\rho }_1$ become zero, marking a critical point during the bounce. Subsequently, the matter continues to accelerate outwards and both H and ${\rho }_1$ increase. However, there is still matter falling, causing the acceleration of the expansion to slow down until the (positive) maxima of H and ${\rho }_1$ are reached. Finally, all matter has positive velocities and slows down during expansion. The amplitudes of H and ${\rho }_1$ decrease. The magnitude of the negative minimum, corresponding to the infall velocity, is higher than that of the positive maximum, representing the bounce expansion. Consequently, the first peak of ${\rho }_1$ is larger than the second. At $t\sim t_c$, when ${\rho }_1=0$, $\rho$ is dominated by ${\rho }_2$. This contribution behaves very similar to ${\overline{\rho }}_{\mathrm{newt}}$ and has a maximum around $t\sim t_c$. That ${\rho }_2$ has its maximum right when ${\rho }_1=0$ in between its two maxima creates the two-peaked structure of $\rho$. Note that initially ${\rho }_2$ is larger than the Newtonian equivalents (pressureless solution or ${\overline{\rho }}_{\mathrm{newt}}$), and it increases slower with time. This is caused by our choice of k in the relativistic case. $k\ne 0$ leads to the dependence ${\rho }_1\sim a^{-2}\sim R^{-2}$ in contrast to ${\overline{\rho }}_{\mathrm{newt}}\sim R^{-3}$ for the Newtonian case. ${\rho }_1$ is also initially smaller than ${\rho }_2$. They intersect for the first time when $r=r_S$. At this point ${\rho }_1={\rho }_2$, which implies that $k={\dot{a}}^2$ and consequently $\dot{r}=\dot{a}\chi =c$, which marks the formation of a BH in GR and shows that the bounce happens inside the BH.

The different dependencies of ${\overline{\rho }}_{\mathrm{newt}}$ and ${\rho }_2$ lead to very different asymptotic behavior for very large a (or R). This difference arises from the nonvanishing k in the GR solution, which is a purely geometrical term absent in the Newtonian ($E=0$) solution. Only with this term, alongside the ground-state pressure (such as neutron degeneracy pressure), which causes the bounce in the Newtonian framework, can matter rebound in GR as well. To see this, note how $P=-\rho$ in GR results in $\ddot{a}>0$ (see Equation (19)). This is a necessary but not sufficient condition for a bounce to occur. We also need $k>0$ (in the Equation (21)) to satisfy the additional bounce condition $\dot{a}=0$ while $\rho >0$.

In Figure 8, we also plot the absolute value of the relativistic pressure. P is equal to the negative density ($P=-\rho$) on the two peaks of $\rho$. This happens because $\dot{\rho }=0$ at these points and, consequently, $P\left(\rho \right)=-(\rho +\dot{\rho }/3H)=-\rho$. Furthermore, when $H(t\sim t_c)=0$, the pressure diverges and ${\rho }_2$ and ${\overline{\rho }}_{\mathrm{newt}}$ peak. This is what we expect from a speculated quantum degeneracy pressure: a bounce caused by a mechanism similar to the Pauli exclusion principle.

We may also interpret this behavior with the help of the analogy to a relativistic scalar field $\psi$ in curved spacetime (see the Appendix B).

When such a system is trapped in a ground state, which is defined by $V\gg K={\dot{\psi }}^2$, the Equations (A14) and (A15) read $P=K-V$ and $\rho =K+V$. Here, $K\left(\psi \right)$ and $V\left(\psi \right)$ are the kinetic and potential energies of the scalar field, respectively. When $V\gg K={\dot{\psi }}^2$, the relativistic EoS can be expressed as $P=-\rho$. The pressure equal to the negative energy density is exactly what we have found in Figure 8. When $\rho$ has its two maxima, we find that $P=-\rho$. After the first peak of $\rho$, P increases further to negative infinity when $H=0$. This huge negative pressure stops gravitational collapse and causes a rebound or bounce. Once the matter expands, the pressure reduces and tends towards $P=0$. During this phase, the curve passes again through $P=-\rho$ at the second peak in the relativistic density plot. The analogy with the scalar field cannot explain how to obtain $P<-\rho$ in between the two density maxima. However, we can see the similarity and consider the case where $P\le -\rho$ is the ground state of our more complex system.

4.2. BHU vs Cosmological Coupling

This paper builds upon the theoretical foundations of the Black Hole Universe (BHU) model, and our results support its predictions. However, the BHU model is not the only alternative explanation for the accelerated expansion of the universe. One compelling proposal is the concept of "cosmological coupling," which posits that certain objects, such as black holes (BHs), are coupled to the large-scale dynamics of the universe [54,55,56,57].

Recent work by [58] presents observational evidence that black holes can grow in mass without consuming gas or stars. This suggests that black holes might be a significant source of dark energy. Specifically, the study found that the mass of supermassive black holes (SMBHs) at the centers of elliptical galaxies increases in a manner consistent with cosmological coupling. As the universe expands, the mass of these black holes grows proportionately, contributing to the overall energy density of the universe. Furthermore, [59] provides theoretical support for cosmological coupling by demonstrating that an exactly static BH event horizon cannot be embedded in a time-dependent geometry like the FLRW metric, inherently linking BHs to the universe’s expansion.

However, if black holes form through the spherical collapse process depicted in Figure 1, they create an empty transition region between the black hole (static only outside the Schwarzschild radius, $r_s$) and the surrounding FLRW metric. According to the corollary of Birkhoff’s theorem, this region would follow the Minkowski metric. Such a transition would effectively decouple the black hole from the FLRW metric, analogous to how regular perturbations, such as dark matter halos, become virialized or collapse without contributing to dark energy coupling or global background dynamics. Additionally, [60] argues that non-singular black holes or horizonless compact objects in General Relativity (GR) can exhibit cosmological coupling but are unlikely to act as a source of dark energy. This conclusion was supported by analyzing the redshift dependence of mass growth in supermassive black holes within a sample of elliptical galaxies and comparing their theoretical predictions with higher redshift mass measurements obtained from the James Webb Space Telescope (JWST).

While both the BHU and cosmological coupling models challenge conventional theories and provide fresh perspectives on the role of black holes in the universe’s expansion, they offer distinct approaches and implications. The BHU model proposes that the universe itself can be described as a black hole, obviating the need for dark energy by suggesting that the expansion of the universe, driven by a non-singular bounce, behaves like a cosmological constant. This framework is consistent with an asymptotically de Sitter interior metric that smoothly transitions into an exterior Schwarzschild metric without discontinuities. In contrast, the cosmological coupling model attributes the mass growth of black holes to their interaction with the universe’s expansion, suggesting that this growth contributes to dark energy, which drives accelerated cosmic expansion.

Interestingly, our BHU could be one of many black holes within a larger manifold, potentially embedded in an expanding or contracting FLRW metric. This external metric could possess its own $\mathsf{\Lambda }$ or dark energy term, different from the one measured inside our BHU, which is given by $\mathsf{\Lambda }\simeq 3/r_s^2$. If the cosmological coupling between our BHU and the external FLRW metric were a significant effect, it could influence our BHU mass and, consequently, the associated local term $\mathsf{\Lambda }\simeq 3/r_s^2$. In this scenario, our local $\mathsf{\Lambda }$ would evolve over time, either decreasing or increasing (corresponding to an effective ${\omega }_0<-1$), depending on whether the external FLRW metric expands or contracts. Such a framework provides a potential explanation for the recent evidence of ${\omega }_0<-1$ reported by DESI2024 measurements (DESI Collaboration et al. [61]), without necessitating the introduction of ghost-like dark energy models.

4.3. Potential Signatures in the CMB

Although the background equations for the Big Bang model and the Black Hole Universe (BHU) are identical—since the ${\mathsf{\Lambda }}_{EH}$ (or event horizon boundary) lies beyond the observable Universe—differences may arise in the origin and evolution of perturbations, as detailed in [62]. Specifically, the size of super-horizon perturbations could be constrained by the BHU boundary during collapse. This constraint might explain observed anomalies in the Cosmic Microwave Background (CMB), such as the lower quadrupole and octopole moments, or the absence of perturbations on angular scales larger than 60 degrees.

It is important to note that this effect is distinct from the parity asymmetry observed at $\ell <20$, which hints at quantum effects on super-horizon scales and could also provide an explanation for the anomalies in the low quadrupole and hemispherical power asymmetry (Gaztañaga and Sravan Kumar [63]).

In the standard model of inflation, super-horizon perturbations are explained as being "frozen" when they cross the horizon during the inflationary phase, only to re-enter at a later stage of the universe’s evolution. Although the BHU model solves the horizon problem through a similar causal mechanism, it introduces a key difference: a cutoff exists for scales larger than $\lambda >2R(k<\pi /R)$ [28], unlike the scale-invariant spectrum predicted by inflation. This distinction results in the anomalous absence of the largest structures in the CMB sky temperature (on top of the quantum parity asymmetry mentioned before), which potentially explains the observed deviations at these scales.

Furthermore, perturbations generated during the bounce phase of the BHU model could lead to notable consequences. For example (see Gaztañaga [28]):

-

A portion or even the entirety of Dark Matter (DM) might consist of primordial black holes or primordial neutron stars formed during the bounce.

-

The bounce mechanism could offer insights into the measured entropy of the Universe, particularly the observed baryon-to-photon ratio.

These implications underscore the need for further investigation, which forms the central motivation for the research presented in this paper.

5. Conclusions

This work is motivated by finding a scenario which has the potential to explain various features of the behavior of our Universe (e.g. the singularity, flatness, dark energy, and horizon problems), which cannot be solved easily within the standard model of cosmology. Following the BHU model [15,16], our Universe is inside a black hole and is currently in a state of expansion. Before this expansion, the universe was thought to have collapsed from a cold, tenuous cloud of matter. Once a critical density had been reached, this matter had to bounce back. As a first step to understanding the processes and necessary conditions of such a bouncing universe model, we investigated the collapse of a spherically symmetric matter distribution and how a potential bounce can occur for a simplistic setup.

We started by discussing the general relativistic equations that describe the collapse of a uniform spherical cloud. For a cloud with zero pressure and zero total energy, these equations are mathematically equivalent to a Newtonian description. This equivalence is a well-known fact in GR. A numerical simulation of one such cloud in Newtonian hydrodynamics with a mass of $8M_{\odot }$ and $P=0$ shows that the radius and density evolve approximately as predicted by the analytical GR solution. Moreover, this cloud remains approximately FLRW (uniform density) at all times. In our approach, we were able to simplify the model such that the formation of a BH, a genuine GR effect, can be inferred from a Newtonian simulation, significantly easing the numerical implementation. This simplification is valid only under the assumption of spherical symmetry and when vorticity and tidal fields are neglected (see Section 2.3 [43]). Proceeding with the Newtonian code and abandoning the approximation of zero pressure, we simulate the collapse of FLRW clouds of masses 5M_⊙, 20M_⊙ and 1000M_⊙ with a polytropic EoS. This is a simple EoS that creates a bounce and is motivated by its frequent use to describe neutron stars that result from a collapse of stars with masses from 8 to 40M_⊙. Since this is a qualitative study, we do not intend to solve this collapse in all detail, but rather, we want to study the basic features on the route to understand the collapse of much more massive configurations. Our scenario of a cold collapse is significantly simpler than a stellar collapse, where all kinds of nuclear reactions and neutrino-matter interaction play important roles, and a much simpler barotropic EoS is justified here. However, the exact equivalent GR solution for an EoS with non-zero pressure is not straightforward and is left to be done in future work.

For the cases with non-zero pressures, initially, the central densities increase in our simulations as in the pressureless solution. However, because of the equally increasing pressure, the matter bounces once it reaches a certain maximum density, and subsequently, the central density starts to decrease. Apart from the lightest model with the stiffest EoS, this happens when $95\%$ of the total mass is already inside its Schwarzschild radius. This means that in the GR framework, a BH has formed, and matter collapses and bounces inside its event horizon. For the masses we simulated here, this repulsion leading to the bounce could be caused by neutron or even quark degeneracy. However, for much larger masses comparable to the mass of our universe, which we finally aim at collapsing, none of the known physical effects could stop the collapse. In the BHU model, we thus assume that at some point, some supranuclear saturation densities will be reached, and an exotic form of quantum matter will instigate the bounce due to a quantum mechanism similar to Pauli’s exclusion principle. Crucially, this quantum mechanism violates the strong energy condition in classical GR and, together with $k>0$, sidesteps the singularity GR theorems proposed by Hawking and Penrose [1], thus presenting a novel solution to a pivotal issue in cosmological theory.

Figure 8 illustrates our key findings. When estimating the properties of the relativistic bounce with our Newtonian simulations with a polytropic EoS, we find conditions where $P=-\rho$. This usually occurs when a ground state of the matter is reached (see the Appendix B for an example of a scalar field). As stated before, such a ground state could be characterized by the Pauli Exclusion Principle for an unknown form of quantum matter.

We started our investigation from the equivalence of the pressureless collapse in GR and Newtonian descriptions and used a Newtonian simulation with pressure to estimate what one can expect in the full GR case. We have focused on distributions containing the order of magnitude of stellar size masses up to 1000M_⊙, where the pressure can come from the neutron degeneracy pressure. For larger masses, the particular mechanism that creates the required extreme pressures is highly speculative. However, from quantum principles, it seems reasonable that matter cannot be compressed arbitrarily. As a next step, we intend to increase the simulated masses and generalize our results by treating the entire collapse within GR. This includes, in particular, the requirement to drop the condition of homogeneous density and pressure in the GR case.