"""AC78 linear benchmark: where do warm and fuzzy dark matter cut the power? Inputs : reference-grid.json (frozen CLASS v3.4.2 output, see generate_frozen_reference_grid.py). No CLASS needed here. Needs : Python 3.9+, NumPy, SciPy. Run : python3 benchmark_linear_cutoffs.py reference-grid.json linear-output.json Three measurements, each compared with something it was not fitted to: 1. Warm dark matter. The Viel et al. (2005, arXiv:astro-ph/0501562v2, eqs. 6-7) fitting function T(k) = [1 + (alpha k)^(2 nu)]^(-5/nu) is compared with the CLASS thermal-relic transfer functions stored in the frozen grid. We locate the half-mode wavenumber, where the amplitude ratio T = 1/2 (power ratio 1/4, the convention of Schneider et al. 2012 and Nadler et al. 2021), and the half-mode mass M_hm = (4 pi / 3) rho_m (pi / k_hm)^3. The 3 keV mass is also compared with the published value M_hm = 5e8 Msun (m / 3 keV)^(-10/3) of Nadler et al. (2021, eq. 5). 2. Fuzzy dark matter. The Hu, Barkana & Gruzinov (2000, arXiv:astro-ph/0003365v2, eq. 8) fit T = cos(x^3) / (1 + x^8) is compared with our own integration of the linear growth equation that includes the quantum-pressure term hbar^2 k^4 / (4 m^2 a^4) of the Madelung fluid, on a radiation + matter + Lambda background. The comparison point is the wavenumber where the power ratio T^2 falls to 1/2, the convention of their eq. 9. 3. Failure cases that the reader should be able to reproduce. All wavenumbers in the output are in h/Mpc unless a key ends in PerMpc. """ import json import sys import numpy as np from scipy.integrate import solve_ivp from scipy.optimize import brentq HBAR_OVER_M22_MPC_KMS = 0.019171 # hbar / (1e-22 eV / c^2) in Mpc km/s RHO_CRIT_H2 = 2.77536627e11 # Msun / Mpc^3 for h = 1 OMEGA_RADIATION_H2 = 2.473e-5 * (1 + 0.2271 * 3.046) RECOMBINATION_SCALE_FACTOR = 1 / 1100.0 def viel_alpha(mass_kev, omega_warm, hubble): """Break scale alpha in Mpc/h, Viel et al. (2005) eq. 7, thermal relics.""" return 0.049 * mass_kev ** -1.11 * (omega_warm / 0.25) ** 0.11 * (hubble / 0.7) ** 1.22 def viel_transfer(k, alpha, nu=1.12): return (1 + (alpha * k) ** (2 * nu)) ** (-5 / nu) def viel_half_mode(alpha, nu=1.12): return (2 ** (nu / 5) - 1) ** (1 / (2 * nu)) / alpha def crossing(k, values, level): """Log-log interpolated first wavenumber where values drops below level.""" index = int(np.argmax(np.asarray(values) < level)) if index == 0: raise ValueError('no crossing') k0, k1 = np.log(k[index - 1]), np.log(k[index]) v0, v1 = np.log(values[index - 1]), np.log(values[index]) return float(np.exp(k0 + (np.log(level) - v0) * (k1 - k0) / (v1 - v0))) def half_mode_mass(k_hm, omega_matter, hubble): """(4 pi / 3) rho_m (pi / k_hm)^3 with k_hm in h/Mpc; result in Msun.""" radius_mpc_over_h = np.pi / k_hm mass_over_h = 4 * np.pi / 3 * omega_matter * RHO_CRIT_H2 * radius_mpc_over_h ** 3 return mass_over_h / hubble def hu_transfer(k_per_mpc, m22): x = 1.61 * m22 ** (1 / 18) * k_per_mpc / (9 * m22 ** 0.5) return np.cos(x ** 3) / (1 + x ** 8) class Background: def __init__(self, hubble, omega_matter, omega_cold): self.hubble, self.omega_matter, self.omega_cold = hubble, omega_matter, omega_cold self.omega_radiation = OMEGA_RADIATION_H2 / hubble ** 2 self.omega_lambda = 1 - omega_matter - self.omega_radiation def expansion_squared(self, a): return self.omega_radiation / a ** 4 + self.omega_matter / a ** 3 + self.omega_lambda def hubble_kms_mpc(self, a): return 100 * self.hubble * np.sqrt(self.expansion_squared(a)) def log_derivative(self, a): e2 = self.expansion_squared(a) return (-4 * self.omega_radiation / a ** 4 - 3 * self.omega_matter / a ** 3) / (2 * e2) def clustering_fraction(self, a): """Only dark matter clusters before recombination; all matter afterwards.""" omega = self.omega_cold if a < RECOMBINATION_SCALE_FACTOR else self.omega_matter return omega / a ** 3 / self.expansion_squared(a) def growth(k_per_mpc, m22, background, method='DOP853', fixed_steps=None, initial_slope=1.0, start_factor=1.0): """delta(a = 1) for comoving k, starting at horizon entry with delta = 1, d delta / d ln a = 1. m22 = inf gives cold dark matter. The quantum-pressure term is Q = (hbar k^2 / (2 m a^2 H))^2, the Madelung sound speed c^2 = hbar^2 k^2 / (4 a^2 m^2). """ entry = brentq(lambda a: k_per_mpc - a * background.hubble_kms_mpc(a) / 299792.458, 1e-12, 1.0) start = np.log(max(entry * start_factor, 1e-9)) def right_hand_side(log_a, state): a = np.exp(log_a) h = background.hubble_kms_mpc(a) quantum = 0.0 if np.isinf(m22) else (HBAR_OVER_M22_MPC_KMS / m22 * k_per_mpc ** 2 / (2 * a * a * h)) ** 2 source = 1.5 * background.clustering_fraction(a) return [state[1], -(2 + background.log_derivative(a)) * state[1] + (source - quantum) * state[0]] if fixed_steps: grid = np.linspace(start, 0.0, fixed_steps) state = np.array([1.0, initial_slope]) for left, right in zip(grid[:-1], grid[1:]): step = right - left k1 = np.array(right_hand_side(left, state)) k2 = np.array(right_hand_side(left + step / 2, state + step / 2 * k1)) k3 = np.array(right_hand_side(left + step / 2, state + step / 2 * k2)) k4 = np.array(right_hand_side(right, state + step * k3)) state = state + step / 6 * (k1 + 2 * k2 + 2 * k3 + k4) return float(state[0]) solution = solve_ivp(right_hand_side, (start, 0.0), [1.0, initial_slope], method=method, rtol=1e-9, atol=1e-12) if not solution.success: raise RuntimeError(solution.message) return float(solution.y[0, -1]) def fuzzy_benchmark(m22, background, wavenumbers_per_mpc): ratios = [] for k in wavenumbers_per_mpc: ratios.append(growth(k, m22, background) / growth(k, np.inf, background)) ratios = np.array(ratios) power_ratio = ratios ** 2 ode_half = crossing(wavenumbers_per_mpc, power_ratio, 0.5) fit_half = brentq(lambda k: hu_transfer(k, m22) ** 2 - 0.5, 0.1, 9 * m22 ** 0.5 / (1.61 * m22 ** (1 / 18))) return { 'm22': m22, 'odePowerHalfPerMpc': ode_half, 'fitPowerHalfPerMpc': fit_half, 'paperEquation9PerMpc': 4.5 * m22 ** (4 / 9), 'odeOverFit': ode_half / fit_half, 'sampledPerMpc': [round(float(k), 5) for k in wavenumbers_per_mpc], 'odeAmplitudeRatio': [round(float(r), 6) for r in ratios], 'fitAmplitudeRatio': [round(float(hu_transfer(k, m22)), 6) for k in wavenumbers_per_mpc], } def main(grid_path, output_path): with open(grid_path, encoding='utf-8') as handle: grid = json.load(handle) hubble = grid['parameters']['h'] omega_matter = grid['omegaMatter'] omega_cold = grid['omegaCdm'] / hubble ** 2 k = np.array(grid['k']) warm_rows = [] for warm in grid['warm']: alpha = viel_alpha(warm['massKeV'], omega_cold, hubble) fit_half = viel_half_mode(alpha) class_half = crossing(k, warm['transfer'], 0.5) fit_values = viel_transfer(k, alpha) well_resolved = np.array(warm['transfer']) > 0.25 worst = float(np.max(np.abs(fit_values[well_resolved] / np.array(warm['transfer'])[well_resolved] - 1))) warm_rows.append({ 'massKeV': warm['massKeV'], 'alphaMpcOverH': alpha, 'fitHalfModeK': fit_half, 'classHalfModeK': class_half, 'fitOverClass': fit_half / class_half, 'fitHalfModeMassMsun': half_mode_mass(fit_half, omega_matter, hubble), 'classHalfModeMassMsun': half_mode_mass(class_half, omega_matter, hubble), 'worstFitErrorWhileTAbove0p25': worst, }) # Published reference value, in the authors' own cosmology (h = 0.7, Omega_m = 0.286). nadler_alpha = viel_alpha(3.0, 0.286 - 0.046, 0.7) nadler_check = { 'massKeV': 3.0, 'publishedMsun': 5e8, 'ourMsunInTheirCosmology': half_mode_mass(viel_half_mode(nadler_alpha), 0.286, 0.7), 'assumedOmegaWarm': 0.24, } background = Background(hubble, omega_matter, omega_cold) fuzzy = [fuzzy_benchmark(m22, background, np.logspace(np.log10(0.6), np.log10(4.0), 23) * 4.5 * m22 ** (4 / 9)) for m22 in (1.0, 4.0, 20.0)] # Failure 1: the FDM fit changes sign; T^2 then grows again although the ODE stays suppressed. x_zero = (np.pi / 2) ** (1 / 3) k_zero = x_zero * 9 / 1.61 fit_rebound = float(hu_transfer(1.6 * k_zero, 1.0) ** 2) # Failure 2: a coarse fixed-step integrator cannot follow quantum-pressure oscillations. k_hard = 30.0 adaptive = growth(k_hard, 1.0, background) / growth(k_hard, np.inf, background) coarse = growth(k_hard, 1.0, background, fixed_steps=40) / growth(k_hard, np.inf, background, fixed_steps=40) # Sensitivity: the crude start (delta = 1 at horizon entry) moves the cutoff; record how much. sensitivity = [] for slope, factor in ((1.0, 1.0), (0.0, 1.0), (1.0, 0.1)): ratios = np.array([growth(kk, 1.0, background, initial_slope=slope, start_factor=factor) / growth(kk, np.inf, background, initial_slope=slope, start_factor=factor) for kk in fuzzy[0]['sampledPerMpc']]) sensitivity.append({'initialSlope': slope, 'startFactor': factor, 'powerHalfPerMpc': crossing(np.array(fuzzy[0]['sampledPerMpc']), ratios ** 2, 0.5)}) # Failure 3: the Viel fit was calibrated below 5 h/Mpc for light relics; its deep tail is not CLASS. tail = [] for warm in grid['warm']: alpha = viel_alpha(warm['massKeV'], omega_cold, hubble) index = int(np.argmax(np.array(warm['transfer']) < 0.02)) tail.append({'massKeV': warm['massKeV'], 'k': float(k[index]), 'classT': warm['transfer'][index], 'fitT': float(viel_transfer(k[index], alpha))}) output = { 'inputGrid': grid_path, 'generator': grid['generator'], 'cosmology': {'h': hubble, 'omegaMatter': omega_matter, 'omegaCold': omega_cold}, 'warmHalfModes': warm_rows, 'nadlerReference': nadler_check, 'fuzzyCutoffs': fuzzy, 'failures': { 'fuzzyFitSignChange': {'kPerMpcAtZero': k_zero, 'fitPowerAt1p6TimesZero': fit_rebound, 'note': 'cos(x^3) turns negative above x = (pi/2)^(1/3); the fit is only a cutoff shape.'}, 'coarseIntegrator': {'kPerMpc': k_hard, 'm22': 1.0, 'adaptiveRatio': adaptive, 'fixed40StepRatio': coarse}, 'fuzzyStartSensitivity': sensitivity, 'warmFitDeepTail': tail, }, 'versions': {'numpy': np.__version__, 'scipy': __import__('scipy').__version__, 'python': sys.version.split()[0]}, } with open(output_path, 'w', encoding='utf-8') as handle: json.dump(output, handle, indent=1) handle.write('\n') for row in warm_rows: print('WDM %.1f keV k_hm fit %.2f CLASS %.2f ratio %.3f M_hm %.3g Msun' % ( row['massKeV'], row['fitHalfModeK'], row['classHalfModeK'], row['fitOverClass'], row['classHalfModeMassMsun'])) print('Nadler 3 keV: published 5e8, ours %.3g' % nadler_check['ourMsunInTheirCosmology']) for row in fuzzy: print('FDM m22=%g power-half ODE %.3f fit %.3f eq9 %.3f per Mpc ratio %.3f' % ( row['m22'], row['odePowerHalfPerMpc'], row['fitPowerHalfPerMpc'], row['paperEquation9PerMpc'], row['odeOverFit'])) print('coarse integrator at k=30/Mpc: adaptive %.4g vs 40 fixed steps %.4g' % (adaptive, coarse)) for row in sensitivity: print('FDM m22=1 start slope %g factor %g: power-half %.3f per Mpc' % (row['initialSlope'], row['startFactor'], row['powerHalfPerMpc'])) if __name__ == '__main__': main(sys.argv[1] if len(sys.argv) > 1 else 'reference-grid.json', sys.argv[2] if len(sys.argv) > 2 else 'linear-output.json')