Helioseismology of the Shell: The Campaign for the Acoustic Peaks
1. The Enemy: What the Peaks Are
The pattern of warm and cool spots on the microwave sky prefers certain sizes, in a harmonic series — ℓ = 220, 537, 810, a drumbeat — and this series’ static cosmology has owed an account of it since the Four Calculations note declared it the standing obstacle. This note chronicles the campaign waged against that obstacle: the one mechanism class ever known to produce spectral peaks from incoherent driving (the resonant cavity, with the Sun as existence proof) was formulated, computed, and honestly defeated at the measured scales; the wall of the photosphere was then derived from atomic physics; a proposed engine (the Cepheid valve) was tested and found quenched; the decisive over-constrained fit was mounted and lost with cause stated; and the final reconnaissance uncovered something larger than the battle — an apparent 10⁴ inconsistency in the cosmology’s deep interior, which the framework’s own principles resolved into a foundational commitment. The campaign ends with the peaks still unexplained, three genuine discoveries banked, one master calculation named, and every step reproducible from the scripts beside this note. Defeats in public: that is the standard, and this note is its fullest exercise.
Strip the sky to its 2.725-kelvin glow and measure the mottling — warm and cool patches at one part in a hundred thousand. The patches are random; the statistics of their sizes are not. Decompose the sky into angular tones (the multipole ℓ, roughly 180° divided by patch size) and the power against ℓ is a drumbeat: peaks at ℓ ≈ 220, 537, 810, evenly spaced near Δℓ ≈ 300, odd peaks slightly enhanced, all fading under a smooth damping beyond ℓ ~ 1000, with a polarization pattern half a beat out of phase. The standard model explains this with a synchronized start: all sound waves in the young plasma released in phase at t ≈ 0, snapshot taken at recombination — a chord frozen mid-song. Models without the synchronized start — continuous, incoherent sources — were computed in 1997 to give one broad hump, and the measured second peak eliminated that entire field within three years. A static, eternal model stands prima facie in the same dock; this series’ own projection calculation (Four Calculations, §4) confirmed the naive expectation: hump, no chord.
2The One Open Door: Resonance
One physical system takes fully incoherent driving and produces needle-sharp spectral peaks anyway: the Sun. Random convection, no synchronization — yet its oscillation spectrum shows thousands of discrete peaks, because the Sun is a cavity, and resonance selects frequencies regardless of phase. Random hammering on a bell still sounds the bell’s note. The shell of this cosmology — stratified by the Tolman gradient, dense below, opaque at the bottom, thinning to transparency above — is structurally a stellar envelope. The campaign’s founding question: does the bell’s note become the drum’s spots? Stated risk, stated first: solar peaks are discrete in frequency, sustained in time; the sky’s snapshot needs discreteness in angular scale, frozen in space.
3The Campaign, Engagement by Engagement
The wall derived,
4The Discovery Beneath the Battlefield
Walking the ground before the assault, the recon checked whether the metric can carry its own bath — and found that it cannot, under standard bookkeeping: the equilibrium bath’s naive energy density grows as (1+z)⁴ against the metric’s derived source at ~(1+z)², crossing near z ≈ 150 and exceeding it by 10⁴ at the wall. A crisis larger than the peaks — resolved not by tuning but by the corpus’s own tension bridge: the bath is the thermodynamic face of the tension, one field read twice, and counting its aT⁴ as a second gravitating fluid is double-entry bookkeeping. The bath is booked once — the single-booking commitment, now installed in the Allgemeine Feldtheorie with its price stated (a departure from standard semiclassical bookkeeping; untested, not contradicted; forced by the postulates). Equivalently: gravity — the pre-tension — ends at the horizon; the mollusk lives strictly inside its sphere.
5The Standing of the Front
The peaks remain unexplained: that is the plain sentence, and it stays in the flagship’s caveats at full width. What the campaign banked: the photosphere derived from atomic physics; the wall-thickness/first-peak coincidence; the lapse-as-floor-mirror; the quarter-offset fingerprint matching the boundary anatomy; the closure of the resonance route with stated reasons; the quenching of the Cepheid engine; the dilution no-go; and the single-booking commitment, which redefines every deep-interior computation. All fronts now converge on one weapon that does not yet exist: the perturbation theory of the tension medium — with its zeroth-order principle fixed (bath fluctuations are not a second gravitating fluid), its battlefield mapped, and three measured numbers (ℓ = 220 with quarter offset; damping at ℓ ≈ 1400; the polarization phase) waiting as its judges. Scripts: cavity.py, solver.py, step3.py, step3b.py, saha.py, battle.py, beside this note.
6Why This Note Exists
Because a theory is its record. The standard model earned its authority partly through the defeats its rivals suffered in public — and a challenger earns the right to its victories only by keeping its losses in the same ledger. This campaign lost its stated objective and found, in losing, a derivation, a fingerprint, a quenched engine, a no-go, and a foundational commitment. The drum still belongs to the enemy. The war does not.
References
P. J. E. Peebles and J. T. Yu (1970); R. A. Sunyaev and Ya. B. Zeldovich (1970) — the acoustic-peak prediction; U.-L. Pen, U. Seljak and N. Turok, Phys. Rev. Lett. 79, 1611 (1997) — the incoherent-source no-go; Planck Collaboration results; R. B. Leighton et al. (1962) and solar p-mode literature — helioseismology; M. Saha (1920); A. S. Eddington, The Internal Constitution of the Stars (1926) — the valve; and the documents of this series (Four Calculations; the Metrics of the Living Spaces; the Allgemeine Feldtheorie, caveats (ii) and (xi); war diary and scripts, Cosmology folder). (Citations from memory; the literature-verification pass applies.)
7Verification
The companion scripts, with their recorded output. Each script's docstring states what it establishes and what it does not; the Source tab shows the file itself, unedited.
battle.py — battle
=== RECON: can the metric carry its own bath? ===
z= 10: bath rho_gamma=6.80e-27 metric-D source=5.84e-26 ratio bath/source=1.2e-01
z= 100: bath rho_gamma=4.83e-23 metric-D source=1.37e-24 ratio bath/source=3.5e+01
z= 300: bath rho_gamma=3.81e-21 metric-D source=7.97e-24 ratio bath/source=4.8e+02
z= 1100: bath rho_gamma=6.82e-19 metric-D source=7.08e-23 ratio bath/source=9.6e+03
-> beyond z~150 the equilibrium bath OUTWEIGHS the metric's entire source budget;
at the wall by 1e4. Consistency requires the tension medium's active mass to CANCEL
the bath's gravity to one part in 1e4 — or the bath must not gravitate conventionally.
=== THE ASSAULT: one dial (local n_b), three observables ===
n_b0 [1/m^3] z_vis ell_1 ell_D R_load rho_b/rho_LCDMrec rho_b/rho_metric
2.0e+03 1939 111 1 0.00 0.00 1.5e-01
1.0e+06 1973 113 inf 0.00 0.01 7.8e+01
1.0e+08 1973 120 inf 0.06 0.94 7.8e+03
5.0e+08 1973 150 inf 0.29 4.72 3.9e+04
1.5e+09 1973 225 inf 0.86 14.16 1.2e+05
3.0e+09 1973 337 inf 1.72 28.32 2.3e+05
1.0e+10 1973 859 inf 5.72 94.41 7.8e+05
CLOSEST FIT: n_b0 ~ 1.5e+09: ell_1=225, z_vis=1973, ell_D=inf
COST: local baryon density = 8.1e-18 kg/m^3 = 14.2x LCDM's recombination
and 1e+05x the metric's own source budget at the wall.
import numpy as np
kB=1.381e-23; me=9.109e-31; h=6.626e-34; eV=1.602e-19; chi=13.6*eV
c=2.998e8; H=2.27e-18; Mpc=3.086e22; sigT=6.652e-29; a_rad=7.566e-16
G=6.674e-11; mp=1.673e-27; rho_c=8.6e-27
kk=H/c*Mpc; r_sh=4280*np.log(1101.)
print("=== RECON: can the metric carry its own bath? ===")
for z in [10,100,300,1100]:
rho_gam=a_rad*(2.725*(1+z))**4/c**2
# metric D source (corpus profile): rho = rho_c[(e^{2kr}+6kr-3k^2r^2-1)/(3k^2r^2)] ~ leading (1+z)^2/(3 ln^2)
lnz=np.log(1+z)
rho_src=rho_c*((1+z)**2+6*lnz-3*lnz**2-1)/(3*lnz**2)
print(f" z={z:5d}: bath rho_gamma={rho_gam:.2e} metric-D source={rho_src:.2e} ratio bath/source={rho_gam/rho_src:.1e}")
print(" -> beyond z~150 the equilibrium bath OUTWEIGHS the metric's entire source budget;")
print(" at the wall by 1e4. Consistency requires the tension medium's active mass to CANCEL")
print(" the bath's gravity to one part in 1e4 — or the bath must not gravitate conventionally.")
print()
print("=== THE ASSAULT: one dial (local n_b), three observables ===")
def saha_profile(nb0):
x=np.linspace(-2500,2500,3000) # Mpc around T=3000K
T=3000*np.exp(kk*x); nb=nb0*np.exp(2*kk*x)
A=(2*np.pi*me*kB*T/h**2)**1.5*np.exp(-chi/(kB*T))/nb
xe=(-A+np.sqrt(A**2+4*A))/2
ne=xe*nb
tau=np.cumsum((ne*sigT*np.gradient(x)*Mpc)[::-1])[::-1]*(-1) # integrate from outside in
tau=np.cumsum(ne[::-1]*sigT*np.abs(np.gradient(x))[::-1]*Mpc)[::-1]
return x,T,xe,ne,tau
def observables(nb0):
x,T,xe,ne,tau=saha_profile(nb0)
i1=np.argmin(np.abs(tau-1.0))
z_vis=T[i1]/2.725*1.0-1 # T=T0(1+z)
z_vis=T[i1]/2.725-1
r_vis=4280*np.log(1+z_vis)
rho_b=ne[i1]/max(xe[i1],1e-9)*mp # total baryons there
rho_g=a_rad*T[i1]**4/c**2
R=3*rho_b/(4*rho_g); cs=c/np.sqrt(3*(1+R))
rho_fl=rho_g+rho_b
lamJ=cs*np.sqrt(np.pi/(G*rho_fl))/Mpc # local proper Mpc
ell1=np.pi*r_vis/(lamJ*(1+z_vis))
mfp=1/(ne[i1]*sigT)/Mpc
# visibility thickness: between tau=0.1 and tau=3
ia=np.argmin(np.abs(tau-0.1)); ib=np.argmin(np.abs(tau-3.0))
Dvis=abs(x[ia]-x[ib])
lamD=np.sqrt(mfp*Dvis/3)
ellD=np.pi*r_vis/(lamD*(1+z_vis))
return z_vis,ell1,ellD,R,rho_b,mfp,Dvis
print(f"{'n_b0 [1/m^3]':>12s} {'z_vis':>6s} {'ell_1':>7s} {'ell_D':>7s} {'R_load':>7s} {'rho_b/rho_LCDMrec':>17s} {'rho_b/rho_metric':>16s}")
rho_lcdm_rec=5.7e-19; rho_metric=6.9e-23
best=None
for nb0 in [2e3,1e6,1e8,5e8,1.5e9,3e9,1e10]:
z,l1,lD,R,rb,mfp,Dv=observables(nb0)
print(f"{nb0:12.1e} {z:6.0f} {l1:7.0f} {lD:7.0f} {R:7.2f} {rb/rho_lcdm_rec:17.2f} {rb/rho_metric:16.1e}")
if best is None or abs(l1-220)<abs(best[1]-220): best=(nb0,l1,z,lD,rb)
print()
nb=best[0]
print(f"CLOSEST FIT: n_b0 ~ {nb:.1e}: ell_1={best[1]:.0f}, z_vis={best[2]:.0f}, ell_D={best[3]:.0f}")
print(f"COST: local baryon density = {best[4]:.1e} kg/m^3 = {best[4]/rho_lcdm_rec:.1f}x LCDM's recombination")
print(f" and {best[4]/rho_metric:.0e}x the metric's own source budget at the wall.")
cavity.py — cavity
Hubble radius D_H = 4280 Mpc; shell at r = 29975 Mpc (z=1100) rho_gamma=6.82e-19, rho_b=3.44e-24 kg/m3 -> baryon loading R_b=3.78e-06 sound speed c_s = 0.577 c (pure radiation fluid; NO baryon loading -> honest flag: no odd/even asymmetry available) pressure scale height H_p = c_s^2/(cH) = 1427 Mpc = D_H/3 acoustic cutoff length 2H_p = 2853 Mpc Thomson wall thickness (0.1 e-fold) = 428 Mpc <- the visibility width corpus mapping: transverse 428 Mpc <-> l=220; so l = 220*(428/L) === candidate cavities and their fundamental scales === visibility wall D= 428 Mpc L=2D: L= 856 Mpc -> l = 110 visibility wall D= 428 Mpc L=D: L= 428 Mpc -> l = 220 pressure scale height D= 1427 Mpc L=2D: L= 2853 Mpc -> l = 33 pressure scale height D= 1427 Mpc L=D: L= 1427 Mpc -> l = 66 2H_p (cutoff length) D= 2853 Mpc L=2D: L= 5706 Mpc -> l = 17 2H_p (cutoff length) D= 2853 Mpc L=D: L= 2853 Mpc -> l = 33 === WKB trapping map: where do modes live? === omega_ac = 1.966e-18 /s (period 101464.0 Myr); N = 1.966e-18 /s equal-time P(l) from the naive sealed slab: monotonic smooth? max|d log P| step = 1.28e-02 (no peaks) === the honest finding and the one remaining lever === Flat sealed slab + steady incoherent driving -> P(l) smooth (risk confirmed analytically). Lever: the VISIBILITY WEIGHTING. We do not see the whole cavity; we see through a wall of thickness ~428 Mpc. A mode of vertical wavelength ~ 2D/n has surface signature weighted by integral of its eigenfunction across the visibility layer: modes with n*wall/2D ~ integer wash out. That filter DOES depend on the ratio (vertical wavelength / wall thickness) -> imposes structure in the n-sum; whether it transfers structure to kh depends on the w-integration. Step 2 computes it.
import numpy as np
# === The shell cavity of metric D, from the corpus's own numbers ===
c=2.998e8; H=2.27e-18 # s^-1 (70 km/s/Mpc)
Mpc=3.086e22
DH=c/H/Mpc # Hubble radius in Mpc
k=1/DH # e-fold per Mpc
zsh=1100.0; rsh=DH*np.log(1+zsh) # shell radius (proper), Mpc
print(f"Hubble radius D_H = {DH:.0f} Mpc; shell at r = {rsh:.0f} Mpc (z=1100)")
# stratification near the shell
# temperature T = T0 e^{kr}; local T=3000K; scale height of T: 1/k = DH (huge)
# gravity felt by static fluid: g = cH (uniform!)
g=c*H
# sound speed: photon-dominated fluid (baryons negligible in metric D at shell: check)
rho_crit=8.6e-27
rho_shell=8e3*rho_crit # tension medium at z=1100 (corpus)
rho_b=0.05*rho_shell # baryons
arad=7.566e-16; Tloc=3000.0
rho_g=arad*Tloc**4/c**2
Rb=3*rho_b/(4*rho_g)
cs=c/np.sqrt(3*(1+Rb))
print(f"rho_gamma={rho_g:.2e}, rho_b={rho_b:.2e} kg/m3 -> baryon loading R_b={Rb:.2e}")
print(f"sound speed c_s = {cs/c:.3f} c (pure radiation fluid; NO baryon loading -> honest flag: no odd/even asymmetry available)")
# pressure scale height of the radiation fluid held static in g=cH:
# H_p = c_s^2/g (isothermal-atmosphere estimate)
Hp=cs**2/g/Mpc/1.0
print(f"pressure scale height H_p = c_s^2/(cH) = {Hp:.0f} Mpc = D_H/3")
# acoustic cutoff frequency (isothermal): w_ac = c_s/(2 H_p); as a LENGTH: c_s/w_ac = 2H_p
print(f"acoustic cutoff length 2H_p = {2*Hp:.0f} Mpc")
# opacity wall: tau=2 within 0.1 e-folds (corpus) -> wall thickness
Dwall=0.1/k
print(f"Thomson wall thickness (0.1 e-fold) = {Dwall:.0f} Mpc <- the visibility width")
# corpus mapping: l=220 <-> 428 Mpc transverse at shell
L220=428.0
print(f"corpus mapping: transverse 428 Mpc <-> l=220; so l = 220*(428/L)")
print()
print("=== candidate cavities and their fundamental scales ===")
for name,D in [("visibility wall",Dwall),("pressure scale height",Hp),("2H_p (cutoff length)",2*Hp)]:
# fundamental standing wave: L1 ~ 2D; overtone spacing dl ~ l1
for mode,L in [("L=2D",2*D),("L=D",D)]:
l=220*428/L
print(f" {name:22s} D={D:6.0f} Mpc {mode}: L={L:6.0f} Mpc -> l = {l:6.0f}")
print()
print("=== WKB trapping map: where do modes live? ===")
# isothermal slab, vertical wavenumber: kv^2 = (w^2 - wac^2)/cs^2 - kh^2*(1 - N^2/w^2)
# radiation fluid, adiabatic index 4/3; N (buoyancy) in isothermal atm: N^2 = g/Hp*(1-1/Gamma)...
Gam=4/3
wac=cs/(2*Hp*Mpc)
N2=(g/(Hp*Mpc))*(Gam-1)/Gam
print(f"omega_ac = {wac:.3e} /s (period {2*np.pi/wac/3.15e13:.1f} Myr); N = {np.sqrt(N2):.3e} /s")
# for each kh, the mode ladder in a cavity of depth D sealed below (velocity node) and cutoff above:
D=Dwall*Mpc
kh_grid=np.linspace(1e-3,30,600)/ (428*Mpc/ (2*np.pi)) # kh in rad/m, spanning l ~ few..~4000
ell=lambda kh: 220*428/(2*np.pi/kh/Mpc)
nmax=[]
for kh in kh_grid:
# trapped acoustic branch: w^2 > wac^2 + cs^2 kh^2 (approx, ignoring N-branch)
# vertical quantization: kv_n = n pi / D -> w_n^2 = wac^2 + cs^2(kh^2 + (n pi/D)^2)
nmax.append(5) # ladder exists for all kh: no kh-window from trapping in this slab
# equal-time spectrum: white noise driving, damping gamma; P(kh) ~ sum_n 1/w_n^2 (velocity^2 ~ T fluctuation power)
Pk=[]
for kh in kh_grid:
wn2=[wac**2+cs**2*(kh**2+(n*np.pi/D)**2) for n in range(1,40)]
Pk.append(sum(1/w2 for w2 in wn2))
Pk=np.array(Pk); Pk/=Pk.max()
# check for structure
d=np.diff(np.log(Pk))
print(f"equal-time P(l) from the naive sealed slab: monotonic smooth? max|d log P| step = {np.max(np.abs(np.diff(d))):.2e} (no peaks)")
print()
print("=== the honest finding and the one remaining lever ===")
print("Flat sealed slab + steady incoherent driving -> P(l) smooth (risk confirmed analytically).")
print("Lever: the VISIBILITY WEIGHTING. We do not see the whole cavity; we see through a wall of")
print("thickness ~428 Mpc. A mode of vertical wavelength ~ 2D/n has surface signature weighted by")
print("integral of its eigenfunction across the visibility layer: modes with n*wall/2D ~ integer wash out.")
print("That filter DOES depend on the ratio (vertical wavelength / wall thickness) -> imposes structure")
print("in the n-sum; whether it transfers structure to kh depends on the w-integration. Step 2 computes it.")
saha.py — saha
Saha wall: x_e=0.5 at T=2902 K (x=-142 Mpc); 10%-90% width = 519 Mpc photon mean free path at full ionization: 204 Mpc; at x_e=0.5: 521 Mpc tau=1 at x=+88 Mpc (T=3062 K); tau=2 at x=+317 Mpc -> visibility depth ~1086 Mpc energy reservoirs at the wall: ionization 4.1e-15 J/m^3 vs radiation 5.37e-02 J/m^3 ratio = 7.6e-14 -> Gamma_1 dip ~ 1e-13: the gas cannot bend the fluid's spring and the deeper quench: kappa-mechanism throttles a background FLUX; Tolman equilibrium has net flux ~ 0 (the bath is maintained, not powered). VERDICT: the Cepheid engine is QUENCHED — the wall breathes across the door (Saha two-way traffic) but coherent self-excitation has no free-energy river to tap. Driving must be stirring/noise; the duct only SELECTS. diffusion scale sqrt(mfp*D/3) = 306 Mpc -> damping onset ell_D ~ 307 measured damping tail onset: ~1000-1400. ratio ell_D/ell_1 predicted 1.4, measured ~4.5-6
import numpy as np
# === PART 1: the wall's true anatomy from Saha ===
kB=1.381e-23; me=9.109e-31; h=6.626e-34; eV=1.602e-19; chi=13.6*eV
c=2.998e8; H=2.27e-18; Mpc=3.086e22; sigT=6.652e-29; a_rad=7.566e-16
kk=H/c*Mpc # e-fold per Mpc
nb0=2.0e3 # baryons per m^3 at the 3000K layer (corpus)
x=np.linspace(-3000,3000,4000) # Mpc, 0 at T=3000K, increasing inward
T=3000.0*np.exp(kk*x)
nb=nb0*np.exp(2*kk*x) # mild density gradient
S=(2*np.pi*me*kB*T/h**2)**1.5*np.exp(-chi/(kB*T))
# Saha: xe^2/(1-xe) = S/nb
A=S/nb
xe=(-A+np.sqrt(A**2+4*A*nb/nb))/2 # solve xe^2 + A xe - A =0 per unit: xe=(-A+sqrt(A^2+4A))/2
xe=(-A+np.sqrt(A**2+4*A))/2
ne=xe*nb
mfp=1/(ne*sigT)/Mpc # photon mean free path in Mpc
# wall location & width
i50=np.argmin(np.abs(xe-0.5)); i10=np.argmin(np.abs(xe-0.1)); i90=np.argmin(np.abs(xe-0.9))
print(f"Saha wall: x_e=0.5 at T={T[i50]:.0f} K (x={x[i50]:+.0f} Mpc); 10%-90% width = {abs(x[i90]-x[i10]):.0f} Mpc")
print(f"photon mean free path at full ionization: {mfp[i90+200]:.0f} Mpc; at x_e=0.5: {mfp[i50]:.0f} Mpc")
# optical depth inward from transparency
tau=np.cumsum(ne*sigT*np.gradient(x)*Mpc)
tau-=tau[np.argmin(np.abs(xe-0.01))]
itau1=np.argmin(np.abs(tau-1)); itau2=np.argmin(np.abs(tau-2))
print(f"tau=1 at x={x[itau1]:+.0f} Mpc (T={T[itau1]:.0f} K); tau=2 at x={x[itau2]:+.0f} Mpc -> visibility depth ~{x[itau2]-x[np.argmin(np.abs(xe-0.01))]:.0f} Mpc")
print()
# === PART 2: is the kappa-engine alive? ===
rho_g=a_rad*T**4/c**2; rho_b=nb*1.673e-27
ion_res=nb*chi # ionization energy reservoir J/m^3
rad_res=a_rad*T**4
i0=i50
print(f"energy reservoirs at the wall: ionization {ion_res[i0]:.1e} J/m^3 vs radiation {rad_res[i0]:.2e} J/m^3")
print(f"ratio = {ion_res[i0]/rad_res[i0]:.1e} -> Gamma_1 dip ~ 1e-13: the gas cannot bend the fluid's spring")
print(f"and the deeper quench: kappa-mechanism throttles a background FLUX; Tolman equilibrium has")
print(f"net flux ~ 0 (the bath is maintained, not powered). VERDICT: the Cepheid engine is QUENCHED —")
print(f"the wall breathes across the door (Saha two-way traffic) but coherent self-excitation has no")
print(f"free-energy river to tap. Driving must be stirring/noise; the duct only SELECTS.")
print()
# === PART 3: damping tail consistency ===
D_vis=abs(x[itau2]-x[np.argmin(np.abs(xe-0.01))])
lam_mfp=mfp[itau1]
lam_D=np.sqrt(lam_mfp*D_vis/3)
ellD=220*428/lam_D
print(f"diffusion scale sqrt(mfp*D/3) = {lam_D:.0f} Mpc -> damping onset ell_D ~ {ellD:.0f}")
print(f"measured damping tail onset: ~1000-1400. ratio ell_D/ell_1 predicted {ellD/220:.1f}, measured ~4.5-6")
np.save("wall.npy",np.vstack([x,T,xe,ne,mfp]))
solver.py — solver
cutoff at wall: 1.97e-18/s; at 4000 Mpc depth: 5.01e-18/s (trapping window x2.5) P(l) computed. structure diagnostics: white-driving max at l=40; local maxima at: [1155 1201 1262 1323 1354] flow-driving max at l=40; local maxima: [ 987 1017 1063 1155 1201 1262 1323 1354]
import numpy as np
# ==== Vertical wave solver for the metric-D shell cavity ====
# depth coordinate x [Mpc], increasing INWARD (deeper, hotter); x=0 at the opacity wall center
c=2.998e8; H=2.27e-18; Mpc=3.086e22
cs=c/np.sqrt(3.0)
k=H/c*Mpc # e-fold per Mpc (1/4280)
g=c*H # static gravity
w_wall=428.0 # opacity transition width [Mpc]
# profiles (coordinate/observer units):
def v_of(x): return cs*np.exp(-k*x) # coordinate sound speed: slower deeper (lapse)
def wac_of(x): return g/(2*v_of(x)) # acoustic cutoff: RISES with depth (the floor mirror)
def N2_of(x): return wac_of(x)**2 # buoyancy ~ cutoff (Gamma=4/3 isothermal-like)
def fluid_frac(x): # 1 deep (coupled fluid) -> 0 above wall (free streaming)
return 0.5*(1+np.tanh(x/(w_wall/2)))
# grid: from deep interior to above the wall
xg=np.linspace(-1500,6000,3000)[::-1] # integrate downward->? we'll go bottom->top: reorder
xg=np.linspace(6000,-1500,3000) # bottom (deep) -> top
dx=(xg[1]-xg[0])*Mpc # negative (moving up)
# horizontal scale mapping: corpus: transverse L Mpc <-> ell = 220*428/L; k_h = 2pi/L
ells=np.linspace(40,1400,90)
Ls=220*428/ells # Mpc
khs=2*np.pi/(Ls*Mpc)
# frequency band: around the trapped window [wac(top fluid), wac(deep)]
wac_top=wac_of(0.0); wac_deep=wac_of(4000.0)
oms=np.linspace(0.3*wac_top,1.8*wac_deep,240)
print(f"cutoff at wall: {wac_top:.2e}/s; at 4000 Mpc depth: {wac_deep:.2e}/s (trapping window x{wac_deep/wac_top:.1f})")
# transfer across the grid for all omegas at once, per ell
def response(kh):
om=oms
psi=np.ones_like(om,dtype=complex); dpsi=np.zeros_like(psi)
# bottom BC: unit up-going wave in deep region: psi=e^{iKx}: dpsi/dx=iK psi
x0=xg[0]; ff=fluid_frac(x0)
K2=( (om**2-wac_of(x0)**2)/v_of(x0)**2 - kh**2*(1-N2_of(x0)/om**2) )*ff + (1-ff)*(-(kh**2))
K=np.sqrt(np.abs(K2))*np.where(K2>0,1,1) # magnitude
dpsi=1j*np.sqrt(K2.astype(complex))*psi
O=np.zeros_like(om,dtype=complex)
for x in xg[1:]:
ff=fluid_frac(x)
K2=(( (om**2-wac_of(x)**2)/v_of(x)**2 - kh**2*(1-N2_of(x)/om**2) )*ff + (1-ff)*(-(kh**2))).astype(complex)
# RK2 step for psi''=-K2 psi
ddpsi=-K2*psi
psi_m=psi+dpsi*dx/2; dpsi_m=dpsi+ddpsi*dx/2
ddpsi_m=-K2*psi_m
psi=psi+dpsi_m*dx; dpsi=dpsi+ddpsi_m*dx
# visibility weighting: we observe within the wall (gaussian at x=0, width w_wall/2)
gvis=np.exp(-x**2/(2*(w_wall/2)**2))
O+= gvis*psi*abs(dx)/Mpc
# normalize occasionally to avoid overflow in evanescent zones
m=np.max(np.abs(psi));
if m>1e100: psi/=m; dpsi/=m; O/=m
# normalize response by field amplitude in fluid (per unit driving flux):
return np.abs(O)**2/np.maximum(np.abs(psi)**2,1e-300)
P_white=[]; P_flow=[]
Uflow=6e5 # 600 km/s FIRAS bound: flow-band driving centered w*=U/w_wall
w_star=Uflow/(w_wall*Mpc)
Sflow=np.exp(-(np.log(oms/max(w_star,1e-22)))**2/2) # broad log-band at flow timescale
for kh in khs:
R=response(kh)
P_white.append(np.trapz(R,oms))
P_flow.append(np.trapz(R*Sflow,oms))
P_white=np.array(P_white); P_flow=np.array(P_flow)
np.save('ells.npy',ells); np.save('Pw.npy',P_white); np.save('Pf.npy',P_flow)
import matplotlib
matplotlib.use("Agg"); import matplotlib.pyplot as plt
fig,ax=plt.subplots(figsize=(9,5.5),dpi=140)
ax.plot(ells,P_white/P_white.max(),'b-',lw=2,label='white incoherent driving (eternal noise)')
ax.plot(ells,P_flow/P_flow.max(),'g--',lw=2,label='flow-band driving (the wind, 600 km/s)')
for l0 in [220,537,810]: ax.axvline(l0,color='r',ls=':',alpha=0.7)
ax.set_xlabel('multipole ℓ'); ax.set_ylabel('P(ℓ) (normalized)')
ax.set_title('Shell-cavity equal-time spectrum vs measured peak positions (dotted red)')
ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout(); plt.savefig('Pl.png')
print("P(l) computed. structure diagnostics:")
lm=ells[np.argmax(P_white)]
print(f"white-driving max at l={lm:.0f}; local maxima at:", ells[1:-1][(P_white[1:-1]>P_white[2:])&(P_white[1:-1]>P_white[:-2])].astype(int))
print(f"flow-driving max at l={ells[np.argmax(P_flow)]:.0f}; local maxima:", ells[1:-1][(P_flow[1:-1]>P_flow[2:])&(P_flow[1:-1]>P_flow[:-2])].astype(int))
step3.py — step3
6 trapped modes found across 60 ell values quality factors: median Q = 4.6e+05, min 9.5e+04, max 2.5e+11 vertical orders n present: [np.int64(0), np.int64(1)]... white-driving local maxima at ell = [] flow-driving local maxima at ell = []
import numpy as np
c=2.998e8; H=2.27e-18; Mpc=3.086e22
cs=c/np.sqrt(3.0); k=H/c*Mpc; g=c*H; w_wall=428.0
x=np.linspace(-1200,6000,1400) # Mpc, ascending: above wall -> deep
dxm=(x[1]-x[0])*Mpc
v=cs*np.exp(-k*x); wac=g/(2*v); N2=wac**2
ff=0.5*(1+np.tanh(x/(w_wall/2))) # fluid fraction
gvis=np.exp(-x**2/(2*(w_wall/2)**2)) # visibility layer
ifl=np.argmax(ff>0.05) # fluid edge index
i0=np.argmin(np.abs(x)) # wall center
ells=np.linspace(40,1200,60); khs=ells/(14986*Mpc) # 1/m
oms=np.linspace(0.3*wac[i0],2.5*wac[-1],900)
w_star=6e5/(w_wall*Mpc)
def K2(om,kh): return ff*((om**2-wac**2)/v**2 - kh**2*(1-N2/om**2)) + (1-ff)*(-kh**2)
modes=[] # (ell, om_n, n, B_top, B_vis, Phi)
for iL,kh in enumerate(khs):
Phis=np.zeros(len(oms)); Bts=np.zeros(len(oms)); Bvs=np.zeros(len(oms)); ok=np.zeros(len(oms),bool)
for io,om in enumerate(oms):
q=K2(om,kh)
m=q>0
m[:ifl]=False
if not m.any(): continue
# connected run nearest the wall
idx=np.where(m)[0]
splits=np.where(np.diff(idx)>1)[0]
runs=np.split(idx,splits+1)
run=runs[0]
if run[-1]>=len(x)-2: continue # leaks into the deep: not trapped
Kv=np.sqrt(q[run])
Phis[io]=np.sum(Kv)*dxm
# barrier above: from fluid edge to run start
bar=slice(ifl,run[0])
Bts[io]=np.sum(np.sqrt(np.maximum(-q[bar],0)))*dxm
# tail to visibility center
barv=slice(min(i0,run[0]),run[0])
Bvs[io]=np.sum(np.sqrt(np.maximum(-q[barv],0)))*dxm
ok[io]=True
# quantization crossings Phi = (n+1/2)pi
for n in range(0,60):
tgt=(n+0.5)*np.pi
d=Phis-tgt
s=np.where(ok[:-1]&ok[1:]&(d[:-1]*d[1:]<0))[0]
for i in s:
f=d[i]/(d[i]-d[i+1]); om_n=oms[i]+(oms[i+1]-oms[i])*f
Bt=Bts[i]+(Bts[i+1]-Bts[i])*f; Bv=Bvs[i]+(Bvs[i+1]-Bvs[i])*f
modes.append((ells[iL],om_n,n,Bt,Bv,tgt))
modes=np.array(modes)
np.save("modes.npy",modes)
print(f"{len(modes)} trapped modes found across {len(ells)} ell values")
if len(modes):
Q=np.exp(2*modes[:,3])*(modes[:,5]/np.pi+0.5)
print(f"quality factors: median Q = {np.median(Q):.1e}, min {Q.min():.1e}, max {Q.max():.1e}")
print(f"vertical orders n present: {sorted(set(modes[:,2].astype(int)))[:12]}...")
# spectrum: P(ell) = sum_n W_n * S(om)/ (Gam*om^2); W ~ e^{-2Bv}; Gam = om e^{-2Bt}/(2(Phi/pi))
P=np.zeros(len(ells)); Pf=np.zeros(len(ells))
for (l,om,n,Bt,Bv,Phi) in modes:
iL=np.argmin(np.abs(ells-l))
Gam=om*np.exp(-2*Bt)/(2*(Phi/np.pi))
W=np.exp(-2*Bv)
Sf=np.exp(-(np.log(om/w_star))**2/2)
P[iL]+= W/( (Gam if Gam>0 else 1e-99) *om**2)
Pf[iL]+= W*Sf/((Gam if Gam>0 else 1e-99)*om**2)
np.save("Pl3.npy",np.vstack([ells,P,Pf]))
pk=ells[1:-1][(P[1:-1]>P[2:])&(P[1:-1]>P[:-2])]
print("white-driving local maxima at ell =",pk.astype(int))
pkf=ells[1:-1][(Pf[1:-1]>Pf[2:])&(Pf[1:-1]>Pf[:-2])]
print("flow-driving local maxima at ell =",pkf.astype(int))
step3b.py — step3b
mode census per ell (min/median/max): 0 0 4 flow target omega* = 4.54e-20/s vs cutoff at wall 1.97e-18/s (ratio 0.023) white P(l) local maxima: [] flow P(l) local maxima: []
import numpy as np
c=2.998e8; H=2.27e-18; Mpc=3.086e22
cs=c/np.sqrt(3.0); k=H/c*Mpc; g=c*H; w_wall=428.0
x=np.linspace(-1200,6000,1200)
dxm=(x[1]-x[0])*Mpc
v=cs*np.exp(-k*x); wac=g/(2*v); N2=wac**2
ff=0.5*(1+np.tanh(x/(w_wall/2)))
ifl=np.argmax(ff>0.05); i0=np.argmin(np.abs(x))
ells=np.linspace(40,1200,60); khs=ells/(14986*Mpc)
oms=np.linspace(0.25*wac[i0],2.6*wac[-1],1200)
w_star=6e5/(w_wall*Mpc); sig=0.15
def K2(om,kh): return ff*((om**2-wac**2)/v**2 - kh**2*(1-N2/om**2)) + (1-ff)*(-kh**2)
P=np.zeros(len(ells)); Pf=np.zeros(len(ells)); counts=np.zeros(len(ells),int)
ridge_l=[]; ridge_o=[]
for iL,kh in enumerate(khs):
Phis=np.full(len(oms),np.nan)
for io,om in enumerate(oms):
q=K2(om,kh); m=q>0; m[:ifl]=False
if not m.any(): continue
idx=np.where(m)[0]; splits=np.where(np.diff(idx)>1)[0]; run=np.split(idx,splits+1)[0]
if run[-1]>=len(x)-2: continue
Phis[io]=np.sum(np.sqrt(q[run]))*dxm
okm=~np.isnan(Phis)
if not okm.any(): continue
# crossings of half-integer multiples of pi
Nf=(Phis-np.pi/2)/np.pi
for io in range(len(oms)-1):
if not (okm[io] and okm[io+1]): continue
n1,n2=np.floor(Nf[io]),np.floor(Nf[io+1])
for n in np.arange(min(n1,n2)+1,max(n1,n2)+1):
f=(n-Nf[io])/(Nf[io+1]-Nf[io]); om_n=oms[io]+(oms[io+1]-oms[io])*f
counts[iL]+=1
E=1.0/om_n**2 # equal-time variance per mode, white noise, ~Q-independent overlap sum
P[iL]+=E
Sf=np.exp(-(np.log(om_n/w_star))**2/(2*sig**2))
Pf[iL]+=E*Sf
if counts[iL]<400: ridge_l.append(ells[iL]); ridge_o.append(om_n)
print("mode census per ell (min/median/max):",counts.min(),int(np.median(counts)),counts.max())
np.save("Pl3b.npy",np.vstack([ells,P,Pf]))
print(f"flow target omega* = {w_star:.2e}/s vs cutoff at wall {wac[i0]:.2e}/s (ratio {w_star/wac[i0]:.3f})")
def peaks(A): return ells[1:-1][(A[1:-1]>A[2:])&(A[1:-1]>A[:-2])].astype(int)
print("white P(l) local maxima:",peaks(P))
print("flow P(l) local maxima:",peaks(Pf))
import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(13,5),dpi=130)
axs[0].plot(ridge_l,np.array(ridge_o)/wac[i0],'k.',ms=1); axs[0].axhline(w_star/wac[i0],color='g',ls='--',label='the wind ω*')
axs[0].set_xlabel('ℓ'); axs[0].set_ylabel('ω / ω_ac(wall)'); axs[0].set_title('mode ridges ω_n(ℓ)'); axs[0].legend()
axs[1].plot(ells,P/P.max(),'b-',label='white driving')
axs[1].plot(ells,Pf/max(Pf.max(),1e-300),'g--',label='narrowband wind')
for l0 in [220,537,810]: axs[1].axvline(l0,color='r',ls=':')
axs[1].set_xlabel('ℓ'); axs[1].set_title('P(ℓ) vs measured peaks'); axs[1].legend()
plt.tight_layout(); plt.savefig('Pl3b.png')