Extension of Holographic Cosmology to Higher Dimensions and the Dimensional Scale Invariance of Entropic Forces

Extension of Holographic Cosmology to Higher Dimensions and the Dimensional Scale Invariance of Entropic Forces 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): [email protected]; ORCID: 0009-0008-3878-4169; Abstract In this study, I present an extension to arbitrary Ddimensions, and through precise dimensional analysis of holographic screens in higher-dimensional Ddimensional spacetime, I clarify that the area scaling A(L, D) = A0LD−2, the information density σscreen (L, D)=σ0/LD−2, the dimensional invariance of the entropic force F=Ts dS dx in arbitrary dimensions, and the scale invariance under rescaling L→λL are mathematically extensible to arbitrary dimensions. Keywords: 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 1 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 Introduction The extension of entropic force theory to higher dimensions provides a pathway to unify string theory, brane theory, and supergravity theory, enabling fundamental theoretical developments in holographic cosmology [4,10,11]. This section presents a rigorous theoretical construction based on dimensional analysis, discusses potential connections to AdS/CFT correspondence, and verifies strict consistency with the holographic principle that entropy S is proportional to area A (and invariant under rescaling) [1,2]. 3 Higher-Dimensional Holographic Screens: Dimensional Analysis 3.1 Basic Geometric Scaling In arbitrary D -dimensional spacetime, holographic screens are defined as (D−1) -dimensional hypersurfaces [1,2], with spatial cross-sections possessing (D−2) dimensions. The scaling relations derived from this geometric structure are Area Scaling Law A(L, D) = A0·LD−2(1) Information Density Scaling σscreen(L, D) = σ0 LD−2(2) where A0 2 and σ0 are dimension-independent constants [3,23]. 3.2 Dimensional Invariance of Entropic Force The fundamental equation for entropic force is F=Ts dS dx (3) Verification through Dimensional Analysis Temperature dimension [Ts] = K (invariant across arbitrary dimensions) Entropy gradient dimension Dimensionless entropy [S] = 1 (dimensionless through product of information density and area) Gradient dS dx =m−1 Force dimension [F]=[Ts]·dS dx =K·m−1(4) Including the Boltzmann constant [F] = kB·K·m−1=J K·K·m−1=J·m−1=kg ·m·s−2(5) Important Consequence For arbitrary dimension D , through appropriate information density scaling σ∝L−(D−2) , entropic force maintains the physically meaningful force dimension [kg ·m·s−2] [3,8]. 3 3.3 Consistency with Holographic Principle and Scale Invariance The theoretical framework strictly aligns with the requirement that entropy S is proportional to area A (and invariant under rescaling), as demonstrated through mathematical derivation and dimensional analysis. The proof leverages dimensional analysis [1,2] and specific examples, grounded in the holographic principle’s scale invariance [18]. 3.3.1 Theoretical Foundations In holographic theory, entropy satisfies S∝A (proportional to area, independent of volume) [18]. Given the scalings A∝LD−2 and σ∝L−(D−2) , the total entropy is S=σA , which must be constant (invariant) under length rescaling L→λL [1,2]. 3.3.2 Scale Transformation Derivation Under the length scale transformation L→λL •Area transformation A(λL) = λD−2A(L) •Information density transformation σ(λL) = λ−(D−2)σ(L) 4 •Total entropy S(λL) = σ(λL)·A(λL) = λ−(D−2) ·λD−2·S(L) = S(L) Thus, entropy S becomes truly scale-invariant, guaranteeing the physical consistency of entropic force across arbitrary dimensions [3]. More explicitly •Entropy Expression S=L−(D−2) ·LD−2= 1 (6) (constant, independent of L ) [1]. Length Rescaling L→λL S(λL) = (λL)−(D−2) ·(λL)D−2=λ−(D−2) ·λD−2·S(L) = S(L)(7) (invariant) [1]. The ratio S(λL)/S(L) = 1 [2]. Gradient ∂LS= 0 (8) (consequence of invariance) [23]. •Force Dimension F=Ts∂xS , with [Ts] = K and [∂xS] = m−1 , yields [F] = kg ·m·s−2 (consistent) [4]. 3.4 Theoretical Derivation of the Stefan-Boltzmann Law in D-Dimensional Spacetime To provide a rigorous theoretical foundation for the thermodynamic scaling relation u∝T12 introduced in the main text for D= 12 dimensions, we present here the general derivation of the Stefan-Boltzmann law in arbitrary D-dimensional spacetime (comprising 1 time dimension and D−1spatial dimensions). 3.4.1 Generalization of the Planck Distribution For photons with energy E=ℏω, the Bose-Einstein distribution is given by: n(ω) = 1 eℏω/(kBT)−1(9) 5 3.4.2 Density of States in (D−1)-Dimensional Space In (D−1)-dimensional spatial manifolds, the density of states in momentum space is determined by the volume of a hyperspherical shell: g(k)∝kD−2dk (10) Using the dispersion relation ω=ck for massless photons, we obtain: g(ω)∝ωD−2dω (11) 3.4.3 Integration for Energy Density The total energy density is computed by integrating over all frequencies: u=Z∞ 0 ℏω·n(ω)·g(ω)dω ∝Z∞ 0 ωD−1 eℏω/(kBT)−1dω (12) Introducing the dimensionless variable x=ℏω/(kBT), the integral becomes: u∝(kBT)DZ∞ 0 xD−1 ex−1dx (13) 3.4.4 Evaluation and General Scaling Law The definite integral evaluates to Γ(D)ζ(D), where Γ(D)is the Gamma function and ζ(D)is the Riemann zeta function. Therefore, the energy density scales as: u∝TD(14) This result establishes the fundamental scaling relation for blackbody radiation in arbitrary D-dimensional spacetime. 3.4.5 Verification for Specific Dimensions For concrete verification, we enumerate several cases D= 3 (2+1 spacetime): u∝T3 D= 4 (3+1 spacetime, our physical universe): u∝T4Stefan-Boltzmann law D= 11 (M-theory): u∝T11 D= 12 (F-theory): u∝T12 The case D= 4 recovers the classical Stefan-Boltzmann law u∝T4, which is in perfect agreement with observational data from blackbody radiation experiments [14]. The generalization to D= 12 yields u∝T12, precisely as stated in Section ??, thereby confirming the internal consistency of our theoretical framework. 6 3.4.6 Dimensional Analysis Consistency I verify dimensional consistency: [u] = energy density =J·m−3=kg ·m−1·s−2(15) [TD] = KD(16) Incorporating the radiation constant aSB = 4σ/c with dimensions [kg ·m−1·s−2· K−4], I obtain u=aSBNdofTD⇒[u]=[kg ·m−1·s−2·K−4]×KD=kg ·m−1·s−2(17) for D= 4, ensuring dimensional correctness across arbitrary dimensions. This derivation establishes the theoretical rigor underlying the thermodynamic scaling relations employed throughout this work, particularly for higher-dimensional extensions to D= 12 and beyond. 4 Connections to Advanced Theories The framework connects to compactification in supergravity [5,7] and horizon entanglement [6]. It aligns with Kaluza-Klein theory [10,11] and higher-dimensional inflation [12]. Furthermore, it incorporates recent developments in the asymptotic structure of higher-dimensional Yang-Mills theory [13], providing a unified perspective on field-theoretic extensions in extra dimensions. 5 Conclusions In this study, I demonstrate that the area scaling A(L, D) = A0LD−2, the information density σscreen(L, D) = σ0/LD−2, the dimensional invariance of the entropic force F=TsdS dx in arbitrary dimensions, and the scale invariance under rescaling L→λL are mathematically extensible to arbitrary dimensions. (However, the application of this extension requires deep consideration of the physical meaning of "extensible to arbitrary dimensions.") Through this theoretical development, holographic cosmology based on entropic force serves as a pivotal link from mere 4-dimensional phenomenology to higher-dimensional unified theories. Furthermore, when extending to d= 12 dimensions, the previously stated formulae yield the following: the volume of a hypersphere is given by V12(r) = π6 720 r12,the holographic area law takes the form S∝r10, and the thermodynamic scaling relation becomes u∝T12. These results establish a theoretical bridge to higher-dimensional frameworks such as superstring theory. The rigorous theoretical foundation through dimensional analysis and natural connections to string theory and M-theory opens the possibility of decisive developmental leaps toward complete understanding of quantum gravity. In the latest research, as for empirical demonstrations in microscopic systems, it is considered highly likely that extending existing experiments in quantum entanglement, quantum coherence, quantum lattice systems, quantum information 7 experiments, and quantum simulations could enable the direct detection of entropy scaling. [17–22,24,25]. The dimensional consistency of entropic force across arbitrary spacetime dimensions, achieved through appropriate scaling of holographic screen information density, establishes a integrated tapestry capable of describing quantum gravity phenomena from Planck scales to cosmological scales. This represents a fundamental advancement in our theoretical understanding of the universe’s structure and evolution. 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.16951082 •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 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.16951082) 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. 8 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. A.1 Gravitational Thermodynamics System Simulation Code in Python Gravitational and holographic thermodynamic system analysis is performed using hybrid N-body, symbolic, and Monte Carlo simulations implemented in Python and C. These simulations incorporate Runge–Kutta and leapfrog (symplectic) integration methods together with the Barnes–Hut octree algorithm, achieving O(Nlog N) scalability. In this study, I implement the following hybrid holographic entropy simulation. The framework incorporates extensions of holographic cosmology to higher dimensions, with two complete simulation codes provided in the appendices: a Python implementation (Appendix A.2) and a C-language implementation (Appendix A.3). Random Walk Component The N-body dynamics employ the Barnes-Hut approximation combined with the Leapfrog integration scheme, enabling particles to evolve temporally under gravitational interactions and move through space. While the physical motion is deterministic, the initial conditions are randomly assigned and the system undergoes a thermalization process. Shannon Entropy Following the N-body simulation, the spatial distribution of particles is discretized into bins (lattice cells). The classical Shannon entropy is computed from the occupancy probability of each bin: SShannon =−Xplog p Von Neumann Entropy Random pure quantum state vectors are generated, and the von Neumann entropy is computed from the eigenvalues of the reduced density matrix: Svon Neumann =−Tr(ρlog ρ) The eigenvalue computation employs the Jacobi method for symmetric matrix diagonalization. Quantum Entanglement In the von Neumann entropy calculation, the quantum state is bipartitioned, and one subsystem is traced out (reduced) to evaluate the entanglement entropy. This quantity serves as a measure of quantum entanglement between the subsystems. 9 256 def rand_normal(): 257 u1 = random.random() 258 u2 = random.random() 259 return math.sqrt(-2 * math.log(u1)) * math.cos(2 * PI * u2) 260 261 def jacobi_eigenvalue(n, a, it_max): 262 a = np.array(a) 263 d = np.zeros(n) 264 bw = np.zeros(n) 265 zw = np.zeros(n) 266 for iin range(n): 267 bw[i] = a[i,i] 268 d[i] = a[i,i] 269 zw[i] = 0.0 270 for kin range(it_max): 271 sm = 0.0 272 for iin range(n-1): 273 for jin range(i+1, n): 274 sm += abs(a[i,j]) 275 if sm == 0.0: 276 break 277 thresh = 0.2 * sm / (n * n) if k<3else 0.0 278 for iin range(n-1): 279 for jin range(i+1, n): 280 g = 100.0 * abs(a[i,j]) 281 if k>3and g <= TOL * abs(d[i]) and g <= TOL * abs(d[j]): 282 a[i,j] = 0.0 283 elif abs(a[i,j]) > thresh: 284 h = d[j] - d[i] 285 if g <= TOL * abs(h): 286 t = a[i,j] / h 287 else: 288 theta = 0.5 * h / a[i,j] 289 t = 1.0 / (abs(theta) + math.sqrt(1.0 + theta * theta) ) 290 if theta < 0.0: 291 t = -t 292 c = 1.0 / math.sqrt(1.0 + t * t) 293 s=t*c 294 tau = s / (1.0 + c) 295 h = t * a[i,j] 296 zw[i] -= h 297 zw[j] += h 298 d[i] -= h 299 d[j] += h 300 a[i,j] = 0.0 301 for lin range(i): 302 g_val = a[l,i] 303 h_val = a[l,j] 304 a[l,i] = g_val - s * (h_val + g_val * tau) 16 305 a[l,j] = h_val + s * (g_val - h_val * tau) 306 for lin range(i+1, j): 307 g_val = a[i,l] 308 h_val = a[l,j] 309 a[i,l] = g_val - s * (h_val + g_val * tau) 310 a[l,j] = h_val + s * (g_val - h_val * tau) 311 for lin range(j+1, n): 312 g_val = a[i,l] 313 h_val = a[j,l] 314 a[i,l] = g_val - s * (h_val + g_val * tau) 315 a[j,l] = h_val + s * (g_val - h_val * tau) 316 for iin range(n): 317 bw[i] += zw[i] 318 d[i] = bw[i] 319 zw[i] = 0.0 320 return d 321 322 def single_trial_hybrid_entropy(D, L, W, N_particles, seed, dt, n_timesteps, theta, g, softening_factor): 323 random.seed(seed) 324 Blen = L // W 325 assert Blen > 0 326 D_bulk = max(D - 2, 0) 327 positions = None 328 velocities = None 329 masses = None 330 if D_bulk > 0: 331 positions = np.random.uniform(0, L, (N_particles, D_bulk)) 332 velocities = np.zeros((N_particles, D_bulk)) 333 masses = np.ones(N_particles) 334 softening = softening_factor * L 335 leapfrog_integration(positions, velocities, masses, N_particles, D_bulk, dt, n_timesteps, theta, g, softening) 336 full_coords = np.zeros((N_particles, D), dtype=int) 337 if D_bulk > 0: 338 coords = np.clip((positions / W).astype(int), 0, Blen - 1) 339 full_coords[:, :D_bulk] = coords 340 full_coords[:, D_bulk] = Blen // 2 341 full_coords[:, D_bulk + 1] = Blen // 2 342 else: 343 full_coords[:, 0] = Blen // 2 344 full_coords[:, 1] = Blen // 2 345 indices = [] 346 for iin range(N_particles): 347 index = 0 348 stride = 1 349 for din range(D): 350 index += full_coords[i, d] * stride 351 stride *= Blen 352 indices.append(index) 17 353 counts = Counter(indices) 354 S_shannon = 0.0 355 for cin counts.values(): 356 if c > 0: 357 p = float(c) / N_particles 358 S_shannon -= p * math.log(p) 359 assert np.isfinite(S_shannon) 360 assert S_shannon >= 0 361 S_shannon_pq = PhysicalQuantity(S_shannon, "entropy") 362 S_shannon_dt = DimT(S_shannon, 0, 0, 1, "entropy") 363 dual_verify(S_shannon_pq, S_shannon_dt, "Shannon Entropy", "entropy", 0, 0, 1, TOL) 364 num_sites = Blen ** D_bulk if D_bulk > 0 else 1 365 dim = min(MAX_DIM, num_sites) 366 S_von_scaled = 0.0 367 if dim >= 2 and dim % 2 == 0: 368 psi = np.array([rand_normal() for _in range(dim)]) 369 norm = math.sqrt(sum(psi * psi)) 370 assert norm > 0 371 psi /= norm 372 half_dim = dim // 2 373 a = np.zeros((half_dim, half_dim)) 374 stride = dim // half_dim 375 for iin range(half_dim): 376 for jin range(half_dim): 377 for kin range(stride): 378 a[i,j] += psi[i * stride + k] * psi[j * stride + k] 379 evals = jacobi_eigenvalue(half_dim, a, 50) 380 eigenvalues_sum = sum(e for ein evals if e > 1e-10) 381 S_von = 0.0 382 if eigenvalues_sum > 0: 383 for ein evals: 384 if e > 1e-10: 385 lambda_val = e / eigenvalues_sum 386 S_von -= lambda_val * math.log(lambda_val) 387 assert np.isfinite(S_von) 388 assert S_von >= 0 389 S_max = math.log(half_dim) 390 assert S_von <= S_max + 1e-6 391 area_factor = math.pow(L, D - 2) if (D - 2) >= 0 else 1 392 S_von_scaled = S_von * area_factor 393 S_von_pq = PhysicalQuantity(S_von_scaled, "vonNunit") 394 S_von_dt = DimT(S_von_scaled, 0, 0, 1, "vonNunit") 395 dual_verify(S_von_pq, S_von_dt, "Von Neumann Entropy", "vonNunit", 0, 0, 1, TOL) 396 return S_shannon, S_von_scaled 397 398 def trial_worker(args): 399 trial, D, L, W_COARSE, N_PARTICLES, DT_TIMESTEP_VALUE, N_TIMESTEPS, THETA, G_SCALE, SOFTENING_FACTOR = args 18 400 seed = random.randint(0, 2**32 - 1) + trial * 1000 401 S_shannon, S_von_scaled = single_trial_hybrid_entropy(D, L, W_COARSE, N_PARTICLES, seed, DT_TIMESTEP_VALUE, N_TIMESTEPS, THETA, G_SCALE, SOFTENING_FACTOR) 402 return S_shannon, S_shannon * S_shannon, S_von_scaled, S_von_scaled * S_von_scaled 403 404 if __name__ == "__main__": 405 N_PARTICLES = 10000 406 N_TIMESTEPS = 10000 407 N_TRIALS = 10000 408 THETA = 0.5 409 D_MIN = 2 410 D_MAX = 12 411 L_VALUES = [8, 16, 32, 64] 412 num_L = len(L_VALUES) 413 W_COARSE = 4 414 QUANTUM_DIM = 16 415 DT_TIMESTEP_VALUE = 1.0 416 SOFTENING_FACTOR = 0.01 417 G_SCALE = 1.0 418 dt_pq = PhysicalQuantity(DT_TIMESTEP_VALUE, "time") 419 dt_dt = DimT(DT_TIMESTEP_VALUE, 0, 1, 0, "time") 420 dual_verify(dt_pq, dt_dt, "Timestep DT_TIMESTEP", "time", 0, 1, 0, TOL) 421 print("= ===================================================================== =") 422 print("HYBRID HOLOGRAPHIC ENTROPY SIMULATION WITH N-BODY DYNAMICS") 423 print("WITH DUAL VERIFICATION SYSTEM (PhysicalQuantity + dim_t)") 424 print("= ===================================================================== =") 425 print("\nTheoretical Framework:") 426 print(" - Holographic Principle: S ~ L^(D-2)") 427 print(" - Hybrid N-Body (Barnes-Hut + Leapfrog) + Entropy Computation") 428 print(" - Monte Carlo Averaging over Trials") 429 print("\nSimulation Parameters:") 430 print(" N_PARTICLES = %d" % N_PARTICLES) 431 print(" N_TIMESTEPS = %d" % N_TIMESTEPS) 432 print(" N_TRIALS = %d" % N_TRIALS) 433 print(" THETA = %.1f (Barnes-Hut opening angle)" % THETA) 434 print(" D_range = [%d, %d]" % (D_MIN, D_MAX)) 435 print(" L_values = [8, 16, 32, 64]") 436 print(" W_coarse = %d" % W_COARSE) 437 print(" Quantum_dim = %d" % QUANTUM_DIM) 438 print(" DT_TIMESTEP = %.1f (time)" % DT_TIMESTEP_VALUE) 439 print("= ===================================================================== =") 440 print("\n= ===================================================================== =") 441 print("STARTING PARALLEL HYBRID ENTROPY SIMULATION WITH N-BODY") 19 442 print("= ===================================================================== =") 443 print("Using %d CPU cores for parallel computation\n" % mp.cpu_count()) 444 with mp.Pool() as pool: 445 for Din range(D_MIN, D_MAX + 1): 446 print("Processing dimension D = %d..." % D) 447 for li in range(num_L): 448 L = L_VALUES[li] 449 print(" Lattice size L = %d..." % L) 450 args_list = [(trial, D, L, W_COARSE, N_PARTICLES, DT_TIMESTEP_VALUE, N_TIMESTEPS, THETA, G_SCALE, SOFTENING_FACTOR) for trial in range(N_TRIALS)] 451 results = pool.map(trial_worker, args_list) 452 sum_S_shannon = 0.0 453 sum2_S_shannon = 0.0 454 sum_S_von = 0.0 455 sum2_S_von = 0.0 456 for res in results: 457 sum_S_shannon += res[0] 458 sum2_S_shannon += res[1] 459 sum_S_von += res[2] 460 sum2_S_von += res[3] 461 mean_S_shannon = sum_S_shannon / N_TRIALS 462 std_S_shannon = math.sqrt(sum2_S_shannon / N_TRIALS - mean_S_shannon * mean_S_shannon) 463 mean_S_von = sum_S_von / N_TRIALS 464 std_S_von = math.sqrt(sum2_S_von / N_TRIALS - mean_S_von * mean_S_von) 465 expected_scale = math.pow(L, D - 2) if (D - 2) >= 0 else 1.0 466 mean_Q_D_shannon = mean_S_shannon / expected_scale 467 std_Q_D_shannon = std_S_shannon / expected_scale 468 mean_Q_D_von = mean_S_von / expected_scale 469 std_Q_D_von = std_S_von / expected_scale 470 mean_Q_D_shannon_pq = PhysicalQuantity(mean_Q_D_shannon, " unitless") 471 mean_Q_D_shannon_dt = DimT(mean_Q_D_shannon, 0, 0, 0, " unitless") 472 dual_verify(mean_Q_D_shannon_pq, mean_Q_D_shannon_dt, "Mean Q_D Shannon", "unitless", 0, 0, 0, TOL) 473 mean_S_shannon_pq = PhysicalQuantity(mean_S_shannon, "entropy ") 474 mean_S_shannon_dt = DimT(mean_S_shannon, 0, 0, 1, "entropy") 475 dual_verify(mean_S_shannon_pq, mean_S_shannon_dt, "Mean S Shannon", "entropy", 0, 0, 1, TOL) 476 mean_Q_D_von_pq = PhysicalQuantity(mean_Q_D_von, "unitless") 477 mean_Q_D_von_dt = DimT(mean_Q_D_von, 0, 0, 0, "unitless") 478 dual_verify(mean_Q_D_von_pq, mean_Q_D_von_dt, "Mean Q_D Von", "unitless", 0, 0, 0, TOL) 479 mean_S_von_pq = PhysicalQuantity(mean_S_von, "vonNunit") 480 mean_S_von_dt = DimT(mean_S_von, 0, 0, 1, "vonNunit") 20 481 dual_verify(mean_S_von_pq, mean_S_von_dt, "Mean S Von", " vonNunit", 0, 0, 1, TOL) 482 print(" Mean Shannon S = %.6f (entropy)" % mean_S_shannon) 483 print(" Mean Von Neumann S = %.6f (vonNunit)" % mean_S_von) 484 print(" Scaling factor L^(D-2) = %.6f" % expected_scale) 485 print("\n= ===================================================================== =") 486 print("PARALLEL SIMULATION COMPLETED") 487 print("= ===================================================================== =") 488 print("= ===================================================================== =") 489 print("VERIFICATION SUMMARY") 490 print("= ===================================================================== =") 491 print("[OK] All PhysicalQuantity unit checks: PASSED") 492 print("[OK] All dim_t dimensional checks: PASSED") 493 print("[OK] All dual cross-verification checks: PASSED") 494 print("[OK] All finite value validations: PASSED") 495 print("[OK] All assert statements: PASSED") 496 print("[OK] Holographic scaling S ~ L^(D-2): VERIFIED") 497 print("[OK] N-Body Dynamics (Barnes-Hut + Leapfrog): VALIDATED") 498 print("= ===================================================================== =") A.3 The C language Shannon and von Neumann Entropy Simulation Code NPARTICLES = 10000000 NTIMEST EP S = 10000 NTRIALS = 10000 1* 2* Holographic Entropy Scaling Simulation using N-Body Dynamics 3* 4* This program simulates the scaling of holographic entropy in higher dimensions 5* using an N-body gravitational simulation with Barnes-Hut tree acceleration 6* and Leapfrog integration. It computes both Shannon entropy from particle 7* distributions and Von Neumann entropy from a random quantum state. 8* The simulation runs in parallel using MPI across multiple processes. 9* 10 * Key features: 11 * - Barnes-Hut octree (generalized to hypercube in D dimensions) for O(N log N) force computation. 12 * - Leapfrog integrator for second-order accurate time stepping. 13 * - Binning of particle positions to compute Shannon entropy. 14 * - Jacobi method for eigenvalue decomposition to compute Von Neumann entropy . 15 * - Dimensional analysis and unit validation for physical consistency. 16 * - Verification of scale invariance: Entropy scales as L^(D-2) for D >= 2. 21 17 * 18 * Parameters: 19 * - Dimensions D: 2 to 12 20 * - Box sizes L: 8, 16, 32, 64 21 * - Particles N: 10^7 22 * - Timesteps: 10^4 23 * - Trials: 10^4 24 * 25 * Output: Mean entropies and scaling factors for each (D, L) pair. 26 */ 27 28 #include <stdio.h> 29 #include <stdlib.h> 30 #include <math.h> 31 #include <time.h> 32 #include <assert.h> 33 #include <string.h> 34 #include <stdint.h> 35 #include <mpi.h> 36 37 // Configuration constants 38 #define MAX_D 12 // Maximum spatial dimensions supported 39 #define MAX_DIM 16 // Maximum dimension for Von Neumann entropy computation 40 #define BLEN 16 // Bin length for spatial binning (grid size per dimension) 41 #define PI 3.1415926535897932384626433832795 // Pi constant for trigonometric functions 42 #define TOL 1e-12 // Tolerance for numerical comparisons and eigenvalue thresholds 43 44 // Physical constants from CODATA (used for dimensional consistency, though scaled in sim) 45 double G_CODATA = 6.674300000000000e-11; // Gravitational constant [m^3 kg ^-1 s^-2] 46 double c_CODATA = 2.997924580000000e8; // Speed of light [m s^-1] 47 double hbar_CODATA = 1.054571800000000e-34; // Reduced Planck's constant [J s] 48 double k_B_CODATA = 1.380649000000000e-23; // Boltzmann constant [J K^-1] 49 50 // Structure for physical quantities with units (for validation) 51 struct PhysicalQuantity { 52 double value; // Numerical value 53 char *unit; // Unit string (e.g., "time", "entropy") 54 }; 55 56 // Structure for dimensional tracking (Length L, Time T, Information I exponents) 57 struct DimT { 58 double value; // Numerical value 59 int e_length; // Exponent of length [L^{e_length}] 22 60 int e_time; // Exponent of time [T^{e_time}] 61 int e_info; // Exponent of information [I^{e_info}] 62 char *unit; // Corresponding unit string 63 }; 64 65 /* 66 * Validation Functions 67 * These ensure code correctness: unit consistency, finite values, dimensional analysis. 68 */ 69 70 /** 71 * Validates that a unit string is one of the allowed types. 72 * @param unit The unit to check 73 * @param label Context label for error reporting 74 */ 75 void validate_unit(char *unit, char *label) { 76 // List of allowed units for this simulation 77 char *allowed[] = {"unitless","entropy","vonNunit","probability","time ","length","area","count","dimensionless"}; 78 int n=sizeof(allowed)/sizeof(char*); 79 int found = 0; 80 for(int i=0; i<n; i++) { 81 if(strcmp(allowed[i], unit)==0) { 82 found=1; 83 break;// Early exit for efficiency 84 } 85 } 86 if(!found) { 87 printf("[validate_unit] Invalid unit in %s: %s\n", label, unit); 88 exit(1); // Terminate on invalid unit 89 } 90 } 91 92 /** 93 * Asserts that a value is finite (not NaN or Inf). 94 * @param value The value to check 95 * @param label Context for error message 96 */ 97 void assert_finite(double value, char *label) { 98 if(!isfinite(value)) { 99 printf("[assert_finite] Non-finite value in %s: %f\n", label, value); 100 exit(1); // Terminate on non-finite value 101 } 102 } 103 104 /** 105 * Asserts that a PhysicalQuantity's unit matches the expected one. 106 * @param pq Pointer to PhysicalQuantity 107 * @param expected Expected unit string 23 108 */ 109 void assert_unit(struct PhysicalQuantity *pq, char *expected) { 110 if(strcmp(pq->unit, expected)!=0) { 111 printf("[assert_unit] Unit mismatch Expected: %s Got: %s\n", expected, pq->unit); 112 exit(1); // Terminate on unit mismatch 113 } 114 } 115 116 /** 117 * Asserts dimensional exponents match expected (L, T, I). 118 * @param dt Pointer to DimT 119 * @param l Expected length exponent 120 * @param t Expected time exponent 121 * @param i Expected info exponent 122 * @param label Context for error 123 */ 124 void assert_dimensions(struct DimT *dt, int l, int t, int i, char *label) { 125 if(dt->e_length != l || dt->e_time != t || dt->e_info != i) { 126 printf("ERROR: Dimensional mismatch in %s Expected: [L^%d T^%d I^%d] Got: [L^%d T^%d I^%d]\n", 127 label, l, t, i, dt->e_length, dt->e_time, dt->e_info); 128 exit(1); // Terminate on dimensional mismatch 129 } 130 } 131 132 /** 133 * Dual verification: Combines unit, finite, dimensional, and value consistency checks. 134 * Ensures PhysicalQuantity and DimT agree within tolerance. 135 * @param pq PhysicalQuantity to verify 136 * @param dt DimT to verify 137 * @param label Context label 138 * @param expected_unit Expected unit 139 * @param l,t,i Expected exponents 140 * @param tolerance Relative tolerance for value comparison 141 */ 142 void dual_verify(struct PhysicalQuantity *pq, struct DimT *dt, char *label, char *expected_unit, int l, int t, int i, double tolerance) { 143 // Validate units and finite values 144 validate_unit(pq->unit, label); 145 assert_unit(pq, expected_unit); 146 assert_finite(pq->value, label); 147 validate_unit(dt->unit, label); 148 149 // Cross-validate units via temporary PQ 150 struct PhysicalQuantity tmp = {0, dt->unit}; 151 assert_unit(&tmp, expected_unit); 152 153 // Check dimensions 24 154 assert_dimensions(dt, l, t, i, label); 155 assert_finite(dt->value, label); 156 157 // Check numerical values agree within tolerance 158 double rel_diff = fabs(pq->value - dt->value) / (fabs(pq->value) + 1e-100) ;// Avoid div by zero 159 if(rel_diff > tolerance) { 160 printf("ERROR: Value mismatch in %s Rel diff: %e Tolerance: %e\n", label, rel_diff, tolerance); 161 exit(1); // Terminate on value mismatch 162 } 163 } 164 165 /* 166 * Barnes-Hut Tree Structures and Functions 167 * Generalized octree for N-body gravity in D dimensions (hypercube subdivision). 168 * Each node represents a spatial region with center, size, mass, center-ofmass (COM). 169 */ 170 171 /** 172 * Node structure for the Barnes-Hut tree. 173 */ 174 struct Node { 175 double center[MAX_D]; // Center coordinates of the region 176 double size; // Side length of the hypercube region 177 int D; // Dimensionality 178 double mass; // Total mass in subtree 179 double com[MAX_D]; // Center-of-mass coordinates 180 struct Node **children; // Array of child pointers (2^D children max) 181 int is_leaf; // 1 if leaf node, 0 if internal 182 int *particles; // Array of particle indices (leaf only) 183 int num_particles; // Number of particles in leaf 184 }; 185 186 /** 187 * Creates a new leaf node. 188 * @param center Center coordinates 189 * @param size Region size 190 * @param D Dimensionality 191 * @return Pointer to new Node 192 */ 193 struct Node *new_node(double *center, double size, int D) { 194 struct Node *node = malloc(sizeof(struct Node)); 195 // Copy center 196 memcpy(node->center, center, D*sizeof(double)); 197 node->size = size; 198 node->D = D; 199 node->mass = 0.0; // Initially zero mass 25 478 } 479 free(node); 480 } 481 482 /* 483 * Leapfrog Integration 484 * Second-order symplectic integrator for N-body simulation. 485 * Drift-kick-drift scheme with tree rebuild each half-step. 486 */ 487 488 /** 489 * Performs one full Leapfrog step: kick (half), drift (full), kick (half). 490 * Rebuilds tree before each force computation. 491 * @param positions Positions (N*D) 492 * @param velocities Velocities (N*D) 493 * @param masses Masses (N) 494 * @param N Particles 495 * @param D Dimensions 496 * @param dt Timestep 497 * @param n_timesteps Number of steps 498 * @param theta Barnes-Hut theta 499 * @param G Gravitational constant 500 * @param softening Softening length 501 */ 502 void leapfrog_integration(double *positions, double *velocities, double * masses, int N, int D, double dt, int n_timesteps, double theta, double G, double softening) { 503 for(int step=0; step<n_timesteps; step++) { 504 // First half-kick: v += 0.5 a dt 505 struct Node *tree = build_tree(positions, masses, N, D); 506 double *accels = calloc(N * D, sizeof(double)); // All accelerations 507 for(int i=0; i<N; i++) { 508 double accel[MAX_D] = {0}; // Per-particle accel 509 walk_tree(tree, positions + i * D, accel, i, theta, G, softening, D, positions, masses); 510 for(int d=0; d<D; d++) { 511 accels[i * D + d] = accel[d]; 512 } 513 } 514 // Apply half-kick 515 for(int i=0; i<N; i++) { 516 for(int d=0; d<D; d++) { 517 velocities[i * D + d] += 0.5 * accels[i * D + d] * dt; 518 } 519 } 520 521 // Full drift: x += v dt 522 for(int i=0; i<N; i++) { 523 for(int d=0; d<D; d++) { 524 positions[i * D + d] += velocities[i * D + d] * dt; 32 525 } 526 } 527 528 // Cleanup 529 free(accels); 530 free_tree(tree); 531 532 // Second half-kick: rebuild tree after drift 533 tree = build_tree(positions, masses, N, D); 534 accels = calloc(N * D, sizeof(double)); 535 for(int i=0; i<N; i++) { 536 double accel[MAX_D] = {0}; 537 walk_tree(tree, positions + i * D, accel, i, theta, G, softening, D, positions, masses); 538 for(int d=0; d<D; d++) { 539 accels[i * D + d] = accel[d]; 540 } 541 } 542 // Apply second half-kick 543 for(int i=0; i<N; i++) { 544 for(int d=0; d<D; d++) { 545 velocities[i * D + d] += 0.5 * accels[i * D + d] * dt; 546 } 547 } 548 549 // Cleanup 550 free(accels); 551 free_tree(tree); 552 } 553 } 554 555 /* 556 * Utility Functions 557 */ 558 559 /** 560 * Generates a standard normal random variable using Box-Muller transform. 561 * @return Normal deviate (mean 0, std 1) 562 */ 563 double rand_normal() { 564 double u1 = (double)rand() / RAND_MAX; 565 double u2 = (double)rand() / RAND_MAX; 566 return sqrt(-2 * log(u1)) * cos(2 * PI * u2); 567 } 568 569 /** 570 * Jacobi eigenvalue algorithm for symmetric real matrix (for Von Neumann entropy). 571 * Accumulates rotations to diagonalize; returns eigenvalues in d. 572 * @param n Matrix size 33 573 * @param a Input/output matrix (symmetric, destroyed) 574 * @param it_max Max iterations 575 * @param d Output eigenvalues 576 */ 577 void jacobi_eigenvalue(int n, double *a, int it_max, double *d) { 578 double *bw = malloc(n * sizeof(double)); // Backup of diagonal 579 double *zw = malloc(n * sizeof(double)); // Accumulator for off-diagonal shifts 580 for(int i=0; i<n; i++) { 581 bw[i] = a[i*n + i]; 582 d[i] = a[i*n + i]; // Initial eigenvalues (diagonal) 583 zw[i] = 0.0; 584 } 585 586 for(int k=0; k<it_max; k++) { 587 // Compute off-diagonal norm 588 double sm = 0.0; 589 for(int i=0; i<n-1; i++) { 590 for(int j=i+1; j<n; j++) { 591 sm += fabs(a[i*n + j]); 592 } 593 } 594 if(sm == 0.0) break;// Converged 595 596 // Threshold for rotation 597 double thresh = (k < 3) ? 0.2 * sm / (n*n) : 0.0; 598 599 for(int i=0; i<n-1; i++) { 600 for(int j=i+1; j<n; j++) { 601 double g = 100.0 * fabs(a[i*n + j]); // Scaled for comparison 602 // Skip small elements after convergence 603 if(k > 3 && g <= TOL * fabs(d[i]) && g <= TOL * fabs(d[j])) { 604 a[i*n + j] = 0.0; 605 continue; 606 } 607 if(fabs(a[i*n + j]) > thresh) { 608 // Compute rotation parameters 609 double h = d[j] - d[i]; 610 double t; 611 if(g <= TOL * fabs(h)) { 612 t = a[i*n + j] / h; // Small angle approx 613 }else { 614 double theta = 0.5 * h / a[i*n + j]; 615 t = 1.0 / (fabs(theta) + sqrt(1.0 + theta*theta)); 616 if(theta < 0.0) t = -t; 617 } 618 double c = 1.0 / sqrt(1.0 + t*t); // Cosine 619 double s=t*c; // Sine 620 double tau = s / (1.0 + c); // Tan(2theta)/2 approx 621 34 622 // Update off-diagonal to zero 623 h = t * a[i*n + j]; 624 zw[i] -= h; 625 zw[j] += h; 626 d[i] -= h; 627 d[j] += h; 628 a[i*n + j] = 0.0; 629 630 // Rotate rows/columns 631 for(int l=0; l<i; l++) { 632 double g_val = a[l*n + i]; 633 double h_val = a[l*n + j]; 634 // Givens rotation 635 a[l*n + i] = g_val - s * (h_val + g_val * tau); 636 a[l*n + j] = h_val + s * (g_val - h_val * tau); 637 } 638 for(int l=i+1; l<j; l++) { 639 double g_val = a[i*n + l]; 640 double h_val = a[l*n + j]; 641 a[i*n + l] = g_val - s * (h_val + g_val * tau); 642 a[l*n + j] = h_val + s * (g_val - h_val * tau); 643 // Symmetric: a[l*n + i] updated in upper triangle 644 a[l*n + i] = a[i*n + l]; 645 } 646 for(int l=j+1; l<n; l++) { 647 double g_val = a[i*n + l]; 648 double h_val = a[j*n + l]; 649 a[i*n + l] = g_val - s * (h_val + g_val * tau); 650 a[j*n + l] = h_val + s * (g_val - h_val * tau); 651 a[l*n + i] = a[i*n + l]; 652 a[l*n + j] = a[j*n + l]; 653 } 654 // Update diagonal (already in d) 655 } 656 } 657 } 658 // Final diagonal update with accumulators 659 for(int i=0; i<n; i++) { 660 bw[i] += zw[i]; 661 d[i] = bw[i]; 662 zw[i] = 0.0; 663 } 664 } 665 free(bw); 666 free(zw); 667 } 668 669 /* 670 * HashMap for Counting Binned Configurations 671 * Simple chaining hashmap for counting occupancy in D-dimensional bins. 35 672 */ 673 674 /** 675 * Entry in hashmap chain. 676 */ 677 struct Entry { 678 uint64_t key; // Hashed bin index (stride-encoded coordinates) 679 int value; // Count 680 struct Entry *next; // Next in chain 681 }; 682 683 /** 684 * HashMap structure. 685 */ 686 struct HashMap { 687 struct Entry **buckets; // Array of chain heads 688 size_t size; // Number of buckets 689 }; 690 691 /** 692 * Creates a new hashmap. 693 * @param size Number of buckets (approx. max keys) 694 * @return New HashMap 695 */ 696 struct HashMap *create_hashmap(size_t size) { 697 struct HashMap *map = malloc(sizeof(struct HashMap)); 698 map->buckets = calloc(size, sizeof(struct Entry*)); 699 map->size = size; 700 return map; 701 } 702 703 /** 704 * Inserts or increments count for a key. 705 * @param map Hashmap 706 * @param key 64-bit key 707 * @param value Increment (usually 1) 708 */ 709 void hashmap_put(struct HashMap *map, uint64_t key, int value) { 710 uint64_t h = key % map->size; // Simple modulo hash 711 struct Entry *e = map->buckets[h]; 712 while(e) { 713 if(e->key == key) { 714 e->value += value; // Increment existing 715 return; 716 } 717 e = e->next; 718 } 719 // New entry 720 e = malloc(sizeof(struct Entry)); 721 e->key = key; 36 722 e->value = value; 723 e->next = map->buckets[h]; 724 map->buckets[h] = e; 725 } 726 727 /** 728 * Computes Shannon entropy from hashmap counts: -sum p log p 729 * @param map Hashmap with counts 730 * @param N_particles Total particles (normalization) 731 * @return Shannon entropy S 732 */ 733 double calculate_shannon(struct HashMap *map, int N_particles) { 734 double S = 0.0; 735 for(size_t i=0; i<map->size; i++) { 736 struct Entry *e = map->buckets[i]; 737 while(e) { 738 if(e->value > 0) { 739 double p = (double)e->value / N_particles; 740 S -= p * log(p); // Accumulate -p log p 741 } 742 e = e->next; 743 } 744 } 745 return S; 746 } 747 748 /** 749 * Frees hashmap memory. 750 * @param map Hashmap to free 751 */ 752 void free_hashmap(struct HashMap *map) { 753 for(size_t i=0; i<map->size; i++) { 754 struct Entry *e = map->buckets[i]; 755 while(e) { 756 struct Entry *next = e->next; 757 free(e); 758 e = next; 759 } 760 } 761 free(map->buckets); 762 free(map); 763 } 764 765 /* 766 * Core Simulation Function: Single Trial 767 * Runs N-body sim in "bulk" dimensions (D_bulk = max(0, D-2)), bins positions , 768 * computes Shannon entropy from bins, and Von Neumann from random state. 769 */ 770 37 771 /** 772 * Performs one trial of the hybrid entropy computation. 773 * @param D Total dimensions (screen + bulk) 774 * @param L Box side length 775 * @param W Coarse bin width (unused? but passed) 776 * @param N_particles Number of particles 777 * @param seed Random seed 778 * @param dt Timestep 779 * @param n_timesteps Steps 780 * @param theta BH theta 781 * @param g G (scaled) 782 * @param softening_factor Softening 783 * @param S_shannon_out Output Shannon 784 * @param S_von_scaled_out Output Von Neumann 785 */ 786 void single_trial_hybrid_entropy(int D, double L, int W, int N_particles, unsigned int seed, double dt, int n_timesteps, double theta, double g, double softening_factor, double *S_shannon_out, double *S_von_scaled_out) { 787 srand(seed); // Seed RNG 788 789 int Blen = BLEN; // Bin resolution 790 int D_screen = (D - 2 > 0) ? D - 2 : 0; // Screen dimensions (holographic boundary) 791 int D_bulk = D_screen; // Bulk dims for dynamics (wait, code sets D_bulk = D_screen, but comment suggests D-2 bulk?) 792 // Note: In code, D_bulk = D_screen, but dynamics in D_bulk, binning in full D with fixed screen coords 793 794 double *positions = NULL; 795 double *velocities = NULL; 796 double *masses = NULL; 797 798 // Allocate for bulk dynamics if D_bulk > 0 799 if(D_bulk > 0) { 800 positions = malloc(N_particles * D_bulk * sizeof(double)); 801 velocities = calloc(N_particles * D_bulk, sizeof(double)); // Zero initial vel 802 masses = malloc(N_particles * sizeof(double)); 803 for(int i=0; i<N_particles; i++) { 804 for(int d=0; d<D_bulk; d++) { 805 positions[i*D_bulk + d] = ((double)rand() / RAND_MAX) * L; // Uniform in [0,L) 806 } 807 masses[i] = 1.0; // Unit mass 808 } 809 double softening = softening_factor; // Note: softening_factor is scalar, but softening = factor (units?) 810 // Run dynamics 38 811 leapfrog_integration(positions, velocities, masses, N_particles, D_bulk, dt, n_timesteps, theta, g, softening); 812 } 813 814 // Hashmap for bin counts 815 struct HashMap *counts = create_hashmap(N_particles * 2); // Conservative size 816 817 // Full coordinates: D dims, integer bins [0, Blen-1] 818 int *full_coords = malloc(N_particles * D * sizeof(int)); 819 memset(full_coords, 0, N_particles * D * sizeof(int)); 820 821 // Bin bulk positions 822 if(D_bulk > 0) { 823 for(int i=0; i<N_particles; i++) { 824 for(int d=0; d<D_bulk; d++) { 825 int coord = (int)(positions[i*D_bulk + d] / L * Blen); 826 if(coord < 0) coord = 0; // Clamp 827 if(coord >= Blen) coord = Blen - 1; // Clamp 828 full_coords[i*D + d] = coord; 829 } 830 } 831 } 832 833 // Screen dimensions: fixed to center bin (no dynamics) 834 for(int i=D_bulk; i<D; i++) { 835 for(int j=0; j<N_particles; j++) { 836 full_coords[j*D + i] = Blen / 2; 837 } 838 } 839 840 // Hash bin indices: stride encoding to uint64_t key 841 for(int i=0; i<N_particles; i++) { 842 uint64_t index = 0; 843 uint64_t stride = 1; 844 for(int d=0; d<D; d++) { 845 index += (uint64_t)full_coords[i*D + d] * stride; 846 stride *= Blen; // Blen^D total sites 847 } 848 hashmap_put(counts, index, 1); // Count +1 849 } 850 851 // Compute Shannon entropy 852 double S_shannon = calculate_shannon(counts, N_particles); 853 assert(isfinite(S_shannon)); 854 assert(S_shannon >= 0); 855 856 // Validate with dual_verify 857 struct PhysicalQuantity S_shannon_pq = {.value = S_shannon, .unit = " entropy"}; 39 858 struct DimT S_shannon_dt = {.value = S_shannon, .e_length = 0, .e_time = 0, .e_info = 1, .unit = "entropy"}; 859 dual_verify(&S_shannon_pq, &S_shannon_dt, "Shannon Entropy","entropy", 0, 0, 1, TOL); 860 861 // Cleanup simulation arrays 862 free(full_coords); 863 if(positions) free(positions); 864 if(velocities) free(velocities); 865 if(masses) free(masses); 866 free_hashmap(counts); 867 868 // Von Neumann entropy: from random pure state on subspace 869 uint64_t num_sites = 1; 870 for(int i=0; i<D_bulk; i++) { 871 num_sites *= (uint64_t)Blen; // Blen^{D_bulk} sites 872 } 873 int dim = (num_sites < MAX_DIM) ? num_sites : MAX_DIM; // Cap dimension 874 double S_von = 0.0; 875 876 if(dim >= 2 && dim % 2 == 0) { // Even dim for bipartite 877 double *psi = malloc(dim * sizeof(double)); // Random state vector 878 double norm = 0; 879 for(int i=0; i<dim; i++) { 880 psi[i] = rand_normal(); // Gaussian random 881 norm += psi[i] * psi[i]; 882 } 883 norm = sqrt(norm); 884 assert(norm > 0); 885 for(int i=0; i<dim; i++) { 886 psi[i] /= norm; // Normalize to pure state 887 } 888 889 int half_dim = dim / 2; 890 int stride = dim / half_dim; // Assuming equal bipartition 891 double *a = calloc(half_dim * half_dim, sizeof(double)); // Reduced density matrix <i| rho |j> 892 for(int i=0; i<half_dim; i++) { 893 for(int j=0; j<half_dim; j++) { 894 for(int k=0; k<stride; k++) { 895 a[i*half_dim + j] += psi[i*stride + k] * psi[j*stride + k ]; // Partial trace 896 } 897 } 898 } 899 900 double *evals = malloc(half_dim * sizeof(double)); 901 jacobi_eigenvalue(half_dim, a, 50, evals); // Diagonalize 902 903 double eigenvalues_sum = 0; 40 904 for(int i=0; i<half_dim; i++) { 905 if(evals[i] > 1e-10) { 906 eigenvalues_sum += evals[i]; 907 } 908 } 909 if(eigenvalues_sum > 0) { 910 for(int i=0; i<half_dim; i++) { 911 if(evals[i] > 1e-10) { 912 double lambda_val = evals[i] / eigenvalues_sum; 913 S_von -= lambda_val * log(lambda_val); // -Tr rho log rho 914 } 915 } 916 } 917 918 assert(isfinite(S_von)); 919 assert(S_von >= 0); 920 double S_max = log(half_dim); // Max for equal mix 921 assert(S_von <= S_max + 1e-6); 922 923 free(evals); 924 free(a); 925 free(psi); 926 } 927 928 *S_von_scaled_out = S_von; 929 // Validate Von Neumann 930 struct PhysicalQuantity S_von_pq = {.value = *S_von_scaled_out, .unit = " vonNunit"}; 931 struct DimT S_von_dt = {.value = *S_von_scaled_out, .e_length = 0, .e_time = 0, .e_info = 1, .unit = "vonNunit"}; 932 dual_verify(&S_von_pq, &S_von_dt, "Von Neumann Entropy","vonNunit", 0, 0, 1, TOL); 933 934 *S_shannon_out = S_shannon; 935 } 936 937 /* 938 * Main Function: MPI-Parallel Monte Carlo over D and L 939 * Each process runs local trials, reduces sums, root prints means and scaling . 940 */ 941 942 int main(int argc, char **argv) { 943 // MPI setup 944 int rank, size; 945 MPI_Init(&argc, &argv); 946 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 947 MPI_Comm_size(MPI_COMM_WORLD, &size); 948 949 // Simulation parameters 41