HyRex.hydrogen module

class HyRex.hydrogen.hydrogen_model(xe_4He, lna_4He, lna_end, last_4He_lna, twog_redshift, integration_spacing=0.0005, Nsteps=800, swift=Array([[1.775000e+03, -3.631352e-02, 3.584254e-02, -7.258989e-01, 1.320773e-04], ..., [4.415000e+03, -5.489227e-02, 6.942021e-02, -4.241939e+00, 6.363997e-04]], dtype=float64))

Bases: Module

Hydrogen recombination model implementation.

Computes hydrogen ionization fraction evolution through multiple phases: Saha equilibrium, post-Saha expansion, HYREC-2 EMLA with two-photon processes, and late-time EMLA-only evolution.

Attributes:

• integration_spacing float

  Step size in log scale factor for integration

• swift array

  Tabulated SWIFT correction function for recombination rates

• concrete_axis_size array

  Pre-allocated array for post-Saha phase integration

• xe_4He array

  Helium ionization fraction from previous calculation

• lna_4He array

  Log scale factor grid from helium recombination

• last_4He_lna float

  Final log scale factor from helium phase

• twog_redshift float

  Redshift at which two-photon processes become negligible

• lna_end float

  Final log scale factor for hydrogen evolution

Methods:

• get_hydrogen_history :

  Compute full hydrogen recombination history (units: dimensionless)

• Saha_equilibrium :

  Compute Saha equilibrium phase (units: dimensionless)

• post_Saha_expansion :

  Compute post-Saha expansion phase (units: dimensionless)

• solve_emla_twophoton :

  Solve EMLA with two-photon processes (units: dimensionless)

• solve_emla :

  Solve EMLA-only evolution (units: dimensionless)

• xe_derivative_twophoton :

  Compute ionization fraction derivative with two-photon (units: dimensionless)

• xe_tm_derivative :

  Compute coupled xe and Tm derivatives (units: dimensionless, eV)

• dxe_dlna_twophoton :

  Compute two-photon recombination rate (units: dimensionless)

• get_current_correction_func :

  Interpolate SWIFT correction function (units: dimensionless)

• steady_state_equations :

  Set up steady-state level population equations (units: dimensionless)

• Lyn_esc_rate :

  Compute Lyman escape rate (units: s^{-1})

Lyn_esc_rate(n, H, nH, x1s)

Computes the Lyman-n escape rate, rate at which photons redshift past the Lyman-n line without being absorbed. We use the convention that n=2 is Ly-alpha, n=3 is Ly-beta…

Parameters:

• n float

  Requested Lyman transition level, should be greater than 2.

• H float

  Hubble parameter in s^-1.

• nH float

  Hydrogen number density in cm^-3.

• x1s float

  Fraction of 1s bound hydrogen.

Returns:

• RLyn float

  Rate of escape of Lyman-n level.

TLA_xe_deriv(lna, state, args)

Compute coupled derivatives for ionization fraction and matter temperature using Peebles three-level atom.

Parameters:

• lna float

  Log scale factor

• state array

  Current state [xe, Tm]

• args tuple

  h, omega_b, omega_cdm, Neff, YHe; the Hubble parameter, the baryon denisty Omega_b h^2, the CDM density Omega_cdm h^2, the effecgive number of neutrinos, and the helium fraction

Returns:

array

Time derivatives [dxe/dlna, dTm/dlna] (units: dimensionless, eV)

dxe_dlna_twophoton(xe, TCMB, Tm, H, nH, Delta)

Compute two-photon recombination rate.

Calculates ionization fraction evolution rate including two-photon transitions using effective multilevel atom approximation.

Parameters:

• xe float

  Current ionization fraction

• TCMB float

  CMB temperature (units: eV)

• Tm float

  Matter temperature (units: eV)

• H float

  Hubble parameter (units: s^{-1})

• nH float

  Hydrogen number density (units: cm^{-3})

• Delta float

  Correction function value

Returns:

float

Recombination rate dxe/dlna (units: dimensionless)

get_current_correction_func(TCMB, omega_b, omega_cdm, YHe, Neff)

Interpolate correction function for current cosmology.

Interpolates correction function and applies cosmological parameter derivatives at current CMB temperature for accurate recombination rates.

Parameters:

• TCMB float

  CMB temperature (units: eV)

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

Returns:

float

Correction function value (units: dimensionless)

get_hydrogen_history(h, omega_b, omega_cdm, Neff, YHe, rtol=1e-06, atol=1e-09, solver=Kvaerno3(), max_steps=1024)

Compute complete hydrogen recombination history through all phases.

Sequentially computes hydrogen ionization fraction through post-Saha expansion, HYREC-2 EMLA with two-photon processes, and late-time EMLA-only evolution phases.

Parameters:

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

• rtol float, optional

  Relative tolerance for ODE solver (default: 1e-6)

• atol float, optional

  Absolute tolerance for ODE solver (default: 1e-9)

• solver diffrax.Solver, optional

  ODE solver instance (default: Kvaerno3())

• max_steps int, optional

  Maximum solver steps (default: 1024)

Returns:

tuple

(xe_full, lna_full, Tm, lna_Tm) containing ionization fraction evolution, log scale factor grid, matter temperature, and temperature grid

post_Saha_expansion(starting_lna, h, omega_b, omega_cdm, Neff, YHe, threshold=1e-05)

Compute post-Saha expansion phase with two-photon corrections.

Calculates ionization fraction including two-photon processes as perturbative corrections to Saha equilibrium until deviations exceed threshold.

Parameters:

• starting_lna float

  Initial log scale factor

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

• threshold float, optional

  Threshold for deviation from Saha (default: 1e-5)

Returns:

tuple

(xe_output, lna_output) - ionization fraction and log scale factor arrays

solve_TLA(lna0, xe0, Tm0, h, omega_b, omega_cdm, Neff, YHe, rtol=1e-07, atol=1e-09, solver=Kvaerno3(), max_steps=4096)

Solve late-time TLA evolution.

Integrates hydrogen recombination using Peebles TLA in the region beyond where SWIFT corrections are tabulated.

Parameters:

• lna0 float

  Starting log scale factor

• xe0 float

  Initial ionization fraction

• Tm0 float

  Starting matter temperature

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

• rtol float, optional

  Relative tolerance (default: 1e-7)

• atol float, optional

  Absolute tolerance (default: 1e-9)

• solver diffrax.Solver, optional

  ODE solver (default: Tsit5())

• max_steps int, optional

  Maximum steps (default: 4096)

Returns:

tuple

(xe_output, Tm_output, lna_output) - ionization fraction, matter temperature, and log scale factor arrays

solve_emla(lna0, xe0, h, omega_b, omega_cdm, Neff, YHe, rtol=1e-07, atol=1e-09, solver=Tsit5(), max_steps=4096)

Solve late-time EMLA evolution without two-photon processes.

Integrates hydrogen recombination using effective multilevel atom approximation for late times when two-photon processes are negligible.

Parameters:

• lna0 float

  Log scale factor at which initial xe is given

• xe0 float

  Initial ionization fraction

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

• rtol float, optional

  Relative tolerance (default: 1e-7)

• atol float, optional

  Absolute tolerance (default: 1e-9)

• solver diffrax.Solver, optional

  ODE solver (default: Tsit5())

• max_steps int, optional

  Maximum steps (default: 4096)

Returns:

tuple

(xe_output, Tm_output, lna_output) - ionization fraction, matter temperature, and lna arrays

solve_emla_twophoton(lna_axis_init, lna_axis_final, xe0, h, omega_b, omega_cdm, Neff, YHe, rtol=1e-06, atol=1e-09, solver=Kvaerno3(), max_steps=1024)

Solve HYREC-2 EMLA evolution with two-photon processes.

Integrates hydrogen recombination including effective multilevel atom approximation with two-photon transitions and correction functions.

Parameters:

• lna_axis_init float

  Initial log scale factor

• lna_axis_final float

  Final log scale factor

• xe0 float

  Initial ionization fraction

• h float

  Hubble parameter

• omega_b float

  The baryon density Omega_b h^2

• omega_cdm float

  The density of Cold Dark Matter Omega_cdm h^2

• Neff float

  Effective number of neutrinos

• YHe float

  Helium fraction

• rtol float, optional

  Relative tolerance (default: 1e-6)

• atol float, optional

  Absolute tolerance (default: 1e-9)

• solver diffrax.Solver, optional

  ODE solver (default: Kvaerno3())

• max_steps int, optional

  Maximum steps (default: 1024)

Returns:

tuple

(xe_output, lna_output) - ionization fraction and log scale factor arrays

steady_state_equations(xe, H, nH, TCMB, A2s, A2p, B2s, B2p, R2p2s, R2s2p, Delta)

Set up steady-state level population equations.

Constructs matrix equation for hydrogen level populations in steady-state approximation for 2s and 2p levels.

Parameters:

• xe float

  Current ionization fraction

• H float

  Hubble parameter (units: s^{-1})

• nH float

  Hydrogen number density (units: cm^{-3})

• TCMB float

  CMB temperature (units: eV)

• A2s float

  2s recombination coefficient (units: cm^3 s^{-1})

• A2p float

  2p recombination coefficient (units: cm^3 s^{-1})

• B2s float

  2s photoionization coefficient (units: s^{-1})

• B2p float

  2p photoionization coefficient (units: s^{-1})

• R2p2s float

  2p→2s transition rate (units: s^{-1})

• R2s2p float

  2s→2p transition rate (units: s^{-1})

• Delta float

  Correction function value

Returns:

tuple

(T, S) - transition matrix and source vector for level populations

xe_derivative_twophoton(lna, xe, args)

Compute ionization fraction derivative including two-photon processes.

Derivative function for hydrogen ionization fraction evolution including two-photon transitions and correction functions.

Parameters:

• lna float

  Log scale factor

• xe float

  Current ionization fraction

• args tuple

  h, omega_b, omega_cdm, Neff, YHe; the Hubble parameter, the baryon denisty Omega_b h^2, the CDM density Omega_cdm h^2, the effecgive number of neutrinos, and the helium fraction

Returns:

float

Time derivative dxe/dlna (units: dimensionless)

xe_tm_derivative(lna, state, args)

Compute coupled derivatives for ionization fraction and matter temperature.

Derivative function for simultaneous evolution of hydrogen ionization fraction and matter temperature including Compton heating/cooling.

Parameters:

• lna float

  Log scale factor

• state array

  Current state [xe, Tm]

• args tuple

  h, omega_b, omega_cdm, Neff, YHe; the Hubble parameter, the baryon denisty Omega_b h^2, the CDM density Omega_cdm h^2, the effecgive number of neutrinos, and the helium fraction

Returns:

array

Time derivatives [dxe/dlna, dTm/dlna] (units: dimensionless, eV)