import numpy as np
from scipy.special import spherical_jn, spherical_yn, sph_harm
from scipy.optimize import root
class QuantumBlackHole3D:
"""3D quantum black hole in spherical box model"""
def __init__(self, M, G=6.67430e-11, c=2.99792458e8, hbar=1.054571817e-34):
self.M = M # Black hole mass (kg)
self.G = G
self.c = c
self.hbar = hbar
# Schwarzschild radius
self.R = 2 * G * M / c**2
# For particle in box, we use effective mass (could be M or reduced mass)
self.m = M # Simplified: black hole as single quantum entity
# Bessel function zeros (αₙℓ) for n=1
self.alpha_zeros = {
(1,0): np.pi, # j₀(π)=0
(1,1): 4.4934094579, # j₁(x)=0 first zero
(1,2): 5.7634591969, # j₂(x)=0 first zero
(1,3): 6.9879320005, # j₃(x)=0 first zero
(2,0): 2*np.pi, # j₀(2π)=0
(2,1): 7.7252518369, # j₁(x)=0 second zero
}
def energy_level(self, n, l):
"""Energy of state (n,l) in spherical infinite well"""
if (n,l) in self.alpha_zeros:
alpha = self.alpha_zeros[(n,l)]
else:
# Find zero of spherical Bessel function j_l(x)
alpha = self.find_bessel_zero(n, l)
return (self.hbar**2 * alpha**2) / (2 * self.m * self.R**2)
def find_bessel_zero(self, n, l):
"""Find n-th zero of spherical Bessel function j_l(x)"""
# Define function whose root we want: j_l(x) = 0
def func(x):
return spherical_jn(l, x)
# Initial guess based on asymptotic formula
# Zeros of j_l(x) ~ π(n + l/2) for large n
guess = np.pi * (n + l/2)
# Use root finding (simplified - in practice need careful handling)
# This is a simplified version
result = root(func, guess)
return result.x[0]
def radial_wavefunction(self, n, l, r_points=1000):
"""Calculate radial wavefunction R_nl(r)"""
alpha = self.alpha_zeros.get((n,l), self.find_bessel_zero(n, l))
k = alpha / self.R
r = np.linspace(0, self.R, r_points)
R = spherical_jn(l, k * r)
# Normalize: ∫|R(r)|² r² dr = 1
norm = np.trapz(R**2 * r**2, r)
R_normalized = R / np.sqrt(norm)
return r, R_normalized
def probability_density(self, n, l, m, theta=0, phi=0, r_points=1000):
"""Calculate full probability density |ψ_nlm(r,θ,φ)|²"""
r, R_nl = self.radial_wavefunction(n, l, r_points)
# Spherical harmonic
Y_lm = sph_harm(m, l, phi, theta)
# Full wavefunction ψ_nlm(r,θ,φ) = R_nl(r) * Y_lm(θ,φ)
# Probability density integrated over angles gives radial probability
psi = R_nl * Y_lm
return r, np.abs(psi)**2
def tunneling_probability(self, n, l, V0, a):
"""
Estimate tunneling probability through finite barrier
V0: Barrier height (J)
a: Barrier width (m)
"""
E = self.energy_level(n, l)
if E >= V0:
return 1.0 # Above barrier, not tunneling
# For ℓ=0, simple WKB approximation
kappa = np.sqrt(2 * self.m * (V0 - E)) / self.hbar
# For ℓ>0, additional centrifugal barrier
centrifugal = self.hbar**2 * l*(l+1) / (2 * self.m * self.R**2)
kappa_eff = np.sqrt(2 * self.m * (V0 + centrifugal - E)) / self.hbar
# Transmission probability (simplified)
T = np.exp(-2 * kappa_eff * a)
return T
def hawking_temperature(self):
"""Hawking temperature for this black hole"""
return self.hbar * self.c**3 / (8 * np.pi * self.G * self.M * 1.380649e-23)
def evaporation_time(self):
"""Time for complete evaporation via Hawking radiation"""
return (5120 * np.pi * self.G**2 * self.M**3) / (self.hbar * self.c**4)
# Example calculation for minimum black hole
M_min = 6.5e-9 # kg (from previous calculation)
bh = QuantumBlackHole3D(M_min)
print("3D Quantum Black Hole Model")
print(f"Mass: {M_min:.2e} kg")
print(f"Schwarzschild radius: {bh.R:.2e} m")
print(f"\nEnergy levels (in Joules):")
for (n,l) in [(1,0), (1,1), (1,2), (2,0)]:
E = bh.energy_level(n, l)
print(f" State (n={n}, l={l}): E = {E:.2e} J = {E/1.602e-19:.2e} eV")
print(f"\nGround state energy: {bh.energy_level(1,0):.2e} J")
print(f"Mass energy (Mc²): {M_min * (2.9979e8)**2:.2e} J")
# Calculate tunneling probability
V0 = M_min * bh.c**2 # Rest mass energy as barrier height
a = 1.616e-35 # Planck length as barrier width
T_tunnel = bh.tunneling_probability(1, 0, V0, a)
print(f"\nTunneling probability (ℓ=0, ground state): {T_tunnel:.2e}")
# Compare with Hawking radiation rate
T_H = bh.hawking_temperature()
print(f"Hawking temperature: {T_H:.2e} K")
print(f"Evaporation time: {bh.evaporation_time():.2e} s")
# Generate radial probability distribution
import matplotlib.pyplot as plt
r, prob = bh.probability_density(1, 0, 0, theta=0, phi=0)
plt.figure(figsize=(10, 6))
plt.plot(r/bh.R, prob**2 * r**2, 'b-', linewidth=2, label='Radial probability density')
plt.xlabel('r / R (normalized distance)')
plt.ylabel('Probability density × r²')
plt.title('3D Quantum Black Hole: Ground State (n=1, ℓ=0)')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()
