import math

# ---- constants ----
MP   = 1.2209e19      # GeV (Planck mass)
alpha= 1/137.036
fa   = MP*alpha**5    # GeV, forced decay constant
mchi = 95.6e-6        # GeV (95.6 keV)
s0   = 2891.2         # cm^-3, entropy density today
rho_crit_over_h2 = 1.0537e-5  # GeV cm^-3
Omega_b_h2 = 0.0224

print("fa (M_P alpha^5)      =", f"{fa:.3e}", "GeV   (corpus 5.04e7)")
print("m_chi                 =", mchi*1e6, "keV")

# ---- target yield to give Omega_chi h^2 = 0.12 ----
# Omega_chi h^2 = m_chi * Y_chi * s0 / (rho_crit/h2)
# Y_chi_target:
Omega_target = 0.12
Y_target = Omega_target*rho_crit_over_h2/(mchi*s0)
print("\nY_chi target (for Omega h^2=0.12) =", f"{Y_target:.3e}")

# baryon yield for comparison
Y_b = Omega_b_h2*rho_crit_over_h2/(0.938*s0)
print("Y_b (baryon)                       =", f"{Y_b:.3e}")
print("ratio target Y_chi/Y_b             =", f"{Y_target/Y_b:.2f}", "  (16/3 in *energy* uses m_chi/m_p)")

# energy ratio check: Omega_chi/Omega_b
print("Omega_chi/Omega_b at target        =", f"{Omega_target/Omega_b_h2:.2f}", "  (target 16/3=5.33)")

# ---- INSTANTON-DRIVEN production yield ----
# S3 Chern-Simons instanton: dominant k=1 winding, S_inst = 2*pi*k = 2*pi
S_inst = 2*math.pi
print("\n--- Instanton sector ---")
print("S_inst (k=1, =2*pi*k)  =", f"{S_inst:.4f}")
print("exp(-S_inst)           =", f"{math.exp(-S_inst):.3e}")

# Instanton-induced PRODUCTION yield.
# Generic instanton/non-thermal production: the instanton vertex injects chi quanta
# with a rate per Hubble suppressed by the instanton fugacity e^{-S_inst}.
# The number density produced is n_chi ~ (instanton density) ~ K * (T or H scale)^3 * e^{-S_inst}
# In an instanton-condensate/"tunnel" picture the comoving yield freezes at
#   Y_chi ~ (n_inst/s) ~ g_eff^{-1} * e^{-S_inst}   (order-one topological prefactor, T_RH-independent)
# i.e. one chi quantum per instanton, instanton number density ~ e^{-S_inst} * (thermal phase space)

# Most-optimistic *topological* estimate (one chi per unit-winding instanton, per dof):
# the dilute-instanton-gas number density in comoving units, normalized to entropy:
#   Y_inst = (45/(2 pi^2 g_*S)) * (instanton fugacity)   with g_*S ~ 106.75
g_star = 106.75
Y_inst_topo = (45/(2*math.pi**2*g_star))*math.exp(-S_inst)
print("Y_chi (one-per-instanton, e^{-2pi}) =", f"{Y_inst_topo:.3e}")
Omega_inst = mchi*Y_inst_topo*s0/rho_crit_over_h2
print("Omega_chi h^2 (this estimate)       =", f"{Omega_inst:.3e}")
print("  vs target 0.12: factor", f"{Omega_inst/0.12:.2e}")
print("  vs misalignment 1e-9: factor", f"{Omega_inst/1e-9:.2e}")

# ---- compare exponential suppression structures ----
print("\n--- the core scaling problem ---")
print("misalignment Omega ~ 1e-9 (DILUTION-driven, m_chi huge -> early osc -> redshift kill)")
print("instanton e^{-2pi}  =", f"{math.exp(-2*math.pi):.3e}", " (single power, NOT a dilution)")
# what S_inst would be needed to hit Y_target with one-per-instanton?
# Y_target = (45/(2pi^2 g)) e^{-S_need}
S_need = -math.log(Y_target*2*math.pi**2*g_star/45)
print("S_inst needed to hit target          =", f"{S_need:.3f}")
print("vs S_inst=2pi=", f"{2*math.pi:.3f}", " -> instanton OVERPRODUCES if one-per-instanton at full thermal phase space")

# More careful: instanton number density in 4D ~ integral d rho/rho^5 D(rho) with
# D(rho) ~ rho^{-5} (2pi/alpha_s)^{2N} e^{-2pi/alpha} for gauge instantons; here CS on S3 gives e^{-2pi k}.
