Computational Appendices

from sympy import symbols, Eq, solve, sqrt, N

# Define symbols (all positive for physical quantities)
omega, k, xi, M_Pl, eta, Delta_t, c = symbols('omega k xi M_Pl eta Delta_t c', positive=True)

# Dispersion equation
disp_eq = Eq(omega**2, k**2 * (1 + xi**2 * (k / M_Pl)**2 + eta * k * Delta_t / c))

# Solve for omega (positive root for physical frequency)
solution_omega = solve(disp_eq, omega)[1]  # Index 1 for positive root

# Print symbolic solution
print("Dispersion Equation:", disp_eq)
print("ω(k) solution:", solution_omega)

# Numerical evaluation with LISA-relevant parameters
# k = 1e-10 Hz (LISA band)
# ξ = 1e10 (loop correction)
# M_Pl = 1e19 GeV ~ 1e28 Hz (natural units)
# η = 0.1 (multi-time coupling)
# Δt_ortho = 1e-18 s (compact radius)
# c = 1 (natural units)

params = {k: 1e-10, xi: 1e10, M_Pl: 1e19, eta: 0.1, Delta_t: 1e-18, c: 1}
num_omega = N(solution_omega.subs(params))

print("\nNumerical Results:")
print(f"ω = {num_omega} Hz")

# Calculate the correction terms separately
quadratic_term = N((k * xi / M_Pl)**2).subs(params)
linear_term = N(eta * k * Delta_t / c).subs(params)

print(f"\nQuadratic correction: {quadratic_term}")
print(f"Linear correction (multi-time): {linear_term}")
print(f"Total correction: {quadratic_term + linear_term}")

from qutip import *
import numpy as np

# System size: position basis truncation
N = 256  # Sufficient for convergence

# Position and momentum operators
x = position(N)
p = momentum(N)

# Moduli potential parameters (swampland compliant, 2T framework)
F = 1.0        # F-term (normalized)
a = 1.414      # √(26/13) from D→13D via Sp(2,R); > √(2/3) ≈ 0.816 (swampland)
kappa = 1.0    # Uplift strength
b = 1.0        # Instanton parameter
mu = 0.5       # 2T periodic modulation from t_ortho
R = 1.0        # Ortho time compactification radius

# Construct potential as diagonal matrix in position basis
x_vals = x.diag()  # Extract diagonal elements

# V(φ) = |F|²e^(-aφ) + κe^(-b/φ) + μcos(φ/R)
# Add small offset 1e-10 to avoid division by zero
V_mat = Qobj(np.diag(
    F**2 * np.exp(-a * x_vals) +
    kappa * np.exp(-b / (x_vals + 1e-10)) +
    mu * np.cos(x_vals / R)
))

# Hamiltonian: H = p²/2 + V(φ)
H = p**2 / 2 + V_mat

# Initial state: coherent state (Gaussian wave packet)
# Represents quantum fluctuation around moduli minimum
alpha = 5.0  # Displacement in phase space
psi0 = coherent(N, alpha)

# Time evolution using master equation solver
times = np.linspace(0, 10, 100)  # Evolution from t=0 to t=10
result = mesolve(H, psi0, times)

# Calculate expectation values
expect_x = expect(x, result.states)  # <φ(t)>
expect_p = expect(p, result.states)  # <p(t)>

# Calculate von Neumann entropy for each state
# S = -Tr(ρ ln ρ) should remain ~0 for unitary evolution
entropy = [entropy_vn(state) for state in result.states]

# Output results
print("=== Moduli Potential Quantum Simulation ===")
print(f"\nInitial state (t=0):")
print(f"  <φ> = {expect_x[0]:.3f}")
print(f"  <p> = {expect_p[0]:.3f}")
print(f"  Entropy = {entropy[0]:.2e}")

print(f"\nFinal state (t=10):")
print(f"  <φ> = {expect_x[-1]:.3f}")
print(f"  <p> = {expect_p[-1]:.3f}")
print(f"  Entropy = {entropy[-1]:.2e}")

print(f"\nStability analysis:")
print(f"  φ drift: {abs(expect_x[-1] - expect_x[0]):.3f}")
print(f"  Max entropy: {max(entropy):.2e}")
print(f"  Unitarity check: Entropy < 1e-10? {max(entropy) < 1e-10}")

# Optional: plot evolution
# import matplotlib.pyplot as plt
# plt.figure(figsize=(12, 4))
# plt.subplot(131)
# plt.plot(times, expect_x)
# plt.xlabel('Time'); plt.ylabel('<φ>'); plt.title('Position')
# plt.subplot(132)
# plt.plot(times, expect_p)
# plt.xlabel('Time'); plt.ylabel('<p>'); plt.title('Momentum')
# plt.subplot(133)
# plt.plot(times, entropy)
# plt.xlabel('Time'); plt.ylabel('S'); plt.title('Entropy')
# plt.tight_layout(); plt.show()

from sympy import symbols, exp, pi, solve, N, sqrt

# Define symbols
Delta_V, sigma, M_Pl, G = symbols('Delta_V sigma M_Pl G', positive=True)

# Bubble radius (thin-wall approximation)
r_b = 3 * sigma / (4 * Delta_V)

# Euclidean action (CDL formula)
S_E = 27 * pi**2 * sigma**4 / (2 * Delta_V**3)

# Tunneling rate Γ ~ exp(-S_E)
Gamma = exp(-S_E)

# Print symbolic expressions
print("=== CDL Instanton Tunneling Rate ===")
print(f"\nBubble Radius: r_b = {r_b}")
print(f"Euclidean Action: S_E = {S_E}")
print(f"Tunneling Rate: Γ = {Gamma}")

# Scenario 1: High barrier (Planck-scale landscape)
print("\n--- Scenario 1: Planck-scale barrier ---")
params1 = {sigma: 1e12, Delta_V: 1e76}  # σ ~ TeV³, ΔV ~ M_Pl⁴/1000
r_b1 = N(r_b.subs(params1))
S_E1 = N(S_E.subs(params1))
Gamma1 = N(Gamma.subs(params1))

print(f"Parameters: σ = 10^12 GeV³, ΔV = 10^76 GeV⁴")
print(f"Bubble radius: r_b = {r_b1:.2e} cm (sub-Planck!)")
print(f"Euclidean action: S_E = {S_E1:.2e}")
print(f"Tunneling rate: Γ ≈ {Gamma1} (exponentially suppressed)")
print("Conclusion: No observable bubbles in Hubble volume")

# Scenario 2: TeV-scale barrier (testable)
print("\n--- Scenario 2: TeV-scale barrier (testable) ---")
params2 = {sigma: 1e12, Delta_V: 1e60}  # σ ~ TeV³, ΔV ~ (10 TeV)⁴
r_b2 = N(r_b.subs(params2))
S_E2 = N(S_E.subs(params2))
Gamma2 = N(Gamma.subs(params2))

print(f"Parameters: σ = 10^12 GeV³, ΔV = 10^60 GeV⁴")
print(f"Bubble radius: r_b = {r_b2:.2e} cm (macroscopic!)")
print(f"Euclidean action: S_E = {S_E2:.2e}")
print(f"Tunneling rate: Γ ≈ {Gamma2:.2e}")

# Estimate number of bubbles in observable universe
H_0 = 2.2e-18  # Hubble rate in s^-1
t_universe = 4.4e17  # Age of universe in s
V_Hubble = (3e10 / H_0)**3  # Hubble volume in cm³
N_bubbles = Gamma2 * V_Hubble * t_universe

print(f"\nExpected bubbles in observable universe: N ~ {N_bubbles:.2e}")
if N_bubbles > 0.001:
    print("Conclusion: TESTABLE via CMB bubble collision searches!")
else:
    print("Conclusion: Too rare for current CMB sensitivity")

# Scenario 3: Intermediate barrier
print("\n--- Scenario 3: Intermediate barrier ---")
params3 = {sigma: 1e12, Delta_V: 1e65}  # ΔV ~ (100 TeV)⁴
r_b3 = N(r_b.subs(params3))
S_E3 = N(S_E.subs(params3))
Gamma3 = N(Gamma.subs(params3))
N_bubbles3 = Gamma3 * V_Hubble * t_universe

print(f"Parameters: σ = 10^12 GeV³, ΔV = 10^65 GeV⁴")
print(f"Bubble radius: r_b = {r_b3:.2e} cm")
print(f"Euclidean action: S_E = {S_E3:.2e}")
print(f"Tunneling rate: Γ ≈ {Gamma3:.2e}")
print(f"Expected bubbles: N ~ {N_bubbles3:.2e}")

from sympy import symbols, exp, log, pi, N, sqrt

# Define symbols
delta_0, sigma, n, A = symbols('delta_0 sigma n A', positive=True)

# Probability of cold spot deeper than δ_0 (Gumbel tail)
P_gumbel = exp(-n * A * exp(-delta_0**2 / (2 * sigma**2)))

print("=== Gaussian CMB Cold Spot Statistics ===")
print(f"\nGumbel Probability: P(δ > δ₀) = {P_gumbel}")

# Numerical evaluation: Cold Spot at 3σ
# Parameters from Planck:
#   σ = 3×10^-3 (CMB RMS)
#   θ = 1° = 0.017 rad (spot angular size)
#   n = 3/(2π θ²) ≈ 3/(2π × 0.017²) ≈ 1650 sr^-1
#   A = 4π ≈ 12.566 sr (full sky)
#   δ₀ = 3σ = 9×10^-3

params_gauss = {
    sigma: 3e-3,
    n: 1650,
    A: 4 * pi,
    delta_0: 3 * 3e-3  # 3σ spot
}

P_3sigma = N(P_gumbel.subs(params_gauss))

print(f"\n--- Cold Spot (3σ) Probability ---")
print(f"Parameters: σ = 3×10⁻³, θ = 1°, δ₀ = 3σ")
print(f"Minima density: n ≈ 1650 sr⁻¹")
print(f"Sky area: A = 4π ≈ 12.57 sr")
print(f"Probability: P(δ > 3σ) ≈ {P_3sigma:.4f}")
print(f"Interpretation: {P_3sigma*100:.2f}% chance (not anomalous)")

# For 5σ (true anomaly threshold)
params_5sigma = params_gauss.copy()
params_5sigma[delta_0] = 5 * 3e-3
P_5sigma = N(P_gumbel.subs(params_5sigma))

print(f"\n--- Rare Event (5σ) Probability ---")
print(f"Probability: P(δ > 5σ) ≈ {P_5sigma:.2e}")
print(f"Expected number on full sky: {P_5sigma * 12.566:.4f}")

# Gumbel parameters for extremes
N_spots = 4 * pi / (0.017**2)  # Number of 1° patches
mu_gumbel = -sigma * sqrt(2 * log(N_spots))
beta_gumbel = sigma / sqrt(2 * log(N_spots))

print(f"\n--- Gumbel Fit Parameters ---")
print(f"Number of patches: N ≈ {N_spots:.0f}")
print(f"Location μ = {N(mu_gumbel.subs(sigma, 3e-3)):.2e}")
print(f"Scale β = {N(beta_gumbel.subs(sigma, 3e-3)):.2e}")

from sympy import symbols, exp, factorial, N

# Define symbols
N_spots, lambda_poiss = symbols('N_spots lambda_poiss', positive=True)

# Poisson distribution P(N) = λ^N exp(-λ) / N!
P_poisson = lambda_poiss**N_spots * exp(-lambda_poiss) / factorial(N_spots)

print("\n=== Multiverse Bubble Collision Statistics ===")
print(f"\nPoisson Distribution: P(N) = {P_poisson}")

# Expected rate from CDL (Appendix C, Scenario 2)
# λ = Γ × V_Hubble × t_universe / (Hubble 4-volume)
# For testable regime: λ ~ 10^-3 (one bubble per 1000 Hubble volumes)

lambda_vals = [0.001, 0.01, 0.1, 1.0]

for lam in lambda_vals:
    print(f"\n--- λ = {lam} (avg bubbles per Hubble volume) ---")

    # Probability of 0, 1, 2 bubbles
    P0 = N(P_poisson.subs({lambda_poiss: lam, N_spots: 0}))
    P1 = N(P_poisson.subs({lambda_poiss: lam, N_spots: 1}))
    P2 = N(P_poisson.subs({lambda_poiss: lam, N_spots: 2}))

    print(f"P(0 bubbles) = {P0:.4f}")
    print(f"P(1 bubble) = {P1:.4f}")
    print(f"P(2+ bubbles) = {1 - P0 - P1:.4f}")

    if P1 > 0.01:
        print("→ TESTABLE: >1% chance of observable collision")
    else:
        print("→ Too rare for current experiments")

# Non-Gaussian kurtosis from bubble disks
print("\n=== Kurtosis Calculation (Non-Gaussianity Test) ===")

# Kurtosis κ = <(ΔT)⁴> / <(ΔT)²>²
# Gaussian: κ = 3
# Bubble excess: κ = 3 + Σ(δ_disk⁴) / σ⁴

sigma_CMB = 3e-3  # CMB RMS
delta_disk = 0.1   # Bubble collision ΔT/T ~ -0.1
N_disk = 1         # Number of disk collisions

excess = N_disk * delta_disk**4 / sigma_CMB**4
kurtosis = 3 + excess

print(f"Parameters: σ = {sigma_CMB}, δ_disk = {delta_disk}, N = {N_disk}")
print(f"Gaussian kurtosis: κ = 3")
print(f"Excess from bubbles: {excess:.2e}")
print(f"Total kurtosis: κ = {kurtosis:.2e}")
print(f"Sigma deviation: {(kurtosis - 3) / 0.1:.1e}σ (if error ~ 0.1)")

if kurtosis > 3.5:
    print("→ SMOKING GUN: Huge non-Gaussianity, supports multiverse!")
else:
    print("→ Marginal excess, need better statistics")

from sympy import symbols, Eq, solve, sqrt, N
import matplotlib.pyplot as plt
import numpy as np

# Define all symbols as positive (physical quantities)
omega, k, xi, M_Pl, eta, Delta_t, c = symbols('omega k xi M_Pl eta Delta_t c', positive=True)

# Dispersion equation from modified gravity
disp_eq = Eq(omega**2, k**2 * (1 + xi**2 * (k / M_Pl)**2 + eta * k * Delta_t / c))

# Solve for omega - take index [1] for positive root
solution_omega = solve(disp_eq, omega)[1]  # Positive root for physical ω > 0

# Print symbolic results
print("=== SymPy GW Dispersion Derivation ===")
print(f"\nDispersion Equation: {disp_eq}")
print(f"\nSymbolic ω(k) solution:")
print(f"  ω = {solution_omega}")

# Numerical evaluation with LISA parameters
# k = 1e-10 Hz (representative LISA frequency in low band)
# ξ = 1e10 (from loop corrections)
# M_Pl = 1e19 GeV ~ 1e28 Hz (in natural units, using E=hf conversion)
# η = 0.1 (baseline multi-time coupling from g/E_F)
# Δt_ortho = 1e-18 s (from compact radius R ~ TeV^-1)
# c = 1 (natural units)

params_baseline = {
    k: 1e-10,
    xi: 1e10,
    M_Pl: 1e19,
    eta: 0.1,
    Delta_t: 1e-18,
    c: 1
}

num_omega = N(solution_omega.subs(params_baseline))

print(f"\n--- Numerical Evaluation (Baseline) ---")
print(f"Parameters: k=1e-10 Hz, ξ=1e10, M_Pl=1e19 GeV, η=0.1, Δt=1e-18 s")
print(f"Numerical ω: {num_omega} Hz")

# Calculate individual correction terms for analysis
quadratic_correction = N((k * xi / M_Pl)**2).subs(params_baseline)
linear_correction = N((eta * k * Delta_t / c)).subs(params_baseline)
total_correction = quadratic_correction + linear_correction

print(f"\nCorrection terms:")
print(f"  Quadratic term (standard QG): ξ²(k/M_Pl)² = {quadratic_correction:.2e}")
print(f"  Linear term (multi-time):     η k Δt/c   = {linear_correction:.2e}")
print(f"  Total correction:                         = {total_correction:.2e}")
print(f"  Boost factor: {linear_correction / quadratic_correction:.2e}")

# Test boosted scenario (asymptotic safety with UV fixed point)
params_boosted = params_baseline.copy()
params_boosted[eta] = 1e9  # Boost η from RG UV fixed point

num_omega_boosted = N(solution_omega.subs(params_boosted))
linear_correction_boosted = N((eta * k * Delta_t / c)).subs(params_boosted)

print(f"\n--- Boosted Scenario (η = 10^9) ---")
print(f"Linear correction: {linear_correction_boosted:.2e}")
print(f"ω (boosted): {num_omega_boosted:.2e} Hz")
print(f"Relative deviation: δω/k = {(num_omega_boosted / 1e-10 - 1):.2e}")

if linear_correction_boosted > 1e-20:
    print("→ TESTABLE by LISA! (threshold ~10^-20)")
else:
    print("→ Still below LISA sensitivity")

# Plot ω(k) to visualize deviation from linear dispersion
print(f"\n--- Generating Plot: ω(k) vs k ---")

# k range from 1e-12 to 1e-8 Hz (broad LISA band)
k_vals = np.logspace(-12, -8, 100)

# Calculate ω(k) for baseline and boosted scenarios
omega_baseline = []
omega_boosted = []
deviation_baseline = []
deviation_boosted = []

for k_val in k_vals:
    params_temp = params_baseline.copy()
    params_temp[k] = k_val
    omega_val = float(N(solution_omega.subs(params_temp)))
    omega_baseline.append(omega_val)
    deviation_baseline.append(omega_val / k_val - 1)  # Fractional deviation

    params_temp[eta] = 1e9
    omega_val_boost = float(N(solution_omega.subs(params_temp)))
    omega_boosted.append(omega_val_boost)
    deviation_boosted.append(omega_val_boost / k_val - 1)

# Create figure with subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: ω(k) vs k
ax1.loglog(k_vals, omega_baseline, label='Baseline (η=0.1)', linewidth=2)
ax1.loglog(k_vals, omega_boosted, label='Boosted (η=10⁹)', linewidth=2, linestyle='--')
ax1.loglog(k_vals, k_vals, 'k:', label='Linear ω=k', alpha=0.5)
ax1.set_xlabel('Frequency k (Hz)', fontsize=12)
ax1.set_ylabel('ω (Hz)', fontsize=12)
ax1.set_title('GW Dispersion ω(k)', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Deviation ω/k - 1
ax2.loglog(k_vals, np.abs(deviation_baseline), label='Baseline (η=0.1)', linewidth=2)
ax2.loglog(k_vals, np.abs(deviation_boosted), label='Boosted (η=10⁹)', linewidth=2, linestyle='--')
ax2.axhline(y=1e-20, color='r', linestyle=':', label='LISA threshold', linewidth=2)
ax2.set_xlabel('Frequency k (Hz)', fontsize=12)
ax2.set_ylabel('|ω/k - 1| (fractional deviation)', fontsize=12)
ax2.set_title('Deviation from Linear Dispersion', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('gw_dispersion_analysis.png', dpi=150, bbox_inches='tight')
print("Plot saved as: gw_dispersion_analysis.png")
plt.show()

from qutip import *
import numpy as np
import matplotlib.pyplot as plt

# System size: position basis truncation
# N=256 is sufficient for convergence (smooth potential, localized states)
N = 256

# Position and momentum operators in N-dimensional Hilbert space
x = position(N)  # QuTiP position operator
p = momentum(N)  # QuTiP momentum operator

# Moduli potential parameters (swampland-compliant)
F = 1.0        # F-term (normalized to 1 for simplicity)
a = 1.414      # √(26/13) > √(2/3) ≈ 0.816 (swampland constraint satisfied!)
kappa = 1.0    # Uplift strength (non-perturbative instanton amplitude)
b = 1.0        # Instanton parameter ~ 2π/(n g_s)
mu = 0.5       # Periodic modulation from multi-time coupling
R = 1.0        # Ortho radius (normalized)

# Construct potential V(φ) as diagonal matrix in position basis
# This is the key step: converting field theory potential to quantum operator
x_vals = x.diag()  # Extract diagonal elements (eigenvalues of x)

# V(φ) = |F|²e^(-aφ) + κe^(-b/φ) + μcos(φ/R)
# Add small offset 1e-10 to x_vals in second term to avoid division by zero
V_array = (
    F**2 * np.exp(-a * x_vals) +           # GKP flux stabilization
    kappa * np.exp(-b / (x_vals + 1e-10)) +  # Non-pert uplift (avoid x=0)
    mu * np.cos(x_vals / R)                 # Multi-time periodic term
)

# Convert to QuTiP quantum object (diagonal operator)
V_mat = Qobj(np.diag(V_array))

# Hamiltonian: H = p²/2 + V(φ)
# This is the full quantum Hamiltonian for the moduli field
H = p**2 / 2 + V_mat

# Initial state: coherent state (Gaussian wave packet)
# Coherent state |α⟩ is displaced ground state of harmonic oscillator
# α=5.0 centers packet near potential minimum (from preliminary analysis)
alpha = 5.0
psi0 = coherent(N, alpha)

print("=== QuTiP Moduli Potential Simulation ===")
print(f"\nSystem parameters:")
print(f"  Basis size N = {N}")
print(f"  Swampland parameter a = {a:.3f} (required: a > {np.sqrt(2/3):.3f})")
print(f"  Swampland check: {'PASS' if a > np.sqrt(2/3) else 'FAIL'}")
print(f"  Initial coherent state: α = {alpha}")

# Time evolution using master equation solver (mesolve)
# mesolve handles both unitary and non-unitary evolution
# For pure Hamiltonian (no dissipation), this reduces to Schrödinger equation
times = np.linspace(0, 10, 100)  # Evolution from t=0 to t=10 (normalized time units)

print(f"\nEvolving system from t=0 to t=10...")
result = mesolve(H, psi0, times)

# Calculate expectation values ⟨x(t)⟩ and ⟨p(t)⟩
expect_x = expect(x, result.states)  # Position expectation ⟨φ(t)⟩
expect_p = expect(p, result.states)  # Momentum expectation ⟨p(t)⟩

# Calculate von Neumann entropy for each state
# S = -Tr(ρ ln ρ) should remain ~0 for pure state unitary evolution
# Any increase in S indicates decoherence or numerical errors
entropy = [entropy_vn(state) for state in result.states]

# Analysis: Compare initial and final states
print(f"\n--- Initial State (t=0) ---")
print(f"  ⟨φ⟩ = {expect_x[0]:.3f}")
print(f"  ⟨p⟩ = {expect_p[0]:.3f}")
print(f"  Entropy S = {entropy[0]:.2e}")

print(f"\n--- Final State (t=10) ---")
print(f"  ⟨φ⟩ = {expect_x[-1]:.3f}")
print(f"  ⟨p⟩ = {expect_p[-1]:.3f}")
print(f"  Entropy S = {entropy[-1]:.2e}")

# Stability analysis
phi_drift = abs(expect_x[-1] - expect_x[0])
max_entropy = max(entropy)

print(f"\n--- Stability Analysis ---")
print(f"  φ drift: Δφ = {phi_drift:.3f}")
print(f"  Drift percentage: {100 * phi_drift / expect_x[0]:.2f}%")
print(f"  Max entropy: {max_entropy:.2e}")
print(f"  Unitarity check (S < 10^-10): {'PASS' if max_entropy < 1e-10 else 'FAIL'}")

if phi_drift < 1.0 and max_entropy < 1e-10:
    print(f"\n→ CONCLUSION: Potential is STABLE! Wave function localized, no runaway.")
else:
    print(f"\n→ WARNING: Potential may be unstable or numerical issues present.")

# Visualization: Plot ⟨φ(t)⟩, ⟨p(t)⟩, and S(t)
print(f"\n--- Generating Plots ---")

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Plot 1: ⟨φ(t)⟩ evolution
axes[0].plot(times, expect_x, linewidth=2, color='#8b7fff')
axes[0].set_xlabel('Time t', fontsize=12)
axes[0].set_ylabel('⟨φ(t)⟩', fontsize=12)
axes[0].set_title('Position Expectation', fontsize=14)
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=expect_x[0], color='gray', linestyle=':', label='Initial value')
axes[0].legend()

# Plot 2: ⟨p(t)⟩ evolution
axes[1].plot(times, expect_p, linewidth=2, color='#ff7eb6')
axes[1].set_xlabel('Time t', fontsize=12)
axes[1].set_ylabel('⟨p(t)⟩', fontsize=12)
axes[1].set_title('Momentum Expectation', fontsize=14)
axes[1].grid(True, alpha=0.3)

# Plot 3: Entropy S(t) - key unitarity check
axes[2].plot(times, entropy, linewidth=2, color='#ffd93d')
axes[2].set_xlabel('Time t', fontsize=12)
axes[2].set_ylabel('von Neumann Entropy S', fontsize=12)
axes[2].set_title('Entropy (Unitarity Check)', fontsize=14)
axes[2].set_yscale('log')  # Log scale to see tiny values
axes[2].grid(True, alpha=0.3)
axes[2].axhline(y=1e-10, color='r', linestyle=':', label='Unitarity threshold')
axes[2].legend()

plt.tight_layout()
plt.savefig('moduli_evolution.png', dpi=150, bbox_inches='tight')
print("Plot saved as: moduli_evolution.png")
plt.show()

# Optional: Plot potential V(φ) to visualize minimum
fig2, ax = plt.subplots(figsize=(10, 6))
phi_range = np.linspace(0.1, 15, 500)
V_vals = (
    F**2 * np.exp(-a * phi_range) +
    kappa * np.exp(-b / phi_range) +
    mu * np.cos(phi_range / R)
)

ax.plot(phi_range, V_vals, linewidth=2, color='#8b7fff')
ax.axvline(x=expect_x[0], color='green', linestyle='--', label=f'Initial ⟨φ⟩={expect_x[0]:.2f}')
ax.axvline(x=expect_x[-1], color='orange', linestyle='--', label=f'Final ⟨φ⟩={expect_x[-1]:.2f}')
ax.set_xlabel('φ (modulus field)', fontsize=12)
ax.set_ylabel('V(φ)', fontsize=12)
ax.set_title('Moduli Potential V(φ) = |F|²e^(-aφ) + κe^(-b/φ) + μcos(φ/R)', fontsize=14)
ax.grid(True, alpha=0.3)
ax.legend()
ax.set_ylim([0, 5])  # Focus on relevant region

plt.tight_layout()
plt.savefig('moduli_potential.png', dpi=150, bbox_inches='tight')
print("Potential plot saved as: moduli_potential.png")
plt.show()