The Origin of Cosmic Structures Part 3 — Supermassive Black Holes and Galaxy Cluster Evolution

J. C. Botke

Nogales, Arizona, USA.

DOI: 10.4236/jhepgc.2022.82028                

Abstract

In Part 1 of this work, we showed that our new model of cosmology can account for the origin of all cosmic structures ranging in size from stars up to superclusters. In this model, at the time of nucleosynthesis, an imprint embedded in the vacuum regulated the creation of the protons (and electrons) that later made up the structures. Immediately after nucleosynthesis and for a considerable period afterward, the evolution was completely determined by the expansion of the universe. Gradually, however, gravitational influences became more important until finally, the expansion of the structures-to-be ceased at their zero velocity points. Stars, galaxies, and galaxy clusters all reached their zero velocity points more or less simultaneously at the usually accepted time of the beginning of galaxy formation. From that point onward, the evolution gravitation came to dominate the evolution although the expansion still exerted its influence. In this paper, we examine the subsequent cluster evolution in some detail. We establish the conditions required to prevent a free-fall collapse of the clusters and then show that galaxies with quasar-like active nuclei located within the cluster were the sources of the necessary radiation. We also show that the required galactic supermassive black holes were a consequence of the initial free-fall collapse of all galaxies.

Keywords

Galaxy Cluster Evolution, Supermassive Black Holes, Early Universe, Time-Varying Curvature, Nucleosynthesis

1. Introduction

In Part 1 of this series of papers [1], we presented a model detailing the initial evolution of all cosmic structures which covered the period between their initial definition at the time of nucleosynthesis through to the point at which the nascent structures reached their zero-velocity point (ZVP) defined to be the point at which they ceased to expand. It was at this point that gravitational interactions began to dominate over the expansion of the universe. We found that all structures reached their ZVPs more or less simultaneously at the time $t_G=3.2\times {10}^{16}\text{s}$ and that they did so with more or less their final masses. We presented some preliminary ideas about the subsequent evolution of galaxy clusters based on an equilibrium model solution and concluded that intense radiation heating must have occurred very soon after the ZVP was reached because otherwise, the clusters would have undergone free-fall collapse.

In this paper, we examine the post-ZVP evolution in more detail. We base our discussion on the equations of motion of a gas cloud and find that the solutions support our earlier contention. The determining factors include first, the initial density distribution set at the time of nucleosynthesis. Obviously, this encompasses the material particles but in addition, it includes a vacuum energy contribution which is the reality of the dark matter of conventional models. The second factor was the intensity and temporal profile of the radiation and the third factor was the simultaneity of the ZVPs which we found was of critical importance for the formation of stable clusters.

The two standard approaches for studying the dynamics of the gas are to treat it as an N-body problem based on Boltzmann’s equation or to treat it as a continuous fluid described by Euler’s equations [2] [3] [4] [5]. We have chosen to follow the Euler equation approach and will present some justifications in the first part of Section 3.

Having established criteria necessary for the stability of the clusters, we then show that galaxies with quasar-like active nuclei located within the cluster itself were the sources of the radiation and that the necessary supermassive black holes at the center of the galaxies developed very early during the initial free-fall collapse of the galaxies.

2. Starting Point

The proton gas making up the inter-cluster medium (ICM) came into existence during nucleosynthesis. Initially, this gas would have been in thermal equilibrium with the CMB and would have remained so up until the time of recombination. By then, the temperature had dropped to a value on the order of 4000 K which was cool enough to allow the protons and electrons to combine into neutral hydrogen. This resulted in their decoupling from the CMB and from that point on the kinetic energies of the particles decreased proportionally to $a{\left(t\right)}^{-2}$, where $a\left(t\right)$ is the scaling of the expanding universe, and by $t=t_G$, their energies would have dropped to an equivalent temperature of less than 1 K. This low temperature, incidentally, was an essential condition for the formation of the first stars. Our model asserts that there was no bulk motion of the particles at the time they were created so that, aside from their thermal motion, the particles were at rest. At the ZVP, with the effective temperature near absolute zero, the particles were very nearly motionless. Also, at the ZVP, the galaxy cluster to-be had reached its present-day size so the ICM particle density was then equal to its present-day value. In what follows, we will use the Virgo cluster as our example. Its particle density is ${\rho }_{avg}={10}^3{\text{m}}^{-3}$. The total mass, initial size, temperature, and velocity profile of the cluster are thus all known.

We now have two problems to solve. The first is to determine what combinations of the initial density and radiation result in a cluster that neither collapses nor evaporates. The second problem is to justify our contention that quasars-like active galaxy nuclei were the sources of the radiation.

3. Equations

We will begin by reviewing a few basic parameters such as mean free paths (MFP) that characterized the dynamics of the cluster. Initially, the gas was composed of neutral hydrogen (or H₂ molecules) so, to estimate their MFP, it is reasonable to assume a hard-sphere model with an effective diameter on the order of $d={10}^{-10}\text{m}$. The MFP is then given by

$\lambda =\frac{1}{\sqrt{2}\pi n_H{d}^2}$. (3-1)

With an ICM particle density of 10³ m⁻³, the MFP is $\lambda =2.3\times {10}^{16}\text{m}$ which is quite small compared to the radius of the cluster ( $R_C=6.8\times {10}^{22}\text{m}$ ). With a nominal temperature of 1 K, the corresponding kinetic particle velocity would have been 160 m∙s⁻¹ so the average time between collisions would have been 1.5 × 10¹⁴ s which is quite small compared with either $t_G$ or the free-fall time. As we will see when we examine the simulation results in Section 8, very soon after the ZVP, the temperature of the gas must have increased dramatically to prevent a collapse so the gas would have been rapidly ionized with a consequent change in the MFP. From [6], the MFP for an ionized gas is given by

$\lambda \approx 7.1\times {10}^{20}{\left(\frac{T}{{10}^8}\right)}^2\sqrt{\frac{{10}^3}{n_p}}\text{m}$. (3-2)

(The previous hard-sphere MFP was defined in terms of any collision between particles. This latter MFT, on the other hand, is defined as the distance that a charged particle must travel for its accumulated scatterings to change its direction by 90˚. This explains why it increases with temperature.) With a temperature of $T=5.7\times {10}^5\text{K}$, we have the same MFP as for neutral hydrogen and for a value of $T=5\times {10}^7\text{K}$, we have $\lambda =1.8\times {10}^{20}\text{m}$ so with temperatures high enough to prevent collapse, the MFP would have been on the order of galaxy dimensions. Although it doesn’t concern us in this paper, inside the galaxies where the interstellar gas density is on the order of 10⁶ m⁻³, the MFP would be reduced by a factor of about 30 but it would still have been on the order of the radius of all but the largest galaxies.

Using these results, we can compute the Knudsen number,

$K_n=\frac{\lambda }{R_C}\le 2.6\times {10}^{-3}$. (3-3)

The rule of thumb is that a continuum model such as Euler’s equations is appropriate whenever the Knudsen number for the flow is less than 0.1 which is the case here.

Starting with this same result, it is shown in [2] that the time scale for equilibration from a non-equilibrium starting point is on the order of 1.9 × 10¹⁶ s which is smaller than but still comparable to both $t_G$ and free-fall time. At temperatures on the order of 5 × 10⁷ K, the electrons would have been nearly relativistic.

We will next consider the escape velocity. The relevant equation is

$\frac{G{M}_c{m}_p}{R_c}=\frac{1}{2}{m}_p{v}_e^2=\frac{3}{2}{k}_b{T}_e$. (3-4)

After substituting values for the Virgo cluster, we find, $T_e=8.0\times {10}^7\text{K}$ which sets an upper limit on the possible temperatures in the outer regions of the ICM. For temperatures higher than this value, the cluster would begin to evaporate. For temperatures lower than this value, the cluster could still expand (positive collective motion) but it would not have evaporated. For comparison, we note that Jean’s critical temperature is given by

$\frac{3}{5}\frac{G{M}_c^2}{R_c}=N_p\frac{3}{2}{k}_b{T}_c$ (3-5)

where $N_p$ is the number of protons in the ICM. But $M_c/N_p=m_p$ so, aside from the factor of 3/5, this is the same result.

We can add one more Jean’s metric. Jean’s length is the product of the sound speed and the free-fall time. We are getting ahead of ourselves a little but we will see shortly that $t_{ff}\approx 1.1t_s$ where $t_s=6.83\times {10}^{16}\text{s}$ is a scale time introduced below. The sound velocity scale value, which is also introduced below, is $v_s=1.3\times {10}^6\text{m}\cdot {\text{s}}^{-1}$ so Jean’s length is $L=t_{ff}{v}_s=9\times {10}^{22}\text{m}$. Jean’s rule is that a system will collapse if its size is greater than L so, according to this rule, a stable cluster is possible sinceL is larger than $R_C$ but the margin isn’t large.

We will now formulate the equations of motion of the system. One issue that we must address is that the proton gas makes up only a portion of the total mass of the cluster and the vacuum energy, which makes up the remainder, must be taken into account. The procedure we will follow will be to develop the equations pretending that only the gas exists and afterward, to point out the adjustments needed to include the vacuum energy.

Clusters are certainly not spherically symmetric and they don’t have a precisely defined center as does, for example, a star. Nevertheless, to finish with a set of equations that can be solved on a laptop, we will need to assume spherical symmetry. This simplifies the work because it reduces the problem to one with a single spatial dimension but at the same time, it introduces two significant, non-physical, issues that complicate the analysis. The first problem is that such models tend to develop troublesome non-physical singularities at the origin. The second issue that the model exhibits even more troublesome, non-physical oscillatory behavior. This will be discussed further after we have developed the equations.

We now imagine the cluster to consist of a series of concentric spherical shells of radius $r\left(t\right)$, density $\rho \left(t,r\right)$, pressure $p\left(t,r\right)$, and temperature $T\left(t,r\right)$. We will make the additional assumption that we can adopt the Lagrangian viewpoint in which one follows the motion of a specific chunk of the gas as it moves about. We emphasize that this is an assumption because later when we consider the oscillations, we will find that this assumption cannot be satisfied because the mixing that must occur in the real ICM makes it impossible to identify a specific chunk of the gas.

Conservation of momentum gives the usual Newton’s law equation,

$\frac{{\partial }^{2}r\left(t\right)}{\partial t^2}=-\frac{Gm\left(r\right)}{r{\left(t\right)}^2}-\frac{1}{\rho \left(t,r\right)}\frac{\partial p\left(t,r\right)}{\partial r}+\frac{r\left(t\right)}{t^2}\left(-{\gamma }_{\ast }+{\left({\gamma }_{\ast }+c_1\frac{t}{t_0}\right)}^2\right)$ (3-6)

The last term, which is explained in [1], is the acceleration resulting from the expansion of the universe. This term is not considered in most simulations but the fact is that at the time of ZVP, its value, which is negative, is equal to that of the gravitational term and it can have a noticeable effect on the early phases of the evolution. The quantity $m\left(r\right)$ is the total mass inside the shell with radius $r\left(t\right)$,

$m\left(r\right)=4\pi {\displaystyle {\int }_0^{r\left(t\right)}\text{d}{r}^{\prime }}{r^{\prime }}^2\rho \left(t,r^{\prime }\right)$. (3-7)

With no mixing between the shells, this quantity will have a constant value for each initial choice of radius.

The next equation expresses the conservation of the internal energy of each shell. This can vary as a result of a change in its temperature, of work done by the chunk of gas, and by the net addition or loss of energy due to radiation. The internal energy per unit mass of an ideal gas is given by $U=\left(3/2\right)p\left(t,r\right)/\rho \left(t,r\right)$ and the work done by the pressure is $p\left(r,t\right)\delta V=p\left(r,t\right)\delta \left(1/\rho \left(t,r\right)\right)$, again per unit mass so the rate of change of the internal energy is

$\frac{3}{2}\frac{\partial }{\partial t}\left(\frac{p\left(t,r\right)}{\rho \left(t,r\right)}\right)+p\left(t,r\right)\frac{\partial }{\partial t}\left(\frac{1}{\rho \left(t,r\right)}\right)=q_{ext}\left(t,r\right)-q_{rad}\left(t,r\right)$ (3-8)

where $q_{ext}\left(t,r\right)$ is the amount of external energy absorbed by the shell and $q_{rad}\left(t,r\right)$ is the amount of energy radiated away by the shell, both per unit mass.

The LHS can be rearranged into a more convenient form. First, we expand the derivatives to find

$\frac{3}{2}\frac{\partial }{\partial t}\left(\frac{p\left(t,r\right)}{\rho \left(t,r\right)}\right)+p\left(t,r\right)\frac{\partial }{\partial t}\left(\frac{1}{\rho \left(t,r\right)}\right)=\frac{3}{2\rho \left(t,r\right)}\left(\frac{\partial p\left(t,r\right)}{\partial t}-\frac{5}{3}\frac{p\left(t,r\right)}{\rho \left(t,r\right)}\frac{\partial \rho \left(t,r\right)}{\partial t}\right)$ (3-9)

We now define the entropy function,

$\psi \left(t,r\right)\equiv \frac{p\left(t,r\right)}{\rho {\left(t,r\right)}^{5/3}}$ (3-10)

and note that

$\frac{\partial \psi \left(t,r\right)}{\partial t}=\frac{1}{\rho {\left(t,r\right)}^{5/3}}\left(\frac{\partial p\left(t,r\right)}{\partial t}-\frac{5}{3}\frac{p\left(t,r\right)}{\rho \left(t,r\right)}\frac{\partial \rho \left(t,r\right)}{\partial t}\right)$. (3-11)

After substituting, (3-9) becomes

$\frac{\partial \psi \left(t,r\right)}{\partial t}=\frac{2}{3\rho {\left(t,r\right)}^{2/3}}q\left(t,r\right)$ (3-12)

where $q\left(t,r\right)\equiv q_{ext}\left(t,r\right)-q_{rad}\left(t,r\right)$. Under adiabatic conditions, $q\left(t,r\right)=0$ so $\psi$ is then constant. Eventually, we will remove the pressure dependence from (3-5) in favor of $\psi$ and the density but first, we will work out the rate of change of the density.

In the Lagrangian viewpoint, the amount of matter in each shell is fixed (again by assumption) which means that the size and position of the sample can change but, because no matter flows in or out, the density can only vary because of changes in the volume of the shell. For a unit mass, we have

$\frac{\partial \rho \left(t,r\right)}{\partial t}=\frac{\partial }{\partial t}\left(\frac{1}{V}\right)=-\frac{1}{V^2}\frac{\partial V}{\partial t}$. (3-13)

The volume, V, is that of a shell with a thickness of $\delta$ so

$V=\frac{4\pi }{3}\left(r_2^3-r_1^3\right)$ (3-14)

where $r_1$ and $r_2$ are the inner and outer radii of the shell. The rate of change is thus

$\frac{\partial V}{\partial t}=4\pi \left(r_2^2\frac{\partial r_2}{\partial t}-r_1^2\frac{\partial r_1}{\partial t}\right)=4\pi \left(r_2^2{v}_2-r_1^2{v}_1\right)$. (3-15)

We now write $r_2=r_1+\delta$, $v_2=v_1+\left(\partial v/\partial r\right)\delta$, and expand to first order. The result is

$\frac{\partial V}{\partial t}=\left(4\pi r^2\delta \right)\left(\frac{2v}{r}+\frac{\partial v}{\partial r}\right)$. (3-16)

But $V=4\pi r^2\delta =1/\rho \left(t,r\right)$ so we have

$\frac{\partial \rho \left(t,r\right)}{\partial t}=-\left(\rho \left(t,r\right)\frac{\partial v\left(t,r\right)}{\partial r}+2\frac{\rho \left(t,r\right)v\left(t,r\right)}{r\left(t\right)}\right)$. (3-17)

This is similar to the familiar conservation equation in the Euler formulation but it lacks the density gradient term.

The final step is to assume the ideal gas law equation of state,

$\frac{p}{\rho }=\psi {\rho }^{2/3}=\frac{k_{B}T}{\mu m_p}$ (3-18)

where $\mu =1$ for neutral hydrogen and $\mu =1/2$ ionized protons and electrons.

We now have the necessary number of equations. We will next perform several steps to further simplify these equations. The first involves the radial coordinate. Each shell is characterized by its radius but the radius varies with time and is not unique to any particular shell because different shells can have the same radius at different times. The quantity $m\left(r\right)$, on the other hand, is unique for each shell (again assuming no mixing) so it simplifies matters to choose $m\left(r\right)$ to be the independent spatial coordinate. From this perspective, the coordinate r becomes a dependent variable on the same footing as $\psi$ and the density. From (3-7), we have $\partial m/\partial r=4\pi r^2\rho \left(t,r\right)$ so we can replace the r derivatives with

$\frac{\partial }{\partial r}=\frac{\partial m}{\partial r}\frac{\partial }{\partial m}=4\pi r{\left(t,m\right)}^2\rho \left(t,m\right)\frac{\partial }{\partial m}$. (3-19)

The next change is to replace the various parameters with dimensionless equivalents which are denoted by an overbar. The definitions are given in Table 1. The numerical values shown were calculated using Virgo cluster values. The last entry is the sound speed, $v_s=\sqrt{\gamma p/\rho }$.

The total mass, $M_c$, includes both the material masses and the effective mass of the vacuum energy density (the dark matter of the conventional model).

Table 1. Definitions of dimensionless parameters.

Next, because we will be starting our simulations at $t=t_G$, it will be convenient to define a shifted time $\bar{t}$ where $t=t_s\left(0.469+\bar{t}\right)$. With this definition, $\bar{t}$ ranges from 0 to 5.91 as t ranges from $t_G$ to $t_0$.

The final adjustment involves the density. There are two issues involved. First, in some cases, the equations develop a singularity at the origin which particularly affects the density. The second problem is that we will need to specify the value of the density variable at the origin as a boundary condition and, in general, the value of the density is either unknown or is infinite. We can solve both problems by introducing a new dimensionless density parameter defined by

$\bar{\rho }\equiv \frac{\tilde{\rho }}{{\bar{r}}^2}$. (3-20)

While we don’t know a priori the value of the density at the origin, any reasonable solution will have $\tilde{\rho }\left(t,0\right)=0$.

Putting this all together, the final set of equations become

$\frac{\partial \bar{r}}{\partial \bar{t}}=\bar{v}$ (3-21a)

$\frac{\partial \bar{v}}{\partial \bar{t}}=-\frac{\bar{m}}{{\bar{r}}^2}-4\pi \frac{{\tilde{\rho }}^{5/3}}{{\bar{r}}^{4/3}}\left(\frac{\partial \bar{\psi }}{\partial \bar{m}}+\frac{5}{3}\frac{\bar{\psi }}{\tilde{\rho }}\frac{\partial \tilde{\rho }}{\partial \bar{m}}\right)+\frac{10}{3}\frac{{\tilde{\rho }}^{2/3}\bar{\psi }}{{\bar{r}}^{7/3}}+\frac{\bar{r}}{{\bar{t}}^2}\left(-{\gamma }_{\ast }+{\left({\gamma }_{\ast }+c_1\frac{\bar{t}}{{\bar{t}}_0}\right)}^2\right)$ (3-21b)

$\text{}\frac{\partial \tilde{\rho }}{\partial \bar{t}}=-4\pi {\tilde{\rho }}^2\frac{\partial \bar{v}}{\partial \bar{m}}$ (3-21c)

$\frac{\partial \bar{\psi }}{\partial \bar{t}}=\frac{2}{3}{\left(\frac{{\bar{r}}^2}{\tilde{\rho }}\right)}^{2/3}\bar{q}$. (3-21d)

Notice that the redefinition of the density eliminates the second term in the density Equation (3-17). We also find that there is no cluster-specific dimensionless parameter remaining in the equations so the results will apply to any cluster.

The pressure, temperature, and sound speed are given by

$\bar{p}={\left(\frac{\tilde{\rho }}{{\bar{r}}^2}\right)}^{5/3}\stackrel{\to }{\psi }$ (3-22a)

$\bar{T}=\frac{\bar{p}{\bar{r}}^2}{\tilde{\rho }}={\left(\frac{\tilde{\rho }}{{\bar{r}}^2}\right)}^{2/3}\bar{\psi }$ (3-22b)

${\bar{v}}_s=\sqrt{\bar{T}}$. (3-22c)

Finally, (3-7) becomes

$\bar{m}=4\pi {\displaystyle {\int }_0^{\bar{r}}\text{d}{\bar{r}}^{\prime }}\tilde{\rho }\left(\bar{t},{\bar{r}}^{\prime }\right)$. (3-23)

At the outer boundary, we have the constraint that

$4\pi {\displaystyle {\int }_0^{\bar{r}}\text{d}{\bar{r}}^{\prime }}\tilde{\rho }\left(\bar{t},{\bar{r}}^{\prime }\right)=1$ (3-24)

which provides a useful check on the accuracy of each solution.

We will now consider the energy density of the vacuum. To do this properly, we would need to reformulate the entire problem in terms of the full set of Einstein’s equations with the vacuum energy included. This, however, is out of reach since we are restricted to what can be done on a laptop. Looking at the equations, in (3-6), the density in the pressure gradient term is the total density of the cluster and this is also true in both (3-7) and (3-17). The density in the internal energy and ideal gas law equations, (3-8) and (3-18), on the other hand, is the density of just the gas. The pressure is also a consequence of the gas dynamics. What we are lacking is a vacuum equivalent of those two equations. That being the case, we will just assume that the gas density is everywhere proportional to the total density, or ${\rho }_{gas}=\left(M_{gas}/M_c\right){\rho }_{total}$. Substituting into (3-18) gives

$\frac{p}{\rho \left(M_{gas}/M_c\right)}=\frac{k_{B}T}{\mu m_p}$ (3-25)

If we redefine a “total” pressure by $p_{total}\equiv p_{gas}/\left(M_{gas}/M_{total}\right)$, this becomes

$\frac{p_{total}}{{\rho }_{total}}=\frac{p_{gas}}{{\rho }_{gas}}=\frac{k_{B}T}{\mu m_p}$ (3-26)

Equation (3-8) only depends on the same ratio of the pressure and density so it retains its original form but with the density and pressure replaced by the total values

$\frac{3}{2}\frac{\partial }{\partial t}\left(\frac{p_{total}\left(t,r\right)}{{\rho }_{total}\left(t,r\right)}\right)+p_{total}\left(t,r\right)\frac{\partial }{\partial t}\left(\frac{1}{{\rho }_{total}\left(t,r\right)}\right)=q_{ext}\left(t,r\right)-q_{rad}\left(t,r\right)$. (3-27)

The final step is to replace the gas pressure in (3-6) with the total pressure.

$\frac{{\partial }^{2}r\left(t\right)}{\partial t^2}=-\frac{Gm\left(r\right)}{r{\left(t\right)}^2}-\frac{1}{{\rho }_{total}\left(t,r\right)}\left(\frac{M_{gas}}{M_c}\right)\frac{\partial p_{total}\left(t,r\right)}{\partial r}+\frac{r\left(t\right)}{t^2}\left(-{\gamma }_{\ast }+{\left({\gamma }_{\ast }+c_1\frac{t}{t_0}\right)}^2\right)$. (3-28)

After redefining $\psi$ by

$\psi \left(t,r\right)\equiv \frac{p_{total}\left(t,r\right)}{{\rho }_{total}{\left(t,r\right)}^{5/3}}$ (3-29)

all the equations that follow (3-8) remain unchanged but with the understanding that the density and pressure refer to their total values rather than to just the gas values.

The net result is that the only change is the presence of the mass ratio in (3-28).

4. Dimensionless Initial and Boundary Conditions

Certain conditions are fixed. Certainly $r\left(\bar{t},0\right)=0$ and because we are starting at the ZVP, $\bar{v}\left(0,\bar{m}\right)=0$. We also know that the origin doesn’t move so $\bar{v}\left(\bar{t},0\right)=0$ and we are assuming that for any reasonable solution, $\tilde{\rho }\left(\bar{t},0\right)=0$. The initial temperature is $T\left(0,m\right)=1\text{K}$ and with the density known, $\psi \left(0,\bar{m}\right)$ is given by (3-22b). At the outer boundary, the density and everything else should vanish but trying to include the surface would require an extremely dense mesh which, in our case, is not practical. Instead, we terminate the calculation slightly inside the surface to avoid having to deal with the large gradients. Because the equations are first order in the spatial derivatives, only a single boundary condition is necessary and we have already specified values at the origin.

We next need to specify the gas/total mass ratio. For the Virgo cluster, the total mass of the gas is thought to be on the order of $3\times {10}^{14}{M}_{\odot }=6\times {10}^{44}\text{kg}$. The total mass is less certain but is estimated to be around $1.2\times {10}^{15}{M}_{\odot }=2.4\times {10}^{45}\text{kg}$ which would give a mass ratio of $M_c/M_g=4$. In Part 1, on the other hand, we argued that the total mass was likely to be about 1/2 that value which would give a mass ratio of 2 which is the value we will use in the simulations.

Remaining to be specified are $\tilde{\rho }\left(0,\bar{m}\right)$ and $q\left(\bar{t},\bar{m}\right)$. If we assume that the radiation came from sources distributed more or less uniformly, it becomes reasonable to assume that q varies with time but not position so $q\left(\bar{t},\bar{m}\right)\equiv q\left(\bar{t}\right)$. (Recall that q has the units of energy per unit mass rather than energy per unit volume.) In different words, we are assuming that the amount of radiation absorbed by a chunk of gas is dependent on the mass of the chunk but not its position. It is by varying these two quantities that we can narrow down the possible choices that lead to a solution to the cluster evolution problem.

5. Oscillations

Referring to the Equation (3-21), there is a problem buried in these equations that leads to an endless amount of trouble. To see this, we first take the time derivative of (3-21b) and ignore everything except the spatial derivative of the density term. We then take the spatial derivative of (3-21c) and substitute it into the first equation again ignoring everything except the term coming from the spatial derivative of the velocity. The result is

$\frac{{\partial }^2\bar{v}}{\partial {\bar{t}}^2}={\left(4\pi \right)}^2\frac{{\tilde{\rho }}^{8/3}}{{\bar{r}}^{4/3}}\frac{5}{3}\bar{\psi }\frac{M_{gas}}{M_c}\frac{{\partial }^2\bar{v}}{\partial {\bar{m}}^2}$. (5-1)

This is a wave equation with a velocity of

$u=4\pi \frac{{\tilde{\rho }}^{4/3}}{{\bar{r}}^{2/3}}\sqrt{\frac{5}{3}\bar{\psi }\frac{M_{gas}}{M_c}}$. (5-2)

After substituting for $\bar{\psi }$ and using $\tilde{\rho }={\bar{r}}^2\bar{\rho }$, this simplifies to

$u=\left(4\pi \tilde{\rho }\right)\sqrt{\gamma \frac{\bar{p}}{\bar{\rho }}}\sqrt{\frac{M_{gas}}{M_c}}=\frac{\text{d}\bar{m}}{\text{d}\bar{r}}\left(\sqrt{\frac{M_{gas}}{M_c}}{\bar{v}}_s\right)$. (5-3)

In Table 1, we defined the sound speed scaling in terms of the total mass. The sound speed that appears here is

$v_{s,gas}=\sqrt{\frac{M_{gas}}{M_c}}{\bar{v}}_s=\sqrt{\gamma \frac{G{M}_{gas}}{R_c}}{\bar{v}}_s$. (5-4)

We can go one step further by replacing the $\bar{m}$ derivative with the $\bar{r}$ derivative. Again, ignoring everything but the second derivative, we find

$\frac{{\partial }^2\bar{v}}{\partial {\bar{t}}^2}={\bar{v}}_{s,gas}^2\frac{{\partial }^2\bar{v}}{\partial {\bar{r}}^2}$. (5-5)

Going back to the first step, if we instead take the time derivative of (3-21c) and the spatial derivative of (3-21b), we end up with

$\frac{{\partial }^2\tilde{\rho }}{\partial {\bar{t}}^2}={\bar{v}}_{s,gas}^2\frac{{\partial }^2\tilde{\rho }}{\partial {\bar{r}}^2}$. (5-6)

If we consider solutions at or very close to equilibrium at all times, the time derivatives are small and the oscillations don’t become a problem. Generally, however, we will be considering situations that, at least initially, are far from equilibrium and in those cases, the time derivatives become large and oscillations are inevitable. The immediate numerical problem is that these oscillations destabilize the equations and the solution fails. This is particularly a problem when using $\bar{m}$ as the spatial variable because the corresponding velocity u vanishes at the origin so oscillations with wavelengths approaching zero are possible.

What we will now show is that, even though the potential of oscillations is contained within the equations, such oscillations are suppressed by the dynamics of the ICM. There are a couple of ways to see this. Consider a region of higher than average density with a radius of $R_B={10}^{20}\text{m}$ and a mass excess of 20%. This radius was chosen because it approximates the proton MFP. Jean’s critical temperature for this bulge is given by

$T_{crit}=\frac{2}{5}\frac{G{M}_B{m}_p}{k_b{R}_B}\approx 50\text{K}$ (5-7)

which is extremely small compared to the actual temperature of the gas. This indicates that the self-gravitation is far too small to bind the bulge so it would rapidly dissipate.

Another approach is to consider the motion of a proton located at the center of the bulge. Because we have chosen a bulge radius on the order of the MFP, the proton will only be affected by the gravitation of the bulge. The equation of motion in scaled coordinates is

$\frac{{\partial }^2\bar{r}}{\partial {\bar{t}}^2}=-\frac{\bar{m}}{{\bar{r}}^2}$. (5-8)

If we assume a constant density, $\bar{m}=\left(4\pi /3\right){\bar{r}}^3{\rho }_B$ and at the outer boundary, $\bar{m}=\bar{r}=1$ so ${\rho }_B=3/4\pi$. The equation then becomes

$\frac{{\partial }^2\bar{r}}{\partial {\bar{t}}^2}+\bar{r}=0$. (5-9)

The initial (thermal) velocity of the proton is $v_p=\sqrt{3k_{b}T/m_p}=158\sqrt{T}$ which results in a scaled value of ${\bar{v}}_p=0.16\sqrt{T}$. The solution is $\bar{r}\left(\bar{t}\right)=0.16\sqrt{T}\sin \left(\bar{t}\right)$ so the scaled time for the proton to reach the edge of the bulge is

$\bar{t}={\sin }^{-1}\left(\frac{1}{0.16\sqrt{T}}\right)$. (5-10)

With a temperature of $T={10}^7\text{K}$, the actual time for the proton to exit the bulge is $t=2.0\times {10}^{14}\text{s}$ which is very small compared to any time characteristic of the cluster evolution so again we see that the bulge will dissipate.

Even though we have shown that oscillations would be suppressed in the actual cluster, we are still left with the problem of numerical oscillations in the solutions. One possibility we considered was to add a contribution to (3-21c) that takes into account the bulge dissipation just discussed. If we consider an imaginary surface to the left of a volume of gas with a particle density of $n_{+}{\text{m}}^{-3}$, the rate at which particles will cross a circle of radius r on the surface is $\pi r^2{n}_{+}v$ where v is the average particle velocity (refer to any source that discusses the kinetic theory of pressure). Similarly, particles from the other side of the surface will cross the surface at a rate of $\pi r^2{n}_{0}v$ in the opposite direction. If we now consider a region of width δ between two such surfaces, the net rate of change of its density would be

$\frac{\partial \rho }{\partial t}=\frac{\delta }{8}\sqrt{\frac{3k_{b}T}{m_p}}\frac{{\partial }^2\rho }{\partial r^2}$. (5-11)

We didn’t include this term, however, because first, it violates our assumption of no mixing, and second, and even worse, it becomes very singular at the origin when the density is replaced by $\tilde{\rho }$. This suggests that a better approach to the entire problem would be to develop a 3-dimensional model in Cartesian coordinates based on the Euler viewpoint which does allow for mixing. Such a model would, however, require far more computer capacity than we have available.

To summarize the situation, we have shown that oscillations are inevitable given the non-linear hyperbolic nature of Euler’s equations. We then showed that the problem is purely a numerical issue because such oscillations cannot occur in actual clusters.

6. Solution Method

The numerical solution of the equations was accomplished by using the “method of lines.” In this case, the spatial dimension is divided up into a sequence of shells. To reduce the oscillation problem, central differences were used everywhere to approximate the spatial derivatives. At the ends, where forward or backward differences would have been required, we instead assumed that the end derivatives were the same as the derivatives evaluated at the first interior point, i.e. $\frac{\partial \tilde{\rho }\left(t,0\right)}{\partial m}\approx \frac{\partial \tilde{\rho }\left(t,m_1\right)}{\partial m}$, etc. When the finite difference approximations are substituted into (3-21), the result is a set of coupled time-domain ordinary differential equations that we solve using the adaptive predictor-corrector ODE solver called “Lsoda.” This is available on the internet in both the Fortran and “C” languages. For our use, we ported the solver to vb.net. For the most part, lsoda was treated as a black box. As the solution proceeds, lsoda repeatedly calls a user-supplied “GetRates()” method with a vector of the predicted values of the unknowns. The time derivatives (LHS of (3-21)) are then computed and returned to lsoda.

As discussed in the previous section, this scheme by itself doesn’t work because of the oscillations. Instead of modifying the equations, we tried to deal with the oscillations through the use of numerical filtering but the result was only moderately successful. Several filtering schemes were tried with the best results coming from cubic spline filtering which we applied to first, the predicted values, then to the derivatives, and finally to the rates. There are numerous references to this method on the internet. The particular form we used is discussed in [7]. The basic idea of the filter is to find a smooth function $g\left(x\right)$ that minimizes the following cost function,

$\frac{\partial \rho }{\partial t}=\frac{\delta }{8}\sqrt{\frac{3k_{b}T}{m_p}}\frac{{\partial }^2\rho }{\partial r^2}p{\displaystyle \underset{i=0}{\overset{N}{\sum }}{\left(g\left(x_i\right)-f_i\right)}^2}+{\displaystyle {\int }_{x_0}^{x_N}{g}^{″}\left(x\right)\text{d}x}$. (6-1)

The $f_i$ are the noisy values and p is a parameter that regulates the degree of filtering. When p is large, $g\left(x\right)$ will duplicate the noisy values and when it is small, $g\left(x\right)$ will correspond to a least-squares fit to the noisy values. The solution to the minimization problem is given by

$g=\left(I-\frac{1}{p}{C}^{\text{T}}{\left(A+B\right)}^{-1}C\right)f$. (6-2)

$A$ and $C$ are banded matrices that are dependent only on the spatial step sizes. The matrix $B=\frac{1}{p}C{C}^{\text{T}}$ is also banded but even though $A$ and $B$ are both banded, the inverse of their sum is not banded. Nevertheless, because the values of the elements of the inverse far from the diagonal are much smaller than those near the diagonal, we can turn the inverse into a banded matrix by setting to zero all the elements that are below a specified cutoff. Physically, the smallness of the far-off-diagonal elements is just a reflection of the fact that the averaging at any spatial location is primarily influenced by near neighbors rather than by distant points. An important point is that, because the various matrices are dependent only on the spatial step sizes, the filter matrix can be pre-computed. The application of the filter then only requires a single multiplication between a banded matrix and a vector.

The principal difficulty with this method is that, while it works well in the interior of the spatial range, it is not optimal near the ends which is where the biggest problem lies. Another attempt was made using a Fourier-Butterworth low-pass filter scheme. This wasn’t successful for much the same reason because it introduces Gibbs-type oscillations at the ends of the spatial range.

It might seem that we could avoid this problem by starting the solution in the interior or at the outer boundary but that isn’t possible because it is only at the origin that we know the values of the variables. In many cases the oscillations initiate close to the origin where the velocity (5-3) approaches zero and given that the origin must be included, we found that jumping out of the origin with a few large steps is a scheme that works somewhat better than using filters. To avoid a sudden transition at the point at which the large steps joined onto the smaller steps which covered the remainder of the range, we divided the jump span into several steps with graduated sizes, for example, $10\delta ,8\delta ,6\delta ,\cdots$ with the overall jump typically being on the order of $\Delta \bar{m}=0.05$. This scheme doesn’t suppress the oscillations entirely but it does postpone the eventual crash.

7. Radiation Model

It will become clear that no matter what initial density profile is assumed; the cluster will collapse without a significant input of radiation. To model the source, we developed code that generates a smooth profile based on three adjustable parameters. A sample profile is shown in Figure 1. As noted earlier, we assume that the radiation varies with time but not location.

The peak value of this scaled profile is always unity with the actual peak value given by the value in Table 1, namely 1.4 × 10⁻⁵ j∙kg⁻¹∙s⁻¹. To allow for an arbitrary magnitude, a multiplier is included in the model.

8. Simulations

Our goal is to determine what restrictions on the initial density profile and radiation are required to prevent the cluster from collapsing. Because of the oscillations, our results will be limited but we will be able to establish some general guidelines. We will start by displaying the free fall solution since this sets the time frame within which the cluster must be stabilized. From Figure 2, we see that without radiation, the cluster would have collapsed by a time of $\bar{t}\approx 1.1$ so to prevent this from happening, the radiation must have been in existence before $\bar{t}\le 0.5$ or $\left(t-t_G\right)\le 3.4\times {10}^{16}\text{s}$.

To better understand this point, recall that at the time of nucleosynthesis, the formation of the cluster was directed by a vacuum imprint. At that time, the outer boundary of the cluster had a velocity relative to the center given by the expansion rate of the universe. As time progressed, the expansion deceleration acted to slow the expansion until eventually, it ceased at the ZVP as a result of the combined action of cluster mass gravitation and the expansion deceleration. At a later time, $\bar{t}=2.22$ ( $t=1.83\times {10}^{17}\text{s}$ ), the expansion acceleration did

Figure 1. Adjustable source profile.

become positive after which, the acceleration acted to increase the size of the cluster. Its magnitude, however, varies as ${\bar{t}}^{-2}$ so the effect diminished with time.

We will now consider several initial density and source profiles which will allow us to narrow down the possible range of choices. First, we will consider a uniform density profile. From (3-7), we find that $\bar{\rho }=3/\left(4\pi \right)$ and $\bar{r}={\bar{m}}^{1/3}$. Figure 3 shows the variation of the outer dimension of the cluster with time. In this case, the peak radiation multiplier was 3 and the source profile was the one shown in Figure 1.

We see that the cluster undergoes collapse and we found that increasing the peak radiation multiplier does not change that result. In the next few figures, we show the calculated parameters. Figures 4-7 show $\bar{r}\left(\bar{m}\right)$, $\bar{v}\left(\bar{m}\right)$, $\tilde{\rho }\left(\bar{m}\right)$, and

Figure 2. Free-fall with and without the expansion acceleration.

Figure 3. Outer dimension of the cluster with uniform density and a peak radiation multiplier of 3.

Figure 4. The radial coordinate profiles, $\bar{r}\left(\bar{m}\right)$, for the solution with a uniform initial density. Notice the beginnings of the oscillations in the most recent line (the lowest one in the figure).

$\bar{\rho }\left(\bar{m}\right)$ respectively. In these figures, each curve corresponds to a single value of $\bar{t}$. The initialization is shown in dark red with the successive curves shown in a progressive lighter shade of red. The output time step was $\Delta \bar{t}=0.04$ or $\Delta t=2.73\times {10}^{16}\text{s}$ in all cases. Figure 4 shows the radial coordinate.

The kinks near the origin are a result of the origin jump discussed earlier. Figure 5 shows the velocity.

The next, Figure 6, shows the density, $\tilde{\rho }\left(\bar{m}\right)$.

Finally, Figure 7 shows the actual density.

In all these figures, one can see the solution failing as the oscillations increase in magnitude.

From these results, we find that the initial density could not have been uniform which makes sense because with uniform density, and, as always, a uniform initial temperature profile, the pressure from (3-18) will also be uniform which means that there was no pressure gradient in (3-6) to oppose the gravitational

Figure 5. The velocities $\bar{v}\left(\bar{m}\right)$ for the same solution. The oscillations are much more pronounced.

Figure 6. The densities $\tilde{\rho }\left(\bar{m}\right)$ for the same solution.

Figure 7. The densities $\bar{\rho }\left(\bar{m}\right)$ for the same solution. The density is not well-defined at the origin because of the 0/0 division in (3-20).

acceleration.

Having established that a uniform (zero gradient) profile doesn’t work, we considered power-law profiles of the form $\rho \left(r\right)\propto 1/r^{\sigma }$ which generate large gradients. To satisfy our constraint that $\tilde{\rho }\left(t,0\right)=0$ and also because the pressure gradient must be negative, $\sigma$ must lie in the range, $0<\sigma <2$. What we found was that, except for a very narrow range of values close to $\sigma =0.9$, the cluster underwent collapse. The exception is shown in Figure 8. The peak radiation multiplier was 5.

In this case, the cluster undergoes an expansion instead of a collapse. However, because of the very limited range $\sigma$ values that result in an expansion, it is very unlikely that the initial density profile was anything like a power-law profile.

In the next series of figures, we show the results for a straight-line profile with a slope equal to $\partial \bar{\rho }/\partial \bar{r}=-0.58$ that lies between the limits just discussed. The profile used was straight-line but a modestly curved profile would serve just as well. In Figure 9, we show the initialization curves.

We found that with a peak radiation multiplier of 4 or more, an expansion occurs. Figure 10 shows the outer dimension of the cluster calculated with a peak radiation multiplier of 5.

Figure 11 shows the radial coordinate.

In Figure 12, we show the velocity.

Figure 8. Outer dimension of the cluster with the density profile, $\rho \left(r\right)\propto 1/r^{0.9}$ and a peak radiation multiplier of 5.

Figure 9. Straight-line density profile. The $\tilde{\rho }\left(\bar{r}\right)$ and $\bar{\psi }\left(\bar{r}\right)$ curves are also shown.

The density is shown in Figure 13.

Finally, we show the temperature profiles in Figure 14. Notice that the temperature rises very rapidly. It reached the bottom of the red band within 1.1 × 10¹⁶ s. Notice also that its maximum value is slightly greater than one ( $T\approx 6\times {10}^7\text{K}$ )

Figure 10. Solution with the density profile of Figure 9 and a peak radiation multiplier of 5.

Figure 11. The radial coordinate profiles, $\bar{r}\left(\bar{m}\right)$ for the negative slope solution.

Figure 12. The velocities $\bar{v}\left(\bar{m}\right)$ for the same solution.

Figure 13. The densities $\tilde{\rho }\left(\bar{m}\right)$ for the same solution.

and that it shows a slight increase towards the center of the cluster.

The previous simulation was run without the expansion acceleration. Figure 15 shows the result when it is included.

We see that, while the difference is small, the expansion is slightly suppressed compared to the result without the expansion acceleration.

We have found that a density profile with a modest negative gradient will result in an expansion of the cluster but with the radiation profile of Figure 1, we seem to have overdone it because the expansion shows no sign of slowing within the time range shown. The issue is that the radiation profile of Figure 1 is too broad.

In the next example, we confined the radiation to the much narrower range shown in Figure 16. In this case, the width of the peak is on the order of 4.1 × 10¹⁶ s.

Figure 14. The temperature profile curves.

Figure 15. The same simulation as the previous but with the expansion acceleration included.

Figure 16. A narrow radiation profile.

The resulting outer boundary curve, including the expansion acceleration, is shown in Figure 17. The density profile was the same as in the previous example. The peak radiation multiplier was 6.

In Figure 18, we show the corresponding temperature profiles.

Note that in this case, while the scaled temperature at the center is somewhat larger than unity, in the outer regions, it is well below unity and thus, is well within the escape velocity limit.

In all the cases we have examined, the temperature profiles are nearly flat but always increase towards the center. This could mean that we just haven’t hit upon the right initial density profile but two other possibilities are more likely. First, we have not yet considered the actual distributions of the radiation sources. In the simulations, we assumed that the distribution was uniform but real sources would be randomly distributed which would affect the temperature profiles. A second possibility is that the cold centers arose later as a result of, for example, Bremsstrahlung radiation.

We now have a result that looks realistic. We have a solution that neither undergoes collapse nor exhibits undue expansion. To achieve this result, we found first, that it is necessary that the density profile has a modest negative slope and second, that the radiation profile must be narrow relative to the free-fall time and have a peak radiation multiplier on the order of 5 or 6.

We could continue with further examples but given the limitations of this model, it would not serve any useful purpose. For many reasons which we will discuss in the next section, further progress will require a 3-dimensional simulation.

Figure 17. Solution with the density profile of Figure 16 and a peak radiation multiplier of 6.

Figure 18. The temperature profile curves.

With these simulations, we have also established important constraints on the radiation profile. From the simulations we find that a radiation energy density a few times larger than the scale value shown in Table 1 was required or $q_{peak}=f1.4\times {10}^{-5}\text{j}\cdot {\text{kg}}^{-1}\cdot {\text{s}}^{-1}$ where f is a multiplier on the order of 5 or 6. Multiplying by the mass of the cluster, we find that total on the order of $f\times {10}^{40}\text{watts}$ or $f2.5M_{\odot }{c}^2$ per year. Of course, that is the absorbed total and the source emission total would undoubtedly have been considerably larger. The second constraint is the persistence of the source or sources. In the previous section, we discovered that the period during which the sources were active had to be short compared to the cluster free-fall time scale. Nevertheless, it still comprised a period on the order of 10⁹ yr. The third constraint applies to the spectrum of the radiation; to heat the gas to a temperature on the order of 10⁷ - 10⁸ K, soft x-rays were required.

Assuming a source, the next step is to understand the heating process.

9. Heating

Since x-ray energies are far greater than the binding energy of the neutral hydrogen, the first effect would have been to ionize the neutral gas. Immediately afterward, the radiation heated the ionized gas via Compton scattering. The total Compton cross section is given by

$\sigma =\pi {\left(\frac{h}{mc}\right)}^2$ (9-1)

and the MFP of the photons by

$\lambda =\frac{1}{4\sqrt{2}n\sigma }$ (9-2)

where n is the number density of the particles. From these equations, we can determine that the heating of the protons was a 2-step process. For protons, the cross section is ${\sigma }_p=5.4\times {10}^{-30}{\text{m}}^2$ so the MFP is ${\lambda }_p=3\times {10}^{25}\text{m}$ which is considerably larger than the size of clusters. Clearly then, direct heating of the protons was not possible; they are just too massive. For electrons, on the other hand, ${\sigma }_e=1.9\times {10}^{-23}{\text{m}}^2$ so ${\lambda }_e=9.5\times {10}^{18}\text{m}$ which is small compared to the characteristic dimensions of the cluster. The first step thus consisted of rapid heating of the electrons by the radiation. Considering now the thermal equilibrium of the electrons, from [2], we find that they would have reached equilibrium on a time scale equal to 10¹³ s so we find that the thermal equilibrium of the electrons would have closely tracked the changing output level of the radiation.

In the second step, the electrons would have interacted with the protons to bring the latter also into thermal equilibrium. The time scale for this second step is 1.9 × 10¹⁶ s which is about ½ of the period during which the sources were active so the protons would have also achieved a degree of equilibrium with the radiation.

10. Galaxy Free-Fall and Supermassive Black Holes

The final problem is to identify the sources. We need a huge power output and even more important, a long lifespan. Supernova produce a prodigious amount of power but they are both short-lived and not particularly common so they don’t qualify. Galaxies with luminous active nuclei (AGN), or in other words, quasars are the only real possibility. Power outputs of observed quasars [8] are typically on the order of 10⁴⁰ watts which is in the range needed. The actual process involves the heading of the gas in an accretion disk surrounding a supermassive black hole at the center of the galaxy. Because the life span of a quasar is determined by the amount of available gas, they can live a long time as evidenced by the existence of the nearby quasar. Once the gas is used up, the quasar process ends and the galaxy becomes a normal galaxy.

Quasars are known to have been more common in the distant past with the peak density occurring at about $t=t_G$. This latter fact is considered to be something of a mystery but this is easily explained by this new model which asserts that all quasars, along with everything else, originated at $t=t_G$. The observed progressive extinction distribution is simply the result of quasars of differing sizes and masses running out of fuel at different times. Also, because galaxies only had so much gas, aside from possible formation during galaxy collisions, there are no second-generation quasars.

To understand how quasar-like objects could have come into existence at precisely the time needed to prevent cluster collapse, we turn to Section 9 of [1]. In that section, we showed that, like the clusters, all galaxies originated during nucleosynthesis and subsequently reached their ZVPs at $t=t_G$. We also showed that they were many times larger than their present-day size at that time and that soon after, began to collapse. The free-fall time would depend on the size and mass of each galaxy but in all cases, it would have been considerably shorter than the cluster free-fall period. The Milky way, for example, had a free-fall time of about 4 × 10¹⁶ s which is approximately the same as the width of the radiation peak of Figure 16. We know that the galaxies did not collapse but in order not to do so, they must have evolved into a stable state complete with their central massive black hole within a period much shorter than their free-fall time. This means that a large percentage, if not all galaxies, already had a black hole at their center well before the time at which the radiation reached its peak. Adding an accretion disk turns the galaxy into a quasar and because all galaxies developed at the same time, so did all quasars. The authors of [9] point out that the quasar UM 287 was in existence at least as early as 3 × 10⁹ yrs after the big bang. This corresponds to a dimensionless time of $\bar{t}=0.92$ which fits very well with our model (compare with Figure 2).

It follows then, that we must understand how a galaxy initially in free-fall could achieve stability. We have already seen that without radiation, the clusters would have collapsed and the same is true of galaxies so a source must have developed within the galaxies. To see how this happened, we return to the simulation equations. Because these are dimensionless, they apply equally well to galaxies. The initial and boundary conditions are the same as for the clusters. The fact that the initial temperature was very close to absolute zero means that without some source of radiation, no pressure would have developed during the initial collapse no matter how much the gas was compressed. This allows us to drop the pressure gradient term from (3-6) which eliminates the oscillation problem. In Figure 19, we show the density profiles that result.

Take note that the scale is logarithmic. We find that the density increases extremely rapidly at the origin and does so long before there is any significant reduction in the outer dimension of the galaxy. Recall that for the Milky Way, for example, the final dimension of the galaxy is about 1/5 of its initial size ( $\bar{r}=0.2$ in the lower panel) so the collapse would have continued for a considerable time beyond what is shown. Eventually, the rotation of the material making up the galaxy would also play a role in the stabilization but initially, as explained in [1], it was the vacuum that was rotating and carrying the particles along with it. That being the case, there was no centrifugal acceleration until after the galaxy had undergone a significant degree of compaction.

We can now understand the evolution. The initial collapse caused an extremely rapid increase in the central density which resulted in the formation of a supermassive black hole. As the collapse continued, the in-falling gas would then have formed an accretion disk which, in turn, produced the radiation that heated the intergalactic gas. This, in turn again, resulted in the development of the pressure gradient that stopped the collapse.

The key point is that the creation of supermassive black holes came first.

11. Virgo

We have seen that the output from just a few powerful quasars would have been enough to heat the cluster gas but that idea has problems. First, just a few

Figure 19. Density profiles for galaxy free-fall. The solution fails at the point shown because of the steepness of the curves near the origin.

randomly located quasars would not be able to heat the gas uniformly. Second, such quasars have long lifespans and we now know that the sources, at least collectively, must have been short-lived. What seems to be more likely than a scenario with a few very powerful quasars is one with a multitude of mini-quasars.

Referring to the Virgo cluster specifically, we saw in [1] that ETGs dominate its galaxy population. There are 513 ETGs inside the cluster with about 75% being either ellipticals or dwarf ellipticals. Their distributions are shown in Figure 20. The left-hand panel shows just the ellipticals, of which there are 26. The right-hand panel includes also the dE and dEN galaxies of which there is a total of 361.

The reason this is interesting is that all the known quasar host galaxies [10] are ETGs with ellipticals predominating. It is also shown in that same paper that the fraction of massive spheroids/black holes that produce quasar-level activity is about 0.1 percent at present and was greater than 10 percent at $z\approx 2\text{-}3$. If we assume a half-life decay model and extrapolate back to $t_G$, we find a fraction somewhat greater than 15%. but, based on our galaxy stability argument, the actual number was probably 100%. Clearly, from the figure, the distribution of the dwarfs would have a much better chance of uniformly heating the gas than would the few ellipticals.

For both this and the stability idea to work, however, dwarf galaxies must contain mini-supermassive black holes. That has long been a question but now, according to some very recent work [11], we find that they do. And, not only do they have black holes, but they also exhibit active nuclei.

The conclusion is that the clusters were heated from within by their own galaxies and the strongest argument for this model is that there is no other scenario that could account for both galaxy and cluster stability.

We are now at the point where a 3-dimensional simulation is needed to settle several questions. First, in the simulations presented here, it was assumed that

Figure 20. Distribution of elliptical galaxies in the Virgo cluster.

the radiation absorption was evenly distributed whereas the actual distribution would to some extent reflect the random distribution of the sources. The question is whether or not the resulting uneven heating will spatially thermalize quickly enough to prevent clumping and collapse. Another issue is the matter of the dispersion process discussed earlier needed to prevent clumping of the particle distributions after the initial heating was completed. We don’t know if the cluster underwent an expansion followed by a contraction or just remained at its present size which is another question to be studied. In connection with an expansion, another effect that needs to be considered is that, as the cluster gas undergoes expansion, the outer portions would be moving away from the radiation sources and hence would cool. Lastly, we know that the initial density distribution was fixed at the time of nucleosynthesis but the simple simulation model presented here is far too limited to allow us to say anything more about it than that initial distributions must have had a modest negative slope.

Lastly, on a more basic level, because vacuum energy density plays an important role in the cluster evolution, 3-dimensional simulations based on the full set of Einstein’s equations including vacuum energy are needed and this is, even more, the case for understanding galaxy evolution.

12. Conclusions

In this paper, we discussed the evolution of galaxy clusters during the period immediately following the ZVP of their development. We modeled the cluster using the normal energy and momentum conservation equations for a dilute gas. Our goal was to determine what constraints are imposed on the initial gas density and radiation profiles by the requirement that the cluster neither undergoes free-fall collapse nor evaporation. By solving the gas equations for many possible initial density and radiation profiles, we found that the allowable density distributions had to have a moderate negative slope and that the radiation profile must be narrowly peaked in time with a width considerably less than the free-fall time of the cluster. (Even though we speak of the radiation peak as narrow, its actual width was on the order of 10⁹ yr.)

We calculated that the power required to stabilize the cluster was on the order of a few times 10⁴⁰ watts which must be sustained for the period just mentioned which leads immediately to the conclusion that quasars were the only possible source of the radiation since they are both long-lived and have power outputs of the necessary magnitude.

The Virgo cluster does not contain a quasar at present which leads us to argue that instead of a few massive quasars, the cluster was powered by a large number of mini-quasars that came into existence immediately after ZVP of the cluster. These proceeded to use up their fuel supplies during the same period to become ordinary ETG galaxies we now see populating the clusters.

Having established that mini-quasars were the source of the radiation, we were then in the position of needing to explain the origin of the quasars. In [1], we showed that galaxies reached their ZVPs coincident with the galaxy clusters and that their sizes at that point were many times larger than their present-day sizes. Initially, all galaxies consisted of nothing but very cold gas and like the cluster, they began to undergo gravitational free-fall collapse since there was no pressure gradient within the gas to pervert that from happening. The consequence was that during the early stages of the free-fall, the density of the gas increased dramatically at the center of the galaxy which resulted in the formation of a supermassive or mini-supermassive black hole. Once the black hole formed, an accretion disk formed from the infalling gas setting the conditions for an active galaxy nucleus. The resulting radiation headed the gas which in turn, allowed pressure gradients to develop that stabilized both the galaxies and the clusters in which they happened to reside. This process applies to all galaxies so they all must have black holes at their centers, for if they didn’t, they would have rapidly undergone gravitational free-fall collapse.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Botke, J.C. (2021) The Origin of Cosmic Structures Part 1—Stars to Superclusters. Journal of High Energy Physics, Gravitation and Cosmology, 7, 1373-1409. https://doi.org/10.4236/jhepgc.2021.74085
[2]	Vogelsberger, M., Marinacci, F., Torrey, P. and Puchwein, E. (2019) Cosmic Simulations of Galaxy Formation. arXiv:1909.07976v6. https://arxiv.org/abs/1909.07976
[3]	Norman, M.L. (2010) Numerical Simulations of Gas in Galaxy Clusters. arXiv:1005.1100. https://arxiv.org/pdf/1005.1100.pdf
[4]	Borgani, S. and Kravtsov, A. (2011) Cosmological Simulations of Galaxy Clusters. Advanced Science Letters, 4, 204-227. https://doi.org/10.1166/asl.2011.1209
[5]	Armitage, T. (2019) Simulations of Galaxy Clusters. PhD Thesis, University of Manchester, Manchester. https://www.research.manchester.ac.uk/portal/en/theses/simulations-of-galaxy-clusters(a450bc31-2865-46e7-9816-7e9ca7129e46).html
[6]	Sarazin, C.L. (1988) Mean Free Paths and Equilibration Time Scales. https://ned.ipac.caltech.edu/level5/March02/Sarazin/Sarazin5_4.html
[7]	Feng, G. (1998) Data Smoothing by Cubic Spline Filter. IEEE Transitions on Signal Processing, 46, 2790-2796. https://doi.org/10.1109/78.720380
[8]	Wikipedia (2022) Quasar. https://en.wikipedia.org/wiki/Quasar
[9]	Cantalupo, S., Arrigoni-Battaia, F., Frochaska, J.X., Hennawi, J.F. and Madau, P. (2014) A Cosmic Web Filament Revealed in Lyman-Alpha Emission around a Luminous High-Redshift Quasar. Nature, 506, 63-66. https://www.nature.com/articles/nature12898 https://doi.org/10.1038/nature12898
[10]	Dunlog, J.S., et al. (2003) Quasars, Their Host Galaxies and Their Central Black Holes. Monthly Notices of the Royal Astronomical Society, 340, 1095-1135. https://academic.oup.com/mnras/article/340/4/1095/1322304 https://doi.org/10.1046/j.1365-8711.2003.06333.x
[11]	Hickox, R.C. and Parker, J. (2022) Uncovering as Hidden Mini-Monster: A Heavily Obscured Active Galactic Nucleus in a Dwarf Star-Forming Galaxy. https://aas.org/sites/default/files/2022-01/AAS239_Mon1_RyanHickox.pdf