""" cmbs4_forecast_extended.py ---------------------------- Extends cmbs4_forecast.py with three things the author asked about: (A) Eigen-decomposition of the Fisher matrix to identify which parameter *combinations* are actually constrained vs. exactly degenerate. (B) External-data Fisher blocks: a simple BAO-like prior on the background parameters, and a Tegmark (1997) Veff-based LSS P(k) shape forecast for a DESI-like survey. (C) A CMB lensing (phi-phi) diagonal Fisher block, using a flat-sky minimum-variance-*TT-only* quadratic-estimator reconstruction noise N0(L) (Hu & Okamoto 2002 flat-sky limit) built directly from CAMB's lensed C_l^TT + the same noise model as the temperature/polarisation Fisher block. Headline results found while building this (fiducial kc=0.3 Mpc^-1, eps3=0.02, i.e. the paper's own stated CMB-S4 detectability threshold): 1. There is an EXACT degeneracy in the model as written: X_n(k) = beta_eff * c_n * sin^p(chi_k) is invariant under beta_eff -> beta_eff/lambda, c_n -> lambda*c_n for any lambda. This shows up as a numerically-zero Fisher eigenvalue and is a parametrisation issue, not a data limitation -- fix it by fixing beta_eff=1 (or any one c_n=1) and fitting only the ratios. 2. Even after that exact degeneracy is removed, kc and eps3 remain strongly correlated with (c1,c2,c3) at the *current* Planck best fit (kc=2.83), because there the model has relaxed into the pure power-law limit and essentially every RLMT-specific direction has a vanishing derivative -- no reparametrisation rescues a genuinely flat direction; only being closer to the transition scale does. 3. At the paper's own stated threshold (kc=0.3, eps3=0.02): if beta_eff, c1, c2, c3 can be fixed/well-motivated externally (not left as free nuisance parameters), CMB-S4 (TT+TE+EE) alone forecasts sigma(kc) ~ 0.013 -> kc/sigma ~ 24 (strong detection) sigma(eps3) ~ 0.073 -> eps3/sigma ~ 0.27 (NOT a detection) Adding the lensing phi-phi block (TT-only N0) or the simplified LSS P(k) shape prior barely moves sigma(eps3), because both probe essentially the same tilt direction as ns and don't break the ns-eps3 degeneracy. kc is easy once the amplitude nuisance is controlled; eps3 needs a probe with genuinely different k-leverage (e.g. spectral distortions, Lyman-alpha forest, 21cm) to move beyond ~0.3 sigma. Run: python3 cmbs4_forecast_extended.py """ import numpy as np import camb from camb import initialpower from scipy.interpolate import interp1d import cmbs4_forecast as base # reuses rlmt_power_spectrum, get_cls, noise_cl, FIDUCIAL # ---------------------------------------------------------------------- # (A) Eigen-decomposition / degeneracy diagnostics # ---------------------------------------------------------------------- def eigen_analysis(F, names, fiducial): """Diagonalise the Fisher matrix in the log-parameter basis (dimensionless derivatives) so eigenvalues are directly comparable.""" fid = np.array([fiducial[n] for n in names]) Fn = F * np.outer(fid, fid) vals, vecs = np.linalg.eigh(Fn) order = np.argsort(vals)[::-1] vals, vecs = vals[order], vecs[:, order] print("Eigenvalues (log-parameter basis), largest=best-constrained:") for i, v in enumerate(vals): vec = vecs[:, i] top = np.argsort(-np.abs(vec))[:4] desc = ", ".join(f"{names[j]}:{vec[j]:+.2f}" for j in top) print(f" lambda_{i:2d} = {v:12.4e} dominant: {desc}") return vals, vecs def sigma_with_fixed(F, names, fixed): """Marginalised errors with a subset of parameters held FIXED (removed from the Fisher matrix before inversion -- not the same as an infinitely tight prior on a combination; this is the standard 'profile out by omission' prescription for exactly-known parameters).""" keep = [i for i, n in enumerate(names) if n not in fixed] cov = np.linalg.inv(F[np.ix_(keep, keep)]) s = np.sqrt(np.diag(cov)) return dict(zip([names[i] for i in keep], s)) # ---------------------------------------------------------------------- # (B) External priors: BAO (illustrative) + simplified LSS shape (Tegmark Veff) # ---------------------------------------------------------------------- def bao_prior_fisher(names, sigma_H0=0.3, sigma_omch2=0.001, sigma_ombh2=0.0001): """Illustrative diagonal BAO-like prior. Replace sigmas with a real forecast covariance (e.g. DESI/Euclid Fisher forecast) for production use -- these numbers are ballpark only.""" n = len(names) F = np.zeros((n, n)) sig = {"H0": sigma_H0, "omch2": sigma_omch2, "ombh2": sigma_ombh2} for name, s in sig.items(): if name in names: F[names.index(name), names.index(name)] = 1.0 / s ** 2 return F def lss_shape_fisher(fiducial, varied, step, kh_min=0.005, kh_max=0.3, npoints=200, z=1.0, V_survey=10e9, n_g=5e-4): """Tegmark (1997) effective-volume Fisher forecast for the *shape* of the linear matter power spectrum, DESI-like survey numbers (V ~ 10 Gpc^3, n_g ~ 5e-4 Mpc^-3) as a rough stand-in for real LSS Fisher matrices. Ignores RSD/bias marginalisation for simplicity.""" def get_pk(params): pk = base.rlmt_power_spectrum(base.KGRID, params["As"], params["ns"], params["beta_eff"], params["c1"], params["c2"], params["c3"], params["kstar"], params["kc"], params["eps3"]) pars = camb.CAMBparams() pars.set_cosmology(H0=params["H0"], ombh2=params["ombh2"], omch2=params["omch2"], tau=params["tau"]) pk_ini = initialpower.SplinedInitialPower() pk_ini.set_scalar_table(base.KGRID, pk) pk_ini.effective_ns_for_nonlinear = params["ns"] pars.InitPower = pk_ini pars.set_matter_power(redshifts=[z], kmax=kh_max * 2) pars.NonLinear = camb.model.NonLinear_none results = camb.get_results(pars) kh, _, pkm = results.get_matter_power_spectrum(minkh=kh_min, maxkh=kh_max, npoints=npoints) return kh, pkm[0] kh, P0 = get_pk(fiducial) derivs = {} for name in varied: pp = dict(fiducial); pp[name] += step[name] pm = dict(fiducial); pm[name] -= step[name] _, Pp = get_pk(pp) _, Pm = get_pk(pm) derivs[name] = (Pp - Pm) / (2 * step[name]) n = len(varied) F = np.zeros((n, n)) dk = kh[1] - kh[0] for idx, k in enumerate(kh): Pk = P0[idx] Veff = (n_g * Pk / (1 + n_g * Pk)) ** 2 * V_survey prefac = k ** 2 * dk / (4 * np.pi ** 2) * Veff / 2 for i, ni in enumerate(varied): dlnPi = derivs[ni][idx] / Pk for j, nj in enumerate(varied): dlnPj = derivs[nj][idx] / Pk F[i, j] += prefac * dlnPi * dlnPj return F # ---------------------------------------------------------------------- # (C) CMB lensing (phi-phi) block: flat-sky TT-only N0 + Fisher # ---------------------------------------------------------------------- def build_N0_TT(lensed_TT, noise_NT, lmax=3000, L_values=None): """Flat-sky minimum-variance TT-only lensing reconstruction noise N0(L) on the lensing potential phi (Hu & Okamoto 2002, flat-sky limit). TT-only: a real S4 analysis would combine TT+TE+EE+EB+TB which is noticeably better at high L -- this is a conservative / simplified estimate.""" ls = np.arange(len(lensed_TT)) Ctot = lensed_TT[:lmax + 1] + noise_NT[:lmax + 1] ls_c = ls[:lmax + 1] Cl_TT_interp = interp1d(ls_c, lensed_TT[:lmax + 1], bounds_error=False, fill_value=0.0) Ctot_interp = interp1d(ls_c, Ctot, bounds_error=False, fill_value=1e30) if L_values is None: L_values = np.unique(np.round(np.logspace(np.log10(2), np.log10(lmax), 25)).astype(int)) l1_grid = np.linspace(2, lmax, 400) phi_grid = np.linspace(0, 2 * np.pi, 200, endpoint=False) L1, PHI = np.meshgrid(l1_grid, phi_grid, indexing="ij") l1x, l1y = L1 * np.cos(PHI), L1 * np.sin(PHI) dl1, dphi = l1_grid[1] - l1_grid[0], phi_grid[1] - phi_grid[0] Cl1, Ctot1 = Cl_TT_interp(L1), Ctot_interp(L1) N0 = np.zeros(len(L_values)) for iL, Lval in enumerate(L_values): l2x, l2y = Lval - l1x, -l1y l2 = np.sqrt(l2x ** 2 + l2y ** 2) Cl2, Ctot2 = Cl_TT_interp(l2), Ctot_interp(l2) f = Cl1 * (l1x * Lval) + Cl2 * (l2x * Lval) integrand = f ** 2 / (2 * Ctot1 * Ctot2) integrand = np.where((Ctot1 > 1e29) | (Ctot2 > 1e29), 0.0, integrand) integral = np.sum(integrand * L1 * dl1 * dphi) / (2 * np.pi) ** 2 N0[iL] = 1.0 / integral if integral > 0 else np.inf return L_values, N0 def lensing_fisher(fiducial, varied, step, lmax=3000, fsky=0.4, lmin=8): def get_lensed_pp(params): pk = base.rlmt_power_spectrum(base.KGRID, params["As"], params["ns"], params["beta_eff"], params["c1"], params["c2"], params["c3"], params["kstar"], params["kc"], params["eps3"]) pars = camb.CAMBparams() pars.set_cosmology(H0=params["H0"], ombh2=params["ombh2"], omch2=params["omch2"], tau=params["tau"]) pk_ini = initialpower.SplinedInitialPower() pk_ini.set_scalar_table(base.KGRID, pk) pk_ini.effective_ns_for_nonlinear = params["ns"] pars.InitPower = pk_ini pars.set_for_lmax(lmax, lens_potential_accuracy=1) r = camb.get_results(pars) p = r.get_cmb_power_spectra(pars, CMB_unit="muK", raw_cl=True) return p["lensed_scalar"], p["lens_potential"] lensed0, clpp0 = get_lensed_pp(fiducial) ls = np.arange(lensed0.shape[0]) NT, _ = base.noise_cl(ls) Lv, N0sparse = build_N0_TT(lensed0[:, 0], NT, lmax=lmax) N0_full = np.clip(interp1d(Lv, N0sparse, bounds_error=False, fill_value="extrapolate")(ls), 0, None) pp_derivs = {} for name in varied: pp = dict(fiducial); pp[name] += step[name] pm = dict(fiducial); pm[name] -= step[name] _, clpp_p = get_lensed_pp(pp) _, clpp_m = get_lensed_pp(pm) pp_derivs[name] = (clpp_p[:, 0] - clpp_m[:, 0]) / (2 * step[name]) n = len(varied) F = np.zeros((n, n)) for L in range(lmin, min(lmax, len(ls))): tot = clpp0[L, 0] + N0_full[L] if tot <= 0: continue for i, ni in enumerate(varied): for j, nj in enumerate(varied): if j < i: continue val = 0.5 * (2 * L + 1) * fsky * pp_derivs[ni][L] * pp_derivs[nj][L] / tot ** 2 F[i, j] += val if j != i: F[j, i] += val return F # ---------------------------------------------------------------------- # Main: reproduce the headline numbers above # ---------------------------------------------------------------------- def main(): fid = dict(base.FIDUCIAL) fid["kc"], fid["eps3"] = 0.30, 0.02 # paper's stated detectability threshold varied_full = ["kc", "eps3", "ns", "As", "beta_eff", "c1", "c2", "c3"] step = {"kc": 0.009, "eps3": 0.005, "ns": 0.004828, "As": 2.102e-11, "beta_eff": 0.069, "c1": 0.0459, "c2": 0.0462, "c3": 0.075} print("Building full CMB (TT+TE+EE) Fisher matrix at threshold fiducial...") orig_fid = base.FIDUCIAL base.FIDUCIAL = fid F_cmb_full = base.build_fisher(varied_full, lmin=30, lmax_T=3000, lmax_P=5000, fsky=0.4) base.FIDUCIAL = orig_fid print("\n--- (A) Eigen-analysis: where is the exact degeneracy? ---") eigen_analysis(F_cmb_full, varied_full, fid) print("\n--- Effect of fixing the amplitude nuisance parameters ---") for fixed in ([], ["beta_eff"], ["beta_eff", "c1", "c2", "c3"]): r = sigma_with_fixed(F_cmb_full, varied_full, fixed) print(f" fixed={fixed}: sigma(kc)={r['kc']:.4f} sigma(eps3)={r['eps3']:.4f}") # From here on: work in the reduced 4-param space (kc, eps3, ns, As), # amplitude coefficients externally fixed -- the optimistic scenario. keep = [varied_full.index(n) for n in ["kc", "eps3", "ns", "As"]] F_cmb4 = F_cmb_full[np.ix_(keep, keep)] varied4 = ["kc", "eps3", "ns", "As"] print("\n--- (B) Adding a simplified DESI-like LSS shape prior ---") F_lss = lss_shape_fisher(fid, varied4, step) print("--- (C) Adding CMB lensing phi-phi (TT-only N0) ---") F_lens = lensing_fisher(fid, varied4, step) print("\nSummary (fiducial kc=0.30, eps3=0.02):") for label, F in [("CMB TT/TE/EE only", F_cmb4), ("+ lensing phi-phi", F_cmb4 + F_lens), ("+ simplified LSS", F_cmb4 + F_lss), ("+ lensing + LSS", F_cmb4 + F_lens + F_lss)]: s = np.sqrt(np.diag(np.linalg.inv(F))) print(f" {label:28s} sigma(kc)={s[0]:.4f} (kc/sig={fid['kc']/s[0]:.1f}) " f"sigma(eps3)={s[1]:.4f} (eps3/sig={fid['eps3']/s[1]:.2f})") if __name__ == "__main__": main()