""" tensor_coupling_candidates.py -------------------------------- Extends the RLMT Omega-layer framework with a tensor (gravitational-wave) sector, under three candidate coupling laws, and forecasts how much each lets CMB-S4 B-modes constrain the existing scalar RLMT parameters (kc, eps3) that could not be reached by LSS/21cm/mu-distortion/lensing alone. All three candidates reuse EXISTING RLMT parameters (kc, eps3, beta_eff) wherever possible and introduce only ONE new amplitude parameter, r0 (the usual tensor-to-scalar ratio at k0), to stay in the spirit of "new predictions without new unconstrained parameters" from the paper's own Outlook. CANDIDATE 1 -- "Shared dissipation" (common Lindblad damping scale): Motivation: if gravitational waves are damped by the SAME Lindblad dissipators L_a that produce the scalar mode-2 cutoff, the tensor spectrum should inherit the identical cutoff scale kc: P_t(k) = r0 * As * (k/k0)^{n_t0} * exp(-(k/kc)^2), n_t0 = 0 (fixed) Prediction: kc leaves an imprint on the B-mode spectrum's shape (suppression at ell beyond ~ kc * chi_rec), independent of eps3. CANDIDATE 2 -- "Root-tilt coupling" (tensor tilt tied to eps3): Motivation: if the same dissipative generator responsible for mode 3's extra scalar tilt eps3 also tilts the tensor spectrum by a proportional amount: P_t(k) = r0 * As * (k/k0)^{n_t}, n_t = -eps3 (illustrative O(1) coupling) Prediction: the SHAPE of the B-mode spectrum across ell directly encodes eps3, not just its amplitude via r0. CANDIDATE 3 -- "Null / standard consistency" (no coupling, baseline): Standard single-field slow-roll consistency relation, independent of kc and eps3: P_t(k) = r0 * As * (k/k0)^{-r0/8} Included to reproduce and confirm the earlier (negligible) result when there genuinely is no coupling -- the honest null case. Method: CAMB tensor Cl^BB via SplinedInitialPower.set_tensor_table, combined with a residual-lensing-B-mode + CMB-S4 noise model, Fisher matrix over (kc, eps3, r0) [+ns, As held at their previous CMB TT/TE/EE Fisher information] to see whether B-modes add real constraining power under each candidate. """ import numpy as np import camb from camb import initialpower import cmbs4_forecast as base def rlmt_tensor_spectrum(k, candidate, As, r0, kc, eps3, k0=0.05, gamma=None, kstar=None, ns=None): if candidate == 1: return r0 * As * np.ones_like(k) * np.exp(-(k / kc) ** 2) elif candidate == 2: n_t = -eps3 return r0 * As * (k / k0) ** n_t elif candidate == 3: n_t = -r0 / 8.0 return r0 * As * (k / k0) ** n_t elif candidate == 4: # Principled: spin-2 (Clebsch-Gordan 1x1->2, M=0) projection of the # SAME Omega-layer tension modes, using the SU(2)-representation # derived c_n = kappa*(2-m_n^2) [c1=c2=kappa, c3=2*kappa] so that # w1(k)=w2(k) exactly, and building the tensor sourcing amplitude # from mode 1 (untilted) and mode 3 (eps3-tilted) via the CG # coefficients for the M=0 channel: <1,1;1,-1|2,0>=<1,-1;1,1|2,0>=1/sqrt6, # <1,0;1,0|2,0>=sqrt(2/3). Normalised to r0*As at k0. cg11, cg00 = 1 / np.sqrt(6), np.sqrt(2 / 3) chik = np.arctan(k / kstar); sinp = np.sin(chik) X = np.array([gamma * 1 * sinp, gamma * 1 * sinp, gamma * 2 * sinp]) X = X - X.min(axis=0); w = np.exp(-X); w = w / w.sum(axis=0) P1 = As * (k / k0) ** (ns - 1) P3 = As * (k / k0) ** (ns - 1 + eps3) A20 = 2 * cg11 * w[0] * P1 + cg00 * w[2] * P3 chik0 = np.arctan(k0 / kstar); sinp0 = np.sin(chik0) X0 = np.array([gamma * sinp0, gamma * sinp0, 2 * gamma * sinp0]) X0 = X0 - X0.min(); w0 = np.exp(-X0); w0 = w0 / w0.sum() A20_k0 = 2 * cg11 * w0[0] * As + cg00 * w0[2] * As return r0 * As * (A20 / A20_k0) elif candidate == 5: # M=+-1 channel: cross term sqrt(w_{+1}P_1 * w_0 P_3), CG coeff sqrt(2)*1/sqrt2=1 # gives an EXACT, k-independent eps3/2 tilt (cleaner than candidate 4's # M=0 channel, whose dilution factor drifts with k as w0(k) shrinks). chik = np.arctan(k / kstar); sinp = np.sin(chik) X = np.array([gamma * 1 * sinp, gamma * 1 * sinp, gamma * 2 * sinp]) X = X - X.min(axis=0); w = np.exp(-X); w = w / w.sum(axis=0) P1 = As * (k / k0) ** (ns - 1) P3 = As * (k / k0) ** (ns - 1 + eps3) A21 = np.sqrt(2 * w[0] * P1 * w[2] * P3) chik0 = np.arctan(k0 / kstar); sinp0 = np.sin(chik0) X0 = np.array([gamma * sinp0, gamma * sinp0, 2 * gamma * sinp0]) X0 = X0 - X0.min(); w0v = np.exp(-X0); w0v = w0v / w0v.sum() A21_k0 = np.sqrt(2 * w0v[0] * As * w0v[2] * As) return r0 * As * (A21 / A21_k0) elif candidate == 6: # NEW POSTULATE (linear, not bilinear): the Lindblad-evolved density # matrix rho_Omega, viewed as an operator on the 3-dim (spin-1) # Hilbert space, contains a spin-2 (quadrupole) component beyond its # diagonal populations w_n(k). We postulate this quadrupole directly # and LINEARLY sources the tensor metric perturbation: # Q(k) = sum_n [3 m_n^2 - j(j+1)] w_n(k) P_n(k) # = 2 w_{+1}(k) P_1(k) - 2 w_0(k) P_3(k) (j=1) # unlike candidates 4/5, eps3 enters Q(k) LINEARLY (not diluted by a # sqrt or an additive mix with a much-larger untilted term). chik = np.arctan(k / kstar); sinp = np.sin(chik) X = np.array([gamma * 1 * sinp, gamma * 1 * sinp, gamma * 2 * sinp]) X = X - X.min(axis=0); w = np.exp(-X); w = w / w.sum(axis=0) P1 = As * (k / k0) ** (ns - 1) P3 = As * (k / k0) ** (ns - 1 + eps3) Q = 2 * w[0] * P1 - 2 * w[2] * P3 chik0 = np.arctan(k0 / kstar); sinp0 = np.sin(chik0) X0 = np.array([gamma * sinp0, gamma * sinp0, 2 * gamma * sinp0]) X0 = X0 - X0.min(); w0v = np.exp(-X0); w0v = w0v / w0v.sum() Q_k0 = 2 * w0v[0] * As - 2 * w0v[2] * As return r0 * As * np.abs(Q / Q_k0) raise ValueError(candidate) def get_BB(params, candidate, lmax=1000): 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"]) pt = rlmt_tensor_spectrum(base.KGRID, candidate, params["As"], params["r0"], params["kc"], params["eps3"], gamma=params.get("gamma"), kstar=params["kstar"], ns=params["ns"]) 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.set_tensor_table(base.KGRID, pt) pk_ini.effective_ns_for_nonlinear = params["ns"] pars.InitPower = pk_ini pars.WantTensors = True pars.set_for_lmax(lmax, lens_potential_accuracy=1) results = camb.get_results(pars) powers = results.get_cmb_power_spectra(pars, CMB_unit="muK", raw_cl=True) # IMPORTANT: use the tensor-ONLY BB, not 'total' (which also contains the # lensing-induced E->B leakage that depends on kc/eps3 through the SCALAR # spectrum via ordinary lensing -- an effect already captured by the # separate lensing phi-phi Fisher block, and unrelated to any new # Omega-layer tensor coupling). Mixing the two double-counts the known # lensing effect and swamps the candidate-specific tensor signal. return powers["tensor"][:, 2] # BB, muK^2, primordial tensor-only def get_BB_lensing_only(params, lmax=1000): """Lensed-scalar-only BB (no tensors) -- the lensing-induced residual that primordial B-mode searches must delens against.""" 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) results = camb.get_results(pars) powers = results.get_cmb_power_spectra(pars, CMB_unit="muK", raw_cl=True) return powers["lensed_scalar"][:, 2] _lens_cache = {} def bb_fisher(fiducial, candidate, varied, step, lmax=1000, fsky=0.4, delensing_residual=0.1): ls = np.arange(lmax + 1) NT, NE = base.noise_cl(ls) # NE ~ B-mode instrumental noise too cache_key = lmax if cache_key not in _lens_cache: _lens_cache[cache_key] = get_BB_lensing_only(fiducial, lmax=lmax)[:lmax + 1] BB_lens_res = _lens_cache[cache_key] * delensing_residual BB0 = get_BB(fiducial, candidate, lmax=lmax)[:lmax + 1] derivs = {} for name in varied: pp = dict(fiducial); pp[name] += step[name] pm = dict(fiducial); pm[name] -= step[name] BBp = get_BB(pp, candidate, lmax=lmax)[:lmax + 1] BBm = get_BB(pm, candidate, lmax=lmax)[:lmax + 1] derivs[name] = (BBp - BBm) / (2 * step[name]) Ctot = BB0 + BB_lens_res + NE n = len(varied) F = np.zeros((n, n)) for l in range(30, lmax + 1): if Ctot[l] <= 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 * derivs[ni][l] * derivs[nj][l] / Ctot[l] ** 2 F[i, j] += val if j != i: F[j, i] += val return F def main(): fid = dict(base.FIDUCIAL) fid["kc"], fid["eps3"] = 0.30, 0.02 fid["r0"] = 0.01 # below current upper limits (~0.06), a plausible target value varied = ["kc", "eps3", "r0"] step = {"kc": 0.009, "eps3": 0.005, "r0": 0.001} lmax = 700 # Previous CMB TT/TE/EE Fisher (amplitude coefficients externally fixed), # restricted to kc, eps3 (no r0 information there since scalar-only). F_full = np.load("F_threshold.npy") names8 = ["kc", "eps3", "ns", "As", "beta_eff", "c1", "c2", "c3"] results = {} for cand, label in [(1, "1: shared kc cutoff"), (2, "2: tensor tilt = -eps3"), (3, "3: null/standard consistency")]: print(f"Computing BB Fisher for candidate {cand}...") F_bb = bb_fisher(fid, cand, varied, step, lmax=lmax, fsky=0.4) results[cand] = F_bb # Build a consistent joint 5-parameter Fisher (kc, eps3, ns, As, r0): # TT/TE/EE contributes to the (kc,eps3,ns,As) block (no r0 info), # BB contributes to the (kc,eps3,r0) block (no ns,As info) -- then # marginalise over ns,As jointly, properly propagating their # degeneracy with kc/eps3 (the earlier bug took a raw 2-parameter # submatrix, which is equivalent to incorrectly fixing ns,As exactly). names5 = ["kc", "eps3", "ns", "As", "r0"] keep4 = [names8.index(n) for n in ["kc", "eps3", "ns", "As"]] F_cmb4 = F_full[np.ix_(keep4, keep4)] print(f"\n{'candidate':40s} sigma(kc) sigma(eps3) sigma(r0)") for cand, label in [(1, "1: shared kc cutoff"), (2, "2: tensor tilt = -eps3"), (3, "3: null/standard consistency")]: F_bb = results[cand] # indexed by varied = [kc, eps3, r0] F5 = np.zeros((5, 5)) F5[:4, :4] += F_cmb4 idx_map = [0, 1, 4] # kc,eps3,r0 -> positions in names5 for a, ia in enumerate(idx_map): for b, ib in enumerate(idx_map): F5[ia, ib] += F_bb[a, b] s = np.sqrt(np.diag(np.linalg.inv(F5))) print(f"{label:40s} {s[0]:.5f} {s[1]:.5f} {s[4]:.5f}") print("\n(for reference, BB-only Fisher on kc,eps3,r0 directly " "-- singular for candidates with no kc/eps3 coupling, as expected:)") for cand, label in [(1, "1: shared kc cutoff"), (2, "2: tensor tilt = -eps3"), (3, "3: null/standard consistency")]: F_bb = results[cand] try: s = np.sqrt(np.diag(np.linalg.inv(F_bb))) print(f" {label:38s} sigma(kc)={s[0]:.4f} sigma(eps3)={s[1]:.4f} sigma(r0)={s[2]:.5f}") except np.linalg.LinAlgError: print(f" {label:38s} singular (BB alone carries no info on some parameter -- expected for null coupling)") if __name__ == "__main__": main()