PUH- BrianMartell

PUH_inflation_sim_v11.py

import numpy as np

# Constants
G = 6.67430e-11  # Gravitational constant (m^3 kg^-1 s^-2)
c = 2.99792458e8  # Speed of light (m/s)
M_sun = 1.989e30  # Solar mass (kg)
l_P = 1.616e-35  # Planck length (m)
hbar = 1.0545718e-34  # Reduced Planck constant (J s)
E_gamma = 1e6 * 1.60218e-19  # Gamma-ray energy ~1 MeV (J)
E_nu = 1e6 * 1.60218e-19  # Neutrino energy ~1 MeV (J)
N_gamma_cmb_0 = 1e89  # Initial CMB photon number
N_gamma_stellar_0 = 0  # Initial stellar photon number
g = 7e-4  # Coupling constant
M_imbh = 225 * M_sun  # GW231123 final mass (kg)
M_planck = 1e12 * M_sun  # Super-Planck star mass (kg)
rho_m = 1e-27  # Matter density (kg/m^3)
rho_m0 = 1e-27  # Reference density (kg/m^3)
R_P = l_P * (rho_m / rho_m0)**0.1  # Planck lattice radius (m)
F_photon = 1e-10  # Flux factor (arbitrary)
a_spin = 0.9  # Spin parameter
t_star =we 1e17  # Stellar lifetime (s, ~3 Gyr)
f_leftover = 0.01  # Fraction of photons as leftovers

# Scale factor (matter-dominated)
def scale_factor(t):
    return (t / 4.35e17)**(2/3) if t > 0 else 1

# Photon number evolution
def photon_numbers(t):
    a = scale_factor(t)
    N_cmb = N_gamma_cmb_0 / (a**3)  # CMB dilution
    N_stellar = 0 if t < 3.15e16 else min(1e83, N_gamma_stellar_0 + 7e81 * (t / 1.89e17)**2)
    N_total = N_cmb + N_stellar
    N_leftover = f_leftover * N_total  # Leftover photons for spacetime creation
    return N_cmb, N_stellar, N_total, N_leftover

# Redshift function
def redshift_energy(M, R, E_emit):
    r_s = 2 * G * M / c**2
    return E_emit / np.sqrt(max(1e-10, 1 - 2 * G * M / (R * c**2)))

# Emission rate
def emission_rate(N, g, phi, M, a):
    J = a * M * c
    xi = g**2 / (l_P**2 * (2.176e-8)**2)
    return (F_photon / l_P**3) * (xi * J**2 * phi**2 / hbar) * (N / 2)

# Photon pressure
def photon_pressure(N, E, R):
    return (N * E) / (4/3 * np.pi * R**3)

# Hubble rate
def hubble_rate(N, E, R):
    rho_gamma = (N * E) / (4/3 * np.pi * R**3)
    return np.sqrt(8 * np.pi * G * rho_gamma / 3)

# Particle mass
def particle_mass(E, config_factor):
    return (E / c**2) * config_factor

# Calculate for IMBH and super-Planck star at t = 13.8 Gyr
phi = 1e-10  # Scalar field
config_factor = 0.1  # Geometric efficiency
t = 4.35e17  # 13.8 Gyr (s)
N_cmb, N_stellar, N_total, N_leftover = photon_numbers(t)
results = []
for M in [M_imbh, M_planck]:
    r_s = 2 * G * M / c**2
    R_values = [1.01 * r_s, R_P]
    E_obs = [redshift_energy(M, R, E_gamma) / 1.60218e-19 for R in R_values]
    gamma_rate = emission_rate(N_leftover, g, phi, M, a_spin)
    nu_rate = emission_rate(N_leftover, g, phi, M, a_spin)
    P_photon = photon_pressure(N_leftover, E_gamma, R_P)
    H = hubble_rate(N_leftover, E_gamma, R_P)
    m_particle = particle_mass(E_gamma, config_factor)
    results.append((M / M_sun, E_obs, gamma_rate, nu_rate, P_photon, H, m_particle))

# Output
print(f"Photon Numbers at t = 13.8 Gyr:")
print(f"  CMB: {N_cmb:.2e}, Stellar: {N_stellar:.2e}, Total: {N_total:.2e}, Leftover: {N_leftover:.2e}")
for i, (M_solar, E_obs, g_rate, n_rate, P_photon, H, m_particle) in enumerate(results):
    name = "GW231123 IMBH" if i == 0 else "Super-Planck Star"
    print(f"{name} ({M_solar:.0f} M☉):")
    print(f"  R = 1.01 r_s, Gamma Energy: {E_obs[0]:.2e} eV")
    print(f"  R = R_P, Gamma Energy: {E_obs[1]:.2e} eV")
    print(f"  Gamma Rate: {g_rate:.2e} photons/s")
    print(f"  Neutrino Rate: {n_rate:.2e} neutrinos/s")
    print(f"  Photon Pressure: {P_photon:.2e} Pa")
    print(f"  Hubble Rate: {H:.2e} s^-1")
    print(f"  Particle Mass: {m_particle:.2e} kg")