""" cmbs4_forecast.py ------------------ Fisher-matrix forecast of CMB-S4 sensitivity to the RLMT primordial power-spectrum parameters, with emphasis on kc and eps3 (Sec. 9, Outlook item 1 of the RLMT/Planck 2018 paper). Method ------ 1. Rebuild the RLMT primordial spectrum P(k) exactly as in Eq. (5.1)-(5.5) of the paper and feed it to CAMB via SplinedInitialPower (same numerical setup as the paper's Planck fit). 2. Compute unlensed C_l^{TT,TE,EE} at a fiducial point and at points shifted by +-delta_i for every varied parameter (central finite differences). 3. Build a simplified single-channel CMB-S4 "Wide" noise model (Delta_T ~ 1 uK-arcmin, beam FWHM ~ 1 arcmin, fsky ~ 0.4, lmin=30, lmax_T=3000, lmax_P=5000) -- this is a standard first-pass forecast noise level used in the CMB-S4 Science Book; a production forecast would use the full multi-frequency noise curves + foreground residuals instead of this single effective channel. 4. Fisher matrix: F_ij = sum_l (2l+1)/2 * fsky * Tr[ C_l^-1 dC_l/dtheta_i C_l^-1 dC_l/dtheta_j ] with C_l the 2x2 (TT,EE) matrix including TE off-diagonal and noise. 5. Invert F to get marginalised 1-sigma forecast errors on all parameters, in particular kc and eps3, and compare to the detectability thresholds quoted in the paper (kc <~ 0.3 Mpc^-1 or |eps3| >~ 0.02). Notes / assumptions the user should double check against their own rlmt_primordial.py: - The exponent p in X_n(k) = beta_eff * c_n * sin^p(chi_k) is not given a numeric value in the paper text; it is left as a keyword here (default p=1). - Only TT/TE/EE are used in the Fisher sum (matches "lite" Planck likelihood structure used in the paper). CMB lensing (phi-phi) can be added as an extra diagonal block once a CMB-S4 lensing-noise (N0) curve is supplied. - Step sizes for finite differences are set relative (1%) or absolute for parameters that can cross zero (eps3, kstar-independent). Run: python3 cmbs4_forecast.py """ import numpy as np import camb from camb import initialpower # ---------------------------------------------------------------------- # 1. RLMT primordial power spectrum (paper Eqs. 5.1-5.5) # ---------------------------------------------------------------------- def rlmt_power_spectrum(k, As, ns, beta_eff, c1, c2, c3, kstar, kc, eps3, p=1.0, k0=0.05): """Three-mode softmax mixture primordial spectrum, Eq. (5.5).""" chik = np.arctan(k / kstar) sinp = np.sin(chik) ** p c = (c1, c2, c3) X = np.array([beta_eff * ci * sinp for ci in c]) X = X - X.min(axis=0) # numerically stable softmax w = np.exp(-X) w = w / w.sum(axis=0) P1 = As * (k / k0) ** (ns - 1) P2 = As * (k / k0) ** (ns - 1) * np.exp(-(k / kc) ** 2) P3 = As * (k / k0) ** (ns - 1 + eps3) return w[0] * P1 + w[1] * P2 + w[2] * P3 # ---------------------------------------------------------------------- # 2. CAMB wrapper # ---------------------------------------------------------------------- PARAM_NAMES = ["H0", "ombh2", "omch2", "tau", "As", "ns", "beta_eff", "c1", "c2", "c3", "kstar", "kc", "eps3"] FIDUCIAL = dict(H0=67.20, ombh2=0.02235, omch2=0.12031, tau=0.05425, As=2.102e-9, ns=0.9655, beta_eff=2.30, c1=1.53, c2=1.54, c3=2.50, kstar=0.273, kc=2.83, eps3=-0.0017) # relative (fraction of value) step sizes; absolute overrides for # parameters that can vanish or change sign REL_STEP = {"H0": 0.005, "ombh2": 0.01, "omch2": 0.01, "tau": 0.03, "As": 0.01, "ns": 0.005, "beta_eff": 0.03, "c1": 0.03, "c2": 0.03, "c3": 0.03, "kstar": 0.03, "kc": 0.03} ABS_STEP = {"eps3": 0.005} LMAX = 3000 KGRID = np.logspace(-6, np.log10(50), 2000) def get_cls(params, lmax=LMAX): """Return unlensed TT, EE, TE (muK^2, raw C_l, l=0..lmax).""" pk = rlmt_power_spectrum(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(KGRID, pk) pk_ini.effective_ns_for_nonlinear = params["ns"] pars.InitPower = pk_ini pars.set_for_lmax(lmax, lens_potential_accuracy=0) pars.DoLensing = False results = camb.get_results(pars) powers = results.get_cmb_power_spectra(pars, CMB_unit="muK", raw_cl=True, spectra=["unlensed_scalar"]) cl = powers["unlensed_scalar"] # columns: TT, EE, BB, TE ls = np.arange(cl.shape[0]) return ls, cl[:, 0], cl[:, 1], cl[:, 3] # TT, EE, TE # ---------------------------------------------------------------------- # 3. CMB-S4-like noise model (single effective channel, "Wide" survey) # ---------------------------------------------------------------------- def noise_cl(ls, delta_T_uKarcmin=1.0, beam_fwhm_arcmin=1.0, pol_factor=np.sqrt(2)): """Deconvolved white + beam noise, Knox (1995) formula.""" arcmin_to_rad = np.pi / 180.0 / 60.0 theta_b = beam_fwhm_arcmin * arcmin_to_rad sigma_T = delta_T_uKarcmin * arcmin_to_rad sigma_E = pol_factor * sigma_T beam2 = np.exp(ls * (ls + 1) * theta_b ** 2 / (8 * np.log(2))) NT = sigma_T ** 2 * beam2 NE = sigma_E ** 2 * beam2 return NT, NE # ---------------------------------------------------------------------- # 4. Fisher matrix # ---------------------------------------------------------------------- def build_fisher(varied_params, lmin=30, lmax_T=3000, lmax_P=5000, fsky=0.4): n = len(varied_params) derivs = {} ls, TT0, EE0, TE0 = get_cls(FIDUCIAL, lmax=max(lmax_T, lmax_P)) NT, NE = noise_cl(ls) for name in varied_params: step = FIDUCIAL[name] * REL_STEP.get(name, 0.02) if name not in ABS_STEP \ else ABS_STEP[name] p_plus = dict(FIDUCIAL); p_plus[name] = FIDUCIAL[name] + step p_minus = dict(FIDUCIAL); p_minus[name] = FIDUCIAL[name] - step _, TTp, EEp, TEp = get_cls(p_plus, lmax=max(lmax_T, lmax_P)) _, TTm, EEm, TEm = get_cls(p_minus, lmax=max(lmax_T, lmax_P)) dTT = (TTp - TTm) / (2 * step) dEE = (EEp - EEm) / (2 * step) dTE = (TEp - TEm) / (2 * step) derivs[name] = (dTT, dEE, dTE) print(f" computed derivative for {name} (step={step:.4g})") F = np.zeros((n, n)) for l in range(lmin, len(ls)): use_T = l <= lmax_T use_P = l <= lmax_P if not (use_T or use_P): continue # build 2x2 covariance in (T,E) restricted to available fields Cl = np.array([[TT0[l] + NT[l] if use_T else 1e30, TE0[l] if (use_T and use_P) else 0.0], [TE0[l] if (use_T and use_P) else 0.0, EE0[l] + NE[l] if use_P else 1e30]]) try: Cl_inv = np.linalg.inv(Cl) except np.linalg.LinAlgError: continue for i, ni in enumerate(varied_params): dTTi, dEEi, dTEi = derivs[ni] dCi = np.array([[dTTi[l], dTEi[l]], [dTEi[l], dEEi[l]]]) Mi = Cl_inv @ dCi for j, nj in enumerate(varied_params): if j < i: continue dTTj, dEEj, dTEj = derivs[nj] dCj = np.array([[dTTj[l], dTEj[l]], [dTEj[l], dEEj[l]]]) Mj = Cl_inv @ dCj val = 0.5 * (2 * l + 1) * fsky * np.trace(Mi @ Mj) F[i, j] += val if j != i: F[j, i] += val return F # ---------------------------------------------------------------------- # 5. Main # ---------------------------------------------------------------------- def main(): varied = ["H0", "ombh2", "omch2", "tau", "As", "ns", "beta_eff", "c1", "c2", "c3", "kstar", "kc", "eps3"] print("Building CMB-S4 Fisher matrix for RLMT parameters...") F = build_fisher(varied) cov = np.linalg.inv(F) sigmas = np.sqrt(np.diag(cov)) print("\nForecast 1-sigma marginalised errors (CMB-S4-like, fsky=0.4, " "lmax_T=3000, lmax_P=5000, Delta_T=1 uK-arcmin, beam=1 arcmin):\n") for name, s in zip(varied, sigmas): print(f" sigma({name:9s}) = {s:.4g} (fiducial {FIDUCIAL[name]:.4g})") idx_kc = varied.index("kc") idx_eps3 = varied.index("eps3") print(f"\nkc: fiducial {FIDUCIAL['kc']:.3f} +- {sigmas[idx_kc]:.3f} Mpc^-1") print(f"eps3: fiducial {FIDUCIAL['eps3']:.4f} +- {sigmas[idx_eps3]:.4f}") print("\nDetectability check against paper's Outlook thresholds:") print(f" kc <~ 0.3 Mpc^-1 needed for a cutoff signature within reach -> " f"CMB-S4 can resolve kc at the {sigmas[idx_kc]:.3f} Mpc^-1 level") print(f" |eps3| >~ 0.02 needed for a detectable tilt -> " f"CMB-S4 sigma(eps3) = {sigmas[idx_eps3]:.4f}, i.e. " f"{'sub-threshold (good)' if sigmas[idx_eps3] < 0.02 else 'above threshold (not yet sufficient)'}") np.save("fisher_matrix_cmbs4.npy", F) np.save("fisher_covariance_cmbs4.npy", cov) print("\nSaved fisher_matrix_cmbs4.npy and fisher_covariance_cmbs4.npy") if __name__ == "__main__": main()