Reheating's Imprint: Hunting for Primordial Features in the Late-Time Matter Power Spectrum

The matter power spectrum we measure from galaxy surveys and Lyman-alpha forests encodes more than the standard cosmological model. Subtle wiggles, bumps, or slope changes can trace back to the reheating epoch—the brief but critical phase when the inflaton decayed and set the initial conditions for the hot Big Bang. For cosmologists hunting for new physics, these features are a direct line to energy scales of 10¹² GeV or higher. But finding them requires a careful strategy: you need to know what to look for, how to model it, and how to avoid mistaking systematics for signals. This guide is for researchers who already understand the matter power spectrum and want to incorporate reheating imprints into their analysis pipeline. We skip the basics of perturbation theory and go straight to the practical workflow: from generating template spectra to running MCMC chains with feature parameters.

Who Needs This and What Goes Wrong Without It

If you are analyzing data from DESI, Euclid, or the Rubin Observatory, you are likely testing for deviations from ΛCDM. Many analyses focus on dark energy or neutrino mass, but reheating imprints are often overlooked because they are subtle and degenerate with other effects. Without accounting for them, you risk misinterpreting a primordial feature as a baryon acoustic oscillation (BAO) scale shift or a spurious tilt in the spectral index. Worse, you could miss a genuine detection because your model does not include the right template.

Consider a typical power spectrum measurement at k = 0.1 h/Mpc. A damped oscillatory feature from particle production during reheating can produce a 1% modulation. If your analysis assumes a smooth power spectrum, that modulation gets absorbed into the BAO wiggle amplitude or the broadband shape, biasing your constraints on Ω_m and H₀. Several groups have reported tensions in these parameters; some of that tension could stem from unmodeled reheating physics.

Who exactly needs to care? Anyone fitting a power spectrum with a resolution better than Δk/k ∼ 0.05, or combining CMB and LSS data. If you work with the galaxy power spectrum, the matter power spectrum from simulations, or the 21 cm power spectrum from radio arrays, you are susceptible. The problem is that reheating imprints are not a single template; they depend on the inflaton decay rate, the equation of state during reheating, and possible non-perturbative effects. Without a systematic approach, you will either miss them or overfit noise.

The consequences of ignoring these features go beyond missed discoveries. They can lead to incorrect conclusions about early-universe physics. For example, a feature at k ∼ 0.01 h/Mpc could be misinterpreted as evidence for a primordial non-Gaussianity of local type, when in fact it is a relic of a sudden transition in the reheating dynamics. This guide gives you the framework to avoid such pitfalls.

Prerequisites and Context You Should Settle First

Before diving into the search for reheating features, you need a solid grasp of the standard matter power spectrum and its measurement. We assume you are comfortable with the transfer function, the growth factor, and the effects of baryons, neutrinos, and dark energy on the power spectrum shape. You should have run a Boltzmann code like CAMB or CLASS at least once and understand the output: the linear matter power spectrum P(k) at z=0 and the nonlinear correction from halofit or HMCode.

Equally important is the theoretical background on reheating. You do not need to be an expert in inflationary model-building, but you should know the basics: reheating is the process by which the inflaton field decays into Standard Model particles, raising the temperature of the universe. The duration and dynamics of this epoch affect the primordial power spectrum through two main mechanisms: (1) the inflaton's oscillatory phase can imprint a periodic modulation on the curvature perturbation, and (2) sudden changes in the equation of state can produce a bump or dip at a specific scale.

Data-wise, you need a measured or simulated matter power spectrum with error bars. This could come from a galaxy redshift survey (e.g., BOSS, eBOSS, DESI), a weak lensing survey (KiDS, HSC, LSST), or a cosmological simulation. The redshift range matters: features are more prominent at higher redshift because nonlinear evolution washes them out. For z > 1, the linear power spectrum is still a good approximation on scales k < 0.3 h/Mpc. For lower redshifts, you need to model nonlinear damping carefully.

A final prerequisite is a likelihood analysis framework. We recommend using MontePython or Cobaya for parameter estimation, with the option to add custom likelihoods. You should be able to modify the theory code to include a feature template. If you have not done this before, start by adding a simple oscillatory term to the primordial power spectrum in CLASS or CAMB and see how it changes the χ².

Core Workflow: From Reheating Model to Power Spectrum Constraints

The workflow breaks into four sequential steps: (1) choose a reheating scenario, (2) compute the primordial feature, (3) propagate it through the Boltzmann code, and (4) fit it to data. Let us walk through each.

Step 1: Select a Reheating Scenario

Reheating models fall into broad classes: perturbative decay, parametric resonance (preheating), and instant reheating. For a first search, the simplest is a sudden transition from matter-like to radiation-like expansion, parametrized by the reheating temperature T_reh and the equation-of-state parameter w_reh. This produces a step-like feature in the primordial power spectrum at a scale k_reh ∼ a_reh H_reh. More complex models with oscillatory features arise from particle production during preheating, leading to a damped sinusoidal modulation on top of the nearly scale-invariant spectrum.

Step 2: Compute the Primordial Feature

For a step feature, the primordial curvature power spectrum P_ζ(k) can be written as P_ζ,0(k) × [1 + A_step tanh((ln k - ln k_reh)/Δ)], where A_step is the amplitude and Δ controls the sharpness. For oscillatory features, the template is P_ζ,0(k) × [1 + A_osc sin(ω ln k + φ) exp(-(ln k - ln k₀)²/(2σ²))]. Choose the form that matches your theoretical prior. We recommend starting with a single damped oscillation, as it is the most generic prediction.

Step 3: Propagate Through Boltzmann Code

Modify the primordial power spectrum in CLASS or CAMB. In CLASS, you can use the 'primordial_custom' module to input a tabulated P_ζ(k). Generate a fine grid in k from 10^-4 to 10 Mpc^-1 with the feature superimposed. Then run CLASS to get the linear matter power spectrum at the desired redshift. Be careful with k-grid resolution: the feature oscillations should be sampled with at least 10 points per cycle, or you will get numerical damping.

Step 4: Fit to Data

Use your likelihood code to vary A_step, k_reh, Δ (or A_osc, ω, etc.) together with the standard ΛCDM parameters. Include priors: for instance, a uniform prior on ln k_reh between ln(10^-4) and ln(10), and a Jeffreys prior on A_step to avoid volume effects. Run MCMC chains until convergence (R-1 < 0.01). Check the posterior for degeneracies: A_step often correlates with n_s and A_s. If you see a bimodal distribution, it may indicate that the feature is mimicking the BAO wiggles.

Tools, Setup, and Environment Realities

You will need a Boltzmann solver (CLASS or CAMB), a sampler (MontePython, Cobaya, or emcee), and a way to compute the matter power spectrum from data (e.g., the power spectrum measurement code from the survey team). For this guide, we assume you are using CLASS and Cobaya, as they are well-integrated and allow custom likelihoods.

Setting Up the Custom Primordial Spectrum

In CLASS, create a new .ini file that includes 'primordial custom = yes' and provide a file with k and P_ζ columns. We suggest writing a Python script that generates the file for each set of feature parameters. Then call CLASS from Cobaya via the 'classy' wrapper. This allows you to pass the feature parameters as nuisance parameters that modify the input file before each CLASS run. Be aware that CLASS can be slow for high-resolution k grids; we recommend using the 'k_max = 10' and 'P_k_max_h/Mpc = 10' settings to cover the range needed for LSS.

Handling Nonlinearities

For scales k > 0.1 h/Mpc, nonlinear evolution damps features and introduces mode coupling. Use a nonlinear prescription like HMCode or a tuned emulator. However, note that the damping is scale-dependent and can erase oscillatory features. A common fix is to apply a Gaussian damping envelope exp(-k²/k_nl²) to the feature component, where k_nl is the nonlinear scale at the redshift of interest. This is an approximation; for precision work, run N-body simulations with the feature imprinted in the initial conditions and measure the nonlinear power spectrum directly.

Computational Budget

A typical MCMC run with 5 extra feature parameters and 6 ΛCDM parameters takes about 10⁵ likelihood evaluations. Each CLASS evaluation takes ~1 second, so total runtime is ~30 hours on a single core. Parallelize over chains (e.g., 8 chains on 8 cores) to reduce wall time. If that is too heavy, use an emulator for the power spectrum, such as the one from the Cosmopower project, but ensure it covers the feature parameter space.

Variations for Different Constraints

Not all searches are the same. Your approach will depend on the data quality, redshift, and scale range. Here we outline three common scenarios.

Scenario A: High-Redshift Galaxy Surveys (z > 1)

Surveys like DESI’s ELG sample or future Euclid data are ideal for reheating features because nonlinear damping is mild. Use the linear matter power spectrum up to k = 0.3 h/Mpc. The main challenge is the limited volume, which increases cosmic variance on large scales. For a feature at k ∼ 0.01 h/Mpc, you need a survey volume > 10 Gpc³ to get a >2σ detection. Combine multiple tracers (ELGs, QSOs, LRGs) to beat cosmic variance via multitracer technique.

Scenario B: CMB Lensing Power Spectrum

The CMB lensing convergence power spectrum C_l^κκ probes the matter power spectrum integrated along the line of sight, weighted by the lensing kernel. This is sensitive to features at z ∼ 2, but the projection smears out oscillations. Use the Planck lensing data or SPT-3G for l < 2000. The advantage is that the CMB lensing power spectrum is less affected by nonlinear bias and redshift-space distortions. However, the signal-to-noise is lower than galaxy surveys for small-scale features.

Scenario C: Low-Redshift Weak Lensing (z < 1)

KiDS and HSC data probe the matter power spectrum at z ∼ 0.5, where nonlinear effects are strong. Here, the feature is heavily damped, and you may only detect the broad step-like modulation, not oscillations. Use a nonlinear model that includes the feature as a correction to the initial power spectrum. A practical approach is to run a suite of N-body simulations with different feature amplitudes and train a Gaussian process emulator for the nonlinear power spectrum. This is computationally expensive but necessary for accurate constraints.

Pitfalls, Debugging, and What to Check When It Fails

Even with a careful setup, things can go wrong. Here are the most common issues and how to diagnose them.

Degeneracy with BAO Wiggles

The BAO signal is a series of oscillations in the power spectrum with a period of Δk ∼ 0.04 h/Mpc. A reheating oscillation with a similar period can be degenerate. To break this, include a BAO template as a separate component and marginalize over its amplitude. Alternatively, use the phase of the oscillation: BAO wiggles have a specific phase (peaks at k = nπ/s), while reheating features can have any phase. If your posterior shows a strong correlation between the feature frequency ω and the BAO scale r_s, then degeneracy is present. Fix the BAO scale to the Planck value and test if the feature remains.

Numerical Artifacts from Boltzmann Solvers

CLASS and CAMB use interpolation for the transfer function. If your feature has sharp features (e.g., a step with Δ < 0.1), the interpolation may smooth it out. Check by comparing the output power spectrum with a high-resolution run (k_max = 50, l_max = 5000). If the feature amplitude is reduced by more than 10%, increase the sampling. Also, ensure that the k-grid for the primordial spectrum is fine enough: at least 1000 points between k_min and k_max.

Prior Volume Effects

When you add a feature parameter like A_step with a uniform prior, the posterior may show a peak at zero simply because the volume of parameter space is larger for small amplitudes. Use a Jeffreys prior (uniform in log A) or a prior based on the expected signal-to-noise. Alternatively, run a profile likelihood to check if the best-fit is significantly different from zero.

Data Systematics

Survey systematics like fiber collisions, photometric redshift errors, and nonlinear bias can mimic features. Always test your pipeline on mock data that includes known systematics. The mock should be generated from a smooth power spectrum and then processed through the same analysis. If the pipeline detects a feature in the mock, your systematics model is insufficient.

FAQ: Common Questions from Practitioners

Q: Can I use the CMB temperature power spectrum instead of the matter power spectrum?
A: Yes, but the CMB is less sensitive to small-scale features because of Silk damping and projection. For features at k > 0.01 Mpc^-1, the matter power spectrum is better. However, the CMB can constrain features at large scales (k < 0.01 Mpc^-1) where cosmic variance dominates LSS.

Q: How do I choose between a step and an oscillatory template?
A: The template should be guided by theory. If you have a specific reheating model (e.g., quadratic inflaton potential with perturbative decay), compute the exact primordial spectrum and use that. For a model-agnostic search, start with a damped oscillation, as it is the most generic. You can then compare the Bayesian evidence for different templates.

Q: What if my data shows a feature at the same scale as the BAO?
A: This is a common degeneracy. First, check if the feature is present in the power spectrum after subtracting the best-fit BAO model. If it remains, it is likely a reheating feature. You can also cross-check with the correlation function: BAO appears as a peak at 100 h^-1Mpc, while a reheating feature may not.

Q: How many parameters can I add before overfitting?
A: With current data (DESI, Planck), you can add up to 3–4 feature parameters without significant overfitting, as long as you use a Bayesian evidence criterion. For future surveys, the number may increase. Always do a simulation-based calibration: inject a known feature and check if your pipeline recovers it.

Q: Do I need to include nonlinear effects for k < 0.1 h/Mpc?
A: Not for linear theory, but if you are using the galaxy power spectrum, you need to model bias and redshift-space distortions. Those can also introduce scale-dependent features. We recommend using a perturbation theory model (e.g., EFTofLSS) for the galaxy power spectrum up to k = 0.2 h/Mpc.

What to Do Next: Specific Steps to Advance Your Search

You now have a workflow and know the pitfalls. Here are concrete next actions to move from theory to results:

1. Implement the template in CLASS or CAMB. Start with a damped oscillation. Run the code with a few parameter values and verify that the output matter power spectrum looks as expected (e.g., a 1% oscillation at k = 0.05 h/Mpc).
2. Run a simulation-based null test. Generate a mock matter power spectrum from a smooth ΛCDM model with the same noise and window function as your survey. Run your MCMC pipeline on it. If you detect a feature at >2σ, your pipeline has a problem—fix it before moving to real data.
3. Combine with CMB lensing. Use the Planck 2018 lensing power spectrum (or ACT DR6) to add independent information. The cross-correlation between CMB lensing and galaxy clustering can help break degeneracies between the feature and the bias.
4. Prepare for Stage-IV surveys. For Euclid and Rubin, the sensitivity will be an order of magnitude better. Generate mock data for these surveys with a known reheating feature and test your detection threshold. This will tell you what feature amplitudes are detectable and help you design the analysis pipeline early.
5. Write up your results as a methods paper. Even a null detection is valuable: it places upper limits on the amplitude of reheating features, constraining the physics of the early universe. Share your templates and chains to allow others to combine constraints.