diff --git a/CITATION.cff b/CITATION.cff index 665f5e7..f2306c7 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -4,6 +4,7 @@ title: "CMB Decoupling: hydrogen recombination and photon last scattering" type: software authors: - given-names: Sethulekshmi + family-names: Beena affiliation: "University of Hertfordshire" orcid: "https://orcid.org/0009-0004-4767-2733" website: "https://sethulekshmibeena.github.io/my-webpage/" diff --git a/LICENSE b/LICENSE index 164ad48..b2b1d82 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2026 Sethulekshmi +Copyright (c) 2026 Sethulekshmi Beena Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index e4eca45..00d304e 100644 --- a/README.md +++ b/README.md @@ -152,7 +152,7 @@ The full list is at the end of [`docs/theory.md`](docs/theory.md#9-references). ## About me -I'm Sethulekshmi, an MSc Astrophysics student at the University of Hertfordshire. +I'm Sethulekshmi Beena, an MSc Astrophysics student at the University of Hertfordshire. ORCID [0009-0004-4767-2733](https://orcid.org/0009-0004-4767-2733) · [homepage](https://sethulekshmibeena.github.io/my-webpage/) diff --git a/docs/theory.md b/docs/theory.md index c36bc07..29a361e 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -5,8 +5,9 @@ and states each approximation at the point where it is made. Numbers quoted in t were computed with this package for the Planck 2018 cosmology (Section 1); they are reproduced by `scripts/make_figures.py` and the notebook. -Items marked **TODO: verify** are references whose exact location (equation, section or -table number) still needs to be checked against the original source. +Where an equation or value is taken from a paper, the equation or table number is given and +has been checked against the original. Textbooks are cited for standard material without +equation numbers. **Contents** @@ -29,7 +30,8 @@ table number) still needs to be checked against the original source. $\Omega_b h^2 = 0.02237$, $\Omega_c h^2 = 0.1200$. The CMB temperature is $T_0 = 2.7255$ K (Fixsen 2009), the effective number of neutrino species is $N_{\rm eff} = 3.046$ (Mangano et al. 2005), and the primordial helium mass fraction is -$Y_p = 0.2454$ (BBN-consistent value; **TODO: verify** its location in Planck 2018 VI). +$Y_p = 0.2454$. This is the BBN-predicted value Planck uses for all columns of its Table 2 +("posterior mean $Y_P \approx 0.2454$", Table 2 caption; computed with PArthENoPE, Sect. 2.1). **Approximation 1.1: massless neutrinos.** Planck's baseline has one massive neutrino, $\sum m_\nu = 0.06$ eV, which is why it quotes $\Omega_m = 0.3153$. Here all neutrinos are @@ -44,8 +46,8 @@ $\rho_{\rm crit} = 3H_0^2/(8\pi G)$: $$\Omega_\gamma = \frac{4\sigma_{\rm SB} T_0^4/c^3}{3H_0^2/(8\pi G)}, \qquad \Omega_\gamma h^2 = 2.47\times 10^{-5}.$$ Each neutrino species has $7/8$ of the photon energy density per degree of freedom, at a -temperature lower by $(4/11)^{1/3}$ because of $e^\pm$ annihilation (Dodelson & Schmidt 2020, Ch. 2; -**TODO: verify** section). Therefore +temperature lower by $(4/11)^{1/3}$ because of $e^\pm$ annihilation (Kolb & Turner 1990; +Dodelson & Schmidt 2020). Therefore $$\Omega_r = \Omega_\gamma\left[1 + \frac{7}{8}\left(\frac{4}{11}\right)^{4/3} N_{\rm eff}\right] = 9.22\times 10^{-5}.$$ @@ -82,7 +84,7 @@ down to $z \sim 200$; below that the gas cools faster ($T_{\rm m} \propto (1+z)^ For a non-relativistic species $i$ (mass $m_i$, internal degeneracy $g_i$, chemical potential $\mu_i$) in kinetic equilibrium at temperature $T$, integrating the Maxwell–Boltzmann distribution $f = e^{-(E - \mu)/k_BT}$ with $E = m_ic^2 + p^2/2m_i$ over momentum gives -(Kolb & Turner 1990; Dodelson & Schmidt 2020, Ch. 3, **TODO: verify** equation number): +(Kolb & Turner 1990; Dodelson & Schmidt 2020): $$n_i = g_i \int \frac{d^3p}{(2\pi\hbar)^3}\, e^{-(m_ic^2 - \mu_i)/k_BT}\, e^{-p^2/2m_ik_BT} = g_i\left(\frac{m_i k_B T}{2\pi\hbar^2}\right)^{3/2} e^{(\mu_i - m_ic^2)/k_BT}.$$ @@ -117,7 +119,7 @@ $n_p = x_e n_{\rm H}$ and $n_{\rm H\,I} = (1 - x_e)\,n_{\rm H}$. Substituting: $$\boxed{\;\frac{x_e^2}{1 - x_e} = \frac{1}{n_{\rm H}}\left(\frac{m_e k_B T}{2\pi\hbar^2}\right)^{3/2} e^{-E_{\rm ion}/k_BT} \equiv S(z)\;}$$ -This is Dodelson & Schmidt (2020), Ch. 4 (**TODO: verify** equation number). The ionisation +This is the standard Saha equation for hydrogen recombination (e.g. Dodelson & Schmidt 2020). The ionisation energy is the reduced-mass-corrected Rydberg energy, $E_{\rm ion} = R_\infty hc\, m_p/(m_e + m_p) = 13.598$ eV, from CODATA 2018 constants. @@ -146,7 +148,7 @@ $z = 1369.0$, where $T = 3734$ K and $k_BT = 0.322$ eV, far below 13.6 eV. ### 2.5 Helium in Saha equilibrium The same argument applied to ${\rm He\,I} \leftrightarrow {\rm He\,II} + e^-$ and -${\rm He\,II} \leftrightarrow {\rm He\,III} + e^-$ gives (Seager, Sasselov & Scott 2000; **TODO: verify** equation numbers): +${\rm He\,II} \leftrightarrow {\rm He\,III} + e^-$ gives (Seager, Sasselov & Scott 2000, eqs. 80–81): $$\frac{n_{\rm He\,II}\, n_e}{n_{\rm He\,I}} = 4\left(\frac{m_ek_BT}{2\pi\hbar^2}\right)^{3/2} e^{-\chi_{\rm He\,I}/k_BT}, \qquad @@ -155,6 +157,8 @@ $$\frac{n_{\rm He\,II}\, n_e}{n_{\rm He\,I}} = 4\left(\frac{m_ek_BT}{2\pi\hbar^2 The prefactors are $g_e g_{\rm upper}/g_{\rm lower}$: $2\cdot 2/1 = 4$ (He I ground state $^1S_0$ has $g = 1$; He II 1s has $g = 2$), and $2 \cdot 1/2 = 1$ (bare He III has $g = 1$). The ionisation energies are $\chi_{\rm He\,I} = 24.587$ eV and $\chi_{\rm He\,II} = 54.418$ eV (NIST Atomic Spectra Database). +Seager et al. (2000) approximate $m_{\rm He}/m_{\rm H} = 4$ in $f_{\rm He} = Y_p/[4(1 - Y_p)]$. Here the +exact atomic-mass ratio is used, which changes $f_{\rm He}$ by 0.7%. Each species fraction is then a function of $x_e$ alone, and charge neutrality @@ -210,8 +214,8 @@ It is conventional (Seager et al. 1999) to work with $\beta_B = \alpha_B\left(m_ek_BT/2\pi\hbar^2\right)^{3/2} e^{-E_{\rm ion}/4k_BT}$, which is what the code computes. **(b) Two-photon decay 2s → 1s.** The 2s → 1s transition is forbidden for a single photon. It -proceeds by two-photon emission at rate $\Lambda_{2s1s} = 8.22458\ {\rm s^{-1}}$ (value used in RECFAST, -Seager et al. 1999; from Goldman 1989, **TODO: verify**). Including the inverse process +proceeds by two-photon emission at rate $\Lambda_{2s1s} = 8.22458\ {\rm s^{-1}}$ (Goldman 1989, as quoted by +Seager et al. 1999). Including the inverse process (two-photon absorption from 1s), with $g_{2s} = g_{1s}$, the net downward rate is $$\Lambda_{2s1s}\left(n_{2s} - n_{1s}\,e^{-E_{{\rm Ly}\alpha}/k_BT}\right),\qquad E_{{\rm Ly}\alpha} = \tfrac34 E_{\rm ion}.$$ @@ -278,7 +282,9 @@ Converting to redshift with $dt = -dz/[(1+z)H]$: $$\boxed{\;\frac{dx_e}{dz} = \frac{C_r}{(1+z)H(z)}\left[\alpha_B\, n_{\rm H}\, x_e^2 - \beta_B\,(1 - x_e)\, e^{-E_{{\rm Ly}\alpha}/k_BT}\right]\;}$$ -This is the form in Seager et al. (1999) and Dodelson & Schmidt (2020), Ch. 4 (**TODO: verify** equation numbers). +This is eq. (1) of Seager et al. (1999) with $x_e = x_p$. There the Peebles factor is written as +$(1 + K\Lambda_{2s1s}n_{1s})/[1 + K(\Lambda_{2s1s} + \beta_B)n_{1s}]$ with $K = \lambda_\alpha^3/(8\pi H)$; +dividing numerator and denominator by $Kn_{1s}$ gives the form of $C_r$ above, since $1/(Kn_{1s}) = \Lambda_\alpha$. **Detailed-balance check.** Using $E_{\rm ion}/4 + E_{{\rm Ly}\alpha} = E_{\rm ion}$: @@ -291,12 +297,13 @@ directly, to avoid multiplying a large and a small exponential (`tests/test_peeb ### 4.5 The recombination coefficient -The case-B coefficient uses the fit of Péquignot, Petitjean & Boisson (1991), as adopted in RECFAST -(Seager et al. 1999; **TODO: verify** coefficients against the paper): +The case-B coefficient uses the fit of Péquignot, Petitjean & Boisson (1991) to the rates of +Hummer (1994), as adopted in RECFAST (Seager et al. 1999, eq. 3): $$\alpha_B(T) = F \times 10^{-19}\,\frac{a\,t^b}{1 + c\,t^d}\ {\rm m^3\,s^{-1}},\qquad t = \frac{T}{10^4\ {\rm K}},$$ -with $a = 4.309$, $b = -0.6166$, $c = 0.6703$ and $d = 0.5300$. As a check, $\alpha_B(10^4\ {\rm K}) = 2.58\times 10^{-19}\ {\rm m^3\,s^{-1}}$, +with $a = 4.309$, $b = -0.6166$, $c = 0.6703$ and $d = 0.5300$ (identical in the RECFAST module +distributed with CAMB). As a check, $\alpha_B(10^4\ {\rm K}) = 2.58\times 10^{-19}\ {\rm m^3\,s^{-1}}$, consistent with the standard case-B value $2.59\times 10^{-13}\ {\rm cm^3\,s^{-1}}$ (Osterbrock & Ferland 2006). @@ -342,8 +349,7 @@ electrons, with cross-section $\sigma_T = \frac{8\pi}{3}\left(\frac{e^2}{4\pi\va $$\Gamma_T = n_e\,\sigma_T\, c = x_e\, n_{\rm H}\,\sigma_T\, c.$$ Photons are tightly coupled while $\Gamma_T \gg H$, meaning many scatterings per Hubble time. -**Decoupling** is commonly defined by $\Gamma_T(z_{\rm dec}) = H(z_{\rm dec})$ (Dodelson & Schmidt 2020, -Ch. 4; **TODO: verify**). +**Decoupling** is commonly defined by $\Gamma_T(z_{\rm dec}) = H(z_{\rm dec})$ (e.g. Dodelson & Schmidt 2020). | Model | $z_{\rm dec}$ ($\Gamma_T = H$) | |---|---| @@ -399,6 +405,11 @@ Three common definitions give slightly different numbers for the same $x_e(z)$: | Peak of $g(z)$ per unit redshift | 1078.80 | 1078.37 | | Peak of the visibility per unit conformal time | 1089.11 | 1088.78 | +For context, Planck's own fit gives $z_* = 1089.92 \pm 0.25$ for this data combination +(Planck 2018 VI, Table 2). That number is not a like-for-like comparison: it includes a massive +neutrino and a recombination history that incorporates the corrections computed by CosmoRec and +HyRec (Planck 2018 VI, Sect. 7.7). The matched CAMB run above is the clean test of this code. + The last row needs one more step. Conformal time obeys $d\eta = c\,dt/a = -c\,dz/H$, so the visibility per unit conformal time is @@ -461,11 +472,13 @@ therefore isolate the recombination physics. * Chluba, J. & Thomas, R. M. 2011, *MNRAS* **412**, 748 — CosmoRec. * Dodelson, S. & Schmidt, F. 2020, *Modern Cosmology*, 2nd ed. (Academic Press). * Fixsen, D. J. 2009, *ApJ* **707**, 916 — $T_0 = 2.7255$ K. -* Goldman, S. P. 1989, *Phys. Rev. A* **40**, 1185 — two-photon decay rate (**TODO: verify**). +* Goldman, S. P. 1989, *Phys. Rev. A* **40**, 1185 — two-photon decay rate. +* Hummer, D. G. 1994, *MNRAS* **268**, 109 — hydrogen recombination coefficients. * Kolb, E. W. & Turner, M. S. 1990, *The Early Universe* (Addison-Wesley). * Lewis, A., Challinor, A. & Lasenby, A. 2000, *ApJ* **538**, 473 — CAMB. * Mangano, G. et al. 2005, *Nucl. Phys. B* **729**, 221 — $N_{\rm eff} = 3.046$. * NIST Atomic Spectra Database (Kramida, A. et al.), https://physics.nist.gov/asd — ionisation energies. +* NIST Atomic Weights and Isotopic Compositions, https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses — masses of ¹H and ⁴He. * Osterbrock, D. E. & Ferland, G. J. 2006, *Astrophysics of Gaseous Nebulae and Active Galactic Nuclei*, 2nd ed. * Peebles, P. J. E. 1968, *ApJ* **153**, 1. * Péquignot, D., Petitjean, P. & Boisson, C. 1991, *A&A* **251**, 680 — case-B fit. @@ -475,6 +488,5 @@ therefore isolate the recombination physics. * Seager, S., Sasselov, D. D. & Scott, D. 2000, *ApJS* **128**, 407. * Sobolev, V. V. 1960, *Moving Envelopes of Stars* (Harvard University Press). * Switzer, E. R. & Hirata, C. M. 2008, *Phys. Rev. D* **77**, 083006 — helium recombination. -* Wang, M. et al. 2017, *Chinese Phys. C* **41**, 030003 — AME2016 atomic masses (**TODO: verify**). * Wong, W. Y., Moss, A. & Scott, D. 2008, *MNRAS* **386**, 1023 — RECFAST updates. * Zel'dovich, Ya. B., Kurt, V. G. & Sunyaev, R. A. 1968, *Zh. Eksp. Teor. Fiz.* **55**, 278. diff --git a/pyproject.toml b/pyproject.toml index ea0c74f..c5ded50 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,7 +9,7 @@ description = "Hydrogen recombination (Saha, Peebles) and CMB photon decoupling, readme = "README.md" requires-python = ">=3.10" license = { text = "MIT" } -authors = [{ name = "Sethulekshmi" }] +authors = [{ name = "Sethulekshmi Beena" }] dependencies = ["numpy>=1.24", "scipy>=1.11", "astropy>=5.3", "matplotlib>=3.7"] [project.optional-dependencies] diff --git a/src/cmb_decoupling/constants.py b/src/cmb_decoupling/constants.py index 73a3ff9..151e12b 100644 --- a/src/cmb_decoupling/constants.py +++ b/src/cmb_decoupling/constants.py @@ -37,12 +37,11 @@ E_ION_H_N2: float = 0.25 * E_ION_H #: Lyman-alpha wavelength, lambda = h c / E_Lya [m]. LAMBDA_LYA: float = H_PLANCK * C_LIGHT / E_LYA -#: Two-photon 2s -> 1s decay rate [s^-1]. Value used in RECFAST -#: (Seager, Sasselov & Scott 1999, 2000), from Goldman (1989). -#: TODO: verify the Goldman (1989) citation details. +#: Two-photon 2s -> 1s decay rate [s^-1]: Goldman (1989, Phys. Rev. A 40, 1185), +#: as quoted by Seager, Sasselov & Scott (1999). LAMBDA_2S1S: float = 8.22458 -#: Mass of the 1H atom, 1.00782503223 u (Wang et al. 2017, AME2016). -#: TODO: verify AME citation. +#: Mass of the 1H atom, 1.00782503223 u +#: (NIST Atomic Weights and Isotopic Compositions database). M_H_ATOM: float = 1.00782503223 * AMU # --- Helium atomic data ----------------------------------------------------- @@ -50,5 +49,6 @@ E_ION_HEI: float = 24.587387 * EV #: Ionisation energy of He II, 54.417760 eV (NIST Atomic Spectra Database). E_ION_HEII: float = 54.417760 * EV -#: Mass of the 4He atom, 4.00260325413 u (Wang et al. 2017, AME2016). +#: Mass of the 4He atom, 4.00260325413 u +#: (NIST Atomic Weights and Isotopic Compositions database). M_HE4_ATOM: float = 4.00260325413 * AMU diff --git a/src/cmb_decoupling/cosmology.py b/src/cmb_decoupling/cosmology.py index ec976b0..eca8d7d 100644 --- a/src/cmb_decoupling/cosmology.py +++ b/src/cmb_decoupling/cosmology.py @@ -20,7 +20,7 @@ #: Ratio of neutrino to photon temperature, (4/11)^(1/3), and the fermion #: factor 7/8, giving rho_nu / rho_gamma = N_eff * 7/8 * (4/11)^(4/3) -#: (Dodelson & Schmidt 2020, Modern Cosmology, 2nd ed., Ch. 2). +#: (Kolb & Turner 1990; Dodelson & Schmidt 2020, Modern Cosmology, 2nd ed.). NEUTRINO_DENSITY_FACTOR: float = 7.0 / 8.0 * (4.0 / 11.0) ** (4.0 / 3.0) @@ -121,8 +121,8 @@ def n_H(self, z: ArrayLike) -> NDArray[np.float64]: #: Planck 2018 VI, Table 2, TT,TE,EE+lowE+lensing (mean values). -#: Y_p is the BBN-consistent helium fraction quoted for this fit. -#: TODO: verify the Y_p = 0.2454 value and its location in Planck 2018 VI. +#: Y_p is the BBN prediction Planck uses for all Table 2 columns +#: ("posterior mean Y_P ~ 0.2454", Table 2 caption; PArthENoPE, Sect. 2.1). #: N_eff = 3.046 (standard-model value used by Planck 2018; Mangano et al. 2005). PLANCK2018 = Cosmology( H0_km_s_Mpc=67.36, diff --git a/src/cmb_decoupling/peebles.py b/src/cmb_decoupling/peebles.py index bf0fe2e..b599662 100644 --- a/src/cmb_decoupling/peebles.py +++ b/src/cmb_decoupling/peebles.py @@ -1,7 +1,7 @@ """Peebles three-level-atom model of hydrogen recombination. The free-electron fraction obeys (Peebles 1968; notation of Seager, -Sasselov & Scott 1999 and Dodelson & Schmidt 2020, eq. 4.39) +Sasselov & Scott 1999, eq. 1) dx_e/dz = C_r / [H(z)(1+z)] * [alpha_B n_H x_e^2 - beta_B (1 - x_e) exp(-E_Lya / k T)], @@ -37,9 +37,9 @@ from .cosmology import Cosmology from .saha import thermal_density, x_e_saha_full, x_e_saha_hydrogen -# Case-B recombination coefficient fit (Pequignot et al. 1991, as used in -# Seager et al. 1999): alpha_B = 1e-19 a t^b / (1 + c t^d) m^3 s^-1, -# with t = T / 1e4 K. TODO: verify coefficients against Seager et al. (1999). +# Case-B recombination coefficient fit (Pequignot et al. 1991), as given in +# Seager et al. (1999), eq. 3: alpha_B = 1e-19 a t^b / (1 + c t^d) m^3 s^-1, +# with t = T / 1e4 K. ALPHA_B_NORM: float = 1.0e-19 # [m^3 s^-1] ALPHA_B_A: float = 4.309 ALPHA_B_B: float = -0.6166 diff --git a/src/cmb_decoupling/saha.py b/src/cmb_decoupling/saha.py index 429dcd2..443815d 100644 --- a/src/cmb_decoupling/saha.py +++ b/src/cmb_decoupling/saha.py @@ -2,10 +2,10 @@ References ---------- -Dodelson & Schmidt (2020), Modern Cosmology, 2nd ed., Sec. 4.2 (hydrogen). -Seager, Sasselov & Scott (2000), ApJS 128, 407 -- helium Saha relations, -including the statistical-weight factors 4 (He I <-> He II) and 1 -(He II <-> He III). TODO: verify equation numbers. +Dodelson & Schmidt (2020), Modern Cosmology, 2nd ed. (hydrogen). +Seager, Sasselov & Scott (2000), ApJS 128, 407, eqs. 80-81 -- helium Saha +relations, including the statistical-weight factors 4 (He I <-> He II) and 1 +(He II <-> He III). """ from __future__ import annotations diff --git a/src/cmb_decoupling/visibility.py b/src/cmb_decoupling/visibility.py index 4f3e856..e678cea 100644 --- a/src/cmb_decoupling/visibility.py +++ b/src/cmb_decoupling/visibility.py @@ -2,8 +2,7 @@ References ---------- -Dodelson & Schmidt (2020), Modern Cosmology, 2nd ed., Ch. 4 (TODO: verify -section numbers for the optical depth and visibility function). +Dodelson & Schmidt (2020), Modern Cosmology, 2nd ed. Planck Collaboration (2020), A&A 641, A6 -- definition of z_* as tau(z_*) = 1. """