Add simple reionization heating model

[cmbant]

This PR implements a simple reionization heating model to address issue #181.

Changes

• Add include_heating flag to TReionizationModel (default True)
• Smooth temperature evolution: Baryon temperature T_b is smoothly heated from the recombination value to ~10^4 K following the reionization ionization fraction x_e shape
• Sound speed mapping: Baryon sound speed transitions smoothly from the original pre-reionization formula to ideal gas form cs^2 = γ k_B T / (μ m_p) with γ=5/3
• Preserves original behavior: Before reionization (x_e ≈ x_e_recomb), the code uses the original cs^2 calculation. During reionization, it smoothly interpolates to the heated model.

Implementation Details

The heating weight yheat follows the reionization progression:

yheat = (x_e - x_e_recomb) / (x_e_final - x_e_recomb)

The final quantities are:

T_gas = T_recomb + yheat * (T_reion - T_recomb)
cs^2 = (1 - yheat) * cs^2_orig + yheat * cs^2_ideal

This provides a leading-order correction for reionization effects on the matter power spectrum while maintaining backward compatibility.

Testing

The implementation has been tested with matter power spectrum calculations showing the expected suppression on small scales due to increased baryon pressure during reionization.

Closes #181

---

Pull Request opened by Augment Code with guidance from the PR author

[cmbant]

[ahallcosmo]

Is the k range right in that plot? It looks much too large an effect on those scales.

[ahallcosmo]

I suppose that technically the sound speed equation you had originally is always correct (although it neglects the time evolution of the mean molecular weight, which !=0 during reionization), so the ideal gas formula with γ=5/3 is only true when Tb is evolving adiabatically, i.e. neglecting the Compton cooling which is happening even after reionization has ended. But perhaps that is negligible as it's proportional to the CMB energy density, falling as (1+z)^4? Or maybe it was that term that was causing the issue in the first place.

[cmbant]

The model is certainly physically wrong in various ways, though not obviously wrong order of magnitude?

detailed_heating_analysis.py
This is the agent's plotting code, it looked OK at quick glance but haven't independently checked..

[ahallcosmo]

OK, the plotting looks fine. I was surprised because I was only seeing changes to P(k) at k > 100 h/Mpc, not the k=1 h/Mpc scales you find here. But the Jeans scale for baryons after reionization is of order 10 h/Mpc, so that makes sense. Chatting to reionization people here, your simple model for T(z) is fine at the order of magnitude level.

[cmbant]

This enhancement in the non-linear Pk is presumably not physical.