#!/usr/bin/env python3 """how-hot-is-the-core-v1.py -- the inference chain behind "the core is about 6,000 K", re-run from its inputs. Companion to the solvetheuniverse.com article "How Hot Is the Core, and How Would Anyone Know". Nobody has measured the core's temperature. The number is inferred: 1. Seismology fixes where the outer core freezes (the inner-core boundary, ICB). 2. The density profile fixes the PRESSURE at that boundary. <- this script 3. A lab measures iron's melting point at that pressure. <- inputs, cited 4. An adiabat carries the number out to the core-mantle boundary. <- this script Step 2 uses the PREM density polynomials (Dziewonski & Anderson 1981, Phys. Earth Planet. Inter. 25, 297-356, Table I). The coefficients below were typed by hand; the self-test checks the transcription by CONSEQUENCE, against numbers the paper itself publishes: Earth's mass 5.974e24 kg, the core-mantle boundary pressure 135.8 GPa, the inner-core boundary pressure 328.9 GPa, and the density on each side of both boundaries (12.7636/12.1663 and 9.9035/5.5665 g/cm^3). What these checks can and cannot catch, measured by two independent reviewers on 2026-09-14 by injecting single-digit typos: the integrals alone let most typos through (a local error averages away); the boundary densities catch leading-digit errors in the core and lower mantle and nothing in the seven layers above. The real safeguard is that every coefficient below was checked line by line against Table I of the paper on 2026-09-14. The self-test is a regression guard, not a proof. Run: python3 how-hot-is-the-core-v1.py python3 how-hot-is-the-core-v1.py --selftest No dependencies beyond the standard library. """ import math, sys G = 6.674e-11 # m^3 kg^-1 s^-2 A_KM = 6371.0 # PREM Earth radius, km KM = 1000.0 # (outer radius km, polynomial in x = r/a, density in g/cm^3) -- PREM Table I LAYERS = [ (1221.5, (13.0885, 0.0, -8.8381, 0.0)), # inner core (3480.0, (12.5815, -1.2638, -3.6426, -5.5281)), # outer core (5701.0, (7.9565, -6.4761, 5.5283, -3.0807)), # lower mantle (5771.0, (5.3197, -1.4836, 0.0, 0.0)), # transition zone (5971.0, (11.2494, -8.0298, 0.0, 0.0)), (6151.0, (7.1089, -3.8045, 0.0, 0.0)), (6346.6, (2.6910, 0.6924, 0.0, 0.0)), # LVZ + LID (6356.0, (2.900, 0.0, 0.0, 0.0)), # lower crust (6368.0, (2.600, 0.0, 0.0, 0.0)), # upper crust (6371.0, (1.020, 0.0, 0.0, 0.0)), # ocean ] R_ICB_KM, R_CMB_KM = 1221.5, 3480.0 # Published checks (same paper): mass, P_CMB, P_ICB. PUBLISHED = {"mass_kg": 5.974e24, "P_CMB_GPa": 135.8, "P_ICB_GPa": 328.9} # Iron melting point at the ICB pressure, four independent determinations. # Each is an INPUT to this chain, not an output of it. ICB_ESTIMATES = [ ("Boehler 1993, Nature 363:534, laser-heated DAC to 200 GPa, melting judged by watching the surface move, extrapolated", 4850, 200), ("Alfe, Price & Gillan 2002, PRB 65:165118, ab initio free-energy calculation (+/-300 is statistical only, per the group's 2004 review)", 6350, 300), ("Anzellini et al. 2013, Science 340:464, laser-heated DAC to 200 GPa, fast XRD, extrapolated", 6230, 500), ("Sinmyo et al. 2019, EPSL 510:45, resistance-heated DAC to 290 GPa, extrapolated (5770 +/- 280 with a Simon-Glatzel fit)", 5500, 220), ("Li et al. 2020, GRL 47:e2020GL087758, new shock-temperature measurements to 256 GPa, reduced to 330 GPa", 5950, 400), ] # The sources' own alloy-corrected estimates, for comparison with the pure-iron adiabat table. ALLOY_CORRECTED = [ "Sinmyo et al. 2019: upper bounds ICB 5120 +/- 390 K, CMB 3760 +/- 290 K (their adiabat also uses gamma = 1.5)", "Anzellini et al. 2013: ~700 K depression (from Alfe's group); CMB 4050 +/- 500 K in the supplement, implying ~5500 K at the ICB", ] RHO_CENTRE, RHO_ICB_SOLID = 13.0885, 12.7636 # PREM g/cm^3, inner-core adiabat PUBLISHED_DENSITIES = { # g/cm^3, PREM Table I evaluated at the boundaries "ICB_solid": 12.7636, "ICB_liquid": 12.1663, "CMB_core": 9.9035, "CMB_mantle": 5.5665} GAMMA_DEFAULT = 1.5 # outer-core Gruneisen parameter; 1.3-1.6 spans the literature def density(r_km, layers=LAYERS): x = r_km / A_KM for r_out, (c0, c1, c2, c3) in layers: if r_km <= r_out: return c0 + c1 * x + c2 * x ** 2 + c3 * x ** 3 return 0.0 def profile(layers=LAYERS, n=200000): """Return radii (m), density (kg/m^3), enclosed mass (kg), pressure (Pa).""" r = [A_KM * KM * i / n for i in range(n + 1)] rho = [density(ri / KM, layers) * 1000.0 for ri in r] m = [0.0] * (n + 1) for i in range(1, n + 1): # mass outward, trapezoid f0 = 4 * math.pi * r[i - 1] ** 2 * rho[i - 1] f1 = 4 * math.pi * r[i] ** 2 * rho[i] m[i] = m[i - 1] + 0.5 * (f0 + f1) * (r[i] - r[i - 1]) g = [0.0] + [G * m[i] / r[i] ** 2 for i in range(1, n + 1)] p = [0.0] * (n + 1) for i in range(n - 1, -1, -1): # pressure inward, trapezoid p[i] = p[i + 1] + 0.5 * (rho[i] * g[i] + rho[i + 1] * g[i + 1]) * (r[i + 1] - r[i]) return r, rho, m, p def at(r_km, r, arr): i = min(range(len(r)), key=lambda k: abs(r[k] - r_km * KM)) return arr[i] def adiabat(t_icb, rho_icb, rho_cmb, gamma): """Adiabatic temperature ratio for a convecting liquid: T ~ rho^gamma.""" return t_icb * (rho_cmb / rho_icb) ** gamma def run(verbose=True): r, rho, m, p = profile() mass = m[-1] p_icb = at(R_ICB_KM, r, p) / 1e9 p_cmb = at(R_CMB_KM, r, p) / 1e9 # densities just on the liquid side of each boundary rho_icb = density(R_ICB_KM + 1e-6) * 1000 rho_cmb = density(R_CMB_KM - 1e-6) * 1000 out = {"mass_kg": mass, "P_ICB_GPa": p_icb, "P_CMB_GPa": p_cmb, "rho_icb": rho_icb, "rho_cmb": rho_cmb} if verbose: print("Step 2: pressure at the freezing surface, from the PREM density profile alone") print(f" Earth mass {mass:.4e} kg (PREM publishes {PUBLISHED['mass_kg']:.3e})") print(f" P at core-mantle bdy {p_cmb:7.1f} GPa (PREM publishes {PUBLISHED['P_CMB_GPa']})") print(f" P at inner-core bdy {p_icb:7.1f} GPa (PREM publishes {PUBLISHED['P_ICB_GPa']})") print(f" liquid density: ICB side {rho_icb:.0f} kg/m^3, CMB side {rho_cmb:.0f} kg/m^3") print() print("Step 3 inputs: PURE-IRON melting point at ~330 GPa, five determinations (K)") for name, t, u in ICB_ESTIMATES: print(f" {t:5d} +/- {u:3d} {name}") print() print("Step 4: carry each pure-iron value out to the core-mantle boundary along an adiabat.") print(" These are UNCORRECTED for light elements, so each is an upper bound on the real temperature.") print(f" T_CMB = T_ICB * (rho_CMB/rho_ICB)^gamma; ratio at gamma=1.3/1.5/1.6 = " f"{(rho_cmb/rho_icb)**1.3:.3f}/{(rho_cmb/rho_icb)**1.5:.3f}/{(rho_cmb/rho_icb)**1.6:.3f}") print(f" {'ICB (K)':>8} {'CMB g=1.3':>10} {'CMB g=1.5':>10} {'CMB g=1.6':>10}") for name, t, u in ICB_ESTIMATES: print(f" {t:8d} {adiabat(t, rho_icb, rho_cmb, 1.3):10.0f} " f"{adiabat(t, rho_icb, rho_cmb, 1.5):10.0f} {adiabat(t, rho_icb, rho_cmb, 1.6):10.0f}") print() for s in ALLOY_CORRECTED: print(" For comparison, the sources' own alloy-corrected estimates:", s) print() r_in = (RHO_CENTRE / RHO_ICB_SOLID) ** GAMMA_DEFAULT print("The clip said 'the inner core', not its boundary. Along the inner core's own adiabat the centre") print(f" runs hotter than the boundary by a factor {r_in:.4f} (gamma {GAMMA_DEFAULT}): +{5500*(r_in-1):.0f} K on 5500 K, " f"+{6000*(r_in-1):.0f} K on 6000 K.") print() print("Scale check: deepest hole, Kola SG-3, 12.262 km of 6371 km radius = " f"{12.262/A_KM*100:.2f} % of the way down; the core-mantle boundary is at " f"{A_KM-R_CMB_KM:.0f} km, the inner-core boundary at {A_KM-R_ICB_KM:.1f} km.") return out def selftest(): ok = True o = run(verbose=False) checks = [ ("mass within 0.1 % of PREM's 5.974e24 kg", abs(o["mass_kg"] / PUBLISHED["mass_kg"] - 1) < 0.001), ("P_CMB within 0.1 % of 135.8 GPa", abs(o["P_CMB_GPa"] / PUBLISHED["P_CMB_GPa"] - 1) < 0.001), ("P_ICB within 0.1 % of 328.9 GPa", abs(o["P_ICB_GPa"] / PUBLISHED["P_ICB_GPa"] - 1) < 0.001), ] # Layer-local checks: density on each side of each boundary, against PREM Table I values. for name, r_km, side, key in [("inner-core side of ICB", R_ICB_KM, -1e-6, "ICB_solid"), ("outer-core side of ICB", R_ICB_KM, +1e-6, "ICB_liquid"), ("core side of CMB", R_CMB_KM, -1e-6, "CMB_core"), ("mantle side of CMB", R_CMB_KM, +1e-6, "CMB_mantle")]: got = density(r_km + side) checks.append((f"density at {name} = {got:.4f} vs PREM {PUBLISHED_DENSITIES[key]} (0.001 %)", abs(got / PUBLISHED_DENSITIES[key] - 1) < 0.00001)) # Negative control: a uniform-density Earth of the same mean density must FAIL the # ICB check. If it passed, the density profile would not be load-bearing. uniform = [(A_KM, (5.515, 0.0, 0.0, 0.0))] r, rho, m, p = profile(uniform) p_icb_uniform = at(R_ICB_KM, r, p) / 1e9 checks.append(("uniform-density control FAILS the ICB check (profile is load-bearing)", abs(p_icb_uniform / PUBLISHED["P_ICB_GPa"] - 1) > 0.20)) for name, passed in checks: print(("PASS " if passed else "FAIL ") + name) ok &= passed print(f"(uniform-density control gave P_ICB = {p_icb_uniform:.1f} GPa)") return ok if __name__ == "__main__": if "--selftest" in sys.argv: sys.exit(0 if selftest() else 1) run()