# 21cm signal from Dark Ages collapsing halos with detailed molecular cooling treatment

Hugo Plombat, Denis Puy

► To cite this version:

Hugo Plombat, Denis Puy. 21cm signal from Dark Ages collapsing halos with detailed molecular cooling treatment. Astronomy & Astrophysics - A&A, 2025, 698, pp.A224. 10.1051/0004-6361/202450380. hal-04557864

HAL Id: hal-04557864

https://hal.science/hal-04557864v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 21 cm signal from dark-age collapsing halos with a detailed molecular cooling treatment

Hugo Plombat*, Denis Puy*

Laboratoire Univers et Particules de Montpellier (LUPM), CNRS & IN2P3 et Université de Montpellier (UMR-5299),
Place Eugène Bataillon, F-34095 Montpellier Cedex 05, France

Received 15 April 2024 / Accepted 23 April 2025

## ABSTRACT

**Context.** To understand the formation of the first stars, a detailed description of the thermal and chemical processes in collapsing gas clouds is essential. Molecular cooling, particularly via H₂, plays a significant role in triggering thermal instabilities that lead to star formation. The 21 cm hydrogen line serves as a potential probe of the first collapsing structures during the dark ages of the early Universe and it is affected by the gas temperature evolution.

**Aims.** We aim to investigate the molecular cooling in the gas halos prior to the formation of the first stars, with a particular focus on how the H₂ cooling affects the gas temperature. Additionally, we explore the sensitivity of the 21 cm hydrogen line to these cooling processes during the collapse of the first overdense regions.

**Results.** We introduce the CHEFAST code, which tracks the evolution of chemical abundances and computes the 21 cm neutral hydrogen signal in collapsing halos. Our results show that molecular cooling significantly affects the gas temperature inside collapsing clouds of mass ranging from 10⁶ to 10⁹ M☉, influencing the 21 cm signal. The signal exhibits an emission feature that is distinct from the one predicted in simpler expansion models.

**Conclusions.** The 21 cm brightness temperature inside collapsing clouds displays an emission feature driven by molecular cooling, closely mirroring the gas temperature evolution. This makes the dark-age 21 cm signal a promising probe for studying the thermal processes and structure formation in the early Universe.

**Key words.** molecular processes – dark ages, reionization, first stars

## 1. Introduction

The chemistry of the early Universe has been the target of large studies aimed at determining the abundances of elements and dedicated to the calculation of reaction rates that govern chemical processes, which are contingent on the prevailing environmental conditions. Such investigations explore the primordial species produced by the first nuclear reactions during standard Big Bang nucleosynthesis (SBBN), (see e.g., Workman 2022; Signore & Puy 2009; Cyburt et al. 2016; Pitrou et al. 2021; Schöneberg 2024), up to the interstellar medium (ISM), hosting a variety of heavier elements forged within stars (Tielens 2005; Herbst 1995).

From the SBBN, the Universe was completely ionized, as the temperature of the baryon-photon fluid, maintained by interactions via Compton scattering between free electrons and photons, was too high for any neutral atom to form. However, as the Universe expanded, the effectiveness of baryon-radiation coupling diminished, allowing ion-electron reactions to dominate. Consequently, this led to the progressive recombination of ionized elemental forms into their neutral counterparts (Zeldovich et al. 1968; Peebles 1968). Simultaneously, despite their much lower reaction rates due to the necessarily high amount of energy for these processes to happen, molecules also began to form (Puy et al. 1993).

Species abundance estimation is a central problem of any coupled reaction network. One of the first detailed models of molecular synthesis were proposed by Solomon & Klemperer (1972), who computed molecular abundances in diffuse interstellar clouds, and by Herbst & Klemperer (1973) who studied the formation of molecules in dense dark clouds. The article of Prasad & Huntress (1980), which gave a comprehensive library consisting of over 1400 reactions for 137 species, presented a first model for gas phase chemistry in interstellar clouds. The UMIST Database (Le Teuff et al. 2000; McElroy et al. 2013), contains the updated rate coefficient, temperature ranges and temperature dependence of 6173 gas-phase reactions important in astrophysical environments. However at early epochs, when a total absence of dust grains appears justified, the chemistry is different from the typical interstellar medium chemistry. From the exhaustive work of Bates (1951) on the radiative association rate, Takayanagi & Nishimura (1960), Hirasawa (1969), Takeda et al. (1969), Saslaw & Zipoy (1967) and Matsuda et al. (1969) pointed out the way of H₂ formation through H₂⁺ and H⁻ without grains in the early Universe. Several chemical networks studies including the primordial molecules (such as H₂, HD, and LiH) and ions were carried out by Lepp & Shull (1984), Latter & Black (1991), Puy et al. (1993), Galli & Palla (1998), Stancil & Dalgarno (1998), Lepp et al. (2002), Pfenniger & Puy (2003), Hirata & Padmanabhan (2006), Glover & Abel (2008), Signore & Puy (2009), Coppola et al. (2011), Gay et al. (2011), Coppola et al. (2012, 2013).

The determination of the molecular abundances is essential in the study of the formation of structures. Indeed,

<footnote>* Corresponding authors: hugo.plombat@umontpellier.fr; denis.puy@umontpellier.fr</footnote>
the excitation of the rotational and vibrational levels of the molecules inside collapsing matter overdensities gives rise to a molecular cooling mechanism, which can induce the appearance of instabilities of a thermal nature, highlighted for the first time in Field (1965) and later in Yoshii & Sabano (1980). This excitation is necessary in order to trigger the formation of the first stars Tegmark et al. (1997), Puy & Signore (1997), Yoshida et al. (2006), Trenti & Stiavelli (2009), Galli & Palla (2013), Bromm (2013), Glover (2012), Ripamonti & Abel (2004). Molecular cooling has been extensively studied for numerous molecular species, the most important being H₂, which despite the absence of a dipole moment is the stronger coolant due to its high abundance (Abel et al. 1997; Galli & Palla 1998; Stancil et al. 1998; Glover & Abel 2008). HD is also expected to be an efficient coolant due to its non-zero dipole moment, but it has been shown that it only dominates the molecular cooling at very low temperatures (~150 K; Glover & Abel 2008), which are typically not reached in primordial gas during the collapse of the first overdense regions (Bromm et al. 2002). Other early molecules, such as H₂⁺ (Yoshida et al. 2007; Glover & Savin 2009), H₃⁺ (Glover & Savin 2006, 2009), and LiH (Stancil et al. 1996; Mizusawa et al. 2005), have minor cooling contributions in the context of setting up the conditions for the formation of the first stars.

Probing the formation of the first structures preceding star formation, during the so-called dark ages, is very challenging because almost no signals are emitted from this epoch. Baryonic matter predominantly existed in a neutral state, as the Compton coupling with radiation waned with cosmic expansion. The 21 cm hydrogen line, arising from hyperfine spin-flip transitions in neutral hydrogen atoms, offers a promising avenue for probing this epoch.

Many future experiments are planned in order to measure the 21 cm signal during cosmic dawn and the Epoch of Reionization, such as Experiment to Detect the Global rEionization Signature (EDGES; Bowman et al. 2018), Shaped Antenna measurement of the background RAdio Spectrum (SARAS; Nambissan et al. 2021), Large-Aperture Experiment to Detect the Dark Ages (LEDA; Price et al. 2018), the Radio Experiment for the Analysis of Cosmic Hydrogen (REACH; de Lera Acedo et al. 2022), and Probing Radio Intensity at high-Z from Marion (PRizM; Philip et al. 2019). These studies are aimed at the observation of the 21 cm global signal, whereas ground based interferometers such as the Hydrogen Epoch of Reionization Array (HERA; DeBoer 2017), Low Frequency Array (LOFAR; van Haarlem et al. 2013; Patil et al. 2016), Murchison Widefield Array (MWA; Bowman et al. 2013; Tingay et al. 2013; Jacobs et al. 2016), the Giant Metre-wave Radio Telescope (GMRT; Swarup et al. 1991), the Precision Array for Probing the Epoch of Reionization (PAPER; Parsons et al. 2010), and the future Square Kilometer Array (SKA; Koopmans et al. 2015) are measuring or planning to measure the spatial fluctuations in the signal.

However, the observation of the redshifted 21 cm power spectrum and global signal from the dark ages is extremely challenging and cannot be done from Earth-based experiments, because of the ionosphere which acts as a highly opaque foreground for the frequencies of interest, as well as the radio frequency interferences of human origin. The greatest hopes to fathom the dark ages are placed on space missions, either in orbit around the Moon, with the satellites and nanosatellites. We also refer to Burns et al. (2021a) and Fialkov et al. (2024), DSL (Shi et al. 2022), NCDL (Bentum et al. 2020), and

CosmoCube (Artuc & de Lera Acedo 2024), or directly on the Moon's surface with ROLSES (Burns et al. 2021b), LuSEE-Night (Bale et al. 2023), the Lunar Crater Radio Telescope (Goel et al. 2022), FarView (Polidan et al. 2024), the Astronomical Lunar Observatory (ALO; Klein-Wolt et al. 2024), FAR-SIDE (Burns et al. 2021b), and LARAF (Chen et al. 2024). A mission with similar objectives, known as Primordial Hydrogen Observation with Nanosatellites Explorer (PHONE), is also in co-development by the Laboratoire Univers et Particules de Montpellier (LUPM) and the Spatial Center, dedicated to the development of nanosatellites, at the University of Montpellier.

Our work takes part in the studies of the 21 cm signal from minihalos. However, unlike aforementioned works, which concentrate on estimating the signal from halos that have already been virialized, we focus on the earlier stages, from the turnaround point of an overdense region, to the end of its collapse. In particular, we aim to highlight the impact of molecular cooling on the 21 cm signal arising from the halos.

The originality of this work lies in coupling the evolution of interactions between atomic and molecular ions, neutral atoms, and molecules with the dynamic and thermal evolution of a medium undergoing gravitational collapse. While this context has been extensively studied in the past to understand the fragmentation of a molecular cloud by the process of thermal instability, in this article, we integrated the calculation of the emissivity of the excitation process of neutral hydrogen producing a hyperfine transition and followed "step by step" this emissivity during a collapse. This coupled approach enables us to predict the emissivity of the 21 cm line at high redshifts. This work is part of a large body of research into the production of the 21 cm line during the cosmological dark ages.

Numerous studies have been carried out to characterize the first structures with 21 cm hydrogen line, before the formation of the first stars. Iliev et al. (2002) proposed a model of 21 cm emission in minihalos from collisional excitation, before reionization. Meiksin (2011) followed spherical collapsing density perturbations, estimating the 21 cm signature arising from the fully collapsed minihalos accounting for the effects of molecular hydrogen formation. A statistical 21 cm signal from the halo has been provided, demonstrating how this signature could distinguish predictions of the small scales of the 21 cm power spectrum. Furugori et al. (2020) calculated the 21 cm signal emitted by ultracompact minihalos formed at high redshift. Novosyadlyj et al. (2020) calculated the 21 cm emission in halos between 10⁶ and 10¹⁰ M☉, virialized between z ~ 10–50. Recently, Novosyadlyj et al. (2025) estimated the global signal in the redshifted hyperfine structure line 21 cm of hydrogen atoms formed during the dark ages and cosmic dawn. They show that the profile of this line crucially depends on the temperature and ionization of baryonic matter as well as the spectral energy distribution of radiation from the first sources.

In this paper, we make use of the code CHEMFAST initiated by Puy et al. (1993), Pfenniger & Puy (2003) and Vonlanthen et al. (2009). The original purpose of this code is to compute the abundances of atomic and molecular species in the context of cosmological expansion of the universe. We incorporated the computation of the global 21 cm line signal arising from collisional coupling during the dark ages. Then, by changing the equations of dynamics in the code, we followed a collapse scenario of a 10⁶–10⁹ M☉ overdense region, modeling the region with a simple homogeneous spherical model, including pressure effects (Lahav 1986; Puy & Signore 1996). We emphasized the cooling effects on the gas temperature, arising from molecular excitation during the collapse. We additionally computed
the 21 cm signal from the forming protocloud, highlighting its distinct features compared to the expected sky-averaged 21 cm global signal.

This paper is organized as follows. In Sect. 2, we start by developing the reaction network and equations of dynamics in the context of expanding homogeneous Universe. We present the CHEMFAST code computation of atomic and molecular abundances. In Sect. 3, we introduce the physics of the 21 cm global signal during the dark ages. In Sect. 4, we describe the dynamics of the spherical collapse model and, finally, we present our analysis of the thermal and 21 cm signal evolution in collapsing halos in Sect. 5. We present our conclusions and outline our future research directions in Sect. 6. Throughout this paper, we adopted the ΛCDM cosmology, using the parameters values from Planck2018 Planck Collaboration VI (2020) TT,TE,EE+lowE+lensing best-fit.

## 2. CHEMFAST in homogeneous universe

In this section, we consider the well-studied case of the background evolution of chemistry during the recombination and the dark ages. Our aim is to present how the code works using this example, and to introduce the concepts that will also be employed in this study.

In CHEMFAST, we follow a system of coupled differential equations in order to track the evolution of species abundances. The SBBN model for the Universe predicts the nuclei abundances of mainly hydrogen, helium, and lithium as well as their isotopes. The chemistry of the early Universe is the chemistry of these light elements and their respective isotopic forms. The ongoing physical reactions are numerous after the nucleus recombination. These reactions fall into three distinct categories:

- Collisional reactions $\xi + \xi' \leftrightarrow \xi_1 + \xi_2$ involve processes such as association, associative detachment, mutual neutralization, charge exchange, and reverse reactions;

- Electronic reactions $\xi + e^- \leftrightarrow \xi_1 + \xi_2$ encompass recombination, radiative attachment, and dissociative attachment processes;

- Photo-processes $\xi + \gamma \leftrightarrow \xi_1 + \xi_2$ include dissociation, detachment, ionization, and radiative association reactions, where $\gamma$ is a photon.

The chemical kinetic of a reactant, $\xi$, which leads to the products $\xi_1$ and $\xi_2$, imposes the following evolution of the average number density $\bar{n}_{\xi}$ (number of species $\xi$ per cm$^3$ in the homogeneous Universe, i.e., the mean density):

$$ \left(\frac{d\bar{n}_{\xi}}{dt}\right)_{\text{chem}} = \sum_{\xi_1 \xi_2} k_{\xi_1 \xi_2} \bar{n}_{\xi_1} \bar{n}_{\xi_2} - \sum_{\xi'} k_{\xi \xi'} \bar{n}_{\xi} \bar{n}_{\xi'}. \quad (1) $$

Here, $k_{\xi_1 \xi_2}$ represents the reaction rate of the $\xi$-formation process from $\xi_1$ and $\xi_2$, and $k_{\xi \xi'}$ denotes the reaction rate of $\xi$-destruction by collision with the reactant, $\xi'$. The typical rate of collisional reactions, denoted as $k_{\text{coll}}$ (in cm$^3$ s$^{-1}$), is calculated by averaging cross-sections $\sigma(E)$ over a Maxwellian velocity distribution at temperature, $T_k$:

$$ k_{\xi, \xi'} = \sqrt{\frac{8}{\pi m_r}} [k_B T_k]^{-3/2} \int_0^\infty \sigma(E) e^{-E/k_B T_k} E \, dE, \quad (2) $$

where E is the collision energy, $m_r$ the reduced mass of the collisional system ($\xi, \xi'$), $T_k$ the kinetic baryon temperature, and $k_B$ the Boltzmann constant.

The radiative rate coefficient, $k_{\text{rad}}$ (in s$^{-1}$), which depends on the radiative cross-section, $\sigma_{\text{rad}}$, at the frequency, $\nu$, and on the

CMB, is defined by

$$ k_{\text{rad}} = \frac{8\pi}{c^2} \int_{v_{\text{th}}}^{\infty} \sigma_{\text{rad}}(\nu) \frac{\nu^2 \, d\nu}{\exp(h_P \nu / k_B T_r) - 1}, \quad (3) $$

where $v_{\text{th}}$ is the threshold frequency above which the radiative process is possible, $T_r$ is the radiation temperature, c is the speed of light, and $h_P$ is the Planck's constant.

### 2.1. Reaction network

The chemical network, namely, every reaction between atoms, ions, and molecules that we take into account in the species abundances computation, is described in Table 1, along with a summary of the reactions and their references. The primordial chemistry of heavier element is also implemented in CHEMFAST (Li, C, N, O, and F) Vonlanthen et al. (2009), but we ignored it as it is irrelevant for this work. The fitting formula of the reaction rates were mainly taken from the compilation given by Schleicher et al. (2008), except for the hydrogen and deuterium recombination, which were calculated from the rate given by Novosadylyj et al. (2022).

### 2.2. Density and temperature evolution

The complete system of differential equations allows us to calculate the abundances of species. We also followed the thermal evolution of the Universe (gas and radiation temperature), which plays a role in the reaction rates.

**Baryon number density.** The mean baryon number density, $\bar{n}_b$, is expressed by

$$ \bar{n}_b = \sum_{\xi} \bar{n}_{\xi}, \quad (4) $$

summed over $\xi$, the species involved in the network. Thus, the evolution of the mean number density is given by

$$ \frac{d\bar{n}_b}{dt} = -3H(t)\bar{n}_b + \sum_{\xi} \left(\frac{d\bar{n}_{\xi}}{dt}\right)_{\text{chem}}, \quad (5) $$

where $H(t) = \dot{a}/a$ is the Hubble parameter. The second term on the right-hand side of this equation accounts for the variation of the number of particles due to the chemical reactions, as detailed in Eq. (1).

**Species densities.** For each species, $\xi$, we have

$$ \frac{d\bar{n}_{\xi}}{dt} = -3H(t)\bar{n}_{\xi} + \left(\frac{d\bar{n}_{\xi}}{dt}\right)_{\text{chem}}, \quad (6) $$

where $\bar{n}_{\xi}$ is the number density of the species, $\xi$. The last term of the second member of this equation expresses the contribution of the chemical kinetics (see Eq. (1)).

**Radiation temperature.** In an expanding universe, the energy density of radiation evolves as $\rho_r \propto (1+z)^4$. The CMB can be likened to a blackbody, and Stefan-Boltzmann's law can be used to relate its temperature to its energy density: $\rho_r \propto T_\gamma^4$. This leads to the following equation $T_\gamma(z) = T_{\gamma,0}(1+z)$, where $T_{\gamma,0} = 2.726$ K Mather et al. (1994) is the blackbody temperature of the CMB. We can then deduce the evolution law:

$$ \frac{dT_\gamma}{dt} = -H(t)T_\gamma. \quad (7) $$
<table><caption>**Table 1.** CHEMFAST chemical network.</caption><thead><tr><th></th><th>Reaction</th><th>Reference</th><th></th><th>Reaction</th><th>Reference</th></tr></thead><tbody><tr><td>(H1)</td><td>H<sup>+</sup> + e<sup>-</sup> → H + γ</td><td>NKMS22</td><td>(H2)</td><td>H + γ → H<sup>+</sup>e<sup>-</sup></td><td>AAZ97</td></tr><tr><td>(H3)</td><td>H + e<sup>-</sup> → H<sup>-</sup> + γ</td><td>GP98</td><td>(H4)</td><td>H<sup>-</sup> + γ → H + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(H5)</td><td>H<sup>-</sup> + H → H<sub>2</sub> + e<sup>-</sup></td><td>GP98</td><td>(H6)</td><td>H<sup>-</sup> + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(H7)</td><td>H<sup>-</sup> + H<sup>+</sup> → H + H</td><td>LSD02</td><td>(H8)</td><td>H + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + γ</td><td>GP98</td></tr><tr><td>(H9)</td><td>H<sub>2</sub><sup>+</sup> + γ → H + H<sup>+</sup></td><td>GP98</td><td>(H10)</td><td>H<sub>2</sub><sup>+</sup> + H → H<sub>2</sub> + H<sup>+</sup></td><td>GP98</td></tr><tr><td>(H11)</td><td>H<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + H</td><td>GP98</td><td>(H12)</td><td>H<sub>2</sub><sup>+</sup> + γ → H<sup>+</sup> + H<sup>+</sup> + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(H13)</td><td>H<sub>2</sub> + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + H</td><td>SKHS04</td><td>(H14)</td><td>H<sub>2</sub> + e<sup>-</sup> → H + H<sup>-</sup></td><td>CCDL07</td></tr><tr><td>(H15)</td><td>H<sub>2</sub> + e<sup>-</sup> → H + H + e<sup>-</sup></td><td>TT02</td><td>(H16)</td><td>H<sub>2</sub> + γ → H<sub>2</sub><sup>+</sup> + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(H17)</td><td>H + H + H → H<sub>2</sub> + H</td><td>GL08</td><td>(H18)</td><td>H<sup>-</sup> + H<sub>2</sub><sup>+</sup> → H<sub>2</sub> + H</td><td>AAZN97</td></tr><tr><td>(D1)</td><td>D<sup>+</sup> + e<sup>-</sup> → D + γ</td><td>NKMS22</td><td>(D2)</td><td>D + γ → D<sup>+</sup> + e<sup>-</sup></td><td>AAZ97</td></tr><tr><td>(D3)</td><td>D + H<sup>+</sup> → D<sup>+</sup> + H</td><td>Savin (2002)</td><td>(D4)</td><td>D<sup>+</sup> + H → D + H<sup>+</sup></td><td>Savin (2002)</td></tr><tr><td>(D5)</td><td>H + D → HD + γ</td><td>GP98</td><td>(D6)</td><td>D + H<sub>2</sub> → H + HD</td><td>GP02</td></tr><tr><td>(D7)</td><td>HD<sup>+</sup> + H → H<sup>+</sup> + HD</td><td>SLD98</td><td>(D8)</td><td>D<sup>+</sup> + H<sub>2</sub> → H<sup>+</sup> + HD</td><td>GP02</td></tr><tr><td>(D9)</td><td>HD + H → H<sub>2</sub> + D</td><td>GP02</td><td>(D10)</td><td>HD + H<sup>+</sup> → H<sub>2</sub> + D<sup>+</sup></td><td>GP02</td></tr><tr><td>(D11)</td><td>D + H<sup>+</sup> → HD<sup>+</sup> + γ</td><td>GP98</td><td>(D12)</td><td>H + D<sup>+</sup> → HD<sup>+</sup> + γ</td><td>GP98</td></tr><tr><td>(D13)</td><td>HD<sup>+</sup> + γ → H + D<sup>+</sup></td><td>SGP+08</td><td>(D14)</td><td>HD<sup>+</sup> + γ → H<sup>+</sup> + D</td><td>SGP+08</td></tr><tr><td>(D15)</td><td>HD<sup>+</sup> + e<sup>-</sup> → H + D</td><td>SLD98</td><td></td><td></td><td></td></tr><tr><td>(He1)</td><td>He<sup>2+</sup> + e<sup>-</sup> → He<sup>+</sup> + γ</td><td>GP98</td><td>(He2)</td><td>He<sup>+</sup> + γ → He<sup>2+</sup> + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(He3)</td><td>He<sup>+</sup> + e<sup>-</sup> → He + γ</td><td>GP98</td><td>(He4)</td><td>He + γ → He<sup>+</sup> + e<sup>-</sup></td><td>GP98</td></tr><tr><td>(He5)</td><td>He + H<sup>+</sup> → He<sup>+</sup> + H</td><td>GP98 &amp; GJ07</td><td>(He6)</td><td>He<sup>+</sup> + H → He + H<sup>+</sup></td><td>ZDKL89</td></tr><tr><td>(He7)</td><td>He<sup>+</sup> + H → HeH<sup>+</sup> + γ<sup>(1)</sup></td><td>SLD98</td><td>(He8)</td><td>He<sup>+</sup> + H + γ → HeH<sup>+</sup> + γ<sup>(2)</sup></td><td>ZSD98</td></tr><tr><td>(He9)</td><td>He + H<sub>2</sub><sup>+</sup> → HeH<sup>+</sup> + H</td><td>GP98</td><td>(He10)</td><td>He + H<sup>+</sup> → HeH<sup>+</sup> + γ</td><td>SLD98</td></tr><tr><td>(He11)</td><td>HeH<sup>+</sup> + e<sup>-</sup> → He + H<sub>2</sub><sup>+</sup></td><td>LJB95</td><td>(He12)</td><td>HeH<sup>+</sup> + H<sub>2</sub> → He + H</td><td>SLD98</td></tr><tr><td>(He13)</td><td>HeH<sup>+</sup> + γ → He + H<sup>+</sup></td><td>JSK95</td><td>(He14)</td><td>HeH<sup>+</sup> + γ → He<sup>+</sup> + H</td><td>GP98</td></tr></tbody></table>

**Notes.** AAZ97: Abel et al. (1997), CCDL07: Capitelli et al. (2007), GLO08: Glover & Abel (2008), GP98: Galli & Palla (1998), GP02: Galli & Palla (2002), JSK95: Juřek et al. (1995), LJB95: Linder et al. (1995), LSD02: Lepp et al. (2002), NKMS22: Novosadylj et al. (2022), SGP+08: Schleicher et al. (2008), SKHS04: Savin et al. (2004), SLD98: 1998, ZDKL89: Zygelman et al. (1998), ZSD98: Zygelman et al. (1998).
(1) Radiative association; (2) stimulated radiative association.

**Gas temperature.** The evolution of gas kinetic temperature $T_k$ is driven by three thermal channels:

$$ \frac{dT_k}{dt} = \frac{2}{3\bar{n}_b k_B} (\Psi_{\text{adia}} + \Psi_{\text{com}} + \Psi_{\text{mol}}); \qquad (8) $$

$\Psi_{\text{adia}}$ is the adiabatic cooling due to the expansion is given by

$$ \Psi_{\text{adia}} = -3H(t)n_b k_B T_k; \qquad (9) $$

$\Psi_{\text{com}}$ is the energy transfer between matter and radiation via Compton scattering of CMB photons on free electrons (Weymann 1965; Peebles 1968);

$$ \Psi_{\text{com}} = 4\bar{n}_e \sigma_T a_r k_B T_r^4 \frac{T_r - T_k}{m_e c}. \qquad (10) $$

Here, $\sigma_T$ is the Thomson cross section, $a_r$ is the radiation constant, $m_e$ the electron mass, and $\bar{n}_e$ the mean electron number density. $\Psi_{\text{mol}}$ defines the energy transfer via excitation and de-excitation of molecular rotational transition. This process is discussed in Sect. 2.4.

Finally, the last equation which translates the conversion between time and redshift dependence completes the system of equations:

$$ \frac{dz}{dt} = -H(z)(1 + z). \qquad (11) $$

Equations (5), (6), (7), (8), and (11) are coupled and must be solved simultaneously.

## 2.3. Numerical approach and initial conditions

The search for solutions for systems of ordinary differential equations is a typical stiff problem. Gear methods have excellent stability properties and are widely used for solving chemical kinetic problems and estimating the molecular abundances (Gear 1971; Hindmarsh & Petzold 1995).

At each integration step, we must evaluate summations or differences of very small quantities, which involve a careful analysis. Summation of floating-point numbers is ubiquitous in numerical analysis and has been extensively studied (Higham 1993). In CHEMFAST, we have rearranged the summation to calculate the differential abundances and minimize the numerical summation errors. Following this scheme, we can evaluate the very low abundances of chemical components. We adopted the initial abundances of atoms and ions out of the SBBN theoretical predictions, based on recent determinations of $^4$H$_e$ abundances, as given in Pitrou et al. (2018), namely, $[D/H] \sim 2.46 \times 10^{-5}$, and the helium mass fraction, $Y_p = 0.24709$.

CHEMFAST solves a three-level hydrogen atom model, in the same way as the RECFAST (Seager et al. 1999, 2000) recombination code, correcting the recombination rate to account for the superior excitation levels, but unlike the RECFAST code, we considered “coupling” with all molecular processes. This method has now been supplanted by state-of-the-art codes such as COSMOREC (Chluba & Thomas 2011) and HYREC (Ali-Haïmoud & Hirata 2011; Lee & Ali-Haïmoud 2020), which developed a recombination computation using an effective multilevel atom method in order to meet the precision required for Planck (Planck Collaboration VI 2020) CMB observation analysis.
Compared to the codes mentioned above, CHEMFAST additionally computes the chemical evolution of the homogeneous universe by following a whole chemical network. It proposes a simplified version compared to radiative transfer codes such as ENZO (O'Shea et al. 2004; Bryan et al. 1995; Norman et al. 2007), which has also been used to follow primordial chemistry (see e.g., Schleicher et al. 2008).

We start our calculations in a fully ionized Universe, at redshift $z_i = 10^4$. This redshift corresponds to the age of the Universe at $t_i$, where $x = 1 + z$:

$$t_i = \frac{1}{H_0} \int_{1+z_i}^{\infty} \frac{dx}{x \sqrt{\Omega_{r,0}x^4 + \Omega_{m,0}x^3 + \Omega_{\Lambda,0}}}. \quad (12)$$

We stopped our calculations at a redshift of $z_f = 10$, without taking into account the reionization processes, as they are out of the scope of this work.

The initial number density of baryons $n_{b,i}$ (i.e., number of baryons per cm$^3$ at $z_i$) is computed as:

$$\bar{n}_{b,i} = \frac{3H_0^2 \Omega_{b,i}}{8\pi G \mu_b m_p} = \frac{3H_0^2 \Omega_{b,0}}{8\pi G \mu_b m_p} (1 + z_i)^3, \quad (13)$$

where $m_p$ is the mass of a proton and $\mu_b$ is the molecular weight per baryon particle. The molecular weight depends on the BBN initial conditions and is computed as Seager et al. (1999):

$$\mu_b = \frac{1}{(1 + f_{\text{He}})(1 - Y_p)}, \quad (14)$$

where the helium number fraction, $f_{\text{He}}$, is given by

$$f_{\text{He}} = \frac{n_{\text{He}}}{n_{\text{H}}} = \frac{m_{\text{H}}}{m_{\text{He}}} \frac{Y_p}{1 - Y_p}. \quad (15)$$

Lastly, the initial radiation temperature is $T_{\gamma}(z_i) = T_{\gamma,0}(1 + z_i)$. At the initial redshift we are considering here, the baryons and photons are fully thermally coupled, hence, we can take $T_k(z_i) = T_{\gamma}(z_i)$.

## 2.4. Species abundances

The cosmological recombination process was not instantaneous because the electrons, captured into different atomic energy levels, could not cascade instantaneously down to the ground state. The universe expanded and cooled faster than recombination could be completed and a small fraction of free electrons remained. Recombination processes become dominant when the reactions of photoionization are negligible. We define the redshift of recombination, $z_{\text{rec}}$, as the time at which the abundance of a neutral species is equal to the one of its corresponding ion. In this study, we followed an approximation of the multi-level atom model from Novosyadlyj et al. (2022), which reached the 400-level calculation result from the Cloudy code Ferland et al. (2017) with a 0.01% precision. In Fig. 1, we show the evolution of the relative abundances (i.e., the ratio of species number densities over the total baryon number densities) of all the atomic species considered in our reaction network. We find the successive redshifts of recombination of $He^{2+}$ ($z_{\text{rec,He}^{2+}} \sim 5977$), $He^{+}$ ($z_{\text{rec,He}^{+}} \sim 2556$), $D^{+}$ and $H^{+}$ ($z_{\text{rec,D}^{+}} = z_{\text{rec,H}^{+}} \sim 1387$). The flatness of the atomic relative abundances at $z < 100$ is caused by the inefficiency of collisional reactions due to the expansion of the Universe. It causes a decrease of species densities and freezes the relative abundances.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Relative number densities of atomic and ion species with respect to the total density of baryons, as a function of redshift. The computation goes from $z = 10^4$ to $z = 10$. Hydrogen species are represented as full lines, helium as dashed lines, and deuterium as dash-dotted lines. We indicate the recombination redshifts of primordial helium, deuterium and hydrogen nuclei ($He^{2+}$, $He^{+}$, then $D^{+}$ and $H^{+}$) with vertical dotted line.</figcaption></figure>

### HeH$^+$ molecule.
Helium is the first neutral atom that appeared in the Universe (see Fig. 1). One of the first molecular bonds to be formed is HeH$^+$ via radiative association processes$^1$: $He^+ + H \rightarrow HeH + +v$ (see Roberge & Dalgarno 1982; Kraemer et al. 1995) and $He + H^+ \rightarrow HeH^+ + v$ (see Zygelman et al. 1998). It then mostly contributes to the production of $H_2^+$ via exchange reaction ($HeH^+ + H \rightarrow H_2^+ + He$). Schleicher et al. (2008) also showed that this molecule can play a role in the smearing of primordial fluctuations by scattering CMB photons. Thus, helium is an important component to initiate, at the earliest times, the primordial chemistry. Once a large number of helium atoms are formed through the recombination process, charge transfer with ions becomes possible, initiating the formation of the other neutral species. However, HeH$^+$ appears to be negligible in the molecular cooling processes.

### Molecular hydrogen.
Molecular hydrogen cannot form directly by radiative association between two neutral atoms because $H_2$ does not have a permanent dipole moment. Charge transfer $H_2^+ + H \rightarrow H_2 + H^+$ becomes the only alternative. Once HeH$^+$ was formed, this ion became an important source of $H_2^+$ production through the exchange reaction with some neutral H. Then, $H_2^+$ charge transfer with a neutral species led to $H_2$ molecules. Moreover, as the radiation temperature decreases, $H_2$ can be formed through $H^-$ by radiative attachment $H + e^- \rightarrow H^- + \gamma$, followed by associative detachment $H^- + H \rightarrow H_2 + e^-$ (Field 2000). For densities over $10^8$ cm$^{-3}$, $H_2$ molecule can be effectively formed via three-body reaction $H + H + H \rightarrow H_2 + H$. However, the baryon density after thermal decoupling is $\sim 550$ cm$^{-3}$, which makes this channel ineffective in the context of a homogeneous Universe. The two major routes of $H_2$ formation are visible in Fig. 2 through the two jumps at the redshifts $z \sim 1000-500$ ($H_2^+$ channel) and $z \sim 150$ ($H^-$ channel). After the freeze-out of the relative abundances, we found $n_{\text{rel}}(H_2) = 2 \times 10^{-6}$ for the $H_2$ molecule. We notice that $H_3^+$ appeared in the early Universe. However, $H_3^+$ reactions, which are pivotal to the chemistry of

<footnote>$^1$ He$_2^+$ ions can form earlier but are quickly removed by photodissociation and dissociative recombination.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Relative number densities of molecular species with respect to the total density of baryons, as a function of redshift. The computation goes from $z_i = 10^4$ to $z_f = 10$.</figcaption></figure>

dense interstellar clouds (see Oka 1980; Tennyson 1995; Herbst 2000), do not play a significant role in primordial chemistry.

**Molecular deuterium.** Deuterated-hydrogen molecules have permanent dipole moments, which provide them the capacity to be formed by radiative association (forbidden between two hydrogen atoms). Nevertheless, HD formation is significant when H₂ appeared, as the mechanism of dissociative collision H₂ + D⁺ → H⁺ + HD becomes efficient (see Palla et al. 1995; Stancil & Dalgarno 1998; Signore & Puy 2009). Due to the slight difference between the electronic structure of D and H, H₂ and HD significantly appeared at the same epoch, as shown in Fig. 2. The relative abundance of the HD molecule lies at $n_{rel}(HD) = 6.1 \times 10^{-10}$ after the freeze-out.

## 2.5. Molecular thermal functions

The formation of primordial molecules such as H₂ and HD generates thermal responses due to the excitation of rotational and vibrational levels of molecules. In the range of kinetic and radiation temperature considered, only the rotational levels can be excited. Rotational level populations depend on collisional reactions and radiative processes (CMB absorption or induced and spontaneous emission). Two mechanisms are able to excite or de-excite molecules rotational levels, illustrated in Fig. 3. A molecule can be radiatively excited from a level J to a superior one J', by absorbing a photon external to the medium (for instance from the radiative background), then become de-excited via a collision with another atom or molecule from the medium. The external photon energy in converted into kinetic energy; hence, we need to consider molecular heating. A molecule can be also excited by a collision from level J to J', then emit a photon by spontaneous or induced de-excitation to the J level. If the photon is not reabsorbed by the medium, the associated energy is lost and we would be considering molecular cooling.

To express the energy per volume unit that can be gained, $\Gamma_{mol}$ (heating), or lost, $\Lambda_{mol}$ (cooling), by the medium due to rotational level transitions, we define the thermal molecular function :

$$ \Psi_{mol} = \Gamma_{mol} - \Lambda_{mol} \qquad (16) $$

We detail in Appendices A and B, the computation of the heating and cooling terms, as well as the rotational level populations.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic illustration of the transitions between 2 rotational levels J and J'. Left: Molecule is excited to J' level by collision (double arrow), then de-excited by photon emission (wavy arrow). Assuming an optically thin medium, the photon is not reabsorbed, and the medium is cooled. Right: Molecule is excited by absorption of a photon (single arrow), then de-excited by collisions. The transition energy injected by the photon is then transmitted to the medium in the form of kinetic energy and the medium is heated.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of cooling and heating components of the gas temperature as a function of redshift, in the homogeneous expanding universe scenario. All the components are described in terms of energy density per second. $\Psi_{adia}$ arises from the expansion of the universe. As it is negative (cooling component), we display it as an absolute value. $\Psi_{com}$ is due to Compton coupling between matter and radiation. $\Psi_{mol} = \Psi_{H_2} + \Psi_{HD}$ corresponds to the thermal contribution of molecules H₂ and HD.</figcaption></figure>

### 2.5.1. Molecules thermal influence

H₂ is homonuclear, which implies that it does not possess a permanent dipole moment, which strongly diminish its thermal effect. However, as it is the most abundant molecule, the thermal function due to H₂ plays a role in numerous astrophysical medium (Le Bourlot et al. 1999). In the post-recombination context, H₂ contributes to heat the medium. Nevertheless, in Fig. 4, where we compare the thermal channel contributions in $T_k$ evolution, we can see that the influence of H₂ on the background gas temperature (see Eq. (8)) is negligible (as already shown in Puy et al. 1993 and Pfenniger & Puy 2003).

Infrared emissivities of HD were computed by Dalgarno & Wright (1972). The existence of a permanent dipole moment and its excitation temperature (21 K compared to 112 K for H₂) give this molecule interesting cooling properties. Puy & Signore (1997) showed, in a collapsing molecular protocluster, HD is the main cooling agent around 200 K, confirmed by Flower et al. (2000) and Nakamura & Umemura (2002). Despite being by far less abundant than H₂, the HD molecule can still provide a considerable thermal contribution at the lowest temperatures, around $T_k < 200$ K because of its permanent dipole (Puy et al. 1993; Galli & Palla 1998; Stancil et al. 1998). However, as H₂, HD cooling does not play a role in the evolution
<figure><img src="image_6.png" /><figcaption>Fig. 5. Top: Evolution of kinetic gas temperature, $T_k$ (blue solid line), radiation temperature, $T_γ$ (dashed green line), and 21 cm spin temperature, $T_s$ (orange dash-dotted line), as a function of redshift. Bottom: Evolution of the global 21 cm differential brightness temperature as a function of redshift. At $z_{dec,1\%}(T_k = 99\% T_γ)$, thermal decoupling starts to be effective. At $z(\delta T_{b,min})$, $\delta T_b$ reaches its lowest point, corresponding to the maximum absorption during the dark ages.</figcaption></figure>

of the gas temperature in a homogeneous expanding universe (Flower 2000; Galli & Palla 1998, 2002). Thus, the molecular thermal function, $\Psi_{mol}$, is negligible in the thermochemistry of the homogeneous Universe. We will see its importance in a scenario of gravitational collapse in Sect. 4.

### 3. 21 cm line implementation

The hyperfine transition line of atomic hydrogen, in the ground state, arises due to the interaction between the electron and proton spins. In the excited triplet state, the spins are parallel, whereas the spins in the ground (singlet) state are antiparallel. The 21 cm line is a forbidden line for which the probability for a spontaneous $1 \rightarrow 0$ transition is given by the Einstein coefficient $A_{10} = 2.85 \times 10^{-15} \text{ s}^{-1}$. 21 cm signal observable quantity is the differential brightness temperature, which can be expressed as (see e.g., Madau et al. 1997)

$$ \delta T_b \approx \frac{3hc^3}{32\pi} \frac{A_{10}}{k_B v_{10}^2} \frac{\bar{n}_{HI}(1 + \delta_b)}{(1 + z)H(z)\left(1 + \frac{1+z}{H(z)}\frac{d\Omega_{\parallel}}{dr_{\parallel}}\right)} \left(1 - \frac{T_\gamma}{T_s}\right), \quad (17) $$

where $v_{10} = 1420.406$ MHz is the 21 cm line rest frequency, $\bar{n}_{HI}$ is the mean neutral hydrogen number density, and $\delta_b$ is the relative baryon density perturbation. $H(z)$ is the Hubble rate, and $d\Omega_{\parallel}/dr_{\parallel}$ is the gradient of line of sight peculiar velocity, describing the relative motion of the gas with respect to the Hubble flow. Lastly, $T_γ$ is the radio background temperature, which we set as the CMB one, and $T_s$ is the spin temperature defined from the ratio between the hyperfine population levels.

Hydrogen atoms can be excited by various processes: the absorption and stimulated emission of photons from the CMB redshifted at the 21 cm wavelength, collisions with other hydrogen atoms, free electrons and protons, and Lyman-$\alpha$ radiation from astrophysical sources through Wouthuysen-Field effect (Wouthuysen 1952; Field 1958). During the dark ages, before the formation of the first stars, only the first two mechanisms compete in the excitation of the hyperfine level. In thermal equilibrium, the spin temperature can be written:

$$ T_s^{-1} = \frac{T_\gamma^{-1} + x_c T_k^{-1}}{1 + x_c}. \quad (18) $$

$x_c$ is the total coupling coefficient for collisions given by the expression:

$$ x_c = \frac{(C_H + C_e + C_p)}{A_{10}} \frac{T_{10}}{T_\gamma}, \quad (19) $$

where $T_{10} = h\nu_{10}/k_b$ is the equivalent temperature of energy splitting between the two hyperfine levels. $C_H$, $C_e$, and $C_p$ are the de-excitation rates of the triplet due to collisions with neutral hydrogen, electrons, and protons. The collision term for each $i$ interacting species can be written as $C_i = n_i k_{10}^i(T_k)$, where $k_{10}^i$ is the effective single-atom rate coefficient for the transition 1-0 from collisions with that species. These rate coefficients are obtained through quantum mechanics computation for collisions with other hydrogen atoms (Zygelman 2005), with free electrons (Furlanetto & Furlanetto 2007a), and with protons (Furlanetto & Furlanetto 2007b).

In the context of homogeneous expanding universe, the brightness temperature of Eq. (17) can be approximated as

$$ \delta T_b \approx \frac{3hc^3}{32\pi} \frac{A_{10}}{k_B v_{10}^2} \frac{n_{HI}}{H(z)(1+z)} \left(1 - \frac{T_\gamma}{T_s}\right), \quad (20) $$

$$ \approx 27x_{HI} \left( \frac{\Omega_b h^2}{0.023} \right) \sqrt{\frac{0.15}{\Omega_m h^2} \frac{1+z}{10}} \left( 1 - \frac{T_\gamma}{T_s} \right). \quad (21) $$

As the universe is considered homogeneous, we only keep the mean hydrogen density $\bar{n}_{HI}$, and set $\delta_b = 0$. Moreover, we consider an average peculiar velocity equal to zero. The second expression explicits the dependence of the 21 cm signal to cosmological parameters $\Omega_b$ and $\Omega_m$, the baryon and matter energy density fraction today, and $h$, the Hubble parameter in units of 100 km/s/Mpc. Lastly, $x_{HI}$ is the hydrogen neutral fraction. We show in Fig. 5 (top), the evolution of spin, radiation and kinetic temperatures. The behavior of the spin temperature can be separated in three distinct epochs that determine the observability of the brightness temperature, which is visible in the bottom panel of Fig. 5. These epochs are described as follows.

- Before matter-radiation decoupling, the spin temperature is coupled to the CMB temperature ($T_s = T_γ$). Differential brightness temperature, $\delta T_b$, remains at zero.

- Around $z \sim z_{dec,1\%}$ (defined as the redshift at which $T_k(z_{dec,1\%}) = 99\%T_γ(z_{dec,1\%})$), matter and radiation temperatures are progressively decoupling. The Universe is dense enough for the collision mechanism to dominate in the 21 cm signal production. The spin temperature is thermalized to the gas temperature, $T_s \rightarrow T_k$, and the brightness temperature reaches a minimum of $\delta T_b \sim -44$ mK at $z(\delta T_{b,min} \approx 89)$.

- Due to the expansion of the universe, the baryons get more and more diluted and the collision mechanism progressively becomes ineffective. The spin temperature relaxes to the CMB temperature, $T_s \rightarrow T_γ$, which brings $\delta T_b$ back to zero.

This section has served to present the functioning of CHEFAST within the known framework of an expanding homogeneous universe. The code can be used to solve a system of stiff coupled differential equations in order to obtain the evolution of the abundance of chemical species based on a network of reactions. The code can also be used to calculate the thermal influence of molecules on temperature through the excitation of their rotational levels. In a homogeneous expanding universe, we can see that this effect is negligible. Finally, by following the thermal history of the universe and the chemical reactions, we can estimate the intensity of the global 21 cm signal due to arising from collisional coupling. The 21 cm global signal computed
with CHEMFAST is in accordance with the state-of-the-art literature (see e.g., Furlanetto et al. 2006; Pritchard & Loeb 2012).

## 4. Dynamics of collapsing clouds and 21 cm hydrogen line emission

In the context of homogeneous expansion seen in the first part, the thermal influence of H₂ and HD molecules on the gas temperature is negligible, due to their very low mean abundances. Their contribution is completely drowned out by the adiabatic cooling of the matter-radiation thermal coupling associated with Compton scattering. However, we know that the universe is actually not homogeneous, but exhibited small primordial density fluctuations, which became seeds for the formation of structures.

When a density perturbation grows enough due to gravity in order to reach a density contrast comparable to unity, it can no longer be described by the linear theory of perturbations. We must then consider a model to follow its evolution until a potential collapse. Molecules in overdensities are known to cool the gas via the excitation of their rotational and vibrational levels and permitting the further gravitational collapse of the gas, thereby setting up the right conditions for the formation of the first structures.

In this section we focus on the evolution of an individual overdensity. We implement a collapse model in CHEMFAST by modifying the equations of dynamics. Our goal is then to estimate the 21 cm signal from the collapsing overdensity under the influence of molecular cooling.

### 4.1. Model of collapsing cloud

We consider a spherical pressure-supported collapse model with basic hypothesis. The model we use is based on the work of Lahav (1986), extended by Puy & Signore (1996). Unlike the free-fall model, the spherical collapse we consider cannot be solved analytically and requires a numerical resolution.

We assume that the overdense region we follow is isothermal, spherical, and without rotation. We also assume the matter density to be constant inside the overdense region. This matter sphere, at first following the expansion of the universe, will progressively slow down with respect to the background, due to the density excess, until the point when it drops out of expansion, and starts to collapse back on itself. In the range of redshift we consider, the universe is fully dominated by matter, and we can assume a linear $r^{2/3}$ growth of the fluctuations with the expansion of the Universe. Hence, we consider a spectrum of isothermal perturbations, that are described by the spectrum (Gott & Rees 1975):

$$ \frac{\delta \rho}{\bar{\rho}} = \left( \frac{M}{M_0} \right)^{-1/3} (1+z)^{-1}, \qquad (22) $$

where the mass M of the overdense region is normalized to $M_0$, defined as a characteristic mass-scale of $M_0 = 10^{15} M_\odot$, which is the typical mass of large galaxy clusters today, also used in Lahav (1986).

### 4.2. Equations

As in Sect. 2, here we solve a set of differential equations, this time for a collapsing overdense region.

**Baryon number density.** The matter conservation inside the cloud gives:

$$ \frac{dn_b}{dt} = -3 \frac{n_b}{r} v_r + \sum_{\xi} \left( \frac{dn_{\xi}}{dt} \right)_{\text{chem}}, \quad \text{with } v_r = \frac{dr}{dt}. \qquad (23) $$

The equation is similar to the expansion case of Sect. 2, the scale factor, $a$, is simply replaced by the radius of the cloud, $r$. The last term in the right-hand side of the equation is the contribution from chemical reactions.

**Species densities.** Similarly, for each chemical component $\xi$ we have:

$$ \frac{dn_{\xi}}{dt} = -3 \frac{n_{\xi}}{r} v_r + \left( \frac{dn_{\xi}}{dt} \right)_{\text{chem}}. \qquad (24) $$

**Temperatures.** The radiation temperature evolves just like in the homogeneous case, with the expansion of the universe:

$$ \frac{dT_{\gamma}}{dt} = -H(t) T_{\gamma}. \qquad (25) $$

The kinetic gas temperature, under the assumption of perfect gas, can be computed as:

$$ \frac{dT_k}{dt} = \frac{2}{3n_b k_B} \left( \Psi_{\text{grav}} + \Psi_{\text{mol}} + \Psi_{\text{com}} \right). \qquad (26) $$

$\Psi_{\text{grav}} = -3n_b k_B T_k \frac{v_r}{r}$ accounts for the density variation arising from the gravitational collapse, $\Psi_{\text{com}}$ describes the Compton coupling between matter and radiation (see Eq. (10)), and $\Psi_{\text{mol}}$ expresses the thermal influence of molecules on the gas temperature through the excitation of their rotational levels (see Appendix A).

**Radius and velocity of collapse.** Starting from the energy balance:

$$ E = -\frac{3GM^2}{5r} + \frac{3}{10}Mv_r^2 + \frac{3}{2}N_b k_b T_k, \qquad (27) $$

where the first term corresponds to the potential energy in a spherical system with a central potential. The second one accounts for the associated kinetic energy, and the last one for the internal energy. Then, $N_b$ the total number of baryons in the cloud.

We inject the temperature into Eq. (26), expressing the energy variation in the system as $\frac{dE}{dt} = r^3(\Psi_{\text{com}} + \Psi_{\text{mol}})$, and use $M = N_b \mu_b m_H \frac{\Omega_m}{\Omega_b}^2$ to obtain

$$ a_r = -\frac{GM}{r^2} + \frac{5k_b T_k}{\mu m_H} \frac{1}{r} \frac{\Omega_b}{\Omega_m} \quad \text{with} \quad a_r = \frac{dv_r}{dt}. \qquad (28) $$

### 4.3. Initial conditions

Considering a single cloud of mass, M, we follow the collapse starting from the turnaround point, when the gravity takes over the expansion and the perturbation starts to increase in density, as it was only diluting slower than the background before this moment. In a spherical collapse model, an overdense region

<footnote n="2">² We suppose the balance between baryons and dark matter is identical as the background.</footnote>
reaches the turnaround point when the density contrast reaches $(\frac{3\pi}{4})^2$ (Padmanabhan 2002; Peebles 1980).

Hence, the baryon number density at turnaround is expressed by

$$n_{b,ta} = \left(\frac{3\pi}{4}\right)^2 \bar{n}_b(z_{ta}). \qquad (29)$$

Where $\bar{n}_b(z_{ta})$ is the background baryon number density at turn-around redshift, $z_{ta}$.

Following the approach of Gunn et al. (1972) and Lahav (1986), we compute the redshift of turn-around for a given mass:

$$1 + z_{ta} = \left(\frac{3\pi}{4}\right)^{-2/3} \left(\frac{M}{M_0}\right)^{-1/3}. \qquad (30)$$

This is the redshift at which we start to follow the collapse in CHEMFAST. It is only dependent of the mass, M, of the cloud.

The radius of the overdense region is directly derived from the density

$$r_{ta} = \left( \frac{3M}{4\pi\rho_{b,ta}} \right)^{1/3}, \qquad (31)$$

with $\rho_{b,ta} = \mu_b m_H n_{b,ta}$. The collapse velocity $v_{ta}$ at turn-around is by definition set at zero:

$$v_{ta} = \frac{dr_{ta}}{dt} = 0. \qquad (32)$$

For the gas temperature, the computation is less evident. If the gas and radiation temperature were fully decoupled, we could consider the adiabatic gas temperature evolves as $T_k(z) \propto \rho_b(z)^{2/3}$. Then, the gas temperature inside an overdensity at turnaround redshift could be written

$$T_{k,ta}^{\text{adia}} = \left(\frac{3\pi}{4}\right)^{4/3} \bar{T}_k(z_{ta}), \qquad (33)$$

where the adia label refers to adiabatic, in the sense that the gas temperature only changes due to the expansion or compression of its environment, without any exchange of heat. However, as we will see, our model can provide very high redshift of turn-around, depending on the chosen mass. As a consequence, the gas temperature does not fully evolve adiabatically, and the Compton scattering between free electrons and photons can still be efficient at the turn-around redshifts. Choosing $T_{k,ta}^{\text{adia}}$ would mean ignoring the effect of Compton scattering in the period preceding the turnaround. We therefore need to estimate the temperature at turn-around, which lies in the range $[\bar{T}_k(z_{ta}), (\frac{3\pi}{4})^{4/3} \bar{T}_k(z_{ta})]$, accounting for the residual Compton coupling. We chose the turnaround temperature empirically, such that the Compton thermal channel, $\Psi_{\text{com}}$, is steady in the very first moments of collapse at the beginning of the run. Actually, we have checked that this detail does not significantly impact the behavior of $T_k$ during the rest of the collapse. It does, however, help to avoid non-physical behavior<fnref n="3" /> in the very first moments of the collapse.

To summarize, our model takes the mass of the overdense region as the only free parameter, from which all the initial

conditions are derived. For a given mass, we ran CHEMFAST a first time in expansion mode up to the turnaround redshift and extracted the background values of the species abundances, baryon density, and temperatures. These data were then used to compute the initial conditions in the collapse setup.

# 5. Results

## 5.1. Thermal evolution of $10^8 M_{\odot}$ collapsing cloud

We focus on the analysis of a $10^8 M_{\odot}$ collapsing overdensity, which corresponds to the order of magnitude expected for “massive minihalos”, the least massive halos that are expected to be able to host star formation (see e.g., Tegmark et al. (1997), Yoshida et al. (2003), Glover (2005), Haiman & Bryan (2006), Trenti & Stiavelli (2009), Greif (2015) for discussions on the formation of the first structures). In our model, this overdensity reaches its redshift of turnaround at $z_{ta} = 93$, with a temperature $T_{k,ta} = 2.1 \cdot \bar{T}_k(z_{ta})$. We first discuss the evolution of the gas temperature from its different thermal channels and then propose an estimate of the 21 cm line arising from the collapsing overdensity.

For all the following figures, we display the evolution of the quantities as a function of the time to collapse. Defining the turnaround time:

$$t_{ta} = \sqrt{\frac{3\pi}{32G\rho_{ta}}}, \qquad (34)$$

the time of re-collapse can be expressed as $t_{\text{coll}} = 2t_{ta}$ (Lahav 1986). We normalized time in the figures to evolve between 0 and 1.

In Fig. 6, we follow the intensity of all the thermal channels acting on the gas temperature within the cloud (see Eq. (26)). The heating due to gravitational contraction, $\Psi_{\text{grav}}$ (solid blue line), gets stronger as the radius of the halo decreases and the collapse velocity increases; $\Psi_{\text{com}}$ (dashed curve) arises from the Compton coupling between matter and radiation. Since our initial conditions provide a gas temperature in the halo above the background radiation temperature, the Compton interaction tends to cool the gas temperature towards the radiation temperature, which explains why it presents a negative contribution. The Compton coupling dominates over the other contributions in the first instants and decreases as the collapse progresses. It eventually becomes negligible.

The $\Psi_{\text{mol}}$ term (dotted curve) represents the molecular thermal function summing over the contributions from the two most important molecules $\Psi_{\text{H}_2}$ and $\Psi_{\text{HD}}$. Furthermore, $\Psi_{\text{mol}}$ is almost completely driven by $\text{H}_2$ molecule. The HD contribution becomes important for lower temperatures than the one we are interested in (<200 K), it is negligible in the case we are studying (McGreer & Bryan 2008). In Fig. 6, $\Psi_{\text{mol}}$ contributes to cool the gas. Initially negligible, its importance increases sharply as the density of the molecules increases and the gas temperature grows, until it becomes significant and even dominates the other thermal channels at the end of the collapse. Compared to the expansion scenario, the molecular thermal function has a strong influence on the gas temperature.

We display the evolution of the radiation and gas temperatures inside the $10^8 M_{\odot}$ collapsing overdensity in Fig. 7 (top panel). The radiation temperature evolves exactly in the same way as in the background, and is shown for comparison. The evolution of the gas temperature can be divided in several stages:

<footnote n="3">³ In particular, we avoid having a very high $T_k$ at turn-around, which decreases drastically when the run starts, as if we just added a Compton coupling that was not here before turnaround.</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 6.** Evolution of cooling and heating components of the gas temperature in a $10^8 M_\odot$ collapsing cloud as a function of time, normalized by $t_{\text{coll}}$. All the components are described in terms of energy density per second. $\Psi_{\text{grav}}$ arises from to gravitational collapse of the gas; $\Psi_{\text{com}}$ is due to Compton coupling. Finally, $\Psi_{\text{mol}} = \Psi_{\text{H}_2} + \Psi_{\text{HD}}$ corresponds to the thermal contribution of molecules H$_2$ and HD. The negative components $\Psi_{\text{com}}$, $\Psi_{\text{mol}}$, $\Psi_{\text{H}_2}$, and $\Psi_{\text{HD}}$ are displayed as absolute values.

- During the very first moments of the collapse, Compton coupling is the dominant thermal mechanism and tends to bring $T_k$ (blue solid line) to $T_\gamma$ (green dashed line) acting against gravitational contraction. As a result, the gas temperature inside the overdensity slightly decreases. This stage holds for a very short time; hence, the decrease is not visible in the figure.

- Heating due to the gravitational collapse rapidly dominates over the Compton coupling, starting from $t = 0.06t_{\text{coll}}$. From this moment, the gas temperature inside the cloud grows drastically, under the influence of gravity.

- Molecular cooling rapidly becomes more important due to the increasing gas temperature and H$_2$ density inside the cloud. This thermal channel dominates, starting from $t = 0.89t_{\text{coll}}$, the molecular thermal function cools the medium more than it heats up by gravitational contraction. The gas temperature inside the cloud decreases despite the ongoing collapse.

## 5.2. 21 cm line emission from collapsing halo

At the top of Fig. 7, we can see the spin temperature, $T_s$ (dashed line), closely follows gas temperature, $T_k$, throughout the collapse. Looking again at the spin temperature (Eq. (18)), we can explain it by the following arguments. We are in a collapse scenario, where the matter number density is increased up to a factor ~177 compared to the background, at the end of the collapse. Therefore, the collisional coupling is important and dominates the 21 cm line production.

The differential brightness temperature, $\delta T_b$, is positive, meaning that the signal presents an emission feature, different from the global signal absorption feature from the expansion scenario. Throughout the collapse, $T_s$ is always greater than radiation temperature, $T_\gamma$. Therefore, the term $(1 - \frac{T_\gamma}{T_s})$ in Eq. (21) stays between 0 and 1. Moreover, we observe in $\delta T_b$ the same type behavior than $T_k$, with a turnover caused by molecular cooling. Indeed, since the spin temperature completely follows $T_k$, it is sensitive to the same thermal processes, transmitted to $\delta T_b$. We can note that, even without the activation of molecular cooling,

<figure><img src="image_8.png" /></figure>

**Fig. 7.** Top: Evolution of radiation temperature, $T_\gamma$ (dash-dotted line), gas temperature, $T_k$ (solid line) and 21 cm spin temperature, $T_s$ (dotted line) for a collapsing $10^8 M_\odot$ cloud, as a function of the time to collapse. We highlight on top of the figure the times when the dominant thermal channel changes. Bottom: Evolution of the associated 21 cm brightness temperature from a collapsing $10^8 M_\odot$ cloud.

$\delta T_b$ would eventually saturate and reach a maximum limit when $T_s \gg T_\gamma$. In this regime, the dependence on $T_s$ in Eq. (21) would become unimportant.

We adopted the Eq. (21) for the computation of the 21 cm differential brightness temperature. By doing so, we made several approximations, to which we dedicate the following discussion. Firstly, we assumed a zero peculiar velocity for the halo. This choice is not problematic by itself, but it should be noted that it is arbitrary. When studying a whole population of halos, we have to take into account the effects of peculiar velocities, which induce an enhancement of the fluctuations at all scales (Bharadwaj & Ali 2004; Barkana & Loeb 2005).

Equation (21) also implies two major approximations for calculating the 21 cm signal from a halo. The supposition that the optical depth is thin can no longer technically hold in a halo. For an accurate calculation, it would be more appropriate to consider a density profile (such as a truncated isothermal sphere used in Iliev et al. 2002, or a Moore profile used in Furugori et al. 2020) for the halo; instead of a homogeneous density as we have done. The optical depth would then vary as a function of radius (maximum at the center of the halo where the maximum amount of matter is traversed) and the brightness temperature shall be integrated over the surface of the halo. The other approximation we make is to assume the 21 cm line profile as a Dirac function. A more appropriate approximation, such as a thermal-Doppler-broadened line profile, would allow the velocity dispersion of particles within the halo to be taken into account.

However, accounting for the variable optical depth and the Doppler broadening would require us to get rid of the assumption of homogeneous density and temperature in the halo; it could also strongly modify the 21 cm signal. The results we show are in that sense approximations, presenting our approach, and we leave the improvement of the collapse model to a future work.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Top: Evolution of gas temperature, $T_k$ (solid lines), for clouds of mass between $10^{6.0}$ and $10^{8.0} M_\odot$, as a function of the time to collapse (top). On the bottom of this plot are displayed the 21 cm differential brightness temperatures. Bottom: Same as the top panel, for masses between $10^{8.0}$ and $10^{8.8} M_\odot$, also displaying the 21 cm spin temperature $T_s$ (dotted lines).</figcaption></figure>

## 5.3. Thermal evolution of a range of halo masses

After focusing on the physics at work in a $10^8 M_\odot$ collapsing halo, we present the temperatures evolution for a whole range of different masses in Fig. 8 and display their associated turnaround redshifts and gas temperature in Table 2. Several masses are shown in the same figure as a function of their time to collapse in order to easily compare them. However, it is important to bear in mind that these collapses take place at different redshifts, which are correlated with the initial mass.

In the top panel of Fig. 8, we show the temperatures evolution for masses between $10^6$ and $10^8 M_\odot$. The gas temperature presents the same general behavior for all the masses: a growth provoked by the contraction of the clouds, followed by a decline in the last moments of the collapse, due to the molecular cooling. However, the fall in $T_k$ occurs at different moments of the collapse depending on the mass. It happens later for the smaller masses, such as $10^6 M_\odot$, resulting in a higher maximal gas temperature, and 21 cm brightness temperature.

On the other hand, we notice an inversion of this behavior starting from masses $> 10^{7.0} M_\odot$. We observe that the maximum

gas temperature increases again and the moment of collapse at which the maximum temperature is reached is delayed.

Looking at the top of Fig. 9, where is plotted the ratio between $\Psi_{mol}$ and $\Psi_{grav}$, we observe that above $10^7 M_\odot$, the more we increase the mass of the cloud, the later the molecular cooling dominates. This effect can be understood by looking at the number density of $H_2$ inside the halos, presented in the bottom panel of Fig. 9. We can expect a higher overall baryon density in the halos of smaller mass, as they are the first ones to form. The density contrast between the overdense region and the background is similar for all halos, so the one forming when the background is denser are also expected to show a higher density. However, the $H_2$ density directly depends on the presence of the species leading to its formation, $H_2^+$ via charge transfer reaction $H_2^+ + H \rightarrow H_2 + H$, and $H^-$ via associative detachment $H^- + H \rightarrow H_2 + e^-$. Finally, $H^-$ is actually rapidly destroyed by CMB photons when forming at very high redshift, before it would be able to contribute to $H_2$ formation.

For the lowest masses (e.g., $10^6 M_\odot$, which has a turn-around redshift of $z_{ta} = 436$ in our model), despite a high overall density in the cloud, the formation of $H_2$ molecules by the $H^-$ channel is strongly suppressed. Up to $\sim 10^7 M_\odot$ ($z_{ta} = 202$), the $H^-$ are progressively less destroyed, allowing more $H_2$ to be formed and resulting in stronger molecular cooling. In the homogeneous universe (see Fig. 2), the growth of the $H_2$ population via $H^-$ takes place around $z = 150$, when the temperature of the photons decreases enough such that it does not allow them to destroy $H^-$ anymore, letting it maximally convert into $H_2$. However, as the $> 10^7 M_\odot$ halos form later and later, even if the $H^-$ population is not suppressed, the baryons density (and therefore the $H_2$ density) progressively decreases, explaining the gradual decline of molecular cooling.

Interestingly, the 21 cm brightness temperature does not present the same behavior and the highest maximum $\delta T_b$ is still reached in the lightest halos. The reason lies in the redshift dependence of the signal, which is not visible in the figure, but is explicit in Eq. (21). The collapses of the largest masses occur too late and the gas temperature increase is not enough to overcome the decline in redshift, even though $T_s$ is still strongly coupled to $T_k$ via collisions, as can be seen in the bottom panel of Fig. 8.

# 6. Discussion and perspectives

We explored the evolution of species abundances in detail during and after successive recombinations. To do so, we improved the code CHEMFAST, which solves a system of stiff coupled differential equations. Part of the system describes the dynamics (density, baryon, and radiation temperatures) and another is dedicated to describing the network of collisional, electronic, and photo-processes between the species. Particular attention has been paid to tracking the molecular abundances.

We first applied CHEMFAST to the simple expanding homogeneous universe framework. We computed the successive recombinations of $He^{2+}$, $He^+$, $H^+$, and $D^+$, as well as the relative abundances of all the atoms, ions, and molecules contributing to the reaction network. We also computed the thermal influence of molecules through the excitation of their rotational levels. Furthermore, we implemented the 21 cm differential brightness temperature global signal computation during the dark ages, arising from collisional excitation. It presents an absorption peak with an intensity of -40 mK at $z=89$. This result is largely corroborated by the literature (e.g., Furlanetto et al. 2006; Pritchard & Loeb 2012) and validates the effectiveness of our code.
<table><caption>**Table 2.** Redshifts and temperatures of turnaround.</caption><thead><tr><th>log<sub>10</sub>M [M<sub>☉</sub>]</th><th>6.0</th><th>6.5</th><th>7.0</th><th>7.5</th><th>8.0</th><th>8.2</th><th>8.5</th><th>8.8</th></tr></thead><tbody><tr><td>z<sub>ta</sub></td><td>436</td><td>297</td><td>202</td><td>137</td><td>93</td><td>80</td><td>63</td><td>50</td></tr><tr><td>z<sub>coll</sub></td><td>185</td><td>127</td><td>86</td><td>59</td><td>40</td><td>34</td><td>26</td><td>21</td></tr><tr><td>T<sub>k,ta</sub> (factor of T&#773;<sub>k</sub>(z<sub>ta</sub>))</td><td>1.05</td><td>1.14</td><td>1.30</td><td>1.60</td><td>2.10</td><td>2.40</td><td>2.90</td><td>3.14</td></tr></tbody></table>

**Notes.** Temperature is shown as a factor of $\bar{T}_k(z_{ta})$, the background temperature at $z_{ta}$. $z_{coll}$ is the redshift at which we stop the run, when the time to collapse $t_{coll}$ is reached.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Top: Evolution of ratio between $\Psi_{mol}$ and $\Psi_{grav}$ thermal channels, arising from molecular excitation and gravitational contraction respectively, for clouds between $10^{6.0}$ and $10^{8.0} M_{\odot}$, as a function of the time to collapse. Bottom: H$_2$ number density evolution as a function of the time to collapse, for clouds between $10^{6.0}$ and $10^{8.0} M_{\odot}$.</figcaption></figure>

In a second study, our goal was to investigate the 21 cm signal evolution in a new context and to assess how it could be distinguished from the sky-averaged global signal. To this end, we flexibly integrated a new set of dynamical equations into our abundance calculation code CHEMFAST, to follow the collapse of an overdense region. We based our work on a spherical collapse model developed by Lahav (1986), taking into account the pressure forces, with the cloud mass as the only free parameter. We followed the collapse of a typical cloud of $10^8 M_{\odot}$, expected to be able to host the formation of the first galaxies. Within the

collapse, molecules are present in much greater density than in the homogeneous case. Calculating the excitation of the rotational levels for the H$_2$ and HD molecules shows us that they have a strong thermal impact on the cloud's gas temperature. At this point, the contribution becomes dominant and takes over the gravitational heating in the last moments of the collapse, leading to a cooling of the gas.

Finally, we computed the 21 cm line brightness temperature within the collapsing cloud. The signal shows a very different signature compared to the global signal of the dark ages. It presents an emission feature and is also affected by molecular cooling, in the same way as the gas temperature to which the brightness temperature is coupled via collisions. We also present the thermal evolution for overdense regions of $10^6$ to $10^{8.8} M_{\odot}$. The strongest 21 cm signals are expected from the halos of the lowest masses, which have lower H$_2$ abundances, due to the suppression of the H$^-$ formation channel by CMB radiation. These halos therefore exhibit less efficient molecular cooling. The intensity of molecular cooling increases progressively with the halo masses until redshifts $z \sim 150$ at which H$^-$ is no longer destroyed. Halos whose formation begins after this redshift see a decline in molecular cooling, due to their lower density. Despite a higher $T_k$, this pattern is not followed by the 21 cm brightness temperature because of its redshift dependence. The H$_2$ chemistry is not trivial, however, and further investigation would be required to understand the H$_2$ density evolution in collapse context, depending on the initial mass. We plan to conduct this investigation in a future work.

These particular signatures of the 21 cm signal are promising with respect to probing the thermal history inside halos. They could also present an observational interest at the smallest scales of the 21 cm power spectrum in the context of the forthcoming lunar-based observations of the dark ages. One way to estimate this impact would be to compute the brightness temperature arising from a whole halo population, by associating our single halo results to the computation of a halo mass function (Lukić et al. 2007; Murray et al. 2013). However, the particularly high turnaround redshifts of our collapse model make the determination of the halo mass function particularly challenging.

Indeed, our model predicts cloud formation at particularly high redshifts, compared to what is currently proposed in the literature (e.g., Trenti & Stiavelli 2009; Glover 2012; Greif 2015). State-of-the-art halo formation models start at redshifts of $z \sim 30$ for the halos of the lowest mass, which are still allowed the gas to collapse ($\sim 10^5-10^6 M_{\odot}$). Similar ranges are found by studies of 21 cm signal predictions from halos Iliev et al. (2002), Furugori et al. (2020), Novosyadlyj et al. (2024), which makes it difficult to make comparisons with them, as our hypotheses are very different. This discrepancy can be explained by the simplicity of the model we have chosen. Furthermore, we have seen that the assumption of homogeneous density and temperature inside the cloud could prove problematic for the study of the 21 cm
signal. We intend to study other collapse models with CHEFAST in future work and include the density and temperature profiles in our modeling. We also plan to integrate elliptical geometries, which are more consistent with collapse simulations.

Lastly, molecular cooling can induce the development of thermal instabilities. These effects were first studied by Parker (1953), followed by Field (1965), and later Balbus (1986), who derived instability criteria, leading to a growth of the perturbations in the gas. It is now understood that these thermal instabilities play an important role in the fragmentation of halos and the formation of the first stars (Safranek-Shrader et al. 2010; Yoshida et al. 2006; Bromm 2013; Glover 2012; Greif et al. 2013). The aim of our simple model has been to highlight the crucial thermal influence of molecules on the thermal history and 21 cm signal, but it does not allow us to explore the evolution of thermal instabilities in depth. We would then have to turn to adaptive mesh refining (AMR) type simulations, such as the ENZO (Bryan et al. 2014) or RAMSES (Teyssier 2002) codes. This type of approach makes it possible to spatially monitor instabilities that can lead to fragmentation within the structures (Tang & Chen 2024).

*Acknowledgements.* We thank Alice Faure for her support in the code optimization. We also thank Daniel Pfenniger and Benoît Sémelin for helpful discussions. These results have been made possible thanks to LUPM’s cloud computing infrastructure, founded by PHONE project - Occitanie Regions, Ocevu labex, and France-Grilles.

## References

Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New Astron., 2, 181

Ali-Haïmoud, Y., & Hirata, C. M. 2011, Phys. Rev. D, 83, 043513

Artuc, K., & de Lera Acedo, E. 2024, arXiv e-prints [arXiv:2406.10096]

Balbus, S. A. 1986, ApJ, 303, L79

Bale, S. D., Bassett, N., Burns, J. O., et al. 2023, arXiv e-prints [arXiv:2301.10345]

Barkana, R., & Loeb, A. 2005, ApJ, 626, 1

Bates, D. R. 1951, MNRAS, 111, 303

Bentum, M. J., Verma, M. K., Rajan, R. T., et al. 2020, Adv. Space Res., 65, 856

Bharadwaj, S., & Ali, S. S. 2004, MNRAS, 352, 142

Bowman, J. D., Cairns, I., Kaplan, D. L., et al. 2013, PASA, 30, e031

Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67

Bromm, V. 2013, Rept. Prog. Phys., 76, 112901

Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23

Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., & Ostriker, J. P. 1995, Comp. Phys. Commun., 89, 149

Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19

Burns, J., Bale, S., Bradley, R., et al. 2021a, arXiv e-prints [arXiv:2103.05085]

Burns, J., Hallinan, G., Chang, T. C., et al. 2021b, arXiv e-prints [arXiv:2103.08623]

Capitelli, M., Coppola, C. M., Diomede, P., & Longo, S. 2007, A&A, 470, 811

Chen, X., Gao, F., Wu, F., et al. 2024, Phil. Trans. R. Soc. A: Math. Phys. Eng. Sci., 382, 20230094

Chluba, J., & Thomas, R. M. 2011, MNRAS, 412, 748

Cohen-Tannoudji, C., Dui, B., & Laloe, F. 1973, Mecanique quantique (CNRS Éditions)

Coppola, C. M., Longo, S., Capitelli, M., Palla, F., & Galli, D. 2011, ApJS, 193, 7

Coppola, C. M., D’Introno, R., Galli, D., Tennyson, J., & Longo, S. 2012, ApJS, 199, 16

Coppola, C. M., Galli, D., Palla, F., Longo, S., & Chluba, J. 2013, MNRAS, 434, 114

Cyburt, R. H., Fields, B. D., Olive, K. A., & Yeh, T.-H. 2016, Rev. Mod. Phys., 88, 015004

Dalgarno, A., & Wright, E. L. 1972, ApJ, 174, L49

DeBoer, D. R., et al. 2017, PASP, 129, 045001

de Lera Acedo, E., de Villiers, D. I. L., Razavi-Ghods, N., et al. 2022, Nat. Astron., 6, 984

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385

Fialkov, A., Gessey-Jones, T., & Dhandha, J. 2024, Phil. Trans. R. Soc. A: Math. Phys. Eng. Sci., 382, 2271

Field, G. B. 1958, Proc. IRE, 46, 240

Field, G. B. 1965, ApJ, 142, 531

Field, D. 2000, A&A, 362, 774

Flower, D. R. 2000, MNRAS, 318, 875

Flower, D. R., Le Bourlot, J., Pineau des Forêts, G., & Roueff, E. 2000, MNRAS, 314, 753

Furlanetto, S. R., & Furlanetto, M. R. 2007a, MNRAS, 374, 547

Furlanetto, S. R., & Furlanetto, M. R. 2007b, MNRAS, 379, 130

Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181

Furugori, K., Abe, K. T., Tanaka, T., et al. 2020, MNRAS, 494, 4334

Galli, D., & Palla, F. 1998, A&A, 335, 403

Galli, D., & Palla, F. 2002, Planet. Space Sci., 50, 1197

Galli, D., & Palla, F. 2013, ARA&A, 51, 163

Gay, C. D., Stancil, P. C., Lepp, S., & Dalgarno, A. 2011, ApJ, 737, 44

Gear, C. W. 1971, Numerical initial value problems in ordinary differential equations (Prentice-Hall Series in Automatic Computation)

Glover, S. 2005, Space Sci. Rev., 117, 445

Glover, S. 2012, The First Stars (Springer Berlin Heidelberg), 103

Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627

Glover, S., & Savin, D. W. 2006, Philos. Trans. R. Soc. Lond. Ser. A, 364, 3107

Glover, S. C. O., & Savin, D. W. 2009, MNRAS, 393, 911

Goel, A., Bandyopadhyay, S., Lazio, J., et al. 2022, arXiv e-prints [arXiv:2205.05745]

Gott, J. R., III, & Rees, M. J. 1975, A&A, 45, 365

Greif, T. H. 2015, Comput. Astrophys. Cosmol., 2, 3

Greif, T. H., Springel, V., & Bromm, V. 2013, MNRAS, 434, 3408

Gunn, J. E., Gott, J., & Richard, I. 1972, ApJ, 176, 1

Haiman, Z., & Bryan, G. L. 2006, ApJ, 650, 7

Herbst, E. 1995, Annu. Rev. Phys. Chem., 46, 27

Herbst, E. 2000, in Astronomy, physics and chemistry of H$^{3+}$, 358, 2359

Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505

Higham, N. J. 1993, SIAM J. Sci. Comput., 14, 783

Hindmarsh, A. C., & Petzold, L. R. 1995, Comput. Phys., 9, 148

Hirasawa, T. 1969, Prog. Theor. Phys., 42, 523

Hirata, C. M., & Padmanabhan, N. 2006, MNRAS, 372, 1175

Iliev, I. T., Shapiro, P. R., Ferrara, A., & Martel, H. 2002, ApJ, 572, L123

Jacobs, D. C., Hazelton, B. J., Trott, C. M., et al. 2016, ApJ, 825, 114

Juřek, M., Špirko, V., & Kraemer, W. P. 1995, Chem. Phys., 193, 287

Klein Wolt, M., Falcke, H., & Koopmans, L. 2024, in American Astronomical Society Meeting Abstracts, 243, 264.01

Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, Proceedings of Advancing Astrophysics with the Square Kilometre Array – PoS(AASKA14) (Sissa Medialab), AASKA14

Kraemer, W. P., Špirko, V., & Juřek, M. 1995, Chem. Phys. Lett., 236, 177

Kutner, M. L. 1984, Fund. Cosmic Phys., 9, 233

Lahav, O. 1986, MNRAS, 220, 259

Latter, W. B., & Black, J. H. 1991, ApJ, 372, 161

Le Bourlot, J., Pineau des Forêts, G., & Flower, D. R. 1999, MNRAS, 305, 802

Lee, N., & Ali-Haïmoud, Y. 2020, Phys. Rev. D, 102, 083517

Lepp, S., & Shull, J. M. 1984, ApJ, 280, 465

Lepp, S., Stancil, P. C., & Dalgarno, A. 2002, J. Phys. B At. Mol. Phys., 35, R57

Le Teuff, Y. H., Millar, T. J., & Markwick, A. J. 2000, A&AS, 146, 157

Linder, F., Janev, R. K., & Botero, J. 1995, in Atomic and Molecular Processes in Fusion Edge Plasmas, ed. R. K. Janev, 397

Lukić, Z., Heitmann, K., Habib, S., Bashinsky, S., & Ricker, P. M. 2007, ApJ, 671, 1160

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429

Mather, J. C., Cheng, E. S., Cottingham, D. A., et al. 1994, ApJ, 420, 439

Matsuda, T., Satō, H., & Takeda, H. 1969, Prog. Theor. Phys., 42, 219

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

McGreer, I. D., & Bryan, G. L. 2008, ApJ, 685, 8

Meiksin, A. 2011, MNRAS, 417, 1480

Mizusawa, H., Omukai, K., & Nishi, R. 2005, PASJ, 57, 951

Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astron. Comput., 3, 23

Nakamura, F., & Umemura, M. 2002, ApJ, 569, 549

Nambissan, T., Subrahmanyan, R., Somashekar, R., et al. 2021, Exp. Astron., 51, 193

Norman, M. L., Bryan, G. L., Harkness, R., et al. 2007, arXiv e-prints [arXiv:0705.1556]

Novosyadlyj, B., Shulga, V., Kulinich, Y., & Han, W. 2020, Phys. Dark Universe, 27, 100422

Novosyadlyj, B., Kulinich, Y., Melekh, B., & Shulga, V. 2022, A&A, 663, A120

Novosyadlyj, B., Kulinich, Y., & Konovalenko, O. 2024, J. Phys. Stud., 28, 1901

Novosyadlyj, B., Kulinich, Y., & Koval, D. 2025, Phys. Rev. D, 111, 083514

Oka, T. 1980, Phys. Rev. Lett., 45, 531
O'Shea, B. W., Bryan, G., Bordner, J., et al. 2004, arXiv e-prints
[arXiv:astro-ph/0403044]

Padmanabhan, T. 2002, Theoretical Astrophysics, Volume III: Galaxies and Cosmology (Cambridge University Press)

Palla, F., Galli, D., & Silk, J. 1995, ApJ, 451, 44

Parker, E. N. 1953, ApJ, 117, 431

Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, AJ, 139, 1468

Patil, A. H., Yatawatta, S., Zaroubi, S., et al. 2016, MNRAS, 463, 4317

Peebles, P. J. E. 1968, ApJ, 153, 1

Peebles, P. J. E. 1980, The large-scale structure of the universe (Princeton University Press)

Pfenniger, D., & Puy, D. 2003, A&A, 398, 447

Philip, L., Abdurashidova, Z., Chiang, H. C., et al. 2019, J. Astron. Instrum., 8, 1950004

Pitrou, C., Coc, A., Uzan, J.-P., & Vangioni, E. 2018, Phys. Rep., 754, 1

Pitrou, C., Coc, A., Uzan, J.-P., & Vangioni, E. 2021, MNRAS, 502, 2474

Planck Collaboration VI. 2020, A&A, 641, A6

Polidan, R. S., Burns, J. O., Ignatiev, A., et al. 2024, Adv. Space Res., 74, 528

Prasad, S. S., & Huntress, W. T., Jr. 1980, ApJS, 43, 1

Price, D. C., Greenhill, L. J., Fialkov, A., et al. 2018, MNRAS, 478, 4193

Pritchard, J. R., & Loeb, A. 2012, Rept. Prog. Phys., 75, 086901

Puy, D., & Signore, M. 1996, A&A, 305, 371

Puy, D., & Signore, M. 1997, New Astron., 2, 299

Puy, D., Alecian, G., Le Bourlot, J., Leorat, J. & Pineau Des Forets, G. 1993, A&A, 267, 337

Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019

Roberge, W., & Dalgarno, A. 1982, ApJ, 255, 489

Safranek-Shrader, C., Bromm, V., & Milosavljević, M. 2010, ApJ, 723, 1568

Saslaw, W. C., & Zipoy, D. 1967, Nature, 216, 976

Savin, D. W., Krstić, P. S., Haiman, Z., & Stancil, P. C. 2004, ApJ, 606, L167

Schleicher, D. R. G., Galli, D., Palla, F., et al. 2008, A&A, 490, 521

Schöneberg, N. 2024, J. Cosmol. Astropart. Phys., 2024, 006

Seager, S., Sasselov, D. D., & Scott, D. 1999, ApJ, 523, L1

Seager, S., Sasselov, D. D., & Scott, D. 2000, ApJS, 128, 407

Shi, Y., Deng, F., Xu, Y., et al. 2022, ApJ, 929, 32

Signore, M., & Puy, D. 2009, Eur. Phys. J. C, 59, 117

Solomon, P. M., & Klemperer, W. 1972, ApJ, 178, 389

Stancil, P. C., & Dalgarno, A. 1998, Farad. Disc., 109, 61

Stancil, P. C., Lepp, S., & Dalgarno, A. 1996, ApJ, 458, 401

Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1

Swarup, G., Ananthakrishnan, S., Kapahi, V. K., et al. 1991, Curr. Sci., 60, 95

Takayanagi, K., & Nishimura, S. 1960, PASJ, 12, 77

Takeda, H., Satō, H., & Matsuda, T. 1969, Prog. Theor. Phys., 41, 840

Tang, C.-Y., & Chen, K.-J. 2024, MNRAS, 529, 4248

Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1

Tennyson, J. 1995, Rep. Progr. Phys., 58, 421

Teyssier, R. 2002, A&A, 385, 337

Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium
(Cambridge, UK: Cambridge University Press)

Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007

Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Vonlanthen, P. 2009, Ph.D. Thesis, Université de Montpellier, France

Vonlanthen, P., Rauscher, T., Winteler, C., et al. 2009, A&A, 503, 47

Weymann, R. 1965, Phys. Fluids, 8, 2112

Workman, R. L., & Particle Data Group 2022, Prog. Theor. Exp. Phys., 2022, 083C01

Wouthuysen, S. A. 1952, AJ, 57, 31

Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645

Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6

Yoshida, N., Omukai, K., & Hernquist, L. 2007, ApJ, 667, L117

Yoshii, Y., & Sabano, Y. 1980, PASJ, 32, 229

Zeldovich, Y. B., Kurt, V. G., & Syunyaev, R. A. 1968, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 55, 278

Zygelman, B. 2005, ApJ, 622, 1356

Zygelman, B., Stancil, P. C., & Dalgarno, A. 1998, ApJ, 508, 151
# Appendix A: Details on molecular thermal function

In this appendix, we detail the heating and cooling terms composing the molecular thermal function:

$$ \Psi_{\text{mol}} = \Gamma_{\text{mol}} - \Lambda_{\text{mol}}. \qquad (A.1) $$

## A.1. Molecular heating

The collisional de-excitation probability expresses the ratio between the collisional de-excitation rate and the total de-excitation rate, which includes collisions, as well as induced and spontaneous emissions, as follows:

$$ P_{J',J}^{c} = \frac{n_{\xi} X_{J'} k_{J',J}^{\xi}}{n_{\xi} X_{J'} k_{J',J}^{\xi} + X_{J'} A_{J',J} + X_{J'} B_{J',J} \rho_{J',J}}. \qquad (A.2) $$

Here, $n_\xi$ is the numerical density of the collision partner species, $\xi$, and $X_{J'}$ is the population fraction of molecule in the rotational level $J'$. Then, $k_{J',J}^\xi$ is the collision rate of this reaction, while $A_{J',J}$ and $B_{J',J}$ are the Einstein coefficients for spontaneous and induced emission, respectively, and $\rho_{J',J}$ is the radiation energy density. Each term of the sum is the de-excitation term associated to one of the three mentioned processes.

Taking the product between the molecule density at $J$ level, $n_{mol}X_J$, the induced radiative excitation probability, $B_{J',J'}\rho_{J',J'}$, the collisional de-excitation probability, $P_{J',J}^c$, and the energy difference between the two levels, $\Delta\epsilon_{J',J} = \epsilon_{J'} - \epsilon_J$, we can compute the energy gained by the medium per unit of volume and time through molecular thermal influence:

$$ \Gamma_{\text{mol}} = \sum_{J} \sum_{J'} n_{\text{mol}} X_J B_{J',J'} \rho_{J',J'} P_{J',J}^{c} \Delta\epsilon_{J',J}. \qquad (A.3) $$

We summed over all the $J'$ levels reachable by excitation as well as over $J$ to take into account the contribution from all the rotational levels.

The energy of a rotational level, $J$, can be computed as:

$$ \epsilon_J = hBJ(J+1) \text{ where } B \text{ is the rotation constant } B = \frac{h}{8\pi^2 c I}, \qquad (A.4) $$

with $I$ the momentum of inertia.

## A.2. Molecular cooling

By analogy, we first compute the radiative de-excitation probability:

$$ P_{J',J}^{r} = \frac{X_{J'} A_{J',J} + X_{J'} B_{J',J} \rho_{J',J}}{n_{\xi} X_{J'} k_{J',J}^{\xi} + X_{J'} A_{J',J} + X_{J'} B_{J',J} \rho_{J',J}}, \qquad (A.5) $$

which is the ratio of both induced and spontaneous radiative de-excitation rates, to the total one.

We change the radiative excitation probability from Eq. A.3 to the collisional excitation probability with species $\xi: n_\xi k_{J',J}^\xi$ and change the collisional de-excitation probability to the radiative one (Eq. A.5). We can then compute the energy lost by the medium per unit of volume and time through molecular thermal influence:

$$ \Lambda_{\text{mol}} = \sum_{J} \sum_{J'} n_{\text{mol}} n_{\xi} k_{J',J}^{\xi} X_{J} P_{J',J}^{r} \Delta\epsilon_{J',J}. \qquad (A.6) $$
# Appendix B: Molecular rotational levels

This appendix is mainly inspired by Puy & Signore (1997) and Vonlanthen (2009). The aim is to compute, $X_J$, the population fraction of molecules in each rotational level, $J$.

## B.1. Einstein coefficients

Let us assume a container inside which there is a radiation field of density $\rho(\nu)d\nu$ and molecules capable of passing from a rotational level $J$ to a higher level $J'$ by the absorption of a photon whose energy $h_p\nu_{J,J'}$ corresponds to the energy difference between the two levels (radiative excitation) or by collision (collisional excitation). De-excitation can also be collisional or radiative (spontaneous or induced).

There are three radiative processes: spontaneous emission, induced emission, and absorption. The Einstein coefficients $A_{J',J}$ and $B_{J,J'}$ are defined as follows:

* **Spontaneous emission:** the probability of spontaneous emission $A_{J',J}$ is the probability that a molecule in the rotational state $J'$ spontaneously emits a photon of energy $h_p\nu_{J,J'}$ corresponding to the transition from $J'$ to $J$.

* **Induced emission:** The probability of induced emission is defined by the probability that the atom passes from $J'$ to $J$ when it is subjected to induced radiation with a frequency between $\nu_{J'}$ and $\nu_{J'} + d\nu$. This probability is denoted $B_{J',J}\rho(\nu)$.

* **Absorption:** absorption occurs when a molecule in the $J$ state, exposed to isotropic radiation of density $\rho(\nu)d\nu$ and frequency between $\nu_{J'}$ and $\nu_{J'} + d\nu$, absorbs a photon of energy $h_p\nu_{J,J'}$ and changes to the $J'$ state. We write the probability of absorption $B_{J,J'}\rho(\nu)$.

If we write $I$ the moment of inertia of a molecule, then the energy of the rotational level $J$ is given by Cohen-Tannoudji et al. (1973) and Kutner (1984):

$$ \epsilon_J = \frac{J(J+1)h_p^2}{8\pi^2 I}, \qquad (B.1) $$

with $h_p$ the Planck constant. We can rewrite this as:

$$ \epsilon_J = h_p BJ(J + 1), \qquad (B.2) $$

where $B = h_p/8\pi^2 I$ is called the molecule's rotation constant. We can see that the more massive a molecule is, the greater its moment of inertia will be, and therefore the smaller its rotation constant will be. Transition frequencies are obtained by taking the energy differences between two levels and dividing by $h_p$. For dipole molecules, transitions obey the rule $|\Delta J| = 1$, so that $J' = J + 1$. In this case we have:

$$ \nu_{J+1,J} = \nu_{J+1} - \nu_J = \frac{\epsilon_{J+1} - \epsilon_J}{h_p} = 2B(J+1). \qquad (B.3) $$

For homopolar molecules like $H_2$, the transition rule is $|\Delta J| = 2$. Hence we get

$$ \nu_{J+2,J} = \nu_{J+2} - \nu_J = \frac{\epsilon_{J+2} - \epsilon_J}{h_p} = 2B(2J + 3). \qquad (B.4) $$

The Einstein coefficient, $A_{J+1,J}$, can be related to the dipole moment, $\mu_r$ (Kutner 1984):

$$ A_{J+1,J} = \frac{64\pi^4 \nu_{J+1,J}^3}{3h_p c^3} \mu_r^2 \frac{J+1}{2J+3}. \qquad (B.5) $$

Furthermore, the coefficients $B_{J+1,J}$ and $A_{J+1,J}$ are related by the following formula:

$$ B_{J+1,J} = A_{J+1,J} \frac{c^3}{8\pi h_p \nu_{J+1,J}^3}. \qquad (B.6) $$

Finally, the probabilities of absorption and stimulated emission are also related:

$$ g_J B_{J,J+1} = g_{J+1} B_{J+1,J}, \qquad (B.7) $$

where the $g_i$ are the statistical weights of the rotational levels $i$: $g_i = 2i + 1$.

In the case of collisional transitions, the de-excitation probability $C_{J',J}$ is computed by the product:

$$ C_{J',J} = \tau_{\text{coll}} n_{\text{coll}}, \qquad (B.8) $$
where $\tau_{\text{coll}}$ is the collision rate between the molecule and the collisional specie whose density is denoted by $n_{\text{coll}}$. Assuming a Maxwellian velocity distribution, the collisional excitation probability $C_{J,J'}$ is related to $C_{J',J}$ by

$$C_{J,J'} = C_{J',J} \frac{g_{J'}}{g_J} \exp(h_p \nu_{J+1,J} / k_B T_m). \qquad (B.9)$$

From the above considerations we conclude that two quantities mainly determine the efficiency of the Ψ thermal function of a molecular species: the dipole moment $\mu_r$ and the rotation constant $B$. Since the Einstein coefficients depend on the square of the dipole moment, molecules for which $\mu_r$ is large will be very efficient thermal agents. However, of the two most abundant molecules in the cosmic gas, H$_2$ does not have a permanent dipole moment, and the one of HD is very weak. For H$_2$, rotational constant is $B_{\text{H}_2} = 85.35$ K with no dipole moment $\mu_r(\text{H}_2) = 0$ D, while for HD we have $B_{\text{HD}} = 64.30$ K and the dipole moment $\mu_r(\text{HD}) = 8.3 \times 10^{-4}$ D.

## B.2. Populations of the rotational levels

To calculate the thermal function of a molecule, we need to know the populations of its rotational levels $X_J$. When molecules have permanent electric dipole moments (along the axis of the molecule), the strongest transitions are those that obey the selection rule $|\Delta J| = 1$. In the case of homopolar molecules such as molecular hydrogen, the dipole moment is zero for symmetry reasons. Transitions between rotational levels are therefore quadrupolar. This type of transition imposes the rule $|\Delta J| = 2$. In this case, transitions between even-numbered levels $J = 0, 2, 4, ...$ are called para transitions. Those occurring between odd levels $J = 1, 3, 5, ...$ are called ortho transitions. We consider 20 rotational levels in our calculations.

### H$_2$ Molecule
H$_2$ is the only primordial homopolar molecule. Its rotational constant is equal to $B_{\text{H}_2} = 170.66$ K. The radiative transition probability is given by

$$A_{J,J-2} = 7.52 \cdot 10^{-13} \frac{J(J-1)(2J-1)^4}{2J+1} \text{ s}^{-1}. \qquad (B.10)$$

Let's consider three consecutive rotational levels $J-2, J$ and $J+2$. The variation of the $J$ level population is:

$$\frac{dX_J}{dt} = \underbrace{C_{J-2,J}X_{J-2} + B_{J-2,J}\rho_{J-2,J}X_{J-2} + C_{J+2,J}X_{J+2} + B_{J+2,J}\rho_{J+2,J}X_{J+2} + A_{J+2,J}X_{J+2}}_{\text{Gain}} - \underbrace{(C_{J,J+2}X_J + B_{J,J+2}\rho_{J,J+2}X_J + C_{J,J-2}X_J + B_{J,J-2}\rho_{J,J-2}X_J + A_{J,J-2}X_J)}_{\text{Loss}}, \qquad (B.11)$$

The first two terms of the second member in the RHS represent the transition from level $J$ to level $J+2$. The last three terms express the transition from $J$ to $J-2$. The $J$ level can therefore be depopulated either by transition to the $J+2$ level or by de-excitation to the lower $J-2$ level. The transitions are almost instantaneous, leading to a pseudo-stationary state $\frac{dX_J}{dt} = 0$. After rearranging the terms we get:

$$\begin{aligned} & C_{J+2,J}X_{J+2} + A_{J+2,J}X_{J+2} + B_{J+2,J}\rho_{J+2,J}X_{J+2} - C_{J,J+2}X_J - B_{J,J+2}\rho_{J,J+2}X_J \\ & = C_{J,J-2}X_J + A_{J,J-2}X_J + B_{J,J-2}\rho_{J,J-2}X_J - C_{J-2,J}X_{J-2} - B_{J-2,J}\rho_{J-2,J}X_{J-2}. \end{aligned} \qquad (B.12)$$

The symmetry of this last expression for the coefficients $J-2, J$ and $J+2$ allows us to conclude that the coefficients $J$ are independent for each member of the equality. We then obtain:

$$C_{J+2,J}X_{J+2} + A_{J+2,J}X_{J+2} + B_{J+2,J}\rho_{J+2,J}X_{J+2} - C_{J,J+2}X_J - B_{J,J+2}\rho_{J,J+2}X_J = \text{cst.} \qquad (B.13)$$

Through recurrence, we can get back to the fundamental level $J=0$:

$$C_{2,0}X_2 + A_{2,0}X_2 + B_{2,0}\rho_{2,0}X_2 - C_{0,2}X_0 - B_{0,2}\rho_{0,2}X_0 = \text{cst.} \qquad (B.14)$$

And according to Eq. B.11, we also have the pseudo-stationarity hypothesis:

$$0 = \frac{dX_0}{dt} = C_{2,0}X_2 + A_{2,0}X_2 + B_{2,0}\rho_{2,0}X_2 - C_{0,2}X_0 - B_{0,2}\rho_{0,2}X_0. \qquad (B.15)$$

Comparing these two expressions, we conclude that the constant in Eq. B.14 is null. Equation B.13 leads to a simple relation between the populations $X_J$ and $X_{J+2}$:

$$X_J = X_{J+2} \frac{C_{J+2,J} + A_{J+2,J} + B_{J+2,J}\rho_{J+2,J}}{C_{J,J+2} + B_{J,J+2}\rho_{J,J+2}}. \qquad (B.16)$$

In addition, the energy density of cosmological blackbody radiation at the transition frequency $\nu_{J+2,J}$ is given by

$$\rho_{J+2,J} = \rho_{J,J+2} = \frac{8\pi h_p \nu_{J+2,J}^3}{c^3} \frac{1}{\exp\left(\frac{h_p \nu_{J+2,J}}{k_B T_J}\right) - 1}. \qquad (B.17)$$
This relation helps us to simplify Eq. B.16:

$$X_{J+2} = X_J \frac{C_{J,J+2} \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \frac{2J+5}{2J+1}}{C_{J+2,J} \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right)}, \quad (B.18)$$

Where we introduced the transition temperature $T_{J,J+2}$ such that $k_B T_{J,J+2} = h_p \nu_{J,J+2}$. With the relation (Eq. B.9) between the collision probabilities, we get to the following expression:

$$X_{J+2} = X_J \frac{C_{J,J+2} \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \frac{2J+5}{2J+1}}{C_{J,J+2} \frac{2J+5}{2J+1} \exp\left(\frac{T_{J+2,J}}{T_m}\right) \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right)}, \quad (B.19)$$

We introduce the probability $X_{J+2} = a_{J+2} X_J$, where:

$$a_{J+2} := \frac{C_{J,J+2} \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \frac{2J+5}{2J+1}}{C_{J,J+2} \frac{2J+5}{2J+1} \exp\left(\frac{T_{J+2,J}}{T_m}\right) \left( \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right) - 1 \right) + A_{J+2,J} \exp\left(\frac{T_{J,J+2}}{T_\gamma}\right)}. \quad (B.20)$$

This gives us the formulas for calculating the J level population:

$$X_{2n} = \left( \prod_{i=1}^{n} a_{2i} \right) X_0 \quad \text{for para transitions,} \quad (B.21)$$

$$X_{2n+1} = \left( \prod_{i=1}^{n} a_{2i+1} \right) X_1 \quad \text{for ortho transitions.} \quad (B.21)$$

Considering that the sum over the populations of all the rotational levels is normalized, the final result for the rotational level populations is the following:

$$X_0 = \frac{1}{1 + \sum_{n=1}^{\infty} \left( \prod_{i=1}^{n} a_{2i} \right)}, \quad (B.22)$$

$$X_1 = \frac{1}{1 + \sum_{n=1}^{\infty} \left( \prod_{i=1}^{n} a_{2i+1} \right)}, \quad (B.23)$$

$$X_{2n} = \frac{\prod_{i=1}^{n} a_{2i}}{1 + \sum_{n=1}^{\infty} \left( \prod_{i=1}^{n} a_{2i} \right)}, \quad (B.24)$$

$$X_{2n+1} = \frac{\prod_{i=1}^{n} a_{2i+1}}{1 + \sum_{n=1}^{\infty} \left( \prod_{i=1}^{n} a_{2i+1} \right)}, \quad (B.25)$$

where n is a non-zero integer.

HD molecule. A similar computation to the homopolar case gives the following results for the rotational level populations:

$$X_0 = \frac{1}{1 + \sum_{J=1}^{\infty} \left( \prod_{i=1}^{J} a_i \right)}, \quad (B.26)$$

$$X_J = \frac{\prod_{i=1}^{J} a_i}{1 + \sum_{J=1}^{\infty} \left( \prod_{i=1}^{J} a_i \right)}, \quad J \neq 0. \quad (B.27)$$

The proportionality factors are now:

$$a_{J+1} = \left[ C_{J,J+1} \left( \exp\left(\frac{T_{J,J+1}}{T_\gamma}\right) - 1 \right) + A_{J+1,J} \frac{2J+3}{2J+1} \right] \times \left[ C_{J,J+1} \frac{2J+3}{2J+1} \exp\left(\frac{T_{J+1,J}}{T_m}\right) \left( \exp\left(\frac{T_{J,J+1}}{T_\gamma}\right) - 1 \right) + A_{J+1,J} \exp\left(\frac{T_{J,J+1}}{T_\gamma}\right) \right]^{-1}. \quad (B.28)$$