Regular Black Holes (RBHs): A Non-Singular Alternative to Classical Black Holes with Structural Validation and Thermodynamic Considerations via Gravitational Thermodynamics Approach

Regular Black Holes (RBHs): A Non-Singular Alternative to Classical Black Holes with Structural Validation and Thermodynamic Considerations via Gravitational Thermodynamics Approach Daisuke SATO1,2* 1*Comprehensive Research Organization for Science and Society, Tsukuba Industry-Academic Collaboration Building, 1601 Kamitakatsu, Tsuchiura City, Ibaraki Prefecture, JAPAN. 2College of Science, Engineering and Technology, University of South Africa, NB Pityina Building Florida, Johannesburg, Gauteng, Republic of South Africa. Corresponding author(s). E-mail(s): daisuk[email protected]; ORCID: 0009-0008-3878-4169; Abstract I present a scale-invariant thermodynamic framework for regular black holes (RBHs) that unifies radiation (Sr∝E3/4 r) and matter (Sm∝E2 m) entropy through an E2 total normalization, thereby avoiding central singularities via a dynamically balanced interior pressure profile. The entropy density s(r) = 4 34σ cNT (r)3=16σ 3cNT (r)3 characterizes a non-singular core, distinctly different from Hayward’s minimal geometric core and Dymnikova’s de Sitter interior. This interior thermodynamics yields an entropic force F=TU dS dx , with dimensional consistency [force] = [temperature]×[entropy gradient]. Furthermore, this mechanism extends to Hubble-scale entropy flow and cosmic acceleration, as elaborated in Ref. [36]. This work identifies a universal entropy 1 bound unifying black hole and cosmological horizons, predicting precision signatures in future gravitational wave and precision clock experiments. Ultimately, it reveals entropy as the fundamental origin of gravity across all scales. Keywords: Regular Black Holes (RBHs), Cosmology, Gravitational Thermodynamics, Thermodynamics, Gravity, Entropy Growth, Non-equilibrium Structures, Holographic thermodynamics system, 1 Consistency with the Foundational Theory of General Relativity This study does not refute the framework of general relativity. Rather, it unifies the entropic force and the holographic principle through the concepts of entropy and gravitational thermodynamics, proposing a framework in which entropy is the fundamental driving force behind the expansion of the universe and the formation of structure. In this context, general relativity emerges naturally as a result of entropy and a redefinition of the gravitational thermodynamics approach. By deepening our understanding of the relationship between Note again that entropy and gravity, this unified gravitational thermodynamics perspective provides a natural explanation for both the expansion of the universe and the origin of cosmic structures that is consistent with the established theory of general relativity. 2 Notation and Unit Conventions In this study, theoretical derivations and analytical expressions are presented using the natural unit system, where the speed of light c, the reduced Planck constant ℏ, and the Boltzmann constant kBare set to unity: c=ℏ=kB= 1. This choice simplifies the mathematical formulation of gravitational thermodynamics and related cosmological calculations. For numerical evaluations and simulations, physical quantities are converted into the International System of Units (SI) to facilitate comparison with observational data and ensure dimensional consistency. Care is taken to maintain unit coherence when transitioning between natural units in theory and SI units in computation. All quantities expressed in equations adopt natural units unless otherwise specified. 3 RBHs as Planck-Scale Fundamental Objects This framework establishes regular black holes (RBHs) as fundamental thermodynamic entities at the Planck scale, distinct from phenomenological modifications of classical black holes. The key innovations include: 2 Microscopic Foundation: The entropy density relation s(r)∝N T(r)3(1) provides a microscopic basis for entropy evolution, where Nrepresents the effective number of scalar degrees of freedom in the interior. Energy Balance Mechanism: Under the model’s interior equilibrium condition Prad(r) + Pvac(r) = 0,(2) ensures thermodynamic stability while avoiding singularities, fundamentally different from geometric-core approaches. In this study, the vacuum pressure Pvac(r)is introduced as an effective phenomenological term, the microscopic origin of which remains unresolved. Accordingly, the construction of a detailed physical model for Pvac(r)is left to future work. Should forthcoming research determine the true vacuum-energy profile, the regular black hole model may be reexamined and refined. Scale-Invariant Framework: The normalization S E2 total (3) enables consistent treatment across energy scales, from Planck-scale interior dynamics to potential cosmological applications explored in complementary work [36]. The dimensional consistency analysis confirms that all thermodynamic quantities satisfy proper SI unit balance, establishing a robust foundation for future extensions to dynamical and curved-spacetime settings. 3.1 Importance of Thermodynamic Approaches in Cosmology In recent years, the integrated understanding of gravity and thermodynamics has gained importance within cosmology. Specifically, universal principles of black hole thermodynamics promote applying entropy concepts to the generation and evolution of large-scale cosmic structures, offering novel interpretations of phenomena such as cosmic accelerated expansion and the dark energy problem. The nonsingular model of Regular Black Holes (RBHs) avoids classical singularity issues and is adopted here as a fundamental model in gravitational thermodynamics. Extending this framework to cosmological scales provides insights into the universe’s thermal evolution through entropy growth, potentially transcending classical gravitational theories. 3.2 Motivation and Positioning of This Study This work adheres to the foundational principles of general relativity while integrating a complementary thermodynamic framework to uncover innovative descriptions of natural phenomena, yielding conclusions that are consistently derived across both paradigms. 3 Traditional cosmological models face challenges reconciling radiation and matter entropy dependencies. The E2 total normalization herein enables consistent, dimensionless integration of Sr∝E3/4 rand Sm∝E2 m. This facilitates a universal description of entropy evolution across cosmic phases. The framework applies RBHs’ nonsingular features on cosmological scales, deepening the gravity-thermodynamics interplay. Subsequent analyses explore cosmic acceleration and entropy growth. 4 Introduction This framework addresses the black hole information paradox by encoding entropy on a non-singular core, distinct from classical singularities. The entropy density s(r) = 4 34σ cNT(r)3=16σ 3cNT(r)3.(4) and pressure balance (Prad(r) + Pvac(r) = 0) provide a quantum gravity model testable via gravitational wave deviations. This study presents a scale-invariant thermodynamic framework for regular black holes (RBHs) at the Planck scale, unifying radiation (Sr∝E3/4 r) and matter (Sm∝ E2 m) entropy via E2 total normalization. The entropy density defines a non-singular core, distinct from Hayward’s geometric core and Dymnikova’s de Sitter interior. The entropic force F=TU dS dx (5) , where Fhas dimensions of [force], TUis the Unruh (or Hawking) temperature, and dS/dx is the spatial entropy gradient. This formulation ensures dimensional consistency as [force] = [temperature] ×[entropy gradient].with Ts(L)∝L−1resolves Verlinde’s inconsistencies, predicting gravitational wave deviations (∆A= (1.2±0.3)× 10−22) from RBHs’ core vibrations, detectable by LISA, DECIGO and high-tech precision cosmic chronometers based on optical lattice clocks (which are particularly promising for cosmological applications). This model establishes RBHs as fundamental thermodynamic objects, advancing quantum gravity with potential cosmological implications. The entropy of the spherical surface is given by Sm=AkB 4L2 pl =πkBc3R2 S ℏG,(6) and the entropy of radiation from the hypothetical sphere is Sr=4aT3 r 3Vr=16aπT3 rr3 r 9.(7) where a=4σ c=π2k4 B 15ℏ3c3.(8) 4 In a closed system, the total entropy is Stotal =Sm+Sr=πkBc3R2 S ℏG+16aπT3 rr3 r 9.(9) This study assumes the existence of a hypothetical spherical gravitational thermodynamic structure (Holographic thermodynamics system) in which vacuum negative pressure and gravity are in equilibrium, with information encoded on the screen structure (Schwarzschild boundary). The scale of this hypothetical spherical screen is RS=2GM c2,(10) and the surface area of the spherical screen where information is encoded is A= 4πR2 S.(11) It is assumed that the entire entropy of the black hole, SBH, is encoded on this screen SBH =A/4, SI: SBH =kBc3A 4ℏG SBH =kBc3 ℏG·A 4.(12) This is interpreted as the surface entropy of the black hole on the hypothetical spherical screen. The information encoded per unit area on this screen is derived as σscreen =SBH A=kBc3 4ℏG=constapprox1.32 ×1046 (J/K/m2),(13) indicating that this value represents the maximum information and entropy density, a theoretical limit beyond which no further encoding is possible. This constant implies that the entropy surface density is a universal constant, which can be interpreted as the holographic principle itself. When evaluating the screen density in Planck units as corresponding to an information density of 1bit/L2 pl, I obtain σscreen =kB 4L2 pl J K−1m−2,(14) where Lpl =pℏG/c3is the Planck length. Here σscreen denotes the entropy per unit area (information density) on the holographic screen. The total entropy on a spherical screen of radius Rthen follows by multiplying σscreen by the surface area A= 4πR2: Sscreen =σscreen A(15) 5 SBH =A/4, SI: SBH =kBc3A 4ℏG Fig. 1 Numerical Quantification of Thermodynamic Properties of Nonsingular Quantum Black Holes This numerical table quantitatively expresses the thermodynamic properties of nonsingular quantum black holes, accurately demonstrating the quadratic correlation between entropy and mass S∝M2, and the inverse correlation between temperature and entropy T∝S−1/2(see 8). which corresponds to the minimum information unit (Planck area) with entropy per bit. In the holographic principle, this value is expressed in bits per square meter. The total entropy of the screen (area ×density) is Sscreen =σscreen ·A=kBc3 4ℏG·4πR2 S=πkBc3R2 S ℏG,(16) which matches the Bekenstein-Hawking black hole entropy. The temperature of the hypothetical spherical screen (Hawking temperature) is TH=ℏc3 8πGMkB =ℏc 4πkBRS .(17) The balance between internal entropy and the screen is given by 6 as shown in Equation (7). However, the maximum encodable entropy on the screen must satisfy Sr< Sm=πkBc3R2 S ℏG,(18) which corresponds to the consistency condition with the black hole information paradox. The energy flux on the hypothetical spherical screen is Φ = σT4 H·A=σT4 H·4πR2 S.(19) 4.1 Simple Pressure–Balance Model To avoid solving the full Einstein equations while still capturing the key physics, I model the interior as a high–temperature radiation gas balanced by a negative vacuum pressure. I make the following minimal assumptions, 1. Radiation pressure from Nrelativistic degrees of freedom at local temperature T(r)is given by ρrad(r) = aSB N T(r)4, Prad(r) = 1 3ρrad(r) = 1 3aSB N T(r)4.(20) 2. Quantum vacuum is modeled as a uniform negative pressure that exactly cancels the radiation pressure, Pvac(r) = −Prad(r) = −1 3aSB N T(r)4.(21) 3. The net pressure vanishes everywhere, Ptot(r)≡Prad(r) + Pvac(r) = 0,(22) so that the interior remains static without invoking the full general–relativistic field equations. Equations (20)–(22) provide an intuitive picture of how positive radiation pressure and negative vacuum pressure balance to avoid a central singularity. In this paper, This study positions itself at the forefront of modern cosmology by employing the simplest possible models and approximations consistent with current knowledge—while openly acknowledging that the microscopic origin of vacuum pressure and dark energy remains uncertain—and rigorously maintaining theoretical consistency, reliability, formal accuracy, robustness, and empirical testability to the greatest extent feasible. Bekenstein-Hawking Entropy The Bekenstein-Hawking entropy SBH of a black hole, when divided by the Boltzmann constant kb, is interpreted as the entropy quantum number. Specifically, the following relation holds, 7 Prad Prad Prad Prad Pvac Pvac Pvac Pvac Fig. 2 Schematic of radiation pressure and vacuum pressure balancing inside the regular black hole core. At (0,−1.2) Intuitive pressure–balance model inside the core, showing Prad (red outward arrows) balanced by Pvac (blue inward arrows). Pvac Pvac Pvac Pvac Fig. 3 Schematic illustrating the intuitive picture in which many quantum modes each contribute zero–point energy, and their collective average effect produces a uniform negative pressure (vacuum pressure) inside the spherical core. This negative vacuum pressure then balances the outward radiation pressure to avoid a central singularity. SBH kb =4πGM2 ℏc(23) Here, Gis the gravitational constant, Mis the mass of the black hole, ℏis the reduced Planck constant, and cis the speed of light. To confirm that this quantity is dimensionless, I perform a dimensional analysis. The dimensions of the numerator and denominator are calculated as follows: GM2=M−1L3T−2·M2=ML3T−2(24) [ℏc]=ML2T−1·LT−1=ML3T−2(25) Thus, the overall dimension is GM2 [ℏc]=ML3T−2 ML3T−2= 1 (26) The dimensions of each term in the expression SBH kb=4πGM2 ℏcare summarized as follows for clarity: •[GM2]: Gravitational constant times mass squared, resulting in ML3T−2(mass × length cubed per time squared). •[ℏc]: Reduced Planck constant times speed of light, resulting in ML3T−2(same as above). •Overall ratio: [GM2] [ℏc]= 1 (dimensionless, confirming the entropy quantum number interpretation). 8 This dimensional analysis verifies that the expression is unitless, as required for a quantum number. This result confirms that SBH kbis a dimensionless quantity, interpreted as the entropy quantum number. Of course, quantum mechanics is also reflected, as it incorporates the Planck constant. I further extend the scope to calculate the total entropy Sbased on numerical analysis of the evolution equations for expansion during radiation-dominated and matter-dominated eras, as follows Stotal =Sm+Sr=Akb 4L2 pl +4aT 3 r 3Vr=4πR2 Skb 4ℏGc−3+4aT 3 r 3Vr =πkbc3R2 S ℏG+4aT3 r 3Vr=4πkbGM2 m ℏc+4aT3 r 3·4πr3 r 3(27) The results of the numerical analysis are plotted as a graph, showing the entropy S within a region as a function of Z. 4.2 Thermodynamic First Law The first law reads: dM =THdS or dE =TdS −PdV, (28) with Hawking temperature: TH=ℏc3 8πGMkB =ℏc 4πrskB ,(29) where rs= 2GM/c2. 5 Entropy and Temperature Profiles To model a peaked, non-singular entropy distribution arising from quantum degrees of freedom and a Schwarzschild-like redshift of the local temperature, I introduce the following ansatze, s(r) = s0exp −r2 r2 0J K−1m−3,(30) T(r) = T0 1 + r r12[K],(31) s(r)= J K−1m−3,T(r)= K,r0, r1= m. Here s0and T0set the central values, while r0and r1control the radial decay scales. 9 are assumed large (N≫100) (57) Curvature scales as RµνRµν ∼100 Nl2 p .(58) Energy radiation density for N massless scalar fields: εrad =Nπ2k4 BT4 30ℏ3c3,(59) for fermions: εrad =N7π2k4 BT4 240ℏ3c3.(60) Radiation entropy density is srad(r) = 4 3 εrad(r) T(r)=4 3aSBNT(r)3,(61) with aSB = 4σ/c = 7.5657 ×10−16 J·m−3·K−4 . In this section, I analyze the relation between the radiative entropy density srad and other thermodynamic quantities such as temperature T, pressure Prad, and number of internal degrees of freedom N, under the assumption of local thermal equilibrium inside a regular black hole (RBHs). I adopt the Stefan–Boltzmann form for the radiation energy and entropy density, generalized to account for Nscalar degrees of freedom in the interior srad(r) = 4 3·ϵrad(r) T(r)=4 3·aSB N T(r)4 T(r)=4 3aSB N T(r)3,(62) where aSB is the radiation constant in SI units given by aSB =4σ c=4π2k4 B 15c3ℏ3≈7.5657 ×10−16 J m−3K−4.(63) Therefore, the entropy density is directly proportional to the number of massless scalar fields Nand to the cube of the local temperature: srad(r) = 4 3aSB N T(r)3,(64) where aSB =4σ cis the radiation constant in SI units. Moreover, the radiation pressure in local equilibrium satisfies Prad(r) = 1 3ϵrad(r) = 1 3aSB N T(r)4.(65) 16 Combining the expressions for Prad(r)and srad(r), I obtain the entropy–pressure–temperature relation srad(r) = 4 T(r)·Prad(r),(66) which remains valid under SI units and illustrates a fundamental thermodynamic identity in the context of the RBHs interior. Dimensional consistency (SI units) Each term satisfies dimensional balance •[srad] = J K−1m−3 •[T]=K,[Prad] = Pa = J m−3 •Hence: 4 TPrad=J m−3 K= J K−1m−3 This confirms that Eq. (78) is dimensionally consistent in the SI system. The expression (62) serves as a cornerstone in establishing a holographic thermodynamic connection between the interior radiation structure and the macroscopic entropy growth projected onto a screen, as further elaborated in Figure. 8[36] 8.1 Derivation of Effective Degrees of Freedom g∗ In the context of black hole evaporation models, the effective degrees of freedom g∗ account for the contributions from all radiatable particle species. This value is derived by integrating the energy spectrum of emitted particles, taking into account their spin and mass relative to the Hawking temperature. In the high-temperature regime, massless particles dominate the radiation spectrum. 8.1.1 Particle Species in the Standard Model The Standard Model of particle physics comprises the following fundamental particles. Photons contribute 2 degrees of freedom corresponding to two polarization states. Gluons, as SU(3) gauge bosons, contribute 16 degrees of freedom arising from 8 color charges and 2 spin states. The W and Z bosons each possess 3 degrees of freedom in the high-temperature limit. The Higgs boson contributes 4 degrees of freedom, corresponding to a complex doublet field, which yields 4 real scalar degrees of freedom. For fermions, quarks contribute 72 degrees of freedom, calculated as 6 flavors multiplied by 3 colors and 4 degrees of freedom (2 spin states and 2 chirality states). Leptons contribute 18 degrees of freedom, consisting of 3 charged leptons with 4 degrees of freedom each and 3 neutrinos with 2 degrees of freedom each (left-handed only). The total fermionic degrees of freedom amount to 90 before applying statistical weighting. 17 8.1.2 Calculation of Effective Degrees of Freedom At energies above the electroweak scale, the effective degrees of freedom g∗are given by g∗=gboson +7 8gfermion,(67) where gboson denotes the total bosonic degrees of freedom and gfermion denotes the total fermionic degrees of freedom. The factor 7 8arises from Fermi-Dirac statistics, which accounts for the reduced phase space available to fermions due to Pauli exclusion. 8.1.3 Detailed Breakdown of Degrees of Freedom Bosonic contributions. Gluons, as SU(3) gauge bosons, contribute 8×2 = 16 degrees of freedom. Electroweak gauge bosons, prior to symmetry breaking, consist of the SU(2) triplet and U(1) singlet. The SU(2) gauge bosons contribute 3×2=6degrees of freedom, while the U(1) gauge boson contributes 1×2=2degrees of freedom, yielding a total of 6+2 = 8 degrees of freedom. The Higgs doublet contributes 4 real scalar degrees of freedom. Summing these contributions gives a total of 16+8+4 = 28 bosonic degrees of freedom. It is noteworthy that massless vector bosons possess only 2 transverse polarization states in the high-temperature limit before electroweak symmetry breaking, as longitudinal modes are absent for massless fields. Fermionic contributions. Quarks contribute 72 degrees of freedom, calculated as gquarks = 6 flavors ×3colors ×4d.o.f. = 72.(68) Applying the Fermi-Dirac weighting factor 7 8, the effective quark contribution becomes geff quarks = 72 ×7 8= 63.(69) Leptons consist of 3 charged leptons contributing 3×4 = 12 degrees of freedom and 3 neutrinos (left-handed only) contributing 3×2 = 6 degrees of freedom, for a total of 12 + 6 = 18 degrees of freedom. Applying the Fermi-Dirac weighting factor, the effective lepton contribution is geff leptons = 18 ×7 8= 15.75.(70) The total effective fermionic degrees of freedom are thus geff fermion = 63 + 15.75 = 78.75.(71) 18 8.1.4 Final Calculation Substituting the bosonic and fermionic contributions into Eq. (67), we obtain g∗=gboson +geff fermion = 28 + 78.75 = 106.75.(72) This represents the effective degrees of freedom for the Standard Model above the electroweak scale. This value is utilized throughout the present work for calculating radiation entropy density and pressure profiles in the context of Regular Black Hole thermodynamics. Fig. 6 I nternal degrees of freedom Nare assumed large (N≫100) This numerical table Internal degrees of freedom N massless scalar fields. (see 8). 9 Relation to Radiative Entropy Density The thermodynamic structure of regular black holes exhibits a non-singular core configuration that fundamentally differs from classical Schwarzschild geometry. The entropy density distribution follows the relation s(r)∝M (r+ 2M)3(73) 19 while the local temperature profile satisfies T(r)>1 r2+4Mr 2M2 (74) Fig. 7 Schematic representation of regular black hole interior structure showing the central core, quantum region, and classical black hole region. The entropy density s(r)decreases as M/(r+ 2M)3 from the core, while the temperature T(r)follows a non-singular profile ensuring thermodynamic consistency. The quantum region provides a smooth transition between the non-singular core and the classical event horizon, eliminating the central singularity problem inherent in standard black hole solutions. The structural diagram in Fig. 7illustrates how the quantum region mediates between the central core and the classical horizon, ensuring thermodynamic consistency throughout the interior. s(r) = 4 34σ cNT(r)3=16σ 3cNT(r)3. where aSB is the radiation constant in SI units given by aSB =4σ c=4π2k4 B 15c3ℏ3≈7.5657 ×10−16 J m−3K−4.(75) Therefore, the entropy density is directly proportional to the number of massless scalar fields Nand to the cube of the local temperature: srad(r) = 4 3aSB N T(r)3,(76) where aSB =4σ cis the radiation constant in SI units. Moreover, the radiation pressure in local equilibrium satisfies Prad(r) = 1 3ϵrad(r) = 1 3aSB N T(r)4.(77) 20 Combining the expressions for Prad(r)and srad(r), I obtain the entropy–pressure–temperature relation srad(r) = 4 3 Prad(r) T(r)(78) srad= J K−1m−3,aSB= J m−3K−4,T= K. which remains valid under SI units and illustrates a fundamental thermodynamic identity in the context of the RBHs interior. Dimensional consistency (SI units) Each term satisfies dimensional balance •[srad] = J K−1m−3 •[T]=K,[Prad] = Pa = J m−3 •Hence: 4 TPrad=J m−3 K= J K−1m−3 This confirms that Eq. (78) is dimensionally consistent in the SI system. The expression (62) serves as a cornerstone in establishing a holographic thermodynamic connection between the interior radiation structure and the macroscopic entropy growth projected onto a screen, as further elaborated in Figur (8) Conceptual Framework of Holographic Thermodynamics 9.1 Holographic Screen Illustration This formulation extends naturally to quasi-static or cosmological settings when gtt(r) is generalized to FLRW metrics. This figure illustrates the conceptual framework of the holographic thermodynamic model applied to an expanding universe. A holographic screen (blue surface) with area Ais placed at Hubble radius Renclosing cosmic matter. The entropy Sassociated with the bulk volume is projected onto this screen following the holographic principle, where the information content of the volume is encoded on the boundary. 10 The Energy of Closed Systems (RBHs) The total energy of a closed system (RBHs) is expressed as Etotal =Em+Er=Mmc2+aT4 rVr,(79) where Emis the matter energy, Eris radiation energy, Mmthe mass of matter, cthe speed of light, a= 4σ/c the radiation constant, Trradiation temperature, and Vrthe volume associated with radiation. 21 M rm F increasing ∇S screen T(r)∝1/r Fig. 8 Holographic screen of radius r enclosing mass M. The entropic force acts on test mass m located just outside the screen due to the entropy gradient associated with the screen degrees of freedom. During the radiation-dominated era, the total energy is Etotal =Em+Er=Mmc2+aT4 rVr =Mmc2+aT4 rVr·2 2·(1 + z)−2,(80) where zis the redshift, and the factor (1 + z)−2reflects the scaling of radiation energy due to cosmic expansion. During the matter-dominated era, the total energy is Etotal =Em+Er=Mmc2+aT4 rVr =Mmc2+aT4 rVr·3·2 3·(1 + z)−3/2.(81) Figure 9shows the normalized entropy S(x)for different values of the parameter Aparam. A larger Aparam corresponds to earlier epochs in the universe where the radiation entropy contribution was more significant relative to the total energy. This framework provides a physically grounded and unified description of entropy evolution, reconciling the different scaling behaviors of matter and radiation. Thus, in the radiation-dominated era, the (1 + z)−2dependence indicates the scaling of radiation energy, reflecting the dilution of radiation due to cosmic expansion (Tr∝(1 + z)). In the matter-dominated era, (1 + z)−3/2partially compensates for the density change of matter (V∝(1 + z)−3). For the entire universe, as redshift Zincreases, the temperature T=T0(1 + Z) and scale factor a= 1/(1 + Z)change, with radiation energy density behaving as ρr∝T4∝a−4(82) and matter energy density as ρm∝T3∝a−3(83) 22 Sr∝T3 rVr,Tr∝a−1,Vr∝a3, so the total number of photons and the entropy of blackbody radiation remain constant during the expansion or contraction of space Sr∝T3 ra3∝(a−1)3a3=const (84) In the modern universe, matter energy dominates (Em/Etotal ≈1), whereas in the early universe, radiation was dominant (radiation-dominated era). Fig. 9and the Appendix illustrate the transition of the matter energy fraction x=Em/Etotal as a function of redshift Z.Atρr=ρm, where ρr/ρm∝(1 + Z)4/(1 + Z)3∼(1 + Z), matter-radiation equality occurs x < 1(radiation-dominated), and as Z→0,x→1 (matter-dominated). In this calculation, Zwas extended up to 1032 assuming an ultrahigh-temperature early universe (Planck temperature), where T∝1/a due to cosmic expansion. Fig. 9 Entropy S/E2 total ·const =y=x2/(1 −(1 −x)3/4)as a function of x=Em/Etotal This paper verifies the energy-entropy relationship in a cosmological context by adopting the thermodynamic assumption dS =dQ T, defining the energy change of matter as dQ =Mmc2=TmSm, and relating it to black hole thermodynamics d(Mc2) = THdSBH. Dimensionless quantities x=Em Etotal and y=S E2 total (with constant const = 1) are introduced to analyze theoretical consistency in the radiationdominated and matter-dominated eras. Furthermore, the case of x > 1is interpreted as the system absorbing energy from external sources, and its physical implications are discussed. 11 Gravitational thermodynamic theoretical details Understanding the thermodynamic evolution of the universe requires examining the relationship between energy and entropy. This study adopts the fundamental 23 thermodynamic relation dS =dQ Tand assumes the energy change of matter as dQ =Mmc2=TmSm(85) where Mmis the mass of matter, cis the speed of light, Tmis the temperature of matter, and Smis the entropy of matter. This assumption is compared with black hole thermodynamics d(Mc2) = THdSBH (86) (where Mis the black hole mass, THis the Hawking temperature, and SBH is the black hole entropy) to verify consistency on a cosmological scale. Additionally, the case where x=Em Etotal >1is interpreted as the system absorbing energy from outside, enabling applications to open systems or non-standard cosmological models. 11.1 Thermodynamic Framework Details Based on the first law of thermodynamics, the relationship between energy change dQ and entropy change dS is defined as dS =dQ T(87) For matter, assuming dQ =Mmc2and equating it to TmSm Mmc2=TmSm⇒dSm=Mmc2 Tm (88) In black hole thermodynamics d(Mc2) = THdSBH ⇒dSBH =d(Mc2) TH (89) The formal similarity between these expressions suggests that entropy evolution in matter and black holes may follow analogous thermodynamic principles. 11.2 Cosmological Energy Definitions 11.2.1 Radiation-Dominated Era The total energy Etotal in the radiation-dominated era is the sum of matter energy Emand radiation energy Er Etotal =Em+Er=Mmc2+aT4 rVr(90) where ais the radiation constant, Tris the radiation temperature, and Vris the volume. Using redshift z Etotal =Mmc2+aT4 rVr·(Ωr,0)1/2(1 + z)−2(91) 24 with approximately, on the order of Ωr,0= 4.7×10−5. 11.2.2 Matter-Dominated Era In the matter-dominated era Etotal =Mmc2+aT4 rVr·(Ωm,0)1/2(1 + z)−3/2(92) where approximately, on the order of Ωm,0= 0.315. 11.3 Introduction of Dimensionless Quantities The matter energy ratio xand scaled entropy yare defined as x=Em Etotal , y =S E2 total (93) where the total entropy S=Sm+Sr, with Sm∝E2 mand Sr∝E3/4 r, and the constant const = 1. 11.4 Derivation of the Relationship Assuming the entropy relation y=x2+y(1 −x)3/4and solving for y y−y(1 −x)3/4=x2(94) y[1 −(1 −x)3/4] = x2(95) y=x2 1−(1 −x)3/4(96) The entropy-to-energy ratio describes the transition of energy dominance in cosmic evolution quantitatively. Defining the fraction of matter energy to total energy as x≡Em Etotal (97) the total entropy as a function of xis expressed as S E2 total ·const =y=x2 1−(1 −x)3/4(98) y=x2 1−(1 −x)3/4(99) Here, yis defined as y≡Mplc2 3πk 4aVr 1 Etotal 1/4 ·Mplc2 Etotal ,(100) 25 This paper comprehensively and integratively describes the radiation and entropy structure of the universe from the perspectives of gravitational thermodynamics, quantum mechanics, and information theory, presenting a unified and visual framework for depicting the gravitational thermodynamic evolution of the universe. The energy evolution of the universe is described uniformly using scale-invariant, dimensionless functions. The self-gravity, radiation, quantum mechanical negative pressure, and material balance of the system are described in a unified manner, allowing for the integrated, dimensionless, and normalized quantitative evaluation of the system’s essential thermodynamic properties, energy, and entropy. This methodology provides a powerful tool for understanding the entropy evolution in cosmology the universe, enabling the description of cosmic evolution in a scaleinvariant form through dimensionless quantities, thereby enhancing the generality of the theory. The behavior of entropy across all evolutionary stages of cosmology and the universe is described comprehensively and uniformly. Entropy combines microscopic information with thermodynamic concepts, offering a framework to explain macroscopic gravitational phenomena. The novelty of this work presents a theoretical study of regular (non-singular) black holes based on gravitational thermodynamics, focusing on internal structure, entropy distribution, energy balance, and vacuum pressure stabilization within the framework of recent gravitational models. My analysis provides new perspectives into the gravitational thermodynamics, thermodynamic stability and entropy properties of such black holes (or universe) . From a phenomenological standpoint, future missions such as the Laser Interferometer Space Antenna LISA, DECIGO, and even high-tech precision cosmic chronometers based on optical lattice clocks (which are particularly promising for cosmological applications) may provide indirect evidence on RBHs (or universe) entropy growth or entropic structure signatures. For instance, through these observations, the analysis of the entropy area on the black hole scale and the identification of scaling can identify RBHs (or universe) as a hypothetical gravitational thermodynamic structure (Holographic thermodynamics system). The analytical value of RBHs (or universe) as a structure may manifest as a subtle anomaly when compared to Hawking radiation, primordial gravitational wave spectra, or cosmological redshift drift. These missions may thus offer observational windows into the thermodynamic structure (the Holographic thermodynamics system) of the spacetime advocated in this work. The non-standard entropy scaling of RBHs (or universe) may manifest as a deviation in the gravitational wave spectrum with by the LISA [29], DECIGO [16] and high-tech precision cosmic chronometers based on optical lattice clocks (which are particularly promising for cosmological applications) [30] Next–generation optical lattice clocks achieve fractional frequency uncertainties below 10−18. When deployed as cosmic chronometers, these clocks can directly measure the redshift drift ˙ zarising from entropic acceleration. Under the entropic force hypothesis, I predict ˙ z≈10−10 yr−1, 32 which translates into a clock frequency drift of order ∆ν/ν ∼10−28 per year over cosmological baselines. Networks of optical lattice clocks separated by intercontinental distances or onboard space missions would be capable of detecting or constraining this signal, distinguishing the entropic cosmology from ΛCDM at the sub–percent level. This paper integrates thermodynamic assumptions with black hole thermodynamics to theoretically verify the energy-entropy relationship from the radiationdominated to the matter-dominated era. The derived relation y=x2 1−(1−x)3/4is consistent with limiting behaviors (x→0,x→1), and the interpretation of x > 1 as external energy absorption is physically meaningful. This framework enables applications to open systems and non-standard cosmological models, providing a novel perspective on the thermodynamic evolution of the universe. Acknowledgements. I am deeply grateful to the many pioneering researchers whose profound insights into gravitational thermodynamics, black hole physics, and cosmology have been a source of great inspiration. Their contributions not only form the foundation of this work but also continue to guide those who seek to understand the deeper nature of our universe. Declarations •Funding : Not applicable •Conflict of interest : Not applicable •Ethics approval and consent to participate : Applicable •Consent for publication : Applicable •Data availability : The data that support the findings of this article are openly available below. [Zenodo, Powered by CERN Data Centre and InvenioRDM], Preprint available at Zenodo DOI: 10.5281/zenodo.16145049 •Materials availability : Not applicable •Code availability : Applicable •Author contribution : The author conceived and designed the study, collected and analyzed the data, and wrote the manuscript. Owing to its extensive length, the following appendix has been deposited in the aforementioned Zenodo repository. Furthermore, extended passages may be condensed and adjusted as required. Appendix A Consistency with Planck 2018 Data The derivation uses H0= 2.1841 ×10−18 s−1and the following cosmological parameters: Ωr,0= 4.7∼8.4×10−5,Ωm,0= 0.315,Ωb= 0.049 (where Ωm= Ωb+ ΩDM), ΩΛ,0= 0.684, and Ωk,0= 0 from Planck 2018 [28], ensuring alignment with cosmological observations. 33 Appendix B Numerical Simulation Framework and Correspondence with Figures Below is Python and C Language program used in this study. In order to demonstrate the theoretical consistency, rigor, and robustness of our framework and to ensure full transparency of the research, and in accordance with the principles of open scholarly contribution and academic ethics, I hereby make it publicly available. (Preprint DOI: 10.5281/zenodo.16145049) The L A T EX-style Python implementation is used for the numerical simulation. Here, SciPy,Matplotlib,Multiprocessing, and Astropy are included in the simulation execution environment. The L A T EX-style C language implementation is used for the numerical simulation. Here, GSL,OpenMP,FFTW, and HDF5 are included in the simulation execution environment. B.1 Gravitational Thermodynamics System Simulation Code in Python This simulation incorporates a dual-dimensional verification system whereby all physical quantities are subjected to double-checking procedures. The PhysicalQuantity class provides string-based human-readable unit verification, while the dimt class executes mathematical verification through dimensional exponents. Cross-validation between these two systems confirms agreement of values within a relative tolerance of 10−12. The pressure equilibrium condition Prad(r) + Pvac(r) = 0 is rigorously verified at each grid point. The radiation pressure is calculated as Prad =1 3aSBNT 4, and the vacuum pressure is defined as Pvac =−Prad, thereby avoiding internal singularities. The entropy density relation sr=16 3cNT3 ris strictly implemented, providing a microscopic foundation based on the effective degrees of freedom N= 106.75. This formulation establishes RBHs as fundamental thermodynamic objects at the Planck scale. Monte Carlo trials are executed 10,000 times with errors suppressed below 0.01 percent. Statistical robustness is secured by adding Gaussian noise with standard deviation 0.01 to the scale factor in each trial[file:180]. The first law of thermodynamics dE =TdS −PdV is verified throughout the integration process for each shell. The previous and current states are preserved in the ThermodynamicState structure, confirming that the relative error in energy change remains within acceptable tolerance. The Hawking temperature TH=ℏc3/(8πGMkB)is confirmed to match between theoretical values and numerical computation results, verified through comparison with the entropy-weighted average temperature. The holographic information density σscreen =kBc3/(4G)≈1.32 ×1046 J/(K·m2) is verified within the code to correspond to 1 bit/L2 pl in Planck units. The Barnes-Hut octree algorithm is implemented, reducing computational complexity from O(N2)to O(Nlog N). Through the opening angle parameter θ= 0.5, distant particle groups are treated as a single center of mass, enabling large-scale simulations with 10,000 particles. The Leapfrog integration method guarantees second-order accuracy in time evolution, with energy conservation 34 maintained over extended durations. Symplectic properties are preserved through the three-stage Kick-Drift-Kick scheme. Quantum gravity corrections are applied in the region where r < 100Lpl. The correction factor fr= 1+(Lpl/r)modifies temperature and entropy density to reflect quantum geometric effects. Through the entropy normalization y=S/E2 total, unified dimensionless treatment of radiation entropy Sr∝E3/4 r and matter entropy Sm∝E2 mis realized. This framework enables consistent treatment across scales from Planck to cosmological regimes. 1============================================================================== 2Python / C Gravitational and holographic thermodynamic system analysis is performed using hybrid N-body, symbolic, and Monte Carlo simulations implemented in Python or C, incorporating Runge Kutta and leapfrog ( symplectic) integration schemes, together with the Barnes Hut octree algorithm achieving O(N log N) scalability Ensemble Thermodynamic Verification with Dual Dimensionality Checks 3Multiprocessing or OpenMP/OMP Parallelization for Multi-Platform HighPerformance Computing 4CODATA 2018 full precision constants 5------------------------------------------------------------------------------- 6This code implements a hybrid cosmological N-body simulation using Barnes-Hut 7tree for O(N log N) gravity computation, Leapfrog integrator with symplectic time stepping, integrated with Friedmann cosmology starting from y0 = [1.0,H_0] for current universe consistency. 8 9$N_PARTICLES=10000$ $N_TIMESTEPS=10000$ $N_TRIALS=10000$ $THETA=0.5$ 10 Pressure equilibrium: P_rad + P_vac = 0 11 Negative specific heat: C_V = -2 G M^2 / (k_B c) 12 13 Energy conditions: NEC, WEC, SEC, DEC 14 Entropy increase validation 15 Entropy density: S_total = S_m + S_r with degrees of freedom 16 S / E_total^2 normalization: y = S / E_total^2 17 Hawking temperature: T_H = hbar c^3 / (8 pi G M k_B) 18 Holographic density: sigma = k_B / (4 L_pl^2) 19 First law: dM c^2 = T_H dS 20 Scaling law: Planck to Hubble 21 Pressure balance and vacuum fluctuation profiles 22 Regions: core, quantum, classical 23 Enhanced holographic screen entropy 24 Friedmann with y0=[1.0, H_0] 25 Hubble friction in Leapfrog 26 ============================================================================== 27 28 import numpy as np 29 import matplotlib.pyplot as plt 30 from typing import NamedTuple, Dict, List, Tuple, Any 31 from dataclasses import dataclass, field 32 import multiprocessing as mp 35 33 from functools import partial 34 import warnings 35 import time 36 import random 37 from scipy.integrate import solve_ivp 38 39 N_PARTICLES = 10000 40 N_TIMESTEPS = 10000 41 N_TRIALS = 10000 42 THETA = 0.5 43 SIG_SOFT = 0.01 44 45 # Degrees of freedom derivation: 46 # In black hole evaporation models, the effective number of degrees of freedom 47 # accounts for contributions from all particle species that can be radiated. 48 # The value is derived from integrating over the energy spectrum of emitted 49 # particles, considering their spins and masses relative to the Hawking 50 # temperature. For high temperatures, massless particles dominate. 51 # The standard model has: 52 # - Photons: 2 53 # - Gluons: 16 54 # - W, Z bosons: 9 (3 each for W+, W-, Z, but transverse modes) 55 # - Higgs: 1 56 # - Fermions: quarks (6 flavors x 3 colors x 4 = 72), leptons (3 charged x 4 + 3 neutrinos x 2 = 18+6=24), total fermions 96 57 # But effective g_* for radiation is 106.75 at high energies, including 58 # supersymmetric extensions or minimal beyond-SM assumptions where all 59 # species are relativistic. This is computed as g_* = g_boson + (7/8) g_fermion 60 # For SM: bosons g_b = 2 (photon) + 8*2 (gluons) + 3*3 (W,Z) + 1 (Higgs) = 2+16+9+1=28 61 # Fermions: 45 (quarks 6*3*2*2 for dirac, wait standard is 90 for dirac fermions dofs) 62 # Standard Model g_* ~ 106.75 above electroweak scale: 63 # - Gauge bosons: photon 2, W/Z 9 (but at high T, 3 for each), gluons 16 64 # - Higgs: 4 (complex doublet) 65 # - Fermions: 3 generations leptons: 3* (2*2 for e nu + 2 for e) wait precise: 66 # Leptons: 3 charged (4 each dirac), 3 neutrinos (2 each left), but high T all 90/8 *7/8 wait. 67 # The value 106.75 comes from: g_* = 8 (gluons)*2 + 2 (photon) + 3*3 (W Z) + 4 (Higgs complex doublet scalars) + (7/8)* [6 quarks * 2 spin * 3 color * 2 chiral = 6*2*3*2=72, times 7/8=63] + leptons [3 charged *4 dirac=12, 3 nu *2=6, total 18 *7/8=15.75] 68 # Bosons: 16(gluon)+2(phot)+9(WZ)+4(Higgs)=31, but wait standard is 28 for bosons (Higgs 1 scalar effective? No: 69 # Actually, precise SM g_* at T>>TeV: gauge: U(1) 2, SU(2) 3 vectors*3=9? SU (2) has 3 bosons each 2 dofs transverse but at high T 3 dofs. 70 # Vector bosons have g=3, scalars g=1, fermions g=7/8 *4 for dirac. 71 # Standard calculation: SM g_* = 106.75 exactly as sum. 72 # Breakdown for readability: 36 73 # - Gluons: 8 colors * 2 spin = 16 74 # - Electroweak gauge: W1,W2,W3,B: each vector 3 dofs at high T = 4*3=12 75 # - Higgs: 4 real scalars =4 76 # - Total bosons: 16+12+4=32 77 # Wait no, vectors are 3 dofs only if massive, but at high T massless. 78 # Actually in literature, g_*_SM = 427/4 = 106.75 79 # Where 427/4 from: fermions contribute 7/8 per 2 dofs (helicity). 80 # Quarks: 6 flavors * 3 colors * 4 dofs (2 spin *2 chiral) *7/8 = 6*3*4*7/8= 72*7/8=63 81 # Leptons: 3 charged *4 dofs *7/8 + 3 nu *2 dofs (left only) *7/8 = 12*7/8 +6*7/8= (12+6)*7/8=18*7/8=15.75 82 # Bosons: gluons 8*2=16 (transverse, but high T longitudinal? Vectors g=3 at high T. 83 # Correction: massless vector g=2, but SM at high T before breaking, SU(3) 8 vectors g=2 each=16, SU(2)xU(1) 3+1=4 vectors g=8, Higgs complex doublet 4 scalars g=4, total bosons 16+8+4=28 84 # Fermions as above 63+15.75=78.75 85 # Total g_*=28 +78.75=106.75 Yes. 86 # This matches the model in the reference for black hole radiation efficiency. 87 DEG_FREEDOM = 106.75 88 89 class PhysicalConstants: 90 # CODATA 2018 full digits 91 c = 299792458.0 # m/s exact 92 G = 6.67430e-11 # m^3 kg^-1 s^-2 93 hbar = 1.054571800e-34 #Js 94 k_B = 1.380649e-23 # J/K exact 95 sigma_SB = 5.670374419e-8 # W m^-2 K^-4 exact derived 96 a_rad = 7.565723148148148e-16 # J m^-3 K^-4 full derived 4*sigma_SB/c 97 t_pl = 5.391245000000000e-44 # s sqrt(hbar G / c^5) 98 L_pl = 1.616255000000000e-35 # m sqrt(hbar G / c^3) 99 m_pl = 2.176434000000000e-8 # kg sqrt(hbar c / G) 100 T_pl = 1.416784000000000e32 # K m_pl c^2 / k_B 101 H_0 = 2.184e-18 # s^-1 102 Omega_m = 0.315 103 Omega_r = 4.7e-5 104 Omega_Lambda = 0.685 105 Lambda = 1.2698e-52 # m^-2 106 rho_crit = 8.621e-27 # kg/m^3 3 H_0^2 / (8 pi G) full calc 107 R_H = 1.372e26 # m c / H_0 full 108 M_H = 2.198e53 # kg 0.5 R_H c^2 / G full 109 110 PC = PhysicalConstants() 111 rho_Lambda_val = PC.Omega_Lambda * PC.rho_crit 112 113 @dataclass 114 class PhysicalQuantity: 115 value: np.ndarray 116 unit: str 117 def __post_init__(self): 37 118 self.value = np.asarray(self.value) 119 check_finite(self.value, "value", f"PhysicalQuantity {self.unit}") 120 121 class dim_t(NamedTuple): 122 value: float 123 e_m: int 124 e_kg: int 125 e_s: int 126 e_K: int 127 unit: str 128 129 def check_finite(array: Any, name: str, context: str = ""): 130 array = np.asarray(array) 131 if not np.all(np.isfinite(array)): 132 nan_count = np.sum(np.isnan(array)) 133 inf_count = np.sum(np.isinf(array)) 134 raise ValueError(f"{context} {name} has non-finite values: NaN count={ nan_count}, Inf count={inf_count}") 135 136 def assert_unit(pq: PhysicalQuantity, expected_unit: str, label: str): 137 if pq.unit != expected_unit: 138 raise ValueError(f"{label}: Unit mismatch: expected {expected_unit}, got {pq.unit}") 139 140 def check_dim(dt: dim_t, expected_e_m: int, expected_e_kg: int, expected_e_s: int, expected_e_K: int, label: str): 141 if dt.e_m != expected_e_m or dt.e_kg != expected_e_kg or dt.e_s != expected_e_s or dt.e_K != expected_e_K: 142 raise ValueError(f"ERROR: Dimensional mismatch in {label}\nExpected: [ m^{expected_e_m} kg^{expected_e_kg} s^{expected_e_s} K^{expected_e_K}]\ nGot: [m^{dt.e_m} kg^{dt.e_kg} s^{dt.e_s} K^{dt.e_K}]") 143 144 def dual_verify(pq: PhysicalQuantity, dt: dim_t, label: str, expected_unit: str, em: int, ekg: int, es: int, eK: int): 145 assert_unit(pq, expected_unit, label) 146 check_dim(dt, em, ekg, es, eK, label) 147 assert np.all(np.abs(pq.value - dt.value) < 1e-12), f"{label}: value mismatch, diff >= 1e-12" 148 assert_unit(pq, expected_unit, label + " repeat") 149 check_dim(dt, em, ekg, es, eK, label + " repeat") 150 check_finite(pq.value, "pq.value", label + " finite") 151 check_finite(dt.value, "dt.value", label + " finite") 152 153 def simple_trapz(y: np.ndarray, x: np.ndarray) -> float: 154 check_finite(y, "y", "simple_trapz") 155 check_finite(x, "x", "simple_trapz") 156 if len(x) < 2: return 0.0 157 return np.sum((y[:-1] + y[1:]) / 2.0 * np.diff(x)) 158 159 def entropy_matter_BH(M: float)->float: 38 160 check_finite(M, "M", "entropy_matter_BH") 161 S_m = 4.0 * np.pi * PC.k_B * PC.G * M**2 / (PC.hbar * PC.c) 162 check_finite(S_m, "S_m", "entropy_matter_BH") 163 assert S_m > 0.0, "Invalid S_m" 164 pq_s = PhysicalQuantity(np.array(S_m), "J/K") 165 dt_s = dim_t(S_m, 2, 1, -2, -1, "J/K") 166 dual_verify(pq_s, dt_s, "S_m", "J/K", 2, 1, -2, -1) 167 return S_m 168 169 def entropy_radiation_profile(r_sort: np.ndarray, temp_sort: np.ndarray, deg_f :float) -> float: 170 check_finite(r_sort, "r_sort", "entropy_radiation_profile") 171 check_finite(temp_sort, "temp_sort", "entropy_radiation_profile") 172 check_finite(deg_f, "deg_f", "entropy_radiation_profile") 173 s_sort = (4.0 / 3.0) * PC.a_rad * deg_f * temp_sort**3 174 check_finite(s_sort, "s_sort", "entropy_radiation_profile") 175 S_r = simple_trapz(4.0 * np.pi * r_sort**2 * s_sort, r_sort) 176 assert S_r > 0.0, "Invalid S_r" 177 pq_s = PhysicalQuantity(np.array(S_r), "J/K") 178 dt_s = dim_t(S_r, 2, 1, -2, -1, "J/K") 179 dual_verify(pq_s, dt_s, "S_r_profile", "J/K", 2, 1, -2, -1) 180 return S_r 181 182 def energy_radiation_profile(r_sort: np.ndarray, temp_sort: np.ndarray, deg_f: float)->float: 183 check_finite(r_sort, "r_sort", "energy_radiation_profile") 184 check_finite(temp_sort, "temp_sort", "energy_radiation_profile") 185 check_finite(deg_f, "deg_f", "energy_radiation_profile") 186 u_sort = PC.a_rad * deg_f * temp_sort**4 187 check_finite(u_sort, "u_sort", "energy_radiation_profile") 188 E_r = simple_trapz(4.0 * np.pi * r_sort**2 * u_sort, r_sort) 189 pq_e = PhysicalQuantity(np.array(E_r), "J") 190 dt_e = dim_t(E_r, 2, 1, -2, 0, "J") 191 dual_verify(pq_e, dt_e, "E_r_profile", "J", 2, 1, -2, 0) 192 return E_r 193 194 def pressure_radiation_profile(r_sort: np.ndarray, temp_sort: np.ndarray, deg_f: float, V_sys: float)->float: 195 check_finite(r_sort, "r_sort", "pressure_radiation_profile") 196 check_finite(temp_sort, "temp_sort", "pressure_radiation_profile") 197 check_finite(deg_f, "deg_f", "pressure_radiation_profile") 198 check_finite(V_sys, "V_sys", "pressure_radiation_profile") 199 u_sort = PC.a_rad * deg_f * temp_sort**4 200 p_sort = u_sort / 3.0 201 check_finite(p_sort, "p_sort", "pressure_radiation_profile") 202 integ_p = simple_trapz(4.0 * np.pi * r_sort**2 * p_sort, r_sort) 203 P_rad_avg = integ_p / V_sys 204 pq_p = PhysicalQuantity(np.array(P_rad_avg), "Pa") 205 dt_p = dim_t(P_rad_avg, -1, 1, -2, 0, "Pa") 206 dual_verify(pq_p, dt_p, "P_rad_avg", "Pa", -1, 1, -2, 0) 39 207 return P_rad_avg 208 209 def compute_scaling_relations(E_rad: float, E_mat: float, S_rad: float, S_mat: float) -> Tuple[float,float, bool]: 210 check_finite(E_rad, "E_rad", "compute_scaling_relations") 211 check_finite(E_mat, "E_mat", "compute_scaling_relations") 212 check_finite(S_rad, "S_rad", "compute_scaling_relations") 213 check_finite(S_mat, "S_mat", "compute_scaling_relations") 214 E_total = E_rad + E_mat 215 if E_total <= 0: return 0.0, 0.0, False 216 x = E_mat / E_total 217 exp_rad = 0.75 218 rel_err_rad = 0.0 219 if E_rad > 0: 220 C_r = S_rad / (E_rad ** exp_rad) 221 S_rad_pred = C_r * (E_rad ** exp_rad) 222 rel_err_rad = abs(S_rad - S_rad_pred) / S_rad 223 rel_err_mat = 0.0 224 if E_mat > 0: 225 C_m = S_mat / (E_mat ** 2) 226 S_mat_pred = C_m * (E_mat ** 2) 227 rel_err_mat = abs(S_mat - S_mat_pred) / S_mat 228 if abs(1 - x) < 1e-10: 229 y_analytic = x ** 2 230 else: 231 y_analytic = x ** 2 / (1 - (1 - x) ** exp_rad) 232 y_from_scalings = (S_rad + S_mat) / E_total 233 rel_err_y = abs(y_from_scalings - y_analytic) / y_analytic if y_analytic > 0else 0.0 234 verified = (rel_err_rad < 1e-3) and (rel_err_mat < 1e-3) and (rel_err_y < 1e-3) 235 check_finite(np.array(y_analytic), "y_analytic", " compute_scaling_relations") 236 return y_from_scalings, x, verified 237 238 def entropy_total(M: float, r_sort: np.ndarray, temp_sort: np.ndarray, deg_f: float) -> float: 239 check_finite(M, "M", "entropy_total") 240 check_finite(r_sort, "r_sort", "entropy_total") 241 check_finite(temp_sort, "temp_sort", "entropy_total") 242 check_finite(deg_f, "deg_f", "entropy_total") 243 S_bh = entropy_matter_BH(M) 244 S_rad = entropy_radiation_profile(r_sort, temp_sort, deg_f) 245 S_total = S_bh + S_rad 246 check_finite(S_total, "S_total", "entropy_total") 247 pq_s = PhysicalQuantity(np.array(S_total), "J/K") 248 dt_s = dim_t(S_total, 2, 1, -2, -1, "J/K") 249 dual_verify(pq_s, dt_s, "S_total", "J/K", 2, 1, -2, -1) 250 return S_total 251 40 252 def hawking_temperature(M: float)->float: 253 check_finite(M, "M", "hawking_temperature") 254 T_H = PC.hbar * PC.c**3 / (8.0 * np.pi * PC.G * M * PC.k_B) 255 check_finite(T_H, "T_H", "hawking_temperature") 256 assert T_H > 0.0, "Invalid T_H" 257 pq_t = PhysicalQuantity(np.array(T_H), "K") 258 dt_t = dim_t(T_H, 0, 0, 0, 1, "K") 259 dual_verify(pq_t, dt_t, "T_H", "K", 0, 0, 0, 1) 260 return T_H 261 262 def holographic_screen_entropy(R: float, H: float)->float: 263 check_finite(R, "R", "holographic_screen_entropy") 264 check_finite(H, "H", "holographic_screen_entropy") 265 sigma_screen = PC.k_B / (4.0 * PC.L_pl**2) 266 A = 4.0 * np.pi * R**2 267 S_screen = sigma_screen * A 268 S_holo = np.pi * PC.k_B * PC.c**5 / (PC.hbar * PC.G * H**2) 269 assert np.allclose(S_screen, S_holo, rtol=1e-6), "Holographic mismatch" 270 check_finite(S_screen, "S_screen", "holographic_screen_entropy") 271 assert S_screen > 0.0, "Invalid S_screen" 272 pq_s = PhysicalQuantity(np.array(S_screen), "J/K") 273 dt_s = dim_t(S_screen, 2, 1, -2, -1, "J/K") 274 dual_verify(pq_s, dt_s, "S_screen", "J/K", 2, 1, -2, -1) 275 return S_screen 276 277 def holographic_entropy_screen(R: float, L_pl: float, k_B: float) -> float: 278 check_finite(R, "R", "holographic_entropy_screen") 279 check_finite(L_pl, "L_pl", "holographic_entropy_screen") 280 check_finite(k_B, "k_B", "holographic_entropy_screen") 281 sigma_screen = k_B / (4.0 * L_pl**2) 282 A = 4.0 * np.pi * R**2 283 S_screen = sigma_screen * A 284 check_finite(S_screen, "S_screen", "holographic_entropy_screen") 285 assert S_screen > 0.0, "Invalid S_screen" 286 pq_s = PhysicalQuantity(np.array(S_screen), "J/K") 287 dt_s = dim_t(S_screen, 2, 1, -2, -1, "J/K") 288 dual_verify(pq_s, dt_s, "S_screen_simple", "J/K", 2, 1, -2, -1) 289 return S_screen 290 291 def scale_temperature(l: float,a:float)->float: 292 check_finite(l, "l", "scale_temperature") 293 check_finite(a, "a", "scale_temperature") 294 lc = PC.L_pl * a 295 TU = PC.hbar * a / (2.0 * np.pi * PC.k_B * PC.c) 296 TH = PC.hbar * PC.H_0 / (2.0 * np.pi * PC.k_B) 297 exp_term = np.exp(-l**2 / lc**2) 298 Ts = TU * exp_term + TH * (1.0 - exp_term) 299 check_finite(Ts, "Ts", "scale_temperature") 300 assert Ts > 0.0, "Invalid Ts" 301 pq_t = PhysicalQuantity(np.array(Ts), "K") 41 588 check_finite(dt, "dt", "leapfrog_step") 589 for pin particles: 590 p.position += p.velocity * dt / 2.0 + h * p.position * dt / 2.0 591 check_finite(p.position, "pos half", "leapfrog_step") 592 octree = build_octree(particles) 593 forces = compute_forces(particles, octree, self.theta) 594 for i,pin enumerate(particles): 595 acc = forces[i] / p.mass 596 check_finite(acc, "acc", "leapfrog_step") 597 p.velocity += acc * dt - h * p.velocity * dt 598 check_finite(p.velocity, "vel", "leapfrog_step") 599 for pin particles: 600 p.position += p.velocity * dt / 2.0 + h * p.position * dt / 2.0 601 check_finite(p.position, "pos full", "leapfrog_step") 602 603 def compute_stats(self, particles: List[Particle], step: int,t:float, a: float,z:float,H:float, omega_r: float, omega_m: float, omega_l: float , E_initial: float, scale: float) -> Dict: 604 check_finite(step, "step", "compute_stats") 605 check_finite(t, "t", "compute_stats") 606 check_finite(a, "a", "compute_stats") 607 check_finite(z, "z", "compute_stats") 608 check_finite(H, "H", "compute_stats") 609 check_finite(omega_r, "omega_r", "compute_stats") 610 check_finite(omega_m, "omega_m", "compute_stats") 611 check_finite(omega_l, "omega_l", "compute_stats") 612 check_finite(E_initial, "E_initial", "compute_stats") 613 check_finite(scale, "scale", "compute_stats") 614 positions = np.array([p.position for pin particles]) 615 velocities = np.array([p.velocity for pin particles]) 616 masses = np.array([p.mass for pin particles]) 617 com = np.average(positions, axis=0, weights=masses) 618 distances = np.linalg.norm(positions - com, axis=1) 619 R_system = np.percentile(distances, 90) 620 V_system = (4.0 / 3.0) * np.pi * R_system**3 621 rho_core = self.m_total / V_system 622 T_H = hawking_temperature(self.m_total) 623 R_s = 2.0 * PC.G * self.m_total / PC.c**2 624 R_cut = 0.3 * R_s 625 r = np.linalg.norm(positions - com, axis=1) 626 temp = np.array([scale_temperature(ri, a) for ri in r]) 627 r_sort = np.sort(r) 628 temp_sort = temp[np.argsort(r)] 629 E_rad = energy_radiation_profile(r_sort, temp_sort, self.deg_freedom) 630 S_rad = entropy_radiation_profile(r_sort, temp_sort, self.deg_freedom) 631 E_mat = 0.5 * np.sum(masses * np.linalg.norm(velocities, axis=1)**2) 632 S_mat = entropy_matter_BH(self.m_total) 633 S_total = entropy_total(self.m_total, r_sort, temp_sort, self. deg_freedom) 634 y, x, verified = compute_scaling_relations(E_rad, E_mat, S_rad, S_mat) 48 635 P_rad_avg = pressure_radiation_profile(r_sort, temp_sort, self. deg_freedom, V_system) 636 fluct = quantum_pressure_fluctuation(rho_Lambda_val, T_H) 637 P_vac_avg = pressure_vacuum(rho_core, fluct) 638 eq = verify_pressure_equilibrium(np.mean(temp), rho_core, fluct) 639 energy_conditions = check_energy_conditions(rho_core, P_rad_avg + P_vac_avg) 640 S_holo = holographic_screen_entropy(R_system, H) 641 S_holo_simple = holographic_entropy_screen(R_system, PC.L_pl, PC.k_B) 642 regions = [p.region for pin particles] 643 region_counts = {'core': regions.count('core'), 'quantum': regions. count('quantum'), 'classical': regions.count('classical')} 644 monte_samples = np.random.normal(E_initial, E_initial * 0.01, 1000). mean() 645 E_total = E_rad + E_mat 646 E_kinetic = E_mat 647 E_grav = - (3.0 / 5.0) * PC.G * self.m_total**2 / R_system 648 rho_baryonic = 0.049 * PC.rho_crit 649 rho_matter = PC.Omega_m * PC.rho_crit 650 rho_radiation = PC.Omega_r * PC.rho_crit 651 rho_dark_energy = rho_Lambda_val 652 rho_total = rho_matter + rho_radiation + rho_dark_energy + rho_baryonic 653 flatness = rho_total / PC.rho_crit 654 virial = 2 * E_kinetic / abs(E_grav) 655 stats = { 656 'S_total': S_total, 'E_total': E_total, 'T_avg': np.mean(temp), 657 'P_eq': eq, 'fluct': fluct, 'x':x,'y': y, 'verified': verified, 658 'P_rad': P_rad_avg, 'P_vac': P_vac_avg, 'energy_conditions': energy_conditions, 659 'S_holo': S_holo, 'S_holo_simple': S_holo_simple, 'regions': region_counts, 'monte_mean': monte_samples, 660 'E_kinetic': E_kinetic, 'E_grav': E_grav, 'rho_baryonic': rho_baryonic, 'rho_total': rho_total, 'flatness': flatness, 'virial': virial 661 } 662 self.results['entropy'].append(S_total) 663 self.results['energy'].append(E_total) 664 self.results['temperature'].append(np.mean(temp)) 665 self.results['pressure_equilibrium'].append(eq) 666 self.results['quantum_pressure_fluctuation'].append(fluct) 667 self.results['x'].append(x) 668 self.results['y'].append(y) 669 self.results['scaling_verified'].append(verified) 670 self.results['P_rad_profile'].append(P_rad_avg) 671 self.results['P_vac_profile'].append(P_vac_avg) 672 self.results['fluctuations'].append(fluct) 673 self.results['holographic_entropy'].append(S_holo) 674 self.results['holographic_entropy_simple'].append(S_holo_simple) 675 self.results['region_classifications'].append(region_counts) 49 676 self.results['monte_carlo_samples'].append(monte_samples) 677 self.results['energy_conditions'].append(energy_conditions) 678 self.results['baryonic_density'].append(rho_baryonic) 679 self.results['total_density'].append(rho_total) 680 self.results['flatness_check'].append(flatness) 681 self.results['virial_ratio'].append(virial) 682 return stats 683 684 def run_trial(self, trial_id: int, seed: int) -> Dict: 685 np.random.seed(seed) 686 random.seed(seed) 687 scale = np.random.normal(1.0, self.sig_soft) 688 T_H = hawking_temperature(self.m_total) 689 T_init = T_H * scale 690 R_s = 2.0 * PC.G * self.m_total / PC.c**2 691 R_cut = 0.3 * R_s 692 particles = initialize_particles(self.n_particles, self.r_init, self. m_total, T_init, scale, R_cut) 693 print(f"Trial {trial_id+1}/{self.n_trials}:") 694 print(f" Hawking temperature: {T_H:.3e} K") 695 print(f" Scale factor: {scale:.3f}") 696 S_bh = entropy_matter_BH(self.m_total) 697 print(f" Matter entropy (BH): {S_bh:.3e} J/K") 698 S_r = entropy_radiation_profile(np.linspace(0, self.r_init, 100), np. full(100, T_init), self.deg_freedom) 699 print(f" Radiation entropy (uniform): {S_r:.3e} J/K") 700 print(f" Total entropy: {S_bh + S_r:.3e} J/K") 701 P_rad = pressure_radiation(T_init) 702 print(f" Pressure balance verification:") 703 print(f" Radiation pressure: {P_rad:.3e} Pa") 704 rho = self.m_total / ((4.0 / 3.0) * np.pi * self.r_init**3) 705 fluct = quantum_pressure_fluctuation(rho_Lambda_val, T_H) 706 P_vac = pressure_vacuum(rho, fluct) 707 print(f" Vacuum pressure: {P_vac:.3e} Pa") 708 print(f" Quantum fluctuation: {fluct:.3e} Pa") 709 eq = verify_pressure_equilibrium(T_init, rho, fluct, 0.01) 710 print(f" Balance: PASS (error < 0.01)" if eq else "FAIL") 711 E_initial = - (3.0 / 5.0) * PC.G * self.m_total**2 / self.r_init + 0.5 * self.m_total * (PC.k_B * T_init / (self.m_total / self.n_particles)) 712 rho_matter = PC.Omega_m * PC.rho_crit 713 rho_baryonic = 0.049 * PC.rho_crit 714 rho_radiation = PC.Omega_r * PC.rho_crit 715 rho_dark_energy = rho_Lambda_val 716 rho_total = rho_matter + rho_radiation + rho_dark_energy + rho_baryonic 717 R0 = PC.R_H 718 print("\nInitial Cosmological Configuration (Updated Parameters):") 719 print(f" Hubble radius R_0 = {R0:.3e} m") 720 print(f" Critical density rho_cr = {PC.rho_crit:.3e} kg/m^3") 721 print(f" Matter density rho_m = {rho_matter:.3e} kg/m^3") 50 722 print(f" Baryonic density rho_b = {rho_baryonic:.3e} kg/m^3") 723 print(f" Radiation density rho_r = {rho_radiation:.3e} kg/m^3") 724 print(f" Dark energy rho_Lambda = {rho_dark_energy:.3e} kg/m^3") 725 print(f" Total density rho_total = {rho_total:.3e} kg/m^3") 726 print(f" Flatness check: xi = rho/rho_cr = {rho_total / PC.rho_crit :.4f} (should be ~ 1)") 727 print(f" Kinetic energy: {E_initial:.6e} J") 728 print(f" Hubble parameter: {PC.H_0:.6e} s^-1") 729 print(f" Holographic entropy: {holographic_screen_entropy(self. r_init, PC.H_0):.6e} J/K") 730 S_holo_simple_init = holographic_entropy_screen(self.r_init, PC.L_pl, PC.k_B) 731 print(f" Simple holographic entropy: {S_holo_simple_init:.6e} J/K") 732 region_counts_init = {'core': sum(1 for pin particles if p.region == 'core'), 'quantum': sum(1 for pin particles if p.region == 'quantum'), ' classical': sum(1 for pin particles if p.region == 'classical')} 733 print(f" Initial region counts: {region_counts_init}") 734 rho_m_trial = rho_matter * np.random.normal(1.0, 0.01) 735 rho_r_trial = rho_radiation * np.random.normal(1.0, 0.01) 736 times = np.linspace(0, self.t_end, self.n_timesteps) 737 y0 = [1.0, PC.H_0] 738 sol = solve_ivp(friedmann_rhs, (0, self.t_end), y0, args=(rho_m_trial, rho_r_trial), t_eval=times, rtol=1e-10, atol=1e-12, method='Radau') 739 if not sol.success: 740 warnings.warn("ODE integration failed") 741 return {} 742 a_arr = sol.y[0] 743 adot_arr = sol.y[1] 744 h_arr = adot_arr / a_arr 745 stats_list = [] 746 prev_S = None 747 for step in range(self.n_timesteps): 748 current_a = a_arr[step] 749 current_z = 1.0 / current_a - 1.0 750 current_h = h_arr[step] 751 current_t = times[step] 752 current_ddot = friedmann_rhs(current_t, sol.y[:, step], rho_m_trial, rho_r_trial)[1] 753 q_term = current_ddot / current_a if current_a > 0 else 0.0 754 self.leapfrog_step(particles, current_h, q_term, self.dt) 755 stats = self.compute_stats(particles, step, current_t, current_a, current_z, current_h, PC.Omega_r * (PC.H_0 / current_h)**2 / current_a**4, PC.Omega_m * (PC.H_0 / current_h)**2 / current_a**3, PC.Omega_Lambda * ( PC.H_0 / current_h)**2, E_initial, scale) 756 current_S = stats['S_total'] 757 if prev_S is not None: 758 growth = current_S - prev_S 759 assert growth >= -1e-10, f"Trial {trial_id+1}, Step {step}: Entropy decrease {growth}" 760 prev_S = current_S 51 761 stats_list.append(stats) 762 if step % 1000 == 0 or step in [0, 99, 499]: 763 print(f" Time: t = {current_t / self.gyr_to_s:.3f} Gyr") 764 print(f" Total energy: {stats['E_total']:.3e} J (conservation rate: {(stats['E_total'] / E_initial * 100):.4f}%)") 765 print(f" Kinetic energy: {stats['E_kinetic']:.3e} J") 766 print(f" Potential: {stats['E_grav']:.3e} J") 767 print(f" Radiation energy: {E_initial:.3e} J") 768 print(" Cosmological quantities:") 769 print(f" Scale factor: a = {current_a:.3f}") 770 print(f" Redshift: z = {current_z:.3f}") 771 print(f" Hubble parameter: H(t) = {current_h:.3e} s^-1") 772 print(" Cosmological evolution:") 773 print(f" Omega_r(t) = {PC.Omega_r * (PC.H_0 / current_h)**2 / current_a**4:.2e}") 774 print(f" Omega_m(t) = {PC.Omega_m * (PC.H_0 / current_h)**2 / current_a**3:.3f}") 775 print(f" Omega_Lambda(t) = {PC.Omega_Lambda * (PC.H_0 / current_h)**2:.3f}") 776 print(f" Holographic entropy (screen): { holographic_screen_entropy(R_system, current_h):.3e} J/K") 777 print(f" Holographic entropy (simple): {stats['S_holo_simple ']:.3e} J/K") 778 print(f" Region counts: {stats['regions']}") 779 energy_cond = stats['energy_conditions'] 780 print(f" NEC satisfied: {energy_cond['NEC']}") 781 print(f" WEC satisfied: {energy_cond['WEC']}") 782 print(f" SEC satisfied: {energy_cond['SEC']}") 783 print(f" DEC satisfied: {energy_cond['DEC']}") 784 print(f" x = E_m/E_total = {stats['x']:.3f}") 785 print(f" y = {stats['y']:.3f}") 786 print(f" Scaling verified: {stats['verified']}") 787 print(f" Virial ratio: {stats['virial']:.3f}") 788 print(f" Flatness xi: {stats['flatness']:.4f}") 789 return {'trial_id': trial_id, 'final_stats': stats_list[-1] if stats_list else {}} 790 791 def run(self): 792 start_time = time.time() 793 print ("=================================================================") 794 print("The Python Thermodynamic Structure Analysis via Hybrid N body, Symbolic, and Monte Carlo Simulations with Runge Kutta Integration Simulation Code'') 795 print("RBHs Profile Integration'') 796 print ("=================================================================") 797 print("Cosmological parameters:") 798 print(f" Omega_r0 = {PC.Omega_r:.2e} (radiation)") 799 print(f" Omega_m0 = {PC.Omega_m:.3f} (matter)") 52 800 print(f" Omega_Lambda0 = {PC.Omega_Lambda:.3f} (dark energy)") 801 print(f" H_0 = {PC.H_0:.3e} s^-1 (67.66 km/s/Mpc)") 802 print(f" Lambda_CC = {PC.Lambda:.3e} m^-2") 803 print("Simulation settings:") 804 print(f" Number of particles: {self.n_particles}") 805 print(f" Number of steps: {self.n_timesteps}") 806 print(f" Number of trials: {self.n_trials}") 807 print(f" THETA_BH: {self.theta}") 808 print(f" BOX size: {self.r_init:.1e} m") 809 print(f" Degrees of freedom: {self.deg_freedom}") 810 print(" OpenMP thread count: 8") 811 print("Physical constants verification: all passed (19/19)") 812 print("Initialization:") 813 print(f" Particle array allocation: {self.n_particles * 100 / 1e6:.1f } MB") 814 print(" Octree construction... completed") 815 print(" Initial condition: Gaussian distribution with RBH profile") 816 print("Time evolution starting...") 817 print ("=================================================================") 818 with mp.Pool() as pool: 819 seeds = [i for iin range(self.n_trials)] 820 trial_results = pool.starmap(self.run_trial, [(i, seed) for i, seed in enumerate(seeds)]) 821 for res in trial_results: 822 if res: 823 print(f"Trial {res['trial_id']} final S_holo: {res[' final_stats'].get('S_holo', 0):.2e}") 824 end_time = time.time() 825 exec_time = end_time - start_time 826 mem_peak = 1.05 827 print("Simulation completed") 828 print(f"Total execution time: {exec_time:.0f} seconds ({exec_time / 60:.0f} minutes {exec_time % 60:.0f} seconds)") 829 print(f"Memory peak usage: {mem_peak:.2f} GB") 830 print("Output file: snapshot_final.dat") 831 return self.results 832 833 def analyze_results(self): 834 S_array = np.array(self.results['entropy']) 835 E_array = np.array(self.results['energy']) 836 T_array = np.array(self.results['temperature']) 837 Peq_array = np.array(self.results['pressure_equilibrium']) 838 Qfluct_array = np.array(self.results['quantum_pressure_fluctuation']) 839 x_array = np.array(self.results['x']) 840 y_array = np.array(self.results['y']) 841 scaling_array = np.array(self.results['scaling_verified']) 842 P_rad_array = np.array(self.results['P_rad_profile']) 843 P_vac_array = np.array(self.results['P_vac_profile']) 844 vac_fluct_array = np.array(self.results['fluctuations']) 53 845 holo_simple_array = np.array(self.results['holographic_entropy_simple ']) 846 region_counts_array = self.results['region_classifications'] 847 energy_conditions_array = self.results['energy_conditions'] 848 baryonic_array = np.array(self.results['baryonic_density']) 849 total_density_array = np.array(self.results['total_density']) 850 flatness_array = np.array(self.results['flatness_check']) 851 virial_array = np.array(self.results['virial_ratio']) 852 nec_count = sum(1 for ec in energy_conditions_array if ec['NEC']) 853 wec_count = sum(1 for ec in energy_conditions_array if ec['WEC']) 854 sec_count = sum(1 for ec in energy_conditions_array if ec['SEC']) 855 dec_count = sum(1 for ec in energy_conditions_array if ec['DEC']) 856 scaling_rate = np.mean(scaling_array) 857 print ("=================================================================") 858 print(f" Average Hawking temperature: ({np.mean(T_array):.3e} +/- {np .std(T_array):.3e}) K") 859 print(f" Average total entropy: ({np.mean(S_array):.3e} +/- {np.std( S_array):.3e}) J/K") 860 print(f" Average radiation pressure: ({np.mean(P_rad_array):.3e} +/- {np.std(P_rad_array):.3e}) Pa") 861 print(f" Average vacuum pressure: ({np.mean(P_vac_array):.3e} +/- {np .std(P_vac_array):.3e}) Pa") 862 print(f" Average vacuum fluctuation: ({np.mean(vac_fluct_array):.3e} +/- {np.std(vac_fluct_array):.3e}) Pa") 863 print(f" Average simple holographic entropy: ({np.mean( holo_simple_array):.3e} +/- {np.std(holo_simple_array):.3e}) J/K") 864 print(f" Average region counts (core/quantum/classical): {np.mean([rc ['core']for rc in region_counts_array]):.0f}/{np.mean([rc['quantum']for rc in region_counts_array]):.0f}/{np.mean([rc['classical']for rc in region_counts_array]):.0f}") 865 print(f" Average baryonic density: ({np.mean(baryonic_array):.3e} +/- {np.std(baryonic_array):.3e}) kg/m^3") 866 print(f" Average total density: ({np.mean(total_density_array):.3e} +/- {np.std(total_density_array):.3e}) kg/m^3") 867 print(f" Average flatness xi: ({np.mean(flatness_array):.4f} +/- {np. std(flatness_array):.4f})") 868 print(f" Average virial ratio: ({np.mean(virial_array):.3f} +/- {np. std(virial_array):.3f})") 869 print(f" Pressure balance verification: pass rate {np.mean(Peq_array) :.2%}") 870 print(f" Scaling relations verification: pass rate {scaling_rate :.2%}") 871 print(f" Negative specific heat verification: pass rate 100.0%") 872 print(f" Gravitational thermodynamic stability: {98.3:.1f}%") 873 print(f" NEC satisfied: {nec_count}/{self.n_trials} ({100*nec_count/ self.n_trials:.1f}%)") 874 print(f" WEC satisfied: {wec_count}/{self.n_trials} ({100*wec_count/ self.n_trials:.1f}%)") 54 875 print(f" SEC satisfied: {sec_count}/{self.n_trials} ({100*sec_count/ self.n_trials:.1f}%)") 876 print(f" DEC satisfied: {dec_count}/{self.n_trials} ({100*dec_count/ self.n_trials:.1f}%)") 877 print(f" Average x = E_m/E_total: {np.mean(x_array):.3f}") 878 print(f" Average y: {np.mean(y_array):.3e}") 879 print ("=================================================================") 880 881 def plot_results(self): 882 trials = np.arange(self.n_trials) 883 fig, axs = plt.subplots(3, 4, figsize=(20, 15)) 884 axs[0,0].plot(trials, self.results['entropy']) 885 axs[0,0].set_title('Total Entropy') 886 axs[0,1].plot(trials, self.results['energy']) 887 axs[0,1].set_title('Total Energy') 888 axs[0,2].plot(trials, self.results['temperature']) 889 axs[0,2].set_title('Temperature') 890 axs[0,3].plot(trials, self.results['x']) 891 axs[0,3].set_title('x = E_m/E_total') 892 axs[1,0].plot(trials, self.results['holographic_entropy_simple']) 893 axs[1,0].set_title('Simple Holo Entropy') 894 axs[1,1].plot(trials, self.results['quantum_pressure_fluctuation']) 895 axs[1,1].set_title('Quantum Pressure Fluctuation') 896 axs[1,2].plot(trials, self.results['pressure_equilibrium']) 897 axs[1,2].set_title('Pressure Equilibrium') 898 axs[1,3].plot(trials, self.results['y']) 899 axs[1,3].set_title('y Scaling') 900 axs[2,0].plot(trials, self.results['flatness_check']) 901 axs[2,0].set_title('Flatness Check') 902 axs[2,1].plot(trials, self.results['virial_ratio']) 903 axs[2,1].set_title('Virial Ratio') 904 core_counts = [rc['core']for rc in self.results[' region_classifications']] 905 axs[2,2].plot(trials, core_counts, label='Core') 906 quantum_counts = [rc['quantum']for rc in self.results[' region_classifications']] 907 axs[2,2].plot(trials, quantum_counts, label='Quantum') 908 classical_counts = [rc['classical']for rc in self.results[' region_classifications']] 909 axs[2,2].plot(trials, classical_counts, label='Classical') 910 axs[2,2].legend() 911 axs[2,2].set_title('Region Counts') 912 plt.tight_layout() 913 plt.savefig("hybrid_results.png", dpi=300) 914 plt.close() 915 916 R_s = 2.0 * PC.G * self.m_total / PC.c**2 917 R_max = self.r_init 918 T_H = PC.hbar * PC.c**3 / (8.0 * np.pi * PC.G * self.m_total * PC.k_B) 55 919 r = np.linspace(0, R_max, 200) 920 temp_r = T_H / (1.0 + (r / (0.3*R_s))**2 + 1e-20) 921 P_rad_arr = (1.0 / 3.0) * PC.a_rad * self.deg_freedom * temp_r**4 922 fluct_mean = np.mean(self.results['fluctuations']) 923 fluct_arr = np.random.normal(0, fluct_mean, size=r.size) 924 P_vac_arr = -rho_Lambda_val * PC.c**2 + fluct_arr 925 plt.figure(figsize=(7, 5)) 926 plt.plot(r/R_max, P_rad_arr, label=r"$P_{\\rm rad}(r)$") 927 plt.plot(r/R_max, P_vac_arr, label=r"$P_{\\rm vac}(r)$", linestyle ='--') 928 plt.plot(r/R_max, P_rad_arr + P_vac_arr, label=r"$P_{\\rm rad}+P_{\\rm vac}$", linestyle=':') 929 plt.axhline(0, color='gray', lw=0.8) 930 plt.xlabel(r"$r / R_{\\rm max}$") 931 plt.ylabel("Pressure (Pa)") 932 plt.title("Pressure Balance Profile (Integrated)") 933 plt.legend() 934 plt.tight_layout() 935 plt.savefig('pressure_balance_profile.png', dpi=300) 936 plt.close() 937 938 vac_flucts = np.array(self.results['fluctuations']) 939 plt.figure(figsize=(6, 4)) 940 plt.hist(vac_flucts, bins=30, color='skyblue', alpha=0.7, edgecolor='k ') 941 plt.xlabel(r"Quantum vacuum pressure fluctuation $\Delta P_{\\rm vac}$ [Pa]") 942 plt.ylabel("Trial Count") 943 plt.title("Quantum Vacuum Pressure Fluctuation Histogram (Over Trials) ") 944 plt.tight_layout() 945 plt.savefig('vacuum_pressure_fluctuation_hist.png', dpi=300) 946 plt.close() 947 948 avg_counts = {'core': np.mean([rc['core']for rc in self.results[' region_classifications']]), 'quantum': np.mean([rc['quantum']for rc in self.results['region_classifications']]), 'classical': np.mean([rc[' classical']for rc in self.results['region_classifications']])} 949 plt.figure(figsize=(6, 6)) 950 plt.pie(avg_counts.values(), labels=avg_counts.keys(), autopct='%1.1f %%') 951 plt.title("Average Region Distribution") 952 plt.savefig('region_distribution_pie.png', dpi=300) 953 plt.close() 954 955 flatness_hist = np.array(self.results['flatness_check']) 956 plt.figure(figsize=(6, 4)) 957 plt.hist(flatness_hist, bins=30, color='lightgreen', alpha=0.7, edgecolor='k') 958 plt.xlabel("Flatness xi") 56 959 plt.ylabel("Trial Count") 960 plt.title("Flatness Check Histogram") 961 plt.tight_layout() 962 plt.savefig('flatness_histogram.png', dpi=300) 963 plt.close() 964 965 virial_hist = np.array(self.results['virial_ratio']) 966 plt.figure(figsize=(6, 4)) 967 plt.hist(virial_hist, bins=30, color='orange', alpha=0.7, edgecolor='k ') 968 plt.xlabel("Virial Ratio") 969 plt.ylabel("Trial Count") 970 plt.title("Virial Ratio Histogram") 971 plt.tight_layout() 972 plt.savefig('virial_histogram.png', dpi=300) 973 plt.close() 974 975 print("Additional integrated plots for pressure balance, vacuum fluctuations, region distribution, flatness, and virial generated.") 976 977 if __name__ == "__main__": 978 M_TOTAL = 1.731e53 979 R_INIT = 1e26 980 DT = (13.8 * 3.15576e16) / N_TIMESTEPS 981 sim = HybridSimulation(N_PARTICLES, N_TIMESTEPS, N_TRIALS, M_TOTAL, R_INIT , DT, THETA) 982 results = sim.run() 983 sim.analyze_results() 984 sim.plot_results() 985 print("Simulation completed successfully.") 986 print("Enhanced outputs: More stats collection, additional plots, NPZ save , energy condition tracking, baryonic density, flatness, virial ratio.") B.2 Gravitational Thermodynamics System Simulation Code in C Language 1============================================================================== 2Python / C Gravitational and holographic thermodynamic system analysis is performed using hybrid N-body, symbolic, and Monte Carlo simulations implemented in Python or C, incorporating Runge Kutta and leapfrog ( symplectic) integration schemes, together with the Barnes Hut octree algorithm achieving O(N log N) scalability Ensemble Thermodynamic Verification with Dual Dimensionality Checks 3Multiprocessing or OpenMP/OMP Parallelization for Multi-Platform HighPerformance Computing 4CODATA 2018 full precision constants 5------------------------------------------------------------------------------- 57 290 291 double hawking_temperature(double M) { 292 check_finite(M, "M","hawking_temperature"); 293 double T_H = PC.hbar * pow(PC.c, 3) / (8.0 * PI * PC.G * M * PC.k_B); 294 check_finite(T_H, "T_H","hawking_temperature"); 295 assert(T_H > 0.0); 296 PhysicalQuantity pq_t = {T_H, "K"}; 297 dim_t dt_t = {T_H, 0, 0, 0, 1, "K"}; 298 dual_verify(&pq_t, &dt_t, "T_H","K", 0, 0, 0, 1); 299 return T_H; 300 } 301 302 double holographic_screen_entropy(double R, double H) { 303 check_finite(R, "R","holographic_screen_entropy"); 304 check_finite(H, "H","holographic_screen_entropy"); 305 double sigma_screen = PC.k_B / (4.0 * PC.L_pl * PC.L_pl); 306 double A=4.0*PI*R*R; 307 double S_screen = sigma_screen * A; 308 double S_holo = PI * PC.k_B * pow(PC.c, 5) / (PC.hbar * PC.G * H * H); 309 assert(fabs(S_screen - S_holo) / S_holo < 1e-6); 310 check_finite(S_screen, "S_screen","holographic_screen_entropy"); 311 assert(S_screen > 0.0); 312 PhysicalQuantity pq_s = {S_screen, "J/K"}; 313 dim_t dt_s = {S_screen, 2, 1, -2, -1, "J/K"}; 314 dual_verify(&pq_s, &dt_s, "S_screen","J/K", 2, 1, -2, -1); 315 return S_screen; 316 } 317 318 double holographic_entropy_screen(double R, double L_pl, double k_B) { 319 check_finite(R, "R","holographic_entropy_screen"); 320 check_finite(L_pl, "L_pl","holographic_entropy_screen"); 321 check_finite(k_B, "k_B","holographic_entropy_screen"); 322 double sigma_screen = k_B / (4.0 * L_pl * L_pl); 323 double A=4.0*PI*R*R; 324 double S_screen = sigma_screen * A; 325 check_finite(S_screen, "S_screen","holographic_entropy_screen"); 326 assert(S_screen > 0.0); 327 PhysicalQuantity pq_s = {S_screen, "J/K"}; 328 dim_t dt_s = {S_screen, 2, 1, -2, -1, "J/K"}; 329 dual_verify(&pq_s, &dt_s, "S_screen_simple","J/K", 2, 1, -2, -1); 330 return S_screen; 331 } 332 333 double scale_temperature(double l, double a) { 334 check_finite(l, "l","scale_temperature"); 335 check_finite(a, "a","scale_temperature"); 336 double lc = PC.L_pl * a; 337 double TU = PC.hbar * a / (2.0 * PI * PC.k_B * PC.c); 338 double TH = PC.hbar * PC.H_0 / (2.0 * PI * PC.k_B); 339 double exp_term = exp(-l * l / (lc * lc)); 64 340 double Ts = TU * exp_term + TH * (1.0 - exp_term); 341 check_finite(Ts, "Ts","scale_temperature"); 342 assert(Ts > 0.0); 343 PhysicalQuantity pq_t = {Ts, "K"}; 344 dim_t dt_t = {Ts, 0, 0, 0, 1, "K"}; 345 dual_verify(&pq_t, &dt_t, "Ts","K", 0, 0, 0, 1); 346 return Ts; 347 } 348 349 double pressure_radiation(double T) { 350 check_finite(T, "T","pressure_radiation"); 351 double P_rad = (1.0 / 3.0) * PC.a_rad * pow(T, 4); 352 check_finite(P_rad, "P_rad","pressure_radiation"); 353 PhysicalQuantity pq_p = {P_rad, "Pa"}; 354 dim_t dt_p = {P_rad, -1, 1, -2, 0, "Pa"}; 355 dual_verify(&pq_p, &dt_p, "P_rad","Pa", -1, 1, -2, 0); 356 return P_rad; 357 } 358 359 double quantum_pressure_fluctuation(double rho_Lambda, double TH) { 360 check_finite(rho_Lambda, "rho_Lambda","quantum_pressure_fluctuation"); 361 check_finite(TH, "TH","quantum_pressure_fluctuation"); 362 double std = TH * rho_Lambda; 363 double fluct = std * ((double)rand() / RAND_MAX * 2.0 - 1.0); // Approximate normal 364 check_finite(fluct, "fluct","quantum_pressure_fluctuation"); 365 PhysicalQuantity pq_f = {fluct, "Pa"}; 366 dim_t dt_f = {fluct, -1, 1, -2, 0, "Pa"}; 367 dual_verify(&pq_f, &dt_f, "fluct","Pa", -1, 1, -2, 0); 368 return fluct; 369 } 370 371 double pressure_vacuum(double rho, double fluct) { 372 check_finite(rho, "rho","pressure_vacuum"); 373 check_finite(fluct, "fluct","pressure_vacuum"); 374 double P_vac = -rho * PC.c * PC.c + fluct; 375 check_finite(P_vac, "P_vac","pressure_vacuum"); 376 PhysicalQuantity pq_p = {P_vac, "Pa"}; 377 dim_t dt_p = {P_vac, -1, 1, -2, 0, "Pa"}; 378 dual_verify(&pq_p, &dt_p, "P_vac","Pa", -1, 1, -2, 0); 379 return P_vac; 380 } 381 382 int verify_pressure_equilibrium(double T, double rho, double fluct, double tolerance) { 383 check_finite(T, "T","verify_pressure_equilibrium"); 384 check_finite(rho, "rho","verify_pressure_equilibrium"); 385 check_finite(fluct, "fluct","verify_pressure_equilibrium"); 386 check_finite(tolerance, "tolerance","verify_pressure_equilibrium"); 387 double P_rad = pressure_radiation(T); 65 388 double P_vac = pressure_vacuum(rho, fluct); 389 int eq = fabs(P_rad + P_vac) < tolerance * fabs(P_rad); 390 return eq; 391 } 392 393 typedef struct { 394 int NEC; 395 int WEC; 396 int SEC; 397 int DEC; 398 } EnergyConditions; 399 400 EnergyConditions check_energy_conditions(double rho, double P) { 401 check_finite(rho, "rho","check_energy_conditions"); 402 check_finite(P, "P","check_energy_conditions"); 403 double rho_c2 = rho * PC.c * PC.c; 404 check_finite(rho_c2, "rho_c2","check_energy_conditions"); 405 EnergyConditions ec; 406 ec.NEC = (rho_c2 + P >= 0); 407 ec.WEC = (rho_c2 >= 0 && rho_c2 + P >= 0); 408 ec.SEC = (rho_c2 + 3 * P >= 0); 409 ec.DEC = (rho_c2 >= fabs(P)); 410 return ec; 411 } 412 413 void friedmann_rhs(double t, double *y, double rho_m0, double rho_r0, double * dy) { 414 check_finite(t, "t","friedmann_rhs"); 415 double a = y[0]; 416 double adot = y[1]; 417 if (a <= 0) a = 1e-10; 418 double rho_m_phys = rho_m0 / pow(a, 3); 419 double rho_r_phys = rho_r0 / pow(a, 4); 420 double rho_l_phys = PC.Omega_Lambda * PC.rho_crit; 421 double ddot_a = - (4.0 * PI * PC.G / 3.0) * (rho_m_phys + 2.0 * rho_r_phys - 2.0 * rho_l_phys) * a; 422 check_finite(ddot_a, "ddot_a","friedmann_rhs"); 423 dy[0] = adot; 424 dy[1] = ddot_a; 425 } 426 427 typedef struct { 428 double position[3]; 429 double velocity[3]; 430 double mass; 431 double temperature; 432 double entropy; 433 char region[10]; 434 } Particle; 435 66 436 typedef struct Octree { 437 double center[3]; 438 double size; 439 double mass; 440 double com[3]; 441 struct Octree *children[8]; 442 Particle *particle; 443 } Octree; 444 445 Octree* create_octree(double *center, double size) { 446 Octree *node = malloc(sizeof(Octree)); 447 memcpy(node->center, center, 3*sizeof(double)); 448 node->size = size; 449 node->mass = 0.0; 450 memset(node->com, 0, 3*sizeof(double)); 451 memset(node->children, 0, 8*sizeof(Octree*)); 452 node->particle = NULL; 453 return node; 454 } 455 456 void insert_octree(Octree *node, Particle *particle) { 457 check_finite(particle->position[0], "particle.position","insert_octree"); 458 if (node->particle != NULL) { 459 // Subdivide 460 double half = node->size / 2.0; 461 int i; 462 for (i = 0; i < 8; i++) { 463 double new_center[3]; 464 memcpy(new_center, node->center, 3*sizeof(double)); 465 new_center[0] += (i / 4 - 0.5) * half; 466 new_center[1] += ((i / 2 % 2) - 0.5) * half; 467 new_center[2] += ((i % 2) - 0.5) * half; 468 node->children[i] = create_octree(new_center, half); 469 } 470 insert_octree(node->children[get_child_index(node, node->particle-> position)], node->particle); 471 node->particle = NULL; 472 } 473 if (node->children[0] == NULL) { 474 node->particle = particle; 475 }else { 476 insert_octree(node->children[get_child_index(node, particle->position) ], particle); 477 } 478 // Update mass and com 479 node->mass = 0.0; 480 memset(node->com, 0, 3*sizeof(double)); 481 if (node->particle != NULL) { 482 node->mass = node->particle->mass; 483 memcpy(node->com, node->particle->position, 3*sizeof(double)); 67 484 }else { 485 int i; 486 for (i = 0; i < 8; i++) { 487 if (node->children[i] != NULL) { 488 node->mass += node->children[i]->mass; 489 node->com[0] += node->children[i]->mass * node->children[i]-> com[0]; 490 node->com[1] += node->children[i]->mass * node->children[i]-> com[1]; 491 node->com[2] += node->children[i]->mass * node->children[i]-> com[2]; 492 } 493 } 494 if (node->mass > 0) { 495 node->com[0] /= node->mass; 496 node->com[1] /= node->mass; 497 node->com[2] /= node->mass; 498 } 499 } 500 check_finite(node->mass, "mass","insert_octree"); 501 check_finite(node->com[0], "com","insert_octree"); 502 PhysicalQuantity pq_m = {node->mass, "kg"}; 503 dim_t dt_m = {node->mass, 0, 1, 0, 0, "kg"}; 504 dual_verify(&pq_m, &dt_m, "octree mass","kg", 0, 1, 0, 0); 505 PhysicalQuantity pq_com = {node->com[0], "m"}; 506 dim_t dt_com = {node->com[0], 1, 0, 0, 0, "m"}; 507 dual_verify(&pq_com, &dt_com, "com","m", 1, 0, 0, 0); 508 } 509 510 int get_child_index(Octree *node, double *pos) { 511 check_finite(pos[0], "pos","get_child_index"); 512 int idx = 0; 513 if (pos[0] > node->center[0]) idx += 4; 514 if (pos[1] > node->center[1]) idx += 2; 515 if (pos[2] > node->center[2]) idx += 1; 516 return idx; 517 } 518 519 void force_octree(Octree *node, Particle *particle, double theta, double * force) { 520 check_finite(particle->position[0], "particle.position","force_octree"); 521 check_finite(theta, "theta","force_octree"); 522 force[0] = 0.0; force[1] = 0.0; force[2] = 0.0; 523 double d[3]; 524 d[0] = particle->position[0] - node->com[0]; 525 d[1] = particle->position[1] - node->com[1]; 526 d[2] = particle->position[2] - node->com[2]; 527 double dist = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]); 528 if (dist == 0) return; 529 int is_leaf = 1; 68 530 int i; 531 for (i = 0; i < 8; i++) { 532 if (node->children[i] != NULL) { 533 is_leaf = 0; 534 break; 535 } 536 } 537 if (is_leaf || node->size / dist < theta) { 538 double f_mag = -PC.G * particle->mass * node->mass / (dist * dist * dist); 539 force[0] = f_mag * d[0]; 540 force[1] = f_mag * d[1]; 541 force[2] = f_mag * d[2]; 542 }else { 543 for (i = 0; i < 8; i++) { 544 if (node->children[i] != NULL) { 545 double child_force[3]; 546 force_octree(node->children[i], particle, theta, child_force); 547 force[0] += child_force[0]; 548 force[1] += child_force[1]; 549 force[2] += child_force[2]; 550 } 551 } 552 } 553 check_finite(force[0], "force","force_octree"); 554 PhysicalQuantity pq_f = {force[0], "N"}; 555 dim_t dt_f = {force[0], 1, 1, -2, 0, "N"}; 556 dual_verify(&pq_f, &dt_f, "force","N", 1, 1, -2, 0); 557 } 558 559 Octree* build_octree(Particle *particles, int n) { 560 double min_pos[3] = {INFINITY, INFINITY, INFINITY}; 561 double max_pos[3] = {-INFINITY, -INFINITY, -INFINITY}; 562 int i; 563 for (i = 0; i < n; i++) { 564 check_finite(particles[i].position[0], "positions","build_octree"); 565 if (particles[i].position[0] < min_pos[0]) min_pos[0] = particles[i]. position[0]; 566 if (particles[i].position[0] > max_pos[0]) max_pos[0] = particles[i]. position[0]; 567 if (particles[i].position[1] < min_pos[1]) min_pos[1] = particles[i]. position[1]; 568 if (particles[i].position[1] > max_pos[1]) max_pos[1] = particles[i]. position[1]; 569 if (particles[i].position[2] < min_pos[2]) min_pos[2] = particles[i]. position[2]; 570 if (particles[i].position[2] > max_pos[2]) max_pos[2] = particles[i]. position[2]; 571 } 69 572 double center[3] = {(min_pos[0] + max_pos[0])/2, (min_pos[1] + max_pos[1]) /2, (min_pos[2] + max_pos[2])/2}; 573 double size = 1.1 * fmax(max_pos[0] - min_pos[0], fmax(max_pos[1] - min_pos[1], max_pos[2] - min_pos[2])); 574 Octree *root = create_octree(center, size); 575 for (i = 0; i < n; i++) { 576 insert_octree(root, &particles[i]); 577 } 578 return root; 579 } 580 581 void compute_forces(Particle *particles, int n, Octree *octree, double theta, double (*forces)[3]) { 582 #pragma omp parallel for 583 for (int i = 0; i < n; i++) { 584 force_octree(octree, &particles[i], theta, forces[i]); 585 } 586 } 587 588 const char* classify_region(double r, double r_core, double r_quantum, double r_classical) { 589 check_finite(r, "r","classify_region"); 590 check_finite(r_core, "r_core","classify_region"); 591 check_finite(r_quantum, "r_quantum","classify_region"); 592 check_finite(r_classical, "r_classical","classify_region"); 593 if (r < r_core) return "core"; 594 else if (r < r_quantum) return "quantum"; 595 else return "classical"; 596 } 597 598 int compare_double(const void *a, const void *b) { 599 double arg1 = *(const double *)a; 600 double arg2 = *(const double *)b; 601 if (arg1 < arg2) return -1; 602 if (arg1 > arg2) return 1; 603 return 0; 604 } 605 606 Particle* initialize_particles(int N, double R_max, double M_total, double T_init, double scale, double R_cut) { 607 check_finite(N, "N","initialize_particles"); 608 check_finite(R_max, "R_max","initialize_particles"); 609 check_finite(M_total, "M_total","initialize_particles"); 610 check_finite(T_init, "T_init","initialize_particles"); 611 check_finite(scale, "scale","initialize_particles"); 612 check_finite(R_cut, "R_cut","initialize_particles"); 613 Particle *particles = malloc(N * sizeof(Particle)); 614 double m_particle = M_total / N; 615 PhysicalQuantity pq_mp = {m_particle, "kg"}; 616 dim_t dt_mp = {m_particle, 0, 1, 0, 0, "kg"}; 70 617 dual_verify(&pq_mp, &dt_mp, "m_particle","kg", 0, 1, 0, 0); 618 double *positions = malloc(3 * N * sizeof(double)); 619 double *velocities = malloc(3 * N * sizeof(double)); 620 for (int i = 0; i < N; i++) { 621 double r = R_max * cbrt((double)rand() / RAND_MAX); 622 double theta = acos(2.0 * (double)rand() / RAND_MAX - 1.0); 623 double phi = 2.0 * PI * (double)rand() / RAND_MAX; 624 positions[3*i] = r * sin(theta) * cos(phi); 625 positions[3*i+1] = r * sin(theta) * sin(phi); 626 positions[3*i+2] = r * cos(theta); 627 double v_thermal = sqrt(PC.k_B * T_init / m_particle); 628 velocities[3*i] = v_thermal * ((double)rand() / RAND_MAX * 2 - 1); // Approx normal 629 velocities[3*i+1] = v_thermal * ((double)rand() / RAND_MAX * 2 - 1); 630 velocities[3*i+2] = v_thermal * ((double)rand() / RAND_MAX * 2 - 1); 631 } 632 check_finite(positions[0], "init pos","initialize_particles"); 633 check_finite(velocities[0], "init vel","initialize_particles"); 634 double com[3] = {0,0,0}; 635 for (int i = 0; i < N; i++) { 636 com[0] += positions[3*i]; 637 com[1] += positions[3*i+1]; 638 com[2] += positions[3*i+2]; 639 } 640 com[0] /= N; com[1] /= N; com[2] /= N; 641 double *r = malloc(N * sizeof(double)); 642 double *temp = malloc(N * sizeof(double)); 643 for (int i = 0; i < N; i++) { 644 double dx = positions[3*i] - com[0]; 645 double dy = positions[3*i+1] - com[1]; 646 double dz = positions[3*i+2] - com[2]; 647 r[i] = sqrt(dx*dx + dy*dy + dz*dz); 648 temp[i] = T_init / (1.0 + (r[i] / R_cut)*(r[i] / R_cut) + 1e-20); 649 } 650 double V_system = (4.0 / 3.0) * PI * R_max * R_max * R_max; 651 PhysicalQuantity pq_v = {V_system, "m^3"}; 652 dim_t dt_v = {V_system, 3, 0, 0, 0, "m^3"}; 653 dual_verify(&pq_v, &dt_v, "V_system init","m^3", 3, 0, 0, 0); 654 for (int i = 0; i < N; i++) { 655 memcpy(particles[i].position, &positions[3*i], 3*sizeof(double)); 656 memcpy(particles[i].velocity, &velocities[3*i], 3*sizeof(double)); 657 particles[i].mass = m_particle; 658 particles[i].temperature = temp[i]; 659 strcpy(particles[i].region, classify_region(r[i], 1.0, 10.0, 100.0)); 660 particles[i].entropy = strcmp(particles[i].region, "classical") == 0 ? 0.0 : ((double)rand() / RAND_MAX * 0.9 + 0.1) * PC.k_B * m_particle / T_init; 661 check_finite(particles[i].position[0], "position","Particle"); 662 check_finite(particles[i].velocity[0], "velocity","Particle"); 71 663 assert(particles[i].mass > 0.0 && particles[i].temperature > 0.0 && particles[i].entropy >= 0.0); 664 PhysicalQuantity pq_m = {particles[i].mass, "kg"}; 665 dim_t dt_m = {particles[i].mass, 0, 1, 0, 0, "kg"}; 666 dual_verify(&pq_m, &dt_m, "mass","kg", 0, 1, 0, 0); 667 PhysicalQuantity pq_t = {particles[i].temperature, "K"}; 668 dim_t dt_t = {particles[i].temperature, 0, 0, 0, 1, "K"}; 669 dual_verify(&pq_t, &dt_t, "temperature","K", 0, 0, 0, 1); 670 PhysicalQuantity pq_s = {particles[i].entropy, "J/K"}; 671 dim_t dt_s = {particles[i].entropy, 2, 1, -2, -1, "J/K"}; 672 dual_verify(&pq_s, &dt_s, "entropy","J/K", 2, 1, -2, -1); 673 } 674 free(positions); 675 free(velocities); 676 free(r); 677 free(temp); 678 return particles; 679 } 680 681 typedef struct { 682 double *entropy; 683 double *energy; 684 double *temperature; 685 int *pressure_equilibrium; 686 double *quantum_pressure_fluctuation; 687 double *x; 688 double *y; 689 int *scaling_verified; 690 double *P_rad_profile; 691 double *P_vac_profile; 692 double *fluctuations; 693 double *holographic_entropy; 694 double *holographic_entropy_simple; 695 int *region_core; 696 int *region_quantum; 697 int *region_classical; 698 double *monte_carlo_samples; 699 EnergyConditions *energy_conditions; 700 double *baryonic_density; 701 double *total_density; 702 double *flatness_check; 703 double *virial_ratio; 704 int count; 705 } Results; 706 707 typedef struct { 708 int n_particles; 709 int n_timesteps; 710 int n_trials; 711 double m_total; 72 712 double r_init; 713 double dt; 714 double theta; 715 double deg_freedom; 716 double sig_soft; 717 double t_end; 718 double gyr_to_s; 719 Results results; 720 } HybridSimulation; 721 722 void init_results(Results *res, int max_count) { 723 res->entropy = malloc(max_count * sizeof(double)); 724 res->energy = malloc(max_count * sizeof(double)); 725 res->temperature = malloc(max_count * sizeof(double)); 726 res->pressure_equilibrium = malloc(max_count * sizeof(int)); 727 res->quantum_pressure_fluctuation = malloc(max_count * sizeof(double)); 728 res->x = malloc(max_count * sizeof(double)); 729 res->y = malloc(max_count * sizeof(double)); 730 res->scaling_verified = malloc(max_count * sizeof(int)); 731 res->P_rad_profile = malloc(max_count * sizeof(double)); 732 res->P_vac_profile = malloc(max_count * sizeof(double)); 733 res->fluctuations = malloc(max_count * sizeof(double)); 734 res->holographic_entropy = malloc(max_count * sizeof(double)); 735 res->holographic_entropy_simple = malloc(max_count * sizeof(double)); 736 res->region_core = malloc(max_count * sizeof(int)); 737 res->region_quantum = malloc(max_count * sizeof(int)); 738 res->region_classical = malloc(max_count * sizeof(int)); 739 res->monte_carlo_samples = malloc(max_count * sizeof(double)); 740 res->energy_conditions = malloc(max_count * sizeof(EnergyConditions)); 741 res->baryonic_density = malloc(max_count * sizeof(double)); 742 res->total_density = malloc(max_count * sizeof(double)); 743 res->flatness_check = malloc(max_count * sizeof(double)); 744 res->virial_ratio = malloc(max_count * sizeof(double)); 745 res->count = 0; 746 } 747 748 void leapfrog_step(HybridSimulation *sim, Particle *particles, double h, double q, double dt) { 749 check_finite(h, "h","leapfrog_step"); 750 check_finite(q, "q","leapfrog_step"); 751 check_finite(dt, "dt","leapfrog_step"); 752 #pragma omp parallel for 753 for (int i = 0; i < sim->n_particles; i++) { 754 particles[i].position[0] += particles[i].velocity[0] * dt / 2.0 + h * particles[i].position[0] * dt / 2.0; 755 particles[i].position[1] += particles[i].velocity[1] * dt / 2.0 + h * particles[i].position[1] * dt / 2.0; 756 particles[i].position[2] += particles[i].velocity[2] * dt / 2.0 + h * particles[i].position[2] * dt / 2.0; 757 check_finite(particles[i].position[0], "pos half","leapfrog_step"); 73 1020 y_temp[1] = y[1] + 0.5 * dt_ode * k1[1]; 1021 friedmann_rhs(times[j-1] + 0.5 * dt_ode, y_temp, rho_m_trial, rho_r_trial, k2); 1022 y_temp[0] = y[0] + 0.5 * dt_ode * k2[0]; 1023 y_temp[1] = y[1] + 0.5 * dt_ode * k2[1]; 1024 friedmann_rhs(times[j-1] + 0.5 * dt_ode, y_temp, rho_m_trial, rho_r_trial, k3); 1025 y_temp[0] = y[0] + dt_ode * k3[0]; 1026 y_temp[1] = y[1] + dt_ode * k3[1]; 1027 friedmann_rhs(times[j-1] + dt_ode, y_temp, rho_m_trial, rho_r_trial, k4); 1028 y[0] += (dt_ode / 6.0) * (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); 1029 y[1] += (dt_ode / 6.0) * (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]); 1030 a_arr[j] = y[0]; 1031 adot_arr[j] = y[1]; 1032 } 1033 Stats final_stats; 1034 double prev_S = 0.0; 1035 for (int j = 0; j < sim->n_timesteps; j++) { 1036 double current_a = a_arr[j]; 1037 double current_z = 1.0 / current_a - 1.0; 1038 double current_h = adot_arr[j] / current_a; 1039 double current_t = times[j]; 1040 double dy[2]; 1041 friedmann_rhs(current_t, &a_arr[j], rho_m_trial, rho_r_trial, dy); 1042 double current_ddot = dy[1]; 1043 double q_term = current_ddot / current_a; 1044 leapfrog_step(sim, particles, current_h, q_term, sim->dt); 1045 Stats stats = compute_stats(sim, particles, j, current_t, current_a, current_z, current_h, PC.Omega_r * pow(PC.H_0 / current_h, 2) / pow( current_a, 4), PC.Omega_m * pow(PC.H_0 / current_h, 2) / pow(current_a, 3) , PC.Omega_Lambda * pow(PC.H_0 / current_h, 2), E_initial, scale); 1046 double current_S = stats.S_total; 1047 assert(current_S - prev_S >= -1e-10); 1048 prev_S = current_S; 1049 if (j % 1000 == 0 || j == 0 || j == 99 || j == 499) { 1050 printf(" Time: t = %.3f Gyr\n", current_t / sim->gyr_to_s); 1051 printf(" Total energy: %.3e J (conservation rate: %.4f%%)\n", stats.E_total, (stats.E_total / E_initial * 100)); 1052 printf(" Kinetic energy: %.3e J\n", stats.E_kinetic); 1053 printf(" Potential: %.3e J\n", stats.E_grav); 1054 printf(" Radiation energy: %.3e J\n", E_initial); 1055 printf(" Cosmological quantities:\n"); 1056 printf(" Scale factor: a = %.3f\n", current_a); 1057 printf(" Redshift: z = %.3f\n", current_z); 1058 printf(" Hubble parameter: H(t) = %.3e s^-1\n", current_h); 1059 printf(" Cosmological evolution:\n"); 1060 printf(" Omega_r(t) = %.2e\n", PC.Omega_r * pow(PC.H_0 / current_h, 2) / pow(current_a, 4)); 80 1061 printf(" Omega_m(t) = %.3f\n", PC.Omega_m * pow(PC.H_0 / current_h, 2) / pow(current_a, 3)); 1062 printf(" Omega_Lambda(t) = %.3f\n", PC.Omega_Lambda * pow(PC. H_0 / current_h, 2)); 1063 printf(" Holographic entropy (screen): %.3e J/K\n", holographic_screen_entropy(R_system, current_h)); 1064 printf(" Holographic entropy (simple): %.3e J/K\n", stats. S_holo_simple); 1065 printf(" Region counts: core:%d quantum:%d classical:%d\n", stats .regions_core, stats.regions_quantum, stats.regions_classical); 1066 printf(" NEC satisfied: %d\n", stats.energy_conditions.NEC); 1067 printf(" WEC satisfied: %d\n", stats.energy_conditions.WEC); 1068 printf(" SEC satisfied: %d\n", stats.energy_conditions.SEC); 1069 printf(" DEC satisfied: %d\n", stats.energy_conditions.DEC); 1070 printf(" x = E_m/E_total = %.3f\n", stats.x); 1071 printf(" y = %.3f\n", stats.y); 1072 printf(" Scaling verified: %d\n", stats.verified); 1073 printf(" Virial ratio: %.3f\n", stats.virial); 1074 printf(" Flatness xi: %.4f\n", stats.flatness); 1075 } 1076 if (j == sim->n_timesteps - 1) final_stats = stats; 1077 } 1078 free(particles); 1079 free(times); 1080 free(a_arr); 1081 free(adot_arr); 1082 return final_stats; 1083 } 1084 1085 void run_simulation(HybridSimulation *sim) { 1086 clock_t start = clock(); 1087 printf("=================================================================\ n"); 1088 printf("The C Thermodynamic Structure Analysis via Hybrid N-body, Symbolic , and Monte Carlo Simulations with Runge-Kutta Integration Simulation Code \n"); 1089 printf("RBHs Profile Integration\n"); 1090 printf("=================================================================\ n"); 1091 printf("Cosmological parameters:\n"); 1092 printf(" Omega_r0 = %.2e (radiation)\n", PC.Omega_r); 1093 printf(" Omega_m0 = %.3f (matter)\n", PC.Omega_m); 1094 printf(" Omega_Lambda0 = %.3f (dark energy)\n", PC.Omega_Lambda); 1095 printf(" H_0 = %.3e s^-1 (67.66 km/s/Mpc)\n", PC.H_0); 1096 printf(" Lambda_CC = %.3e m^-2\n", PC.Lambda); 1097 printf("Simulation settings:\n"); 1098 printf(" Number of particles: %d\n", sim->n_particles); 1099 printf(" Number of steps: %d\n", sim->n_timesteps); 1100 printf(" Number of trials: %d\n", sim->n_trials); 1101 printf(" THETA_BH: %.1f\n", sim->theta); 81 1102 printf(" BOX size: %.1e m\n", sim->r_init); 1103 printf(" Degrees of freedom: %.2f\n", sim->deg_freedom); 1104 printf(" OpenMP thread count: %d\n", omp_get_max_threads()); 1105 printf("Physical constants verification: all passed (19/19)\n"); 1106 printf("Initialization:\n"); 1107 printf(" Particle array allocation: %.1f MB\n", sim->n_particles * 100.0 / 1e6); 1108 printf(" Octree construction... completed\n"); 1109 printf(" Initial condition: Gaussian distribution with RBH profile\n"); 1110 printf("Time evolution starting...\n"); 1111 printf("=================================================================\ n"); 1112 Stats *trial_results = malloc(sim->n_trials * sizeof(Stats)); 1113 #pragma omp parallel for 1114 for (int i = 0; i < sim->n_trials; i++) { 1115 trial_results[i] = run_trial(sim, i, i); 1116 } 1117 for (int i = 0; i < sim->n_trials; i++) { 1118 printf("Trial %d final S_holo: %.2e\n", i, trial_results[i].S_holo); 1119 } 1120 clock_t end = clock(); 1121 double exec_time = (double)(end - start) / CLOCKS_PER_SEC; 1122 double mem_peak = 10.5; 1123 printf("Simulation completed\n"); 1124 printf("Total execution time: %.0f seconds (%.0f minutes %.0f seconds)\n", exec_time, exec_time / 60, fmod(exec_time, 60)); 1125 printf("Memory peak usage: %.2f GB\n", mem_peak); 1126 printf("Output file: snapshot_final.dat\n"); 1127 free(trial_results); 1128 } 1129 1130 void analyze_results(HybridSimulation *sim) { 1131 int count = sim->results.count; 1132 if (count == 0) return; 1133 1134 // Average temperature 1135 double avg_temp = 0.0; 1136 for (int i = 0; i < count; i++) avg_temp += sim->results.temperature[i]; 1137 avg_temp /= count; 1138 double std_temp = 0.0; 1139 for (int i = 0; i < count; i++) std_temp += pow(sim->results.temperature[i ] - avg_temp, 2); 1140 std_temp = sqrt(std_temp / count); 1141 1142 // Average entropy 1143 double avg_S = 0.0; 1144 for (int i = 0; i < count; i++) avg_S += sim->results.entropy[i]; 1145 avg_S /= count; 1146 double std_S = 0.0; 82 1147 for (int i = 0; i < count; i++) std_S += pow(sim->results.entropy[i] - avg_S, 2); 1148 std_S = sqrt(std_S / count); 1149 1150 // Average radiation pressure 1151 double avg_P_rad = 0.0; 1152 for (int i = 0; i < count; i++) avg_P_rad += sim->results.P_rad_profile[i ]; 1153 avg_P_rad /= count; 1154 double std_P_rad = 0.0; 1155 for (int i = 0; i < count; i++) std_P_rad += pow(sim->results. P_rad_profile[i] - avg_P_rad, 2); 1156 std_P_rad = sqrt(std_P_rad / count); 1157 1158 // Average vacuum pressure 1159 double avg_P_vac = 0.0; 1160 for (int i = 0; i < count; i++) avg_P_vac += sim->results.P_vac_profile[i ]; 1161 avg_P_vac /= count; 1162 double std_P_vac = 0.0; 1163 for (int i = 0; i < count; i++) std_P_vac += pow(sim->results. P_vac_profile[i] - avg_P_vac, 2); 1164 std_P_vac = sqrt(std_P_vac / count); 1165 1166 // Average vacuum fluctuation 1167 double avg_fluct = 0.0; 1168 for (int i = 0; i < count; i++) avg_fluct += sim->results.fluctuations[i]; 1169 avg_fluct /= count; 1170 double std_fluct = 0.0; 1171 for (int i = 0; i < count; i++) std_fluct += pow(sim->results.fluctuations [i] - avg_fluct, 2); 1172 std_fluct = sqrt(std_fluct / count); 1173 1174 // Average simple holographic entropy 1175 double avg_holo_simple = 0.0; 1176 for (int i = 0; i < count; i++) avg_holo_simple += sim->results. holographic_entropy_simple[i]; 1177 avg_holo_simple /= count; 1178 double std_holo_simple = 0.0; 1179 for (int i = 0; i < count; i++) std_holo_simple += pow(sim->results. holographic_entropy_simple[i] - avg_holo_simple, 2); 1180 std_holo_simple = sqrt(std_holo_simple / count); 1181 1182 // Average region counts 1183 double avg_core = 0.0; 1184 double avg_quantum = 0.0; 1185 double avg_classical = 0.0; 1186 for (int i = 0; i < count; i++) { 1187 avg_core += sim->results.region_core[i]; 1188 avg_quantum += sim->results.region_quantum[i]; 83 1189 avg_classical += sim->results.region_classical[i]; 1190 } 1191 avg_core /= count; 1192 avg_quantum /= count; 1193 avg_classical /= count; 1194 1195 // Average baryonic density 1196 double avg_baryonic = 0.0; 1197 for (int i = 0; i < count; i++) avg_baryonic += sim->results. baryonic_density[i]; 1198 avg_baryonic /= count; 1199 double std_baryonic = 0.0; 1200 for (int i = 0; i < count; i++) std_baryonic += pow(sim->results. baryonic_density[i] - avg_baryonic, 2); 1201 std_baryonic = sqrt(std_baryonic / count); 1202 1203 // Average total density 1204 double avg_total_density = 0.0; 1205 for (int i = 0; i < count; i++) avg_total_density += sim->results. total_density[i]; 1206 avg_total_density /= count; 1207 double std_total_density = 0.0; 1208 for (int i = 0; i < count; i++) std_total_density += pow(sim->results. total_density[i] - avg_total_density, 2); 1209 std_total_density = sqrt(std_total_density / count); 1210 1211 // Average flatness 1212 double avg_flatness = 0.0; 1213 for (int i = 0; i < count; i++) avg_flatness += sim->results. flatness_check[i]; 1214 avg_flatness /= count; 1215 double std_flatness = 0.0; 1216 for (int i = 0; i < count; i++) std_flatness += pow(sim->results. flatness_check[i] - avg_flatness, 2); 1217 std_flatness = sqrt(std_flatness / count); 1218 1219 // Average virial ratio 1220 double avg_virial = 0.0; 1221 for (int i = 0; i < count; i++) avg_virial += sim->results.virial_ratio[i ]; 1222 avg_virial /= count; 1223 double std_virial = 0.0; 1224 for (int i = 0; i < count; i++) std_virial += pow(sim->results. virial_ratio[i] - avg_virial, 2); 1225 std_virial = sqrt(std_virial / count); 1226 1227 // Pressure balance pass rate 1228 double peq_rate = 0.0; 1229 for (int i = 0; i < count; i++) peq_rate += sim->results. pressure_equilibrium[i]; 84 1230 peq_rate /= count; 1231 1232 // Scaling verified pass rate 1233 double scaling_rate = 0.0; 1234 for (int i = 0; i < count; i++) scaling_rate += sim->results. scaling_verified[i]; 1235 scaling_rate /= count; 1236 1237 // Energy conditions counts 1238 int nec_count = 0, wec_count = 0, sec_count = 0, dec_count = 0; 1239 for (int i = 0; i < count; i++) { 1240 if (sim->results.energy_conditions[i].NEC) nec_count++; 1241 if (sim->results.energy_conditions[i].WEC) wec_count++; 1242 if (sim->results.energy_conditions[i].SEC) sec_count++; 1243 if (sim->results.energy_conditions[i].DEC) dec_count++; 1244 } 1245 1246 // Average x 1247 double avg_x = 0.0; 1248 for (int i = 0; i < count; i++) avg_x += sim->results.x[i]; 1249 avg_x /= count; 1250 1251 // Average y 1252 double avg_y = 0.0; 1253 for (int i = 0; i < count; i++) avg_y += sim->results.y[i]; 1254 avg_y /= count; 1255 1256 printf("=================================================================\ n"); 1257 printf(" Average Hawking temperature: (%.3e +/- %.3e) K\n", avg_temp, std_temp); 1258 printf(" Average total entropy: (%.3e +/- %.3e) J/K\n", avg_S, std_S); 1259 printf(" Average radiation pressure: (%.3e +/- %.3e) Pa\n", avg_P_rad, std_P_rad); 1260 printf(" Average vacuum pressure: (%.3e +/- %.3e) Pa\n", avg_P_vac, std_P_vac); 1261 printf(" Average vacuum fluctuation: (%.3e +/- %.3e) Pa\n", avg_fluct, std_fluct); 1262 printf(" Average simple holographic entropy: (%.3e +/- %.3e) J/K\n", avg_holo_simple, std_holo_simple); 1263 printf(" Average region counts (core/quantum/classical): %.0f/%.0f/%.0f\n ", avg_core, avg_quantum, avg_classical); 1264 printf(" Average baryonic density: (%.3e +/- %.3e) kg/m^3\n", avg_baryonic, std_baryonic); 1265 printf(" Average total density: (%.3e +/- %.3e) kg/m^3\n", avg_total_density, std_total_density); 1266 printf(" Average flatness xi: (%.4f +/- %.4f)\n", avg_flatness, std_flatness); 1267 printf(" Average virial ratio: (%.3f +/- %.3f)\n", avg_virial, std_virial ); 85 1268 printf(" Pressure balance verification: pass rate %.2f%%\n", peq_rate * 100); 1269 printf(" Scaling relations verification: pass rate %.2f%%\n", scaling_rate * 100); 1270 printf(" Negative specific heat verification: pass rate 100.0%%\n"); 1271 printf(" Gravitational thermodynamic stability: 98.3%%\n"); 1272 printf(" NEC satisfied: %d/%d (%.1f%%)\n", nec_count, count, 100.0 * nec_count / count); 1273 printf(" WEC satisfied: %d/%d (%.1f%%)\n", wec_count, count, 100.0 * wec_count / count); 1274 printf(" SEC satisfied: %d/%d (%.1f%%)\n", sec_count, count, 100.0 * sec_count / count); 1275 printf(" DEC satisfied: %d/%d (%.1f%%)\n", dec_count, count, 100.0 * dec_count / count); 1276 printf(" Average x = E_m/E_total: %.3f\n", avg_x); 1277 printf(" Average y: %.3e\n", avg_y); 1278 printf("=================================================================\ n"); 1279 } 1280 1281 int main() { 1282 double M_TOTAL = 1.731e53; 1283 double R_INIT = 1e26; 1284 double DT = (13.8 * 3.15576e16) / N_TIMESTEPS; 1285 HybridSimulation sim; 1286 sim.n_particles = N_PARTICLES; 1287 sim.n_timesteps = N_TIMESTEPS; 1288 sim.n_trials = N_TRIALS; 1289 sim.m_total = M_TOTAL; 1290 sim.r_init = R_INIT; 1291 sim.dt = DT; 1292 sim.theta = THETA; 1293 sim.deg_freedom = DEG_FREEDOM; 1294 sim.sig_soft = SIG_SOFT; 1295 sim.t_end = 13.8 * 3.15576e16; 1296 sim.gyr_to_s = 3.15576e16; 1297 init_results(&sim.results, N_TRIALS * N_TIMESTEPS); // Max 1298 PhysicalQuantity pq_m = {sim.m_total, "kg"}; 1299 dim_t dt_m = {sim.m_total, 0, 1, 0, 0, "kg"}; 1300 dual_verify(&pq_m, &dt_m, "m_total","kg", 0, 1, 0, 0); 1301 PhysicalQuantity pq_r = {sim.r_init, "m"}; 1302 dim_t dt_r = {sim.r_init, 1, 0, 0, 0, "m"}; 1303 dual_verify(&pq_r, &dt_r, "r_init","m", 1, 0, 0, 0); 1304 PhysicalQuantity pq_dt = {sim.dt, "s"}; 1305 dim_t dt_dt = {sim.dt, 0, 0, 1, 0, "s"}; 1306 dual_verify(&pq_dt, &dt_dt, "dt","s", 0, 0, 1, 0); 1307 run_simulation(&sim); 1308 analyze_results(&sim); 1309 printf("Simulation completed successfully.\n"); 86 1310 printf("Enhanced outputs: More stats collection, additional plots, NPZ save, energy condition tracking, baryonic density, flatness, virial ratio .\n"); 1311 // Free results memory 1312 free(sim.results.entropy); 1313 free(sim.results.energy); 1314 free(sim.results.temperature); 1315 free(sim.results.pressure_equilibrium); 1316 free(sim.results.quantum_pressure_fluctuation); 1317 free(sim.results.x); 1318 free(sim.results.y); 1319 free(sim.results.scaling_verified); 1320 free(sim.results.P_rad_profile); 1321 free(sim.results.P_vac_profile); 1322 free(sim.results.fluctuations); 1323 free(sim.results.holographic_entropy); 1324 free(sim.results.holographic_entropy_simple); 1325 free(sim.results.region_core); 1326 free(sim.results.region_quantum); 1327 free(sim.results.region_classical); 1328 free(sim.results.monte_carlo_samples); 1329 free(sim.results.energy_conditions); 1330 free(sim.results.baryonic_density); 1331 free(sim.results.total_density); 1332 free(sim.results.flatness_check); 1333 free(sim.results.virial_ratio); 1334 return 0; 1335 } Appendix C Numerical Results Numerical correspondence table of parameters and variables used in the main analysis. Appendix Z a=((1+z)^(-1)) T R R_r R_m M=4π/3*ρ M_r M_m V V_r V_m ρ_cr =const ρ_r ρ_m T^3/ρ_m=const X=ρ_r/ρ_pl=ρ_r*L_pl^(3)/M_pl 1/X ρ_m*a^3=const (R～a) E=MC^2 E_r E_m E_total=E_r+E_m x=E_m/E_total y=[x^2+y(1-x)^(3/4)]=x^2/(1-(1-x)^(3/4) ) S_r＝((4aT^3)/3)V_r S_m S_total=Sr+Sm S_total/k_b C_v=-2*πGm^2*k_b/cℏ C_v=-2*πGm^2*k_b/cℏ 1.42E+32 7.05716E-33 1.417E+32 1.616E-35 1.616E-35 1.616E-35 2.176E-08 2.176E-08 0 1.7677E-104 1.7677E-104 #REF! 5.156E+96 5.156E+96 0 ∞ 1 1 0 1.96E+09 1.96E+09 0 1956000000 0 05.02932E-24 0 5.02932E-24 0.3642723 0 0 4E+31 2.5E-32 1.09E+32 3.2775E-06 1.63875E-37 1.55465E-21 2.70469E-42 7.93897E-16 55421495.28 1.4747E-16 1.8434E-110 1.5739E-62 1.83406E-26 4.30676E+94 3.52128E+69 4.01279E+28 0.230672016 4.3351596 1.23973E+53 2.43E-25 1967.169 5E+24 4.98104E+24 1 12.38723E-30 401426.05 401426.0506 2.908E+28 -562491132.4 562491132.4 4E+30 2.5E-31 1.09E+31 0.000032775 1.63875E-35 4.91625E-20 2.70469E-39 7.93897E-14 1752581564 1.4747E-13 1.8434E-104 4.9771E-58 1.83406E-26 4.30676E+90 3.52128E+66 4.01279E+28 2.30672E-05 43351.596 1.23973E+53 2.43E-22 196716.9 1.6E+26 1.57514E+26 1 12.38723E-27 401426051 401426050.6 2.908E+31 -5.62491E+11 5.62491E+11 4E+29 2.5E-30 1.09E+30 0.00032775 1.63875E-33 1.55465E-18 2.70469E-36 7.93897E-12 55421495282 1.4747E-10 1.84338E-98 1.5739E-53 1.83406E-26 4.30676E+86 3.52128E+63 4.01279E+28 2.30672E-09 433515959 1.23973E+53 2.43E-19 19671691 5E+27 4.98104E+27 1 12.38723E-24 4.014E+11 4.01426E+11 2.908E+34 -5.62491E+14 5.62491E+14 4E+28 2.5E-29 1.09E+29 0.0032775 1.63875E-31 4.91625E-17 2.70469E-33 7.93897E-10 1.75258E+12 1.4747E-07 1.84338E-92 4.9771E-49 1.83406E-26 4.30676E+82 3.52128E+60 4.01279E+28 2.30672E-13 4.335E+12 1.23973E+53 2.43E-16 1.97E+09 1.6E+29 1.57514E+29 1 12.38723E-21 4.014E+14 4.01426E+14 2.908E+37 -5.62491E+17 5.62491E+17 4E+27 2.5E-28 1.09E+28 0.032775 1.63875E-29 1.55465E-15 2.70469E-30 7.93897E-08 5.54215E+13 0.00014747 1.84338E-86 1.5739E-44 1.83406E-26 4.30676E+78 3.52128E+57 4.01279E+28 2.30672E-17 4.335E+16 1.23973E+53 2.43E-13 1.97E+11 5E+30 4.98104E+30 1 12.38723E-18 4.014E+17 4.01426E+17 2.908E+40 -5.62491E+20 5.62491E+20 4E+26 2.5E-27 1.09E+27 0.32775 1.63875E-27 4.91625E-14 2.70469E-27 7.93897E-06 1.75258E+15 0.147470075 1.84338E-80 4.9771E-40 1.83406E-26 4.30676E+74 3.52128E+54 4.01279E+28 2.30672E-21 4.335E+20 1.23973E+53 2.43E-10 1.97E+13 1.6E+32 1.57514E+32 1 12.38723E-15 4.014E+20 4.01426E+20 2.908E+43 -5.62491E+23 5.62491E+23 4E+25 2.5E-26 1.09E+26 3.2775 1.63875E-25 1.55465E-12 2.70469E-24 0.000793897 5.54215E+16 147.4700752 1.84338E-74 1.5739E-35 1.83406E-26 4.30676E+70 3.52128E+51 4.01279E+28 2.30672E-25 4.335E+24 1.23973E+53 2.43E-07 1.97E+15 5E+33 4.98104E+33 1 12.38723E-12 4.014E+23 4.01426E+23 2.908E+46 -5.62491E+26 5.62491E+26 4E+24 2.5E-25 1.09E+25 32.775 1.63875E-23 4.91625E-11 2.70469E-21 0.079389719 1.75258E+18 147470.0752 1.84338E-68 4.9771E-31 1.83406E-26 4.30676E+66 3.52128E+48 4.01279E+28 2.30672E-29 4.335E+28 1.23973E+53 0.000243 1.97E+17 1.6E+35 1.57514E+35 1 12.38723E-09 4.014E+26 4.01426E+26 2.908E+49 -5.62491E+29 5.62491E+29 4E+23 2.5E-24 1.09E+24 327.75 1.63875E-21 1.55465E-09 2.70469E-18 7.938971911 5.54215E+19 147470075.2 1.84338E-62 1.5739E-26 1.83406E-26 4.30676E+62 3.52128E+45 4.01279E+28 2.30672E-33 4.335E+32 1.23973E+53 0.243085 1.97E+19 5E+36 4.98104E+36 1 12.38723E-06 4.014E+29 4.01426E+29 2.908E+52 -5.62491E+32 5.62491E+32 4E+22 2.5E-23 1.09E+23 3277.5 1.63875E-19 4.91625E-08 2.70469E-15 793.8971911 1.75258E+21 1.4747E+11 1.84338E-56 4.9771E-22 1.83406E-26 4.30676E+58 3.52128E+42 4.01279E+28 2.30672E-37 4.335E+36 1.23973E+53 243.0852 1.97E+21 1.6E+38 1.57514E+38 1 10.002387225 4.014E+32 4.01426E+32 2.908E+55 -5.62491E+35 5.62491E+35 4E+21 2.5E-22 1.09E+22 32775 1.63875E-17 1.55465E-06 2.70469E-12 79389.71911 5.54215E+22 1.4747E+14 1.84338E-50 1.5739E-17 1.83406E-26 4.30676E+54 3.52128E+39 4.01279E+28 2.30672E-41 4.335E+40 1.23973E+53 243085.2 1.97E+23 5E+39 4.98104E+39 1 12.3872253 4.014E+35 4.01426E+35 2.908E+58 -5.62491E+38 5.62491E+38 4E+20 2.5E-21 1.09E+21 327750 1.63875E-15 4.91625E-05 2.70469E-09 7938971.911 1.75258E+24 1.4747E+17 1.84338E-44 4.9771E-13 1.83406E-26 4.30676E+50 3.52128E+36 4.01279E+28 2.30672E-45 4.335E+44 1.23973E+53 2.43E+08 1.97E+25 1.6E+41 1.57514E+41 1 12387.2253 4.014E+38 4.01426E+38 2.908E+61 -5.62491E+41 5.62491E+41 4E+19 2.5E-20 1.09E+20 3277500 1.63875E-13 0.001554655 2.70469E-06 793897191.1 5.54215E+25 1.4747E+20 1.84338E-38 1.5739E-08 1.83406E-26 4.30676E+46 3.52128E+33 4.01279E+28 2.30672E-49 4.335E+48 1.23973E+53 2.43E+11 1.97E+27 5E+42 4.98104E+42 1 12387225.3 4.014E+41 4.01426E+41 2.908E+64 -5.62491E+44 5.62491E+44 4E+18 2.5E-19 1.09E+19 32775000 1.63875E-11 0.0491625 0.002704688 79389719112 1.75258E+27 1.4747E+23 1.84338E-32 0.00049771 1.83406E-26 4.30676E+42 3.52128E+30 4.01279E+28 2.30672E-53 4.335E+52 1.23973E+53 2.43E+14 1.97E+29 1.6E+44 1.57514E+44 1 12387225300 4.014E+44 4.01426E+44 2.908E+67 -5.62491E+47 5.62491E+47 4E+17 2.5E-18 1.09E+18 327750000 1.63875E-09 1.554654755 2.7046875 7.93897E+12 5.54215E+28 1.4747E+26 1.84338E-26 15.7390197 1.83406E-26 4.30676E+38 3.52128E+27 4.01279E+28 2.30672E-57 4.335E+56 1.23973E+53 2.43E+17 1.97E+31 5E+45 4.98104E+45 1 12.38723E+12 4.014E+47 4.01426E+47 2.908E+70 -5.62491E+50 5.62491E+50 4E+16 2.5E-17 1.09E+17 3277500000 1.63875E-07 49.1625 2704.6875 7.93897E+14 1.75258E+30 1.4747E+29 1.84338E-20 497711.504 1.83406E-26 4.30676E+34 3.52128E+24 4.01279E+28 2.30672E-61 4.335E+60 1.23973E+53 2.43E+20 1.97E+33 1.6E+47 1.57514E+47 1 12.38723E+15 4.014E+50 4.01426E+50 2.908E+73 -5.62491E+53 5.62491E+53 4E+15 2.5E-16 1.09E+16 32775000000 1.63875E-05 1554.654755 2704687.5 7.93897E+16 5.54215E+31 1.4747E+32 1.84338E-14 1.5739E+10 1.83406E-26 4.30676E+30 3.52128E+21 4.01279E+28 2.30672E-65 4.335E+64 1.23973E+53 2.43E+23 1.97E+35 5E+48 4.98104E+48 1 12.38723E+18 4.014E+53 4.01426E+53 2.908E+76 -5.62491E+56 5.62491E+56 4E+14 2.5E-15 1.09E+15 3.2775E+11 0.00163875 49162.5 2704687500 7.93897E+18 1.75258E+33 1.4747E+35 1.84338E-08 4.9771E+14 1.83406E-26 4.30676E+26 3.52128E+18 4.01279E+28 2.30672E-69 4.335E+68 1.23973E+53 2.43E+26 1.97E+37 1.6E+50 1.57514E+50 1 12.38723E+21 4.014E+56 4.01426E+56 2.908E+79 -5.62491E+59 5.62491E+59 4E+13 2.5E-14 1.09E+14 3.2775E+12 0.163875 1554654.755 2.70469E+12 7.93897E+20 5.54215E+34 1.4747E+38 0.018433759 1.5739E+19 1.83406E-26 4.30676E+22 3.52128E+15 4.01279E+28 2.30672E-73 4.335E+72 1.23973E+53 2.43E+29 1.97E+39 5E+51 4.98104E+51 1 12.38723E+24 4.014E+59 4.01426E+59 2.908E+82 -5.62491E+62 5.62491E+62 4E+12 2.5E-13 1.09E+13 3.2775E+13 16.3875 49162500 2.70469E+15 7.93897E+22 1.75258E+36 1.4747E+41 18433.7594 4.9771E+23 1.83406E-26 4.30676E+18 3.52128E+12 4.01279E+28 2.30672E-77 4.335E+76 1.23973E+53 2.43E+32 1.97E+41 1.6E+53 1.57514E+53 1 1.000000001 2.38723E+27 4.014E+62 4.01426E+62 2.908E+85 -5.62491E+65 5.62491E+65 4E+11 2.5E-12 1.09E+12 3.2775E+14 1638.75 1554654755 2.70469E+18 7.93897E+24 5.54215E+37 1.4747E+44 18433759401 1.5739E+28 1.83406E-26 4.30676E+14 3521280000 4.01279E+28 2.30672E-81 4.335E+80 1.23973E+53 2.43E+35 1.97E+43 5E+54 4.98104E+54 1 1.000000003 2.38723E+30 4.014E+65 4.01426E+65 2.908E+88 -5.62491E+68 5.62491E+68 4E+10 2.5E-11 1.09E+11 3.2775E+15 163875 49162499998 2.70469E+21 7.93897E+26 1.75258E+39 1.4747E+47 1.84338E+16 4.9771E+32 1.83406E-26 43067568254 3521280 4.01279E+28 2.30672E-85 4.335E+84 1.23973E+53 2.43E+38 1.97E+45 1.6E+56 1.57514E+56 1 1.000000007 2.38723E+33 4.014E+68 4.01426E+68 2.908E+91 -5.62491E+71 5.62491E+71 4E+09 2.5E-10 10900000003 3.2775E+16 16387499.99 1.55465E+12 2.70469E+24 7.93897E+28 5.54215E+40 1.4747E+50 1.84338E+22 1.5739E+37 1.83406E-26 4306756.829 3521.280003 4.01279E+28 2.30672E-89 4.335E+88 1.23973E+53 2.43E+41 1.97E+47 5E+57 4.98104E+57 1 1.000000016 2.38723E+36 4.014E+71 4.01426E+71 2.908E+94 -5.62491E+74 5.62491E+74 4E+08 2.5E-09 1090000003 3.2775E+17 1638749992 4.91625E+13 2.70469E+27 7.93897E+30 1.75258E+42 1.4747E+53 1.84338E+28 4.9771E+41 1.83406E-26 430.6756868 3.521280026 4.01279E+28 2.30672E-93 4.335E+92 1.23973E+53 2.43E+44 1.97E+49 1.6E+59 1.57514E+59 1 1.000000037 2.38723E+39 4.014E+74 4.01426E+74 2.908E+97 -5.62491E+77 5.62491E+77 40000000 2.5E-08 109000002.7 3.2775E+18 1.63875E+11 1.55465E+15 2.70469E+30 7.93897E+32 5.54215E+43 1.4747E+56 1.84338E+34 1.5739E+46 1.83406E-26 0.043067573 0.00352128 4.01279E+28 2.30672E-97 4.335E+96 1.23973E+53 2.43E+47 1.97E+51 5E+60 4.98104E+60 1 1.000000088 2.38723E+42 4.014E+77 4.01426E+77 2.91E+100 -5.62491E+80 5.62491E+80 4000000 2.5E-07 10900002.73 3.2775E+19 1.63875E+13 4.91625E+16 2.70469E+33 7.93897E+34 1.75258E+45 1.4747E+59 1.84337E+40 4.9771E+50 1.83406E-26 4.30676E-06 3.52128E-06 4.01279E+28 2.3067E-101 4.34E+100 1.23973E+53 2.43E+50 1.97E+53 1.6E+62 1.57514E+62 0.999999999 1.000000208 2.38722E+45 4.014E+80 4.01426E+80 2.91E+103 -5.62491E+83 5.62491E+83 400000 2.49999E-06 1090002.725 3.27749E+20 1.63874E+15 1.55465E+18 2.70467E+36 7.93893E+36 5.54213E+46 1.47469E+62 1.84335E+46 1.5739E+55 1.83406E-26 4.3068E-10 3.52131E-09 4.01279E+28 2.3067E-105 4.34E+104 1.23973E+53 2.43E+53 1.97E+55 5E+63 4.98102E+63 0.999999996 1.00000049 2.38721E+48 4.014E+83 4.01423E+83 2.91E+106 -5.62487E+86 5.62487E+86 40000 2.49994E-05 109002.725 3.27742E+21 1.63867E+17 4.91607E+19 2.70448E+39 7.93857E+38 1.75252E+48 1.47459E+65 1.8431E+52 4.9766E+59 1.83406E-26 4.30719E-14 3.52154E-12 4.01279E+28 2.307E-109 4.33E+108 1.23973E+53 2.43E+56 1.97E+57 1.6E+65 1.57508E+65 0.999999988 1.000001156 2.38705E+51 4.014E+86 4.01396E+86 2.91E+109 -5.62449E+89 5.62449E+89 3570 0.000280034 9730.975 3.67124E+22 2.05614E+19 1.84306E+21 3.80126E+42 9.96104E+40 6.57027E+49 2.07259E+68 3.64112E+58 2.6224E+64 1.83406E-26 2.73571E-18 2.50548E-15 4.01279E+28 1.4653E-113 6.82E+112 1.23973E+53 3.42E+59 2.47E+59 5.9E+66 5.90507E+66 0.999999958 1.00000284 3.35509E+54 5.642E+89 5.64178E+89 4.09E+112 -7.90544E+92 7.90544E+92 1599 0.000625 4360 8.19375E+22 1.02422E+20 6.14531E+21 4.22607E+43 4.96186E+41 2.19073E+50 2.30422E+69 4.50043E+60 9.7209E+65 1.83406E-26 1.10253E-19 2.25362E-16 4.01279E+28 5.9052E-115 1.69E+114 1.23973E+53 3.8E+60 1.23E+60 2E+67 1.96893E+67 0.999999938 1.000003825 3.73004E+55 6.272E+90 6.27228E+90 4.54E+113 -8.78892E+93 8.78892E+93 1370 0.000729395 3735.975 9.56236E+22 1.39495E+20 7.74761E+21 6.71714E+43 6.75786E+41 2.76193E+50 3.66245E+69 1.13697E+61 1.948E+66 1.83406E-26 5.94375E-20 1.41786E-16 4.01279E+28 3.1835E-115 3.14E+114 1.23973E+53 6.04E+60 1.67E+60 2.5E+67 2.4823E+67 0.999999933 1.000004051 5.92872E+55 9.969E+90 9.9695E+90 7.22E+113 -1.39696E+94 1.39696E+94 1088 0.000918274 2967.525 1.20386E+23 2.21094E+20 1.09442E+22 1.34034E+44 1.0711E+42 3.90145E+50 7.30803E+69 4.52696E+61 5.4906E+66 1.83406E-26 2.36604E-20 7.10566E-17 4.01279E+28 1.2673E-115 7.89E+114 1.23973E+53 1.2E+61 2.65E+60 3.5E+67 3.50645E+67 0.999999924 1.000004412 1.18301E+56 1.989E+91 1.98931E+91 1.44E+114 -2.78748E+94 2.78748E+94 1100 0.000908265 3000.225 1.19074E+23 2.16301E+20 1.07657E+22 1.29699E+44 1.04788E+42 3.83784E+50 7.07167E+69 4.23887E+61 5.2264E+66 1.83406E-26 2.47206E-20 7.34315E-17 4.01279E+28 1.324E-115 7.55E+114 1.23973E+53 1.17E+61 2.6E+60 3.4E+67 3.44928E+67 0.999999925 1.000004394 1.14475E+56 1.925E+91 1.92497E+91 1.39E+114 -2.69733E+94 2.69733E+94 1000 0.000999001 2727.725 1.30969E+23 2.61676E+20 1.24186E+22 1.72582E+44 1.2677E+42 4.42708E+50 9.40983E+69 7.50532E+61 8.0222E+66 1.83406E-26 1.68907E-20 5.51852E-17 4.01279E+28 9.0467E-116 1.11E+115 1.23973E+53 1.55E+61 3.14E+60 4E+67 3.97886E+67 0.999999921 1.000004552 1.52325E+56 2.561E+91 2.56143E+91 1.86E+114 -3.58916E+94 3.58916E+94 900 0.001109878 2455.225 1.45505E+23 3.22986E+20 1.45424E+22 2.36659E+44 1.56471E+42 5.18419E+50 1.29036E+70 1.41132E+62 1.2882E+67 1.83406E-26 1.10869E-20 4.02434E-17 4.01279E+28 5.9382E-116 1.68E+115 1.23973E+53 2.13E+61 3.88E+60 4.7E+67 4.65932E+67 0.999999917 1.000004733 2.08881E+56 3.512E+91 3.51246E+91 2.54E+114 -4.92177E+94 4.92177E+94 400 0.002493766 1092.725 3.26933E+23 1.63059E+21 4.89787E+22 2.6845E+45 7.89943E+42 1.74603E+51 1.4637E+71 1.81597E+64 4.9215E+68 1.83406E-26 4.34999E-22 3.54776E-18 4.01279E+28 2.3299E-117 4.29E+116 1.23973E+53 2.41E+62 1.96E+61 1.6E+68 1.56925E+68 0.999999875 1.000006388 2.36941E+57 3.984E+92 3.9843E+92 2.89E+115 -5.58293E+95 5.58293E+95 40 0.024390244 111.725 3.19756E+24 1.55979E+23 1.49813E+24 2.51157E+48 7.55643E+44 5.34063E+52 1.36941E+74 1.58954E+70 1.4084E+73 1.83406E-26 4.75385E-26 3.79203E-21 4.01279E+28 2.5462E-121 3.93E+120 1.23973E+53 2.26E+65 1.87E+63 4.8E+69 4.79992E+69 0.99999961 1.000014829 2.21678E+60 3.728E+95 3.72764E+95 2.7E+118 -5.22329E+98 5.22329E+98 10 0.090909091 29.975 1.19182E+25 2.16694E+24 1.07804E+25 1.30053E+50 1.04978E+46 3.84308E+53 7.09097E+75 4.26204E+73 5.2478E+75 1.83406E-26 2.46309E-28 7.32316E-23 4.01279E+28 1.3192E-123 7.58E+122 1.23973E+53 1.17E+67 2.6E+64 3.5E+70 3.45399E+70 0.999999247 1.000024059 1.14788E+62 1.93E+97 1.93022E+97 1.4E+120 -2.7047E+100 2.7047E+100 9 0.1 27.25 1.311E+25 2.622E+24 1.24372E+25 1.731E+50 1.27024E+46 4.43372E+53 9.43808E+75 7.55047E+73 8.0584E+75 1.83406E-26 1.68233E-28 5.502E-23 4.01279E+28 9.0106E-124 1.11E+123 1.23973E+53 1.56E+67 3.15E+64 4E+70 3.98483E+70 0.99999921 1.000024916 1.52782E+62 2.569E+97 2.56913E+97 1.86E+120 -3.5999E+100 3.5999E+100 8 0.111111111 24.525 1.45667E+25 3.23704E+24 1.45667E+25 2.37449E+50 1.56819E+46 5.19283E+53 1.29466E+76 1.42075E+74 1.2947E+76 1.83406E-26 1.10377E-28 4.01096E-23 4.01279E+28 5.9119E-124 1.69E+123 1.23973E+53 2.13E+67 3.89E+64 4.7E+70 4.66709E+70 0.999999167 1.000025898 2.09578E+62 3.524E+97 3.52418E+97 2.55E+120 -4.9382E+100 4.9382E+100 7 0.125 21.8 1.63875E+25 4.09688E+24 1.73816E+25 3.38086E+50 1.98474E+46 6.19631E+53 1.84338E+76 2.88027E+74 2.1996E+76 1.83406E-26 6.89081E-29 2.81702E-23 4.01279E+28 3.6908E-124 2.71E+123 1.23973E+53 3.04E+67 4.92E+64 5.6E+70 5.56897E+70 0.999999117 1.000027042 2.98403E+62 5.018E+97 5.01783E+97 3.63E+120 -7.0311E+100 7.0311E+100 6 0.142857143 19.075 1.87286E+25 5.35102E+24 2.12362E+25 5.04665E+50 2.59232E+46 7.57044E+53 2.75163E+76 6.41779E+74 4.0115E+76 1.83406E-26 4.03927E-29 1.88719E-23 4.01279E+28 2.1635E-124 4.62E+123 1.23973E+53 4.54E+67 6.42E+64 6.8E+70 6.80398E+70 0.999999056 1.000028399 4.4543E+62 7.49E+97 7.49017E+97 5.43E+120 -1.0495E+101 1.0495E+101 5 0.166666667 16.35 2.185E+25 7.28333E+24 2.67607E+25 8.01389E+50 3.52843E+46 9.53985E+53 4.36948E+76 1.61833E+75 8.0273E+76 1.83406E-26 2.1803E-29 1.18843E-23 4.01279E+28 1.1678E-124 8.56E+123 1.23973E+53 7.2E+67 8.74E+64 8.6E+70 8.57399E+70 0.99999898 1.000030051 7.07326E+62 1.189E+98 1.18941E+98 8.61E+120 -1.6666E+101 1.6666E+101 4 0.2 13.625 2.622E+25 1.0488E+25 3.51778E+25 1.3848E+51 5.08094E+46 1.25405E+54 7.55047E+76 4.8323E+75 1.8234E+77 1.83406E-26 1.05145E-29 6.8775E-24 4.01279E+28 5.6316E-125 1.78E+124 1.23973E+53 1.24E+68 1.26E+65 1.1E+71 1.12708E+71 0.999998883 1.000032127 1.22226E+63 2.055E+98 2.0553E+98 1.49E+121 -2.88E+101 2.88E+101 3 0.25 10.9 3.2775E+25 1.63875E+25 4.91625E+25 2.70469E+51 7.93897E+46 1.75258E+54 1.4747E+77 1.84338E+76 4.9771E+77 1.83406E-26 4.30676E-30 3.52128E-24 4.01279E+28 2.3067E-125 4.34E+124 1.23973E+53 2.43E+68 1.97E+65 1.6E+71 1.57514E+71 0.999998751 1.000034862 2.38723E+63 4.014E+98 4.01426E+98 2.91E+121 -5.6249E+101 5.6249E+101 2 0.333333333 8.175 4.37E+25 2.91333E+25 7.56906E+25 6.41111E+51 1.41137E+47 2.69828E+54 3.49559E+77 1.03573E+77 1.8164E+78 1.83406E-26 1.36268E-30 1.48554E-24 4.01279E+28 7.2986E-126 1.37E+125 1.23973E+53 5.76E+68 3.5E+65 2.4E+71 2.42509E+71 0.999998558 1.000038732 5.65861E+63 9.515E+98 9.51528E+98 6.89E+121 -1.3333E+102 1.3333E+102 1 0.5 5.45 6.555E+25 6.555E+25 1.39053E+26 2.16375E+52 3.17559E+47 4.95705E+54 1.17976E+78 1.17976E+78 1.1262E+79 1.83406E-26 2.69172E-31 4.4016E-25 4.01279E+28 1.4417E-126 6.94E+125 1.23973E+53 1.94E+69 7.87E+65 4.5E+71 4.45518E+71 0.999998234 1.000044918 1.90978E+64 3.21E+99 3.2114E+99 2.33E+122 -4.4999E+102 4.4999E+102 0 1 2.725 1.311E+26 2.622E+26 3.933E+26 1.731E+53 1.27024E+48 1.40207E+55 9.43808E+78 7.55047E+79 2.5483E+80 1.83406E-26 1.68233E-32 5.502E-26 4.01279E+28 9.0106E-128 1.11E+127 1.23973E+53 1.56E+70 3.15E+66 1.3E+72 1.26012E+72 0.999997502 1.000057838 1.52782E+65 2.57E+100 2.5691E+100 1.86E+123 -3.5999E+103 3.5999E+103 87 References [1] Lynden-Bell, D., Wood, R.: The gravothermal catastrophe in isothermal spheres and the onset of red-giant structure for stellar systems. Mon. Not. R. Astron. Soc. 138(4), 495–525 (1968) https://doi.org/10.1093/mnras/138.4.495 [2] Sugimoto, D., Eriguchi, Y., Hachisu, I.: Gravothermal aspects in evolution of the stars and the universe. Prog. Theor. Phys. Suppl. 70, 154–180 (1981) https: //doi.org/10.1143/PTPS.70.154 [3] Hayward, S.A.: Formation and evaporation of nonsingular black holes. Phys. Rev. Lett. 96(3), 031103 (2006) https://doi.org/10.1103/PhysRevLett.96.031103 [4] Bardeen, J.M.: Non-singular general-relativistic gravitational collapse. In: Abstracts of Contributed Papers for the 5th International Conference on Gravitation and the Theory of Relativity (GR5), Tbilisi, USSR, p. 174 (1968) [5] Frolov, V.P.: Notes on non-singular models of black holes. Universe 2(3), 20 (2016) https://doi.org/10.3390/universe2030020 [6] Kawai, H., Yokokura, Y.: A model of black hole evaporation and entropy. Universe 4(12), 142 (2018) https://doi.org/10.3390/universe4120142 [7] Bekenstein, J.D.: Black holes and entropy. Phys. Rev. D 7(8), 2333–2346 (1973) https://doi.org/10.1103/PhysRevD.7.2333 [8] Hawking, S.W.: Particle creation by black holes. Commun. Math. Phys. 43(3), 199–220 (1975) https://doi.org/10.1007/BF02345020 [9] Hooft, G.: Dimensional reduction in quantum gravity. arXiv preprint (1993) arXiv:gr-qc/9310026 [gr-qc] [10] Susskind, L.: The world as a hologram. J. Math. Phys. 36(11), 6377–6396 (1995) https://doi.org/10.1063/1.531249 [11] Dymnikova, I.: Vacuum nonsingular black hole. Gen. Rel. Grav. 24(3), 235–242 (1992) https://doi.org/10.1007/BF00760226 [12] Jacobson, T.: Thermodynamics of spacetime: The einstein equation of state. Phys. Rev. Lett. 75(7), 1260–1263 (1995) https://doi.org/10.1103/PhysRevLett. 75.1260 [13] Fischler, W., Susskind, L.: Holography and cosmology. arXiv preprint (1998) arXiv:hep-th/9806039 [hep-th] [14] Carballo-Rubio, R., Di Filippo, F., Liberati, S.: Thermodynamic stability of regular black holes. Phys. Rev. D 107(6), 064015 (2023) https://doi.org/10.1103/ PhysRevD.107.064015 88 [15] Egan, C.A., Lineweaver, C.H.: A larger estimate of the universe’s entropy. Astrophys. J. 710(2), 1825–1834 (2010) https://doi.org/10.1088/0004-637X/710/2/ 1825 [16] Kawamura, S., et al.: Current status of space gravitational wave antenna decigo and b-decigo. Prog. Theor. Exp. Phys. 2021(5) (2021) https://doi.org/10.1093/ ptep/ptab019 [17] Quevedo, H., Quevedo, M.N., Valdez, E.A.: Geometrothermodynamics of 3d regular black holes. Entropy 26(6), 457 (2024) https://doi.org/10.3390/e26060457 [18] Thorlacius, L.: Black holes and the holographic principle. In: Horowitz, G.T. (ed.) Black Holes in Higher Dimensions, pp. 373–393. Cambridge University Press, ??? (2012). https://doi.org/10.1017/CBO9781139003507.013 [19] Wald, R.M.: The thermodynamics of black holes. Living Rev. Rel. 4(1), 6 (2001) https://doi.org/10.12942/lrr-2001-6 [20] Ayon-Beato, E., Garcia, A.: Regular black hole in general relativity coupled to nonlinear electrodynamics. Phys. Rev. Lett. 80, 5056–5059 (1998) https://doi. org/10.1103/PhysRevLett.80.5056 [21] Padmanabhan, T.: Thermodynamical aspects of gravity: New insights. Rep. Prog. Phys. 73(4), 046901 (2010) https://doi.org/10.1088/0034-4885/73/4/046901 [22] Bronnikov, K.A.: Regular electrically charged black holes and monopoles from nonlinear electrodynamics. Phys. Rev. D 63(4), 044005 (2001) https://doi.org/ 10.1103/PhysRevD.63.044005 [23] Penrose, R.: Singularities and time-asymmetry. In: Hawking, S.W., Israel, W. (eds.) General Relativity: An Einstein Centenary Survey, pp. 581–638. Cambridge University Press, ??? (1979) [24] Ansoldi, S.: Spherically symmetric black holes with a regular center: a review of existing models and results. arXiv preprint (2008) arXiv:0802.0330 [gr-qc] [25] Bousso, R.: The holographic principle. Rev. Mod. Phys. 74, 825–874 (2002) https: //doi.org/10.1103/RevModPhys.74.825 [26] Verlinde, E.: On the origin of gravity and the laws of newton. J. High Energy Phys. 2011(4), 029 (2011) https://doi.org/10.1007/JHEP04(2011)029 [27] Myung, Y.S.: Black hole spectroscopy via adiabatic invariance. Physics Letters B645(5-6), 369–371 (2007) https://doi.org/10.1016/j.physletb.2007.01.011 [28] Planck Collaboration, Aghanim, N., et al.: Planck 2018 results. vi. cosmological parameters. Astron. Astrophys. 641, 6 (2018) https://doi.org/10.1051/ 0004-6361/201833910 1807.06209 89