Stephen W. Hawking · 1942 – 2018

Black Holes,
Quantum Gravity
& Euler's Formula

How a physicist confined to a wheelchair rewrote the laws of the universe — and how the equation e = cos θ + i sin θ runs through every major discovery he made.

▼   enter the event horizon   ▼

Section 01

Stephen Hawking — Life & Legacy

Born January 8, 1942, in Oxford, England — sharing a birthday with Galileo Galilei — Stephen William Hawking became the most celebrated theoretical physicist of the modern era. Diagnosed with motor neurone disease (ALS) at age 21, given two years to live, he went on to work for more than five more decades, revolutionising our understanding of black holes, cosmology, and the quantum nature of spacetime.

"My goal is simple. It is a complete understanding of the universe, why it is as it is and why it exists at all." — Stephen Hawking

Major Theoretical Contributions

1970

Area Theorem

The total surface area of black hole event horizons never decreases — the first hint of a connection between gravity and thermodynamics.

1971–73

Black Hole Thermodynamics

With Bardeen, Carter, and others, established the four laws of black hole mechanics, mirroring the four laws of thermodynamics exactly.

1974–75

Hawking Radiation

Showed that black holes are not truly black — quantum mechanics forces them to emit thermal radiation with temperature T = ℏc³/8πGMkB.

1971

Penrose–Hawking Singularity Theorems

With Roger Penrose, proved that singularities (points of infinite density) are inevitable inside black holes and at the Big Bang, under GR.

1983

Hartle–Hawking No-Boundary Proposal

With James Hartle, proposed that the universe has no initial boundary in imaginary time — a radical use of Euler's formula and Wick rotation.

1970s–2016

Information Paradox

Identified then spent decades revisiting the question: does information fall into a black hole forever, or does it leak out via Hawking radiation?

Where Euler's Formula Appears

At every major junction of Hawking's work, complex exponentials and the formula e = cos θ + i sin θ play a central role:

TheoryRole of e
Hawking RadiationQuantum field mode expansions e±iωt; thermal spectrum derivation via analytic continuation
Imaginary TimeWick rotation t→−iτ transforms oscillatory eiS/ℏ to damped e−SE/ℏ
No-Boundary ProposalPath integral over Euclidean geometries; wave function Ψ = ∫ D[g] e−SE/ℏ
BH ThermodynamicsPartition function Z = Tr e−βH; Green's functions periodic in imaginary time
Penrose DiagramsConformal transformations use complex coordinates; Kruskal extension via complex rotation
The thread: Euler's formula is the bridge between quantum mechanics (wavefunctions oscillate as e) and statistical mechanics (partition functions decay as e−E/kT). Hawking's deepest insights all live at exactly this bridge.

Section 02

The Penrose–Hawking Singularity Theorems

In 1965, Roger Penrose proved that black holes must contain singularities. Hawking extended this to the Big Bang itself (1966–70). Together, these theorems shattered the hope that singularities were mere artifacts of symmetry — they are mathematically inevitable within Einstein's General Relativity.

What Is a Singularity?

A singularity is a point (or region) where spacetime curvature becomes infinite and the equations of General Relativity break down. At a singularity, density, tidal forces, and spacetime curvature all blow up — the theory predicts its own failure.

The Key Ingredients

1. Geodesic Incompleteness

A spacetime is singular if some freely-falling particle (geodesic) reaches the end of its proper time in finite time — it "hits" the singularity and ceases to exist.

2. Trapped Surfaces

A trapped surface is a 2D surface where even outward-pointing light rays are converging. Once matter collapses past a certain point, trapped surfaces form, and the singularity theorems guarantee a singularity must follow.

3. Energy Conditions

The theorems require that matter obey the strong energy condition: gravity must always attract, never repel. This holds for classical matter (though quantum effects can violate it).

60%
Adjust collapse to see trapped surface formation

The Theorems (Simplified)

Penrose (1965) — Black Hole Singularity: If spacetime contains a trapped surface and energy conditions hold, then at least one geodesic is incomplete — a singularity exists.

Hawking (1966) — Big Bang Singularity: Time-reversing Penrose's argument: if the universe is expanding (as observed), and energy conditions hold, then geodesics are incomplete in the past — the Big Bang is a singularity.

The Implication: GR predicts its own breakdown. A theory of quantum gravity — combining GR with quantum mechanics — is needed to describe the Planck-scale physics near singularities. This quest drove Hawking's entire career.

The Raychaudhuri Equation — Mathematics of Focusing

The mathematical heart of the singularity theorems is the Raychaudhuri equation, which describes how a family of geodesics converges or diverges:

dθ/dτ = −(1/3)θ² − σμνσμν + ωμνωμν − Rμνuμuν

Here θ is the expansion of geodesics, σ is shear, ω is rotation, and Rμν is the Ricci curvature. The energy condition forces Rμνuμuν ≥ 0, making the right side negative — geodesics must converge and eventually cross, signalling a singularity.

Euler's formula enters here via complex null coordinates. The Penrose diagram — the conformal map of all of spacetime onto a finite page — is constructed using complex coordinate transformations of the form u + iv = f(t + x) that involve exactly the rotational structure of e.

Section 03

Black Hole Thermodynamics

In the early 1970s, Hawking, Bardeen, Carter, and Bekenstein uncovered a perfect mathematical analogy between the laws of thermodynamics and the mechanics of black holes. This analogy, initially thought to be mere coincidence, turned out to be deep physical reality once Hawking radiation was discovered.

The Four Laws

LawThermodynamicsBlack Hole Mechanics
Zeroth Temperature T is uniform in thermal equilibrium Surface gravity κ is constant on the event horizon of a stationary black hole
First dE = T dS − P dV (energy conservation) dM = (κ/8π)dA + Ω dJ + Φ dQ (mass-energy conservation)
Second Entropy S never decreases: dS ≥ 0 Total horizon area A never decreases: dA ≥ 0 (Hawking's Area Theorem)
Third Cannot reach T = 0 by finite processes Cannot reduce κ to zero by finite physical processes

The Bekenstein–Hawking Entropy

Jacob Bekenstein (1972) proposed that black holes carry entropy proportional to their horizon area. Hawking initially resisted — if a black hole has entropy, it must have temperature, and something hot must radiate. But black holes absorb everything...

Then Hawking calculated it properly using quantum field theory in 1974 and found that black holes do radiate. The entropy is:

SBH = kB A / 4ℓP²

where A is the horizon area and ℓP = √(ℏG/c³) ≈ 1.6 × 10−35 m is the Planck length. This is one of the most important formulas in theoretical physics — it relates a macroscopic quantity (area) to a microscopic one (entropy), but requires both GR (area) and quantum mechanics (ℏ).

Why is this profound? For any other physical system, entropy is extensive (scales with volume). Black holes are the first objects whose entropy scales with area — the holographic principle. This suggests our 3D universe might be encoded on a 2D surface.
10 M☉

The Area Theorem

Hawking's 1971 Area Theorem states: the total area of all black hole event horizons can never decrease in any physical process. Formally, if Atotal is the total horizon area at late times and early times:

Afinal ≥ Ainitial

This mirrors the second law of thermodynamics (entropy never decreases) and was the first hint that black holes have a thermodynamic nature.

The Partition Function — Euler's Formula Hidden in Temperature

The temperature of a black hole is connected to a deep quantum statistical mechanics result. The thermal partition function is:

Z = Tr[e−βH] = Tr[e−H/kBT]

Compare this to the quantum time-evolution operator:

U(t) = e−iHt/ℏ

They are related by Wick rotation: setting t = −iτ (i.e., rotating time by 90° in the complex plane using Euler's formula), we get e−iHt/ℏ → e−Hτ/ℏ. The partition function is quantum time evolution in imaginary time. For a black hole, periodicity in imaginary time with period β = ℏ/kBT defines the temperature.

The Hawking Temperature: TH = ℏc³ / (8πGMkB). For a solar-mass black hole, TH ≈ 6×10−8 K — colder than the CMB (2.7 K). Only primordial micro–black holes would be hot enough to observe directly.

Section 04

Hawking Radiation — The Quantum Miracle

In 1974, Hawking published what is widely considered the most important result in theoretical physics of the 20th century: black holes are not black. They emit thermal radiation due to quantum effects near the event horizon, with a perfect Planck spectrum at temperature TH. This result combines all three great theories of physics: general relativity (black holes), quantum mechanics (particle creation), and thermodynamics (temperature).

Why this was shocking: Classically, nothing escapes a black hole — not even light. Hawking showed that quantum mechanics changes this fundamentally. The derivation relies entirely on how quantum fields (described by complex exponential modes e±iωt) behave across the event horizon.

The Physical Picture — Pair Creation

Near the horizon, quantum fluctuations continuously create virtual particle–antiparticle pairs. Near the event horizon, one particle of a pair can fall in while the other escapes. The infalling particle carries negative energy (relative to infinity), reducing the black hole's mass, while the escaping particle carries positive energy — this is the Hawking radiation.

50 t=0
Hawking radiation: virtual pair creation near the horizon

The Mathematical Derivation — Euler's Formula at the Core

The full derivation uses quantum field theory (QFT) in curved spacetime. Here are the key steps, each involving complex exponentials:

Mode Expansion — A quantum field φ in flat spacetime is expanded in frequency modes. Each mode is a complex exponential:

φ = ∫dω [aω · e−iωt+ikx + aω · e+iωt−ikx] / √(2ω)

Here aω and aω are annihilation/creation operators, and eiωt is directly Euler's formula. In flat spacetime, this defines the vacuum state |0⟩.

The Bogoliubov Transformation — Near a black hole, the notion of "positive frequency" changes across the horizon. The modes seen by a distant observer (Schwarzschild coordinates) and by a freely-falling observer (Kruskal coordinates) are related by a Bogoliubov transformation:

ãω = ∫dω′ [αωω′ aω′ + βωω′ aω′]

The β coefficients are nonzero — the vacuums are different. Particles seen by one observer are not particles to another. The β coefficients involve analytic continuation of mode functions through the complex plane.

Analytic Continuation Across the Horizon — The key insight: a mode function e−iωu (where u is a retarded time coordinate) must be analytically continued through the horizon. Near the horizon, u = −(1/κ) log(−v) where κ is the surface gravity. Continuing to the other side of the horizon (v → 0⁺) requires:

e−iω·(−1/κ)log(−v) = e−iω·(−1/κ)[log|v| − iπ] = e−πω/κ · eiωlog|v|/κ

The factor e−πω/κ comes from the imaginary part of log(−v) = log|v| − iπ (rotating around the branch cut), which is precisely Euler's formula applied to the complex logarithm.

The Planck Spectrum — The Bogoliubov coefficients satisfy |βωω′|²/|αωω′|² = e−2πω/κ, giving a thermal particle number:

⟨Nω⟩ = 1 / (eℏω/kBTH − 1)

This is exactly the Planck blackbody distribution at temperature TH = ℏκ/2πkBc. The exponential in the denominator traces directly to the phase factor from step 3 — from Euler's formula applied to the complex log.

20 nK
Planck blackbody spectrum ⟨N_ω⟩ = 1/(e^(ℏω/kT)−1) — the Hawking radiation distribution

The Temperature Formula

TH = ℏc³ / (8πGMkB)

For a Schwarzschild black hole of mass M:

ObjectMassTHEvaporation Time
Solar-mass BH2×10³⁰ kg6×10⁻⁸ K~10⁶⁷ years
Stellar BH (10M☉)2×10³¹ kg6×10⁻⁹ K~10⁷⁰ years
Supermassive BH10⁹ M☉6×10⁻¹⁷ K~10⁹⁶ years
Primordial Micro-BH10⁻⁸ kg (Planck)~10³¹ K~10⁻⁴³ s

Section 05

Imaginary Time & Wick Rotation

One of Hawking's most powerful and elegant tools is imaginary time — replacing real time t with imaginary time τ = it. This is not a physical observation but a mathematical technique, enabled entirely by Euler's formula, that converts quantum mechanics problems into statistical mechanics ones.

The Wick Rotation

In quantum mechanics, the time-evolution operator is:

U(t) = e−iHt/ℏ

This is an oscillatory complex exponential — directly Euler's formula. Now perform the substitution t = −iτ (rotate by 90° in the complex time plane):

U(−iτ) = e−iH(−iτ)/ℏ = e−Hτ/ℏ

The oscillatory factor e becomes a real decaying exponential e−θ. This transformation — the Wick rotation — converts between:

Real time (Minkowski)
Path integral: Z = ∫D[φ] eiS[φ]/ℏ
Highly oscillatory, hard to compute
Spacetime metric: ds² = −c²dt² + dx²
Imaginary time (Euclidean)
Path integral: Z = ∫D[φ] e−SE[φ]/ℏ
Rapidly damped, easier to define
Metric: ds² = +c²dτ² + dx² (positive definite!)
Real time axis: e^(−iHt/ℏ) oscillates

Temperature from Periodicity

Hawking's deepest use of imaginary time is to derive the temperature of a black hole from the geometry of spacetime in imaginary time.

Near a black hole horizon, the Schwarzschild metric in imaginary time (τ = it) looks like:

ds²E = (1 − 2GM/rc²) c²dτ² + dr²/(1 − 2GM/rc²) + r²dΩ²

Near the horizon r = rs = 2GM/c², this simplifies to:

ds²E ≈ ρ²d(κτ)² + dρ²

where ρ = √(r − rs) is the proper distance and κ is the surface gravity. This is the metric of a flat plane in polar coordinates — but only if κτ is periodic with period !

κτ ~ κτ + 2π  ⟹  τ ~ τ + 2π/κ  ⟹  β = ℏ/kBT = 2πℏ/κ
The magic: If we demand that the imaginary-time geometry be non-singular (no conical singularity at the horizon), time must be periodic. This forces a specific temperature: TH = ℏκ/2πkB. The periodicity is exactly the ei·2π = 1 property of Euler's formula — going around a full circle in the complex time plane brings you back to the start.

The Conical Singularity

If the periodicity is wrong, the imaginary-time metric has a conical singularity — like a cone with a deficit angle — at the horizon. This would represent a physical singularity (infinite curvature). The Hawking temperature is exactly the value that makes the cone into a smooth flat plane (zero deficit angle).

Exact TH
At exact Hawking temperature: flat plane, no conical singularity

Section 06

The Hartle–Hawking No-Boundary Proposal

In 1983, Stephen Hawking and James Hartle proposed one of the most radical ideas in the history of cosmology: the universe has no boundary in imaginary time, and therefore no initial singularity. The Big Bang is replaced by a smooth, rounded beginning — like the South Pole of a sphere. There is no "before" the Big Bang, not because time began, but because the question is ill-formed.

"The boundary condition of the universe is that it has no boundary." — Stephen Hawking, A Brief History of Time (1988)

The Wave Function of the Universe

Quantum cosmology seeks a wave function Ψ for the universe itself. In the path integral formulation, Ψ is defined by summing over all possible spacetime geometries that satisfy certain boundary conditions:

Ψ[hij] = ∫ D[gμν] D[φ] · eiS[g,φ]/ℏ

Here hij is the induced metric on a 3-surface (what we can observe), and the integral is over all 4D spacetimes with that boundary. This is quantum mechanics — specifically Euler's formula in the path integral — applied to all of spacetime at once.

Wick-Rotating the Universe

The oscillatory path integral eiS/ℏ is formally ill-defined. The Hartle-Hawking proposal is to evaluate it in Euclidean (imaginary) time, where it becomes the much better-behaved:

Ψ[hij] = ∫ D[g] e−SE[g]/ℏ

And crucially: sum only over compact, boundaryless Euclidean geometries — 4D manifolds with no edges. This removes the initial singularity by fiat: there simply is no boundary of the universe in (imaginary) time.

No-Boundary: the universe begins as a smooth 4-sphere in imaginary time

The Saddle-Point Approximation

The Euclidean path integral is dominated by the saddle point — the geometry that minimises the Euclidean action SE. For the no-boundary proposal, this is the de Sitter space (a 4-sphere):

Ψ ≈ e−SE[gsaddle]/ℏ

Computing SE for the de Sitter saddle gives:

Ψ ≈ exp(+3/8ΛℓP²)

where Λ is the cosmological constant. The probability of creating a large universe P ∝ |Ψ|² ≈ exp(+3/4ΛℓP²) — it prefers small Λ, i.e., large universes.

What the Proposal Implies

No Initial Singularity: In Euclidean (imaginary) time, the universe begins smoothly with no boundary — no moment before which the theory breaks down.
Quantum Creation from Nothing: The universe spontaneously nucleates from a quantum fluctuation in nothing. There is no external "creator" — the wave function evaluates to a finite, calculable probability.
Inflation: The no-boundary wave function naturally selects an early inflationary phase, explaining the large-scale homogeneity and isotropy of the universe observed today.
The Role of Euler's Formula: Every step of this argument requires t → −iτ (Wick rotation), which is a rotation in the complex plane via e. Without Euler's formula, imaginary time has no meaning.

Section 07

The Black Hole Information Paradox

Hawking's 1974 calculation of black hole radiation contained a devastating implication that haunted him for the rest of his life: information might be destroyed. If true, this would violate one of the most fundamental principles of quantum mechanics — unitarity.

The Paradox

Quantum Mechanics Says:

The quantum state of a system always evolves unitarily — via a unitary operator U = e−iHt/ℏ. This is Euler's formula applied to quantum time evolution. Unitarity means: information is never lost. If you know the final state exactly, you can always reconstruct the initial state.

|ψ(t)⟩ = e−iHt/ℏ|ψ(0)⟩

U is unitary: U†U = 1. No information can be created or destroyed.

Hawking Radiation Says:

Hawking radiation is thermal — it carries no information about what fell into the black hole. An encyclopedia thrown into a black hole and a chair made of the same mass produce identical radiation. When the black hole evaporates, all that information is gone.

The contradiction: Quantum mechanics says information cannot be destroyed. Hawking's calculation says it is. One of them must be wrong.

Timeline of the Debate

1974–75

Hawking's Original Claim

Information is lost when black holes evaporate. Hawking argued this requires a modification of quantum mechanics — pure states could evolve into mixed states.

1981

The Black Hole War Begins

Leonard Susskind and Gerard 't Hooft argued vigorously that information must be preserved, setting off decades of debate. They bet against Hawking.

1997

Maldacena's AdS/CFT Correspondence

A duality between gravity in Anti-de Sitter space and a conformal field theory on its boundary, where unitarity is manifest. Information is preserved in the boundary theory.

2004

Hawking Concedes

At a conference in Dublin, Hawking publicly conceded that information is preserved — it leaks out via subtle correlations in Hawking radiation. He paid his bets to Preskill.

2012

The Firewall Paradox

AMPS (Almheiri, Marolf, Polchinski, Sully) argued that information preservation leads to a "firewall" — infalling observers would be incinerated at the horizon.

2016–2022

The Island Formula

New calculations using quantum extremal surfaces and "islands" show the Page curve (entropy first rising, then decreasing) can be reproduced, suggesting information is indeed encoded in late Hawking radiation.

Where Euler's formula matters: The information paradox is fundamentally about whether the time-evolution operator e−iHt/ℏ — which embeds Euler's formula — remains unitary throughout black hole evaporation. Every proposed resolution (AdS/CFT, firewalls, islands) can be formulated as a question about the analyticity of complex exponentials describing quantum states across the horizon.
0%
The Page curve: entropy of Hawking radiation over evaporation time

Section 08

Euler's Formula: The Mathematical Spine of Hawking's Work

Let us now trace the precise role of e = cos θ + i sin θ through every major result of Hawking's career.

Mode Expansions — Quantum Fields as Euler Waves

Every quantum field is a collection of harmonic oscillators. Each mode is a complex exponential solution to the wave equation:

φk(x,t) = A · ei(k·x − ωt) = A · eikx · e−iωt

This is Euler's formula twice over: once in space (eikx) and once in time (e−iωt). The vacuum state |0⟩ is defined as the state annihilated by ak for all k, where ak and ak satisfy:

[ak, ak′] = δ(k − k′)

The crucial point for Hawking radiation: near a black hole horizon, the natural modes for an infalling observer differ from those for a distant observer. The two are related by Bogoliubov transformations, and the mismatch (the β coefficients) grows exponentially due to the analytic structure of eiωt near the horizon.

3 0.10
Mode expansion across horizon: blue=flat spacetime, gold=near-horizon modified mode

Wick Rotation — Rotating the Complex Time Plane

The Wick rotation t = −iτ is a rotation in the complex time plane by 90°:

e−iH·t/ℏ → e−iH·(−iτ)/ℏ = e−Hτ/ℏ

This is exactly the rotation e with θ = −π/2. The oscillatory quantum evolution becomes the thermal Boltzmann weight. The implications are profound:

  • Quantum ↔ Statistical: The path integral in imaginary time computes thermodynamic partition functions
  • Geometry ↔ Temperature: Periodicity in imaginary time = inverse temperature β = 1/kBT
  • Euclidean geometry: The Minkowski metric −dt² + dx² becomes +dτ² + dx², making calculations well-defined
t = −iτ : e = ei·(−iτ·ω) = eωτ  (real, no oscillation)

Thermal Green's Functions — KMS Condition

A thermal quantum field theory at temperature T has a Green's function that satisfies the KMS condition (Kubo–Martin–Schwinger). For the two-point function:

G(t) = G(t + iβ)    where β = ℏ/kBT

This periodicity in imaginary time — a shift by iβ in the complex t-plane — is the defining property of a thermal state. It follows directly from the cyclic property of the trace in Z = Tr[e−βH] and the Heisenberg equation:

O(t) = eiHt/ℏ O(0) e−iHt/ℏ

For Hawking radiation, the fact that the quantum field state outside a black hole is thermal means it automatically satisfies the KMS condition with β = 2πℏ/κ. This periodicity is equivalent to saying the state is periodic under t → t + 2πi/κ — a complex shift, underpinned by Euler's formula.

20
KMS condition: G(t) = G(t+iβ) — thermal periodicity in imaginary time

Path Integrals — Euler's Formula Over All Spacetimes

The quantum gravity path integral sums over all possible 4D geometries with given boundary conditions. In Lorentzian signature:

Z = ∫D[g] D[φ] eiS[g,φ]/ℏ

Every geometry contributes a phase factor eiS/ℏ — exactly Euler's formula with the action S as the angle. Geometries that extremise S contribute most (stationary phase approximation, i.e., classically allowed paths).

After Wick rotation to Euclidean time, this becomes:

ZE = ∫D[g] D[φ] e−SE[g,φ]/ℏ

The Hartle-Hawking wave function uses this directly, summing over compact Euclidean 4-manifolds with no boundary. The saddle-point (dominant contribution) is the 4-sphere — de Sitter space in imaginary time — and gives the probability for the universe to nucleate.

Unitarity — Euler's Formula and Information

Quantum mechanical time evolution is governed by:

U(t) = e−iHt/ℏ

Unitarity follows from the Hermiticity of H: U†U = e+iH†t/ℏe−iHt/ℏ = 1 when H† = H. This is directly the property ee−iθ = 1 from Euler's formula — unit complex numbers preserve norms.

The information paradox reduces to: is e−iHt/ℏ genuinely unitary when H includes a black hole? If Hawking radiation is exactly thermal (maximally random), it cannot be unitary. The modern consensus (from AdS/CFT) is that black holes are unitary — but the mechanism by which unitarity is preserved in the radiation is still not fully understood.

|U(t)|² = |e−iHt/ℏ|² = 1  (unitarity = Euler's formula!)

Section 09

GNU Octave Code Examples

Ten complete, runnable scripts demonstrating Hawking's physics computationally. Each script is self-contained with detailed comments.

1. Hawking Temperature and Black Hole Properties

Calculate Hawking temperature, entropy, and evaporation time for black holes of different masses.

% ── Hawking Radiation: Temperature, Entropy, Evaporation ─────────────────
% Physical constants (SI units)
hbar = 1.0546e-34;   % ℏ = h/(2π)  [J·s]
c    = 2.998e8;      % speed of light [m/s]
G    = 6.674e-11;    % gravitational constant [m³/kg/s²]
kB   = 1.381e-23;    % Boltzmann constant [J/K]
Msun = 1.989e30;    % solar mass [kg]
lP   = sqrt(hbar*G/c^3);  % Planck length [m]

% ── Hawking Temperature: T_H = ℏc³/(8πGMk_B) ────────────────────────────
% This comes from the Bogoliubov coefficient calculation.
% The e^(−πω/κ) factor from analytically continuing e^(iωt) through
% the horizon (using Euler's formula) gives the thermal factor.
T_hawking = @(M) hbar * c^3 / (8*pi * G * M * kB);

% ── Schwarzschild Radius: r_s = 2GM/c² ───────────────────────────────────
r_s = @(M) 2 * G * M / c^2;

% ── Bekenstein-Hawking Entropy: S = k_B A / (4 l_P²) ────────────────────
% A = 4π r_s² = 16π G² M² / c⁴
S_bh = @(M) kB * 4*pi * r_s(M)^2 / (4 * lP^2);

% ── Evaporation Time: t_evap ≈ 5120 π G² M³ / (ℏ c⁴) ────────────────────
t_evap = @(M) 5120 * pi * G^2 * M^3 / (hbar * c^4);

% ── Print table for different black hole masses ───────────────────────────
masses_Msun = [1e-6, 1e-3, 1, 10, 1e6, 1e9];
names = {'Micro (μ)', 'Micro (m)', 'Solar', 'Stellar', 'Intermediate', 'Supermassive'};

fprintf('%-15s %12s %12s %12s %15s\n', 'Type', 'M (Msun)', 'T_H (K)', 'r_s (km)', 't_evap (s)');
fprintf('%s\n', repmat('-',1,65));
for i = 1:length(masses_Msun)
  M   = masses_Msun(i) * Msun;
  T   = T_hawking(M);
  rs  = r_s(M);
  S   = S_bh(M);
  tev = t_evap(M);
  fprintf('%-15s %12.3g %12.3g %12.3g %15.3g\n', names{i}, masses_Msun(i), T, rs/1000, tev);
end

% ── Plot T_H vs M ─────────────────────────────────────────────────────────
M_range = logspace(-3, 12, 500) * Msun;   % from milisolar to billion solar
T_range = arrayfun(T_hawking, M_range);

figure('Name', 'Hawking Temperature');
loglog(M_range/Msun, T_range, 'b-', 'LineWidth', 2); grid on;
yline(2.725, 'r--', 'CMB (2.725 K)');    % CMB temperature
xlabel('M / M_{sun}'); ylabel('T_H (K)');
title('Hawking Temperature: T_H = \hbar c^3 / (8\pi G M k_B)');
text(1, T_hawking(Msun)*3, sprintf('T(M_\odot) = %.2e K', T_hawking(Msun)), 'FontSize', 10);

2. Planck Spectrum — Euler's Formula in the Blackbody Distribution

The Hawking radiation spectrum is a perfect Planck distribution. The exponential in it traces directly to Euler's formula via the Bogoliubov analytic continuation.

% ── Planck Spectrum for Hawking Radiation ────────────────────────────────
% ⟨N_ω⟩ = 1 / (exp(ℏω/k_B T_H) − 1)  — Bose-Einstein distribution
% The exp() here comes from the phase factor e^(-πω/κ) in Bogoliubov theory,
% which arises from analytically continuing e^(iωt) around the horizon.

hbar = 1.0546e-34; kB = 1.381e-23; c = 2.998e8;
G = 6.674e-11;    Msun = 1.989e30;

T_hawking = @(M) hbar * c^3 / (8*pi * G * M * kB);

% Planck occupation number for massive/massless bosons
bose_einstein = @(omega, T) 1 ./ (exp(hbar * omega ./ (kB * T)) - 1);
fermi_dirac   = @(omega, T) 1 ./ (exp(hbar * omega ./ (kB * T)) + 1);

% Spectral energy density (Planck formula)
planck_u = @(omega, T) (hbar * omega.^3) ./ (pi^2 * c^3) .* bose_einstein(omega, T);

% ── Compare Hawking radiation at different BH masses ──────────────────────
figure('Name', 'Hawking Radiation Spectra');
BH_masses = [1e10, 1e11, 1e12] * 1e3;  % tiny primordial BHs in kg
colors = {'b-', 'r-', 'g-'};
legends = {};

subplot(1,2,1);
for i = 1:length(BH_masses)
  M = BH_masses(i);
  T = T_hawking(M);
  omega_peak = 2.821 * kB * T / hbar;   % Wien peak
  omega = linspace(0.01, 10) * omega_peak;
  N_omega = bose_einstein(omega, T);
  plot(omega/omega_peak, N_omega, colors{i}, 'LineWidth', 2); hold on;
  legends{end+1} = sprintf('M=%.0e kg, T=%.1f K', M, T);
end
xlabel('\omega / \omega_{peak}'); ylabel('');
title('Occupation number: 1/(e^{\hbar\omega/kT}-1)');
legend(legends, 'Location', 'NorthEast'); grid on;

% ── Show how e^(-πω/κ) factor arises ──────────────────────────────────────
% From analytic continuation: log(-v) = log|v| - iπ → picks up e^(-πω/κ)
subplot(1,2,2);
kappa = 1;  % surface gravity (units where kappa=1)
omega_range = linspace(0.01, 5, 200);
ratio_beta = abs(exp(-pi*omega_range/kappa)).^2;  % |beta/alpha|^2
N_thermal   = 1 ./ (exp(2*pi*omega_range/kappa) - 1);  % thermal result
plot(omega_range, ratio_beta, 'b-', 'LineWidth', 2); hold on;
plot(omega_range, N_thermal, 'r--', 'LineWidth', 2);
xlabel('\omega/\kappa'); ylabel('amplitude');
legend({'|β|² = e^{-2πω/κ}', 'Thermal: 1/(e^{2πω/κ}-1)'});
title('Bogoliubov factor from analytic continuation');
grid on;
sgtitle('Euler\'s formula in Hawking radiation: e^{i\omega t} → e^{-\pi\omega/\kappa}');

3. Black Hole Evaporation — Full Time Evolution

Simulate how a black hole evaporates via Hawking radiation, tracking mass, temperature, and luminosity.

% ── Black Hole Evaporation via Hawking Radiation ─────────────────────────
% Stefan-Boltzmann-like power law: dM/dt = -ℏc⁴/(15360πG²M²)
% Solution: M(t) = M_0 * (1 - t/t_evap)^(1/3)
% Temperature: T_H = ℏc³/(8πGMk_B) — diverges as M→0

hbar=1.0546e-34; c=2.998e8; G=6.674e-11; kB=1.381e-23; Msun=1.989e30;

% Evaporation time: t_evap = 5120π G² M₀³ / (ℏ c⁴)
t_evap_fn = @(M0) 5120*pi * G^2 * M0^3 / (hbar * c^4);

% Use a primordial black hole (observable in our universe's lifetime)
M0_kg = 1e12;                    % kg — about the mass of a mountain
t_total = t_evap_fn(M0_kg);       % should be ~age of universe
fprintf('M0 = %.2e kg, t_evap = %.3e s = %.3e yr\n', M0_kg, t_total, t_total/3.156e7);

% Normalized time 0→1
t_norm = linspace(0, 0.999, 2000);
M_t = M0_kg * (1 - t_norm).^(1/3);
T_t = hbar * c^3 ./ (8*pi * G * M_t * kB);
L_t = hbar * c^6 ./ (15360*pi * G^2 * M_t.^2);  % luminosity
r_t = 2 * G * M_t / c^2;                          % Schwarzschild radius

figure('Name', 'Black Hole Evaporation');

subplot(2,2,1);
  plot(t_norm, M_t/M0_kg, 'b-', 'LineWidth', 2); grid on;
  xlabel('t / t_{evap}'); ylabel('M / M_0');
  title('Mass: M(t) = M_0(1-t/t_{evap})^{1/3}');

subplot(2,2,2);
  semilogy(t_norm, T_t, 'r-', 'LineWidth', 2); grid on;
  xlabel('t / t_{evap}'); ylabel('T_H (K)');
  title('Hawking Temperature: T_H \propto 1/M');

subplot(2,2,3);
  semilogy(t_norm, L_t, 'g-', 'LineWidth', 2); grid on;
  xlabel('t / t_{evap}'); ylabel('Luminosity (W)');
  title('Hawking Luminosity: L \propto 1/M²');

subplot(2,2,4);
  plot(t_norm, r_t, 'm-', 'LineWidth', 2); grid on;
  xlabel('t / t_{evap}'); ylabel('r_s (m)');
  title('Schwarzschild Radius: r_s = 2GM/c²');

sgtitle('Black Hole Evaporation — Hawking Radiation');

4. Wick Rotation — Visualizing the Complex Time Plane

Visualize how rotating time into the imaginary axis (Wick rotation) transforms oscillatory quantum evolution into thermal Boltzmann weights.

% ── Wick Rotation: t → -iτ using Euler's Formula ─────────────────────────
% e^(-iHt/ℏ) oscillates in real time t
% e^(-Hτ/ℏ)  decays in imaginary time τ = it
% The rotation in the complex time plane IS Euler's formula: e^(iθ)

H = 3.0;    % energy eigenvalue (ℏ = 1 units)
hbar = 1;

% ── Real-time evolution: oscillatory ─────────────────────────────────────
t_real = linspace(0, 6*pi/3, 400);
U_real = exp(-1i * H * t_real / hbar);   % e^(-iHt/ℏ) from Euler's formula

% ── Imaginary-time evolution: thermal weight ──────────────────────────────
tau_imag = linspace(0, 5, 300);
U_imag = exp(-H * tau_imag / hbar);       % e^(-Hτ/ℏ): decays to zero

% ── Complex time plane: Wick rotation by angle θ ─────────────────────────
theta_wick = linspace(0, pi/2, 100);   % rotate from real to imaginary axis
t0 = 1.0;                              % fixed magnitude
t_complex = t0 * exp(1i * theta_wick);  % Euler: rotate in complex plane
U_complex = exp(-1i * H * t_complex);   % amplitude as we rotate

figure('Name', 'Wick Rotation');
subplot(2,2,1);
  plot(t_real, real(U_real), 'b-', t_real, imag(U_real), 'r-', 'LineWidth', 2);
  legend({'Re[e^{-iHt}]', 'Im[e^{-iHt}]'}); grid on;
  xlabel('real time t'); title('Real time: oscillatory');

subplot(2,2,2);
  plot(tau_imag, U_imag, 'g-', 'LineWidth', 2); grid on;
  xlabel('imaginary time 	au'); title('Imaginary time: Boltzmann weight');
  ylabel('e^{-H	au}');

subplot(2,2,3);
  plot(real(t_complex), imag(t_complex), 'k-', 'LineWidth', 2);
  hold on; plot(t0, 0, 'go', 0, -t0, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'auto');
  grid on; axis equal; axis([-0.2 1.3 -1.3 0.2]);
  xlabel('Re(t)'); ylabel('Im(t)');
  title('Wick rotation path: t → -i	au (via Euler''s formula)');
  legend({'rotation arc', 'real t=1', 'imag t=-i'});

subplot(2,2,4);
  plot(theta_wick*180/pi, abs(U_complex), 'b-', 'LineWidth', 2); hold on;
  plot(theta_wick*180/pi, real(U_complex), 'r-', 'LineWidth', 2);
  grid on; xlabel('Wick angle (degrees)'); 
  legend({'|e^{-iHt}|', 'Re[e^{-iHt}]'});
  title('Amplitude during Wick rotation');
  xline(90, 'k--', '	au axis');

sgtitle('Wick Rotation: e^{-iHt/\hbar} → e^{-H\tau/\hbar}');

5. Quantum Field Mode Expansion Near a Black Hole

Illustrate how field modes e±iωt behave differently inside vs. outside the horizon.

% ── Field Mode Expansion Near Black Hole Horizon ─────────────────────────
% Outside the horizon (r > r_s): modes ~ e^(-iωt) in Schwarzschild coords
% Inside the horizon: modes modified by log singularity → e^(-πω/κ)
% This factor from Euler's formula is what causes Hawking radiation.

omega = 2.5;   % mode frequency
kappa = 1.0;   % surface gravity (units: ℏ = c = k_B = 1)

% Outside mode: u(t) = e^(-iωt) — standard oscillation from Euler's formula
t_out = linspace(-3, 3, 500);
mode_out = exp(-1i * omega * t_out);

% Tortoise coordinate: r_star = r + r_s * log(|r/r_s - 1|)
% Near the horizon, retarded time u = t - r_star → -∞
% Mode e^(-iω u) needs analytic continuation through u = -∞ horizon

% The key calculation: log(-v) = log|v| - iπ (Euler's formula for complex log)
v_vals = linspace(-2, 2, 1000);
v_vals(abs(v_vals) < 0.001) = NaN;  % avoid zero

% Outside (v > 0): mode = e^(-iω/κ * log(v))
mode_outside  = exp(-1i*omega/kappa * log(abs(v_vals)));
% Inside (v < 0): analytically continue log(-v) = log|v| - iπ
% → e^(-iω/κ * [log|v| - iπ]) = e^(-πω/κ) * e^(-iω/κ * log|v|)
beta_factor = exp(-pi * omega / kappa);  % Bogoliubov β coefficient from Euler!
mode_inside = beta_factor * exp(-1i*omega/kappa * log(abs(v_vals)));

figure('Name', 'Mode Expansion Near Horizon');
subplot(1,2,1);
  v_pos = v_vals(v_vals > 0 & ~isnan(v_vals));
  v_neg = v_vals(v_vals < 0 & ~isnan(v_vals));
  m_out = mode_outside(v_vals > 0 & ~isnan(v_vals));
  m_in  = mode_inside(v_vals < 0 & ~isnan(v_vals));
  plot(v_pos, real(m_out), 'b-', 'LineWidth', 1.5); hold on;
  plot(v_neg, real(m_in),  'r-', 'LineWidth', 1.5);
  xline(0, 'k--', 'Horizon'); grid on;
  xlabel('v (Kruskal coordinate)'); ylabel('Re[\phi]');
  legend({'Outside (v>0)', 'Inside (v<0)'});
  title(sprintf('Mode across horizon; |\beta| = e^{-\pi\omega/\kappa} = %.4f', beta_factor));

subplot(1,2,2);
  % Show the thermal spectrum arising from |β|²/|α|²
  omega_range = linspace(0.1, 5, 200);
  thermal_n = 1 ./ (exp(2*pi*omega_range/kappa) - 1);
  beta_sq = exp(-2*pi*omega_range/kappa);
  semilogy(omega_range, thermal_n, 'b-', omega_range, beta_sq, 'r--', 'LineWidth', 2);
  legend({'⟨N_ω⟩ = 1/(e^{2πω/κ}-1)', '|β|² = e^{-2πω/κ}'}); grid on;
  xlabel('\omega/\kappa'); ylabel('amplitude');
  title('Thermal spectrum from Bogoliubov coefficients');
sgtitle('Euler\'s formula at the horizon: e^{i\omega t} → thermal factor');

6. Bekenstein–Hawking Entropy and the Holographic Principle

% ── Bekenstein-Hawking Entropy: S = k_B A / (4 l_P²) ────────────────────
% This formula unites GR (A = area), QM (ℏ in l_P), and thermo (k_B, S)

hbar=1.0546e-34; c=2.998e8; G=6.674e-11; kB=1.381e-23; Msun=1.989e30;
lP = sqrt(hbar*G/c^3);   % Planck length ≈ 1.616e-35 m
tP = lP/c;                  % Planck time
mP = sqrt(hbar*c/G);        % Planck mass ≈ 2.18e-8 kg

fprintf('Planck units: l_P = %.3e m, m_P = %.3e kg\n', lP, mP);

% Schwarzschild BH: A = 4π r_s² = 16π G² M² / c⁴
A_bh     = @(M) 16*pi * G^2 * M.^2 / c^4;
S_bh     = @(M) kB * A_bh(M) / (4 * lP^2);

% In Planck units: S = 4π M² (dimensionless, with M in Planck masses)
S_planck = @(M_planck) 4*pi * M_planck.^2;

% Compare BH entropy with ordinary entropy ────────────────────────────────
M_sun_in_solar = 1;
S_sun_actual = 1e58 * kB;   % rough estimate: ~10^58 k_B
S_sun_as_BH  = S_bh(Msun);
fprintf('Sun entropy (actual):      %.2e k_B\n', S_sun_actual/kB);
fprintf('Sun as black hole entropy: %.2e k_B\n', S_sun_as_BH/kB);
fprintf('Ratio (BH/Sun): %.2e  (BH is vastly higher entropy!)\n\n', S_sun_as_BH/S_sun_actual);

% ── Plot S vs M ────────────────────────────────────────────────────────────
M_range = logspace(-3, 12, 500) * Msun;
S_range = S_bh(M_range);

figure('Name', 'Bekenstein-Hawking Entropy');
subplot(1,2,1);
  loglog(M_range/Msun, S_range/kB, 'b-', 'LineWidth', 2); grid on;
  xlabel('M/M_{sun}'); ylabel('S/k_B');
  title('Bekenstein-Hawking Entropy S = k_B A/4l_P^2');
  text(1, S_bh(Msun)/kB * 3, sprintf('S(M_\odot) = %.2e', S_bh(Msun)/kB), 'FontSize', 9);

subplot(1,2,2);
  % Holographic: S ∝ A (not Volume!) — plot S vs r_s
  r_range = 2*G*M_range/c^2;  % Schwarzschild radii
  loglog(r_range, S_range/kB, 'r-', 'LineWidth', 2); hold on;
  loglog(r_range, (r_range/lP).^3, 'b--', 'LineWidth', 1);  % Volume scaling comparison
  legend({'S_{BH} ∝ r^2 (area)', 'V ∝ r^3 (volume, for comparison)'}, 'Location', 'NW');
  grid on; xlabel('r_s (m)'); ylabel('S/k_B');
  title('Holographic Principle: S \propto Area, not Volume');

sgtitle('Bekenstein-Hawking Entropy — connecting GR, QM, and Thermodynamics');

7. Euclidean Path Integral & No-Boundary Saddle Point

Compute the Euclidean action for the de Sitter saddle point of the Hartle-Hawking path integral.

% ── Euclidean Path Integral: No-Boundary Wave Function ───────────────────
% Ψ[h] = ∫ D[g] e^(-S_E[g]/ℏ)  (Euclidean path integral over 4-geometries)
% Dominant contribution: de Sitter saddle, S_E = -3/(8Λ l_P²)
% This uses Wick rotation: e^(iS/ℏ) → e^(-S_E/ℏ), powered by Euler's formula

hbar=1.0546e-34; c=2.998e8; G=6.674e-11;
lP = sqrt(hbar*G/c^3);

% ── de Sitter space: the no-boundary saddle ───────────────────────────────
% In imaginary time τ, de Sitter space = 4-sphere of radius L = sqrt(3/Λ)
% ds² = dτ² + L²sin²(τ/L) dΩ³  (τ ∈ [0, πL])

Lambda_range = logspace(-120, -80, 500);  % cosmological constant (m^-2)
% Euclidean action for de Sitter: S_E = -3π/(G Λ)  (units G=ℏ=c=1)
% Wave function: Ψ ∝ e^(-S_E/ℏ) = e^(+3π/GΛ)  (note: negative action!)
% This means larger universes (smaller Λ) are more probable
G_nat = 1; hbar_nat = 1;  % natural units
S_E_dS = @(Lambda) -3*pi ./ (G_nat * Lambda);  % de Sitter Euclidean action

% ── de Sitter geometry in imaginary time ──────────────────────────────────
tau = linspace(0, pi, 200);
L = 1;   % de Sitter radius (normalized)
a_dS = L * sin(tau);  % scale factor: starts at 0, returns to 0

% ── Comparison: singular initial condition vs no-boundary ─────────────────
% Singular: a(τ) starts with singular derivative
a_singular = tau;  % linear start (conical singularity at τ=0)

figure('Name', 'No-Boundary Proposal');
subplot(2,2,1);
  plot(tau/pi, a_dS, 'b-', 'LineWidth', 2); hold on;
  plot([0 pi/2/pi], [0 1], 'r--', 'LineWidth', 2);
  grid on; xlabel('	au / \pi'); ylabel('a(	au)');
  legend({'No-boundary (sin)', 'Singular (linear)'});
  title('Scale factor in imaginary time');

subplot(2,2,2);
  % Show wave function amplitude vs Λ
  Lambda_vals = linspace(0.1, 5, 300);
  Psi_sq = exp(2 * 3*pi ./ Lambda_vals);    % |Ψ|² ∝ e^(+6π/Λ)
  Psi_sq = Psi_sq / max(Psi_sq);   % normalize
  semilogy(Lambda_vals, Psi_sq, 'g-', 'LineWidth', 2); grid on;
  xlabel('\Lambda (cosmological constant)'); ylabel('|\Psi|^2 (normalized)');
  title('Wave function |\Psi|^2 \propto e^{+3\pi/\Lambda}: prefers small \Lambda');

subplot(2,2,[3,4]);
  % Visualize the 4-sphere (de Sitter) in imaginary time
  [theta, phi] = meshgrid(linspace(0,pi,60), linspace(0,2*pi,60));
  X = sin(theta).*cos(phi); Y = sin(theta).*sin(phi); Z = cos(theta);
  surf(X, Y, Z, 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'FaceColor', 'interp');
  colormap(cool); axis equal; grid on;
  title('de Sitter 4-sphere (visualized as S²): no boundary, no singularity');
  xlabel('x'); ylabel('y'); zlabel('z');
  view(35, 25);

sgtitle('Hartle-Hawking No-Boundary Proposal via Euclidean Path Integral');

8. KMS Condition — Thermal Periodicity in Imaginary Time

Verify the KMS condition and demonstrate how periodicity in imaginary time encodes temperature.

% ── KMS Condition: G(t) = G(t + iβ) for thermal states ──────────────────
% This periodicity in imaginary time β = ℏ/k_B T is the fingerprint
% of a thermal state — directly from Euler's formula e^(iω(t+iβ)) = e^(iωt)e^(-ωβ)

kB = 1; hbar = 1;  % natural units
T = 0.5;             % temperature
beta = hbar / (kB * T);  % inverse temperature (= 2π/κ for Hawking)

% ── Simple harmonic oscillator Green's function ───────────────────────────
% G(t) = Tr[e^{-βH} A(t) B(0)] / Z
omega0 = 1.5;       % oscillator frequency
n_avg  = 1/(exp(beta*omega0) - 1);   % thermal average occupation

% Time-ordered propagator at finite temperature:
% G(t) = (n+1)e^(-iω₀t) + n·e^(+iω₀t)  — from Euler's formula
t_range = linspace(-3*beta, 3*beta, 1000);
G_t = (n_avg+1)*exp(-1i*omega0*t_range) + n_avg*exp(1i*omega0*t_range);

% ── Verify KMS: G(t) should equal G(t + iβ) ──────────────────────────────
t_test = 1.2;
G_at_t     = (n_avg+1)*exp(-1i*omega0*t_test) + n_avg*exp(1i*omega0*t_test);
G_at_t_iB  = (n_avg+1)*exp(-1i*omega0*(t_test+1i*beta)) + n_avg*exp(1i*omega0*(t_test+1i*beta));
fprintf('KMS check at t = %.2f:\n', t_test);
fprintf('  G(t)     = %.6f + %.6fi\n', real(G_at_t), imag(G_at_t));
fprintf('  G(t+iβ) = %.6f + %.6fi\n', real(G_at_t_iB), imag(G_at_t_iB));
fprintf('  |difference| = %.2e (should be ~0)\n\n', abs(G_at_t - G_at_t_iB));

% ── Plot G(t) in complex time plane ──────────────────────────────────────
figure('Name', 'KMS Condition');
subplot(1,2,1);
  plot(t_range/beta, real(G_t), 'b-', 'LineWidth', 2); hold on;
  plot(t_range/beta, imag(G_t), 'r-', 'LineWidth', 2);
  xlabel('t/eta'); legend({'Re G(t)', 'Im G(t)'}); grid on;
  title('Thermal Green''s function G(t)');
  xline(-1, 'k--'); xline(0, 'k--'); xline(1, 'k--');

subplot(1,2,2);
  % Periodicity in imaginary time: show G(it) — should be periodic with period β
  tau_range = linspace(0, 3*beta, 500);
  G_imag_t = (n_avg+1)*exp(-1i*omega0*(1i*tau_range)) + n_avg*exp(1i*omega0*(1i*tau_range));
  plot(tau_range/beta, real(G_imag_t), 'g-', 'LineWidth', 2); grid on;
  xlabel('	au / eta'); ylabel('G(i	au) — real');
  title('G(i	au): periodic with period eta = 1/T');
  xline([1 2 3], 'k--');  % period markers

sgtitle('KMS Condition: Thermal Periodicity = Euler\'s formula e^{i(t+i\beta)\omega} = e^{it\omega}e^{-\beta\omega}');

9. Bogoliubov Transformation — Particle Creation

Simulate the Bogoliubov transformation between in- and out-vacua that underlies Hawking radiation.

% ── Bogoliubov Transformation and Particle Creation ──────────────────────
% In curved spacetime, different observers define different vacua.
% Bogoliubov: ã_ω = ∫ [α_{ωω'} a_{ω'} + β_{ωω'} a†_{ω'}] dω'
% For Hawking: β_{ωω'} ≠ 0 → vacuum has particles for Schwarzschild observer
% The β coefficients arise from e^(iωt) having branch cut at horizon

kappa = 1.0;  % surface gravity

% Bogoliubov coefficients for Hawking radiation:
% |α_{ωω'}|² - |β_{ωω'}|² = δ(ω-ω') (normalization)
% Key ratio from analytic continuation: |β/α|² = e^(-2πω/κ)

omega_vals = linspace(0.01, 8, 500);
ratio_sq = exp(-2*pi * omega_vals / kappa);  % |β/α|² from Euler's formula
% From |α|² - |β|² = 1: |α|² = 1/(1-e^(-2πω/κ))
alpha_sq = 1 ./ (1 - exp(-2*pi*omega_vals/kappa));
beta_sq  = alpha_sq - 1;
% Thermal occupation number:
N_omega = beta_sq;  % = 1/(e^(2πω/κ) - 1) = Planck distribution!

fprintf('Bogoliubov coefficients at ω/κ = 1:\n');
idx = find(abs(omega_vals - kappa) == min(abs(omega_vals - kappa)));
fprintf('  |α|² = %.6f\n  |β|² = %.6f\n  N(ω) = %.6f\n', alpha_sq(idx), beta_sq(idx), N_omega(idx));
fprintf('  Check |α|²-|β|² = %.8f (should be 1)\n\n', alpha_sq(idx)-beta_sq(idx));

% Simulate squeezing — Bogoliubov transformation as optical squeezing
% â_out = cosh(r) â_in + sinh(r) â†_in  (squeezing parameter r = πω/κ)
r_vals = pi * linspace(0, 3, 200) / kappa;
n_squeezed = sinh(r_vals).^2;  % ⟨n⟩ = sinh²(r)

figure('Name', 'Bogoliubov Transformation');
subplot(2,2,1);
  semilogy(omega_vals/kappa, alpha_sq, 'b-', omega_vals/kappa, beta_sq, 'r-', 'LineWidth', 2);
  legend({'|α_{ωω}|²', '|β_{ωω}|² = ⟨N_ω⟩'}); grid on;
  xlabel('\omega/\kappa'); title('Bogoliubov coefficients');

subplot(2,2,2);
  plot(omega_vals/kappa, N_omega, 'g-', 'LineWidth', 2); grid on;
  xlabel('\omega/\kappa'); ylabel('⟨N_ω⟩');
  title('Particle content: Planck spectrum');

subplot(2,2,[3,4]);
  % Phase space picture of squeezing (Wigner function)
  [X, P] = meshgrid(linspace(-4,4,100), linspace(-4,4,100));
  r_sq = 1.2;  % squeezing for ω/κ = 1.2/π
  % Squeezed vacuum Wigner function:
  W = (2/pi) * exp(-2*(X.^2 * exp(-2*r_sq) + P.^2 * exp(2*r_sq)));
  contourf(X, P, W, 15, 'LineColor', 'none'); colormap(hot); colorbar;
  xlabel('X (quadrature)'); ylabel('P (quadrature)');
  title('Squeezed vacuum Wigner function (analog of Hawking state)');

sgtitle('Bogoliubov Transformation: vacuum → particles via complex phases');

10. The Page Curve — Information Recovery

Compute and visualize the Page curve showing how entanglement entropy of Hawking radiation should evolve if information is preserved.

% ── The Page Curve: Entanglement Entropy of Hawking Radiation ─────────────
% If information is preserved (unitary evolution via e^(-iHt/ℏ)),
% the entanglement entropy must follow the Page curve:
% - Rises from 0 until the "Page time" (≈ t_evap/2 in coarse approximation)
% - Then decreases back to 0 when BH fully evaporates
% Hawking's original calculation gives entropy that always increases (dashed)
% — the information paradox in a graph.

N_total = 1000;   % total Hilbert space dimension (proportional to BH entropy)
t_page  = 0.5;    % Page time (when BH has radiated half its information)

t_norm = linspace(0, 1, 1000);  % normalized evaporation time

% Remaining BH dimension (decreasing as it evaporates)
n_bh  = round(N_total * max(0, 1 - t_norm));
n_bh(n_bh < 1) = 1;

% Radiation Hilbert space dimension (increasing)
n_rad = max(1, N_total - n_bh);

% Page entropy: S = log(min(n_bh, n_rad)) for maximally entangled subsystems
S_page = log(min(n_bh, n_rad));  % this is the Page curve

% Hawking result: entropy always increases (never decreases)
S_hawking = log(n_rad);  % keeps growing even after Page time

% Max entropy of BH (Bekenstein-Hawking)
S_BH_max = log(N_total) * (1 - t_norm);  % decreasing entropy of BH

figure('Name', 'Page Curve — Information Paradox');
subplot(1,2,1);
  plot(t_norm, S_hawking/log(N_total), 'r-', 'LineWidth', 2); hold on;
  plot(t_norm, S_page/log(N_total), 'b-', 'LineWidth', 2.5);
  plot(t_norm, S_BH_max/log(N_total), 'g--', 'LineWidth', 1.5);
  xline(t_page, 'k:', 'Page time');
  xlabel('t / t_{evap}'); ylabel('S / S_{max}'); grid on;
  legend({'Hawking (info lost)', 'Page curve (unitary)', 'S_{BH} of black hole'});
  title('The Page Curve');

subplot(1,2,2);
  % Mutual information: I = S_A + S_B - S_AB
  % For unitarity: I should grow, encoding information leaking out
  S_total = log(N_total) * ones(size(t_norm));  % pure state: constant
  I_mutual_hawking = S_hawking + S_BH_max - S_total;
  I_mutual_page    = S_page    + S_BH_max - S_total;
  I_mutual_hawking = max(0, I_mutual_hawking);
  I_mutual_page    = max(0, I_mutual_page);
  plot(t_norm, I_mutual_hawking/log(N_total), 'r-', 'LineWidth', 2); hold on;
  plot(t_norm, I_mutual_page/log(N_total), 'b-', 'LineWidth', 2.5);
  xline(t_page, 'k:', 'Page time'); grid on;
  xlabel('t / t_{evap}'); ylabel('Mutual Info / S_{max}');
  legend({'Hawking', 'Page (unitary)'});
  title('Mutual Information: does radiation remember the BH?');

sgtitle('Page Curve: Information preserved (e^{-iHt} unitary) vs Hawking''s paradox');