# Molecular Chemistry for Dark Matter II: Recombination, Molecule Formation,
and Halo Mass Function in Atomic Dark Matter

James Gurian<EMAIL_ADDRESS>Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Department of
Astronomy and Astrophysics, The Pennsylvania State University, University
Park, PA, 16802, USA  Donghui Jeong<EMAIL_ADDRESS>Institute for Gravitation
and the Cosmos, The Pennsylvania State University, University Park, PA 16802,
USA  Department of Astronomy and Astrophysics, The Pennsylvania State
University, University Park, PA, 16802, USA  School of Physics, Korea
Institute for Advanced Study (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul, 02455,
Republic of Korea  Michael Ryan<EMAIL_ADDRESS>Institute for Gravitation and
the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
Department of Physics, The Pennsylvania State University, University Park, PA,
16802, USA  Sarah Shandera<EMAIL_ADDRESS>Institute for Gravitation and the
Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
Department of Physics, The Pennsylvania State University, University Park, PA,
16802, USA

###### Abstract

Dissipative dark matter predicts rich observable phenomena that can be tested
with future large-scale structure surveys. As a specific example, we study
atomic dark matter, consisting of a heavy particle and a light particle
charged under a dark electromagnetism. In particular, we calculate the
cosmological evolution of atomic dark matter focusing on dark recombination
and dark-molecule formation. We have obtained the relevant interaction-rate
coefficients by re-scaling the rates for normal hydrogen, and evolved the
abundances for ionized, atomic, and molecular states using a modified version
of Recfast++ (which we have released publicly
111https://github.com/jamesgurian/RecfastJulia). We also provide an analytical
approximation for the final abundances. We then calculate the effects of the
atomic dark matter on the linear power spectrum, which enter through a dark-
photon diffusion and dark acoustic oscillations. At the formation time, the
atomic dark matter model suppresses halo abundances on scales smaller than the
diffusion scale, just like the warm dark matter models suppress the abundance
below the free-streaming scale. The subsequent evolution with radiative
cooling, however, will alter the halo mass function further.

cosmology: theory – dark matter – molecular processes

††software: Recfast++ (Chluba et al., 2010), Julia (Bezanson et al., 2012),
DifferentialEquations library (Rackauckas & Nie, 2017)

##  1 Introduction

About 84% of the matter in the universe is dark, revealing itself only through
its gravitational charge, which is its mass. Various large-scale structure
observation campaigns (Bennett et al., 2013; Aghanim et al., 2020; Alam et
al., 2021) show that the observed spectrum of large-scale structure is
consistent with the cold dark matter (CDM) model, and detailed study of the
Bullet Cluster collision shows that the dark matter self-interaction is much
weaker than gas particles (Markevitch et al., 2004). There are small-scale
phenomena such as cored density-profiles of dwarf galaxies (Ostriker, 2003; de
Blok, 2010), the missing satellite problem (Diemand et al., 2008; Springel et
al., 2008; Bullock, 2010) and the too-big-to-fail problem (Boylan-Kolchin et
al., 2011, 2012) that might suggest a deviation from prevailing CDM model,
such as self-interacting dark matter (Spergel & Steinhardt, 2000) or warm dark
matter with non-negligible velocity dispersion (Bode et al., 2001; Sommer-
Larsen & Dolgov, 2001; Hogan & Dalcanton, 2000). None of them, however, have
conclusive evidence, and more accurate accounting for baryonic physics can
also yield some of these phenomena (see Brooks (2019) for a comprehensive
review).

On the other hand, despite decades of searching through direct detection
(Agnese et al., 2018; Collaboration et al., 2013; Akerib et al., 2017),
indirect detection (Albert et al., 2017; Bergström et al., 2013) and pair
creation (Boveia & Doglioni, 2018; Mitsou, 2015), there is no definitive
evidence that dark matter interacts with Standard Model particles. It may be
that the dark matter is for all practical purposes completely secluded from
the Standard Model. Although this prospect would doom many current detection
programs, it still leaves open the possibility of a rich dark sector
consisting of one or more new particles interacting via new forces. Because
such a model predicts distinctive dark-matter halo structure, it is possible
to constrain the dark matter self interactions using only the gravitational
interaction which the dark matter must possess (Buckley & Peter, 2018).

Here, we consider the particular possibility often referred to as “atomic dark
matter” (Goldberg & Hall, 1986; Ackerman et al., 2009; Feng et al., 2009;
Kaplan et al., 2010, 2011; Cyr-Racine & Sigurdson, 2013; Cyr-Racine et al.,
2014; Fan et al., 2013a; Cline et al., 2014; Boddy et al., 2016; Foot &
Vagnozzi, 2015, 2016; Randall & Scholtz, 2015; Agrawal et al., 2017a; Ghalsasi
& McQuinn, 2018). In this model, the dark matter consists of a heavy and a
light fermion charged under a “dark electromagnetism”. This dark matter can
form hydrogen-like bound states including molecules. From one perspective,
this can be seen as a toy model which captures some of the essential physics
of a more complexly interacting dark sector. This model can satisfy
constraints on the dark matter self-interaction cross section from, for
example, colliding galaxy clusters and dwarf galaxy density profiles (Agrawal
et al., 2017b). In fact, dark matter models which introduce a mechanism to
suppress the linear matter power spectrum on small scales (including warm dark
matter and atomic dark matter) may alleviate challenges facing the
$\Lambda$CDM model on small scales (Bullock & Boylan-Kolchin, 2017).

The parameters of the atomic-dark-matter model are the mass of the heavy
particle $M$, that of the light particle $m$, the coupling constant $\alpha$
and the ratio of the dark-photon temperature to the cosmic microwave
background (CMB) temperature, $\xi=T_{\gamma,{\rm D}}/T_{\gamma}$.
Observations of Big-Bang Nucleosynthesis (BBN) and the CMB temperature
anisotropies constrain the energy density in relativistic species other than
the three known neutrino generations, parameterized as the change in the
effective number of neutrino species $\Delta N_{\rm eff}$. The most recent
such constraint is $\Delta N_{\rm eff}\lesssim 0.2$ (Ade et al., 2016), which
restricts $\xi\lesssim 0.4$. This upper limit is easy to accommodate since the
dark-photon temperature need not be the same as the CMB temperature. In some
regions of the parameter space, the atomic dark matter model allows formation
of a double disk (Fan et al., 2013b), but this high dark-photon temperature
region is strongly constrained by the lack of dark acoustic oscillations (DAO)
in the observed galaxy power spectrum (Cyr-Racine et al., 2014).

One particularly interesting aspect of atomic dark matter is that the dark
matter is dissipative: atomic dark matter halos can radiate away kinetic
energy. This means that these halos can undergo collapse, fragmentation, and
can eventually form dark compact objects. This raises the prospect of probing
an atomic dark matter model through gravitational wave events from their
compact-binary coalescence (Shandera et al., 2018). Crucially, the
Chandrasekhar mass $M_{\rm Ch}\propto 1/M^{2}$ only depends on the mass of the
dark ‘proton’, and it can be much smaller than the nominal value $M_{\rm
Ch}\simeq 1.4M_{\odot}$ in the Standard Model. That is, the mass spectrum of
merger events, particularly for the sub-solar mass region, can carry
information about the physics of even a fully secluded dark sector. Shandera
et al. (2018) has showed that a range of model parameters can produce merger
rates accessible to the current generation of advanced-LIGO detectors. More
recently, Singh et al. (2020) has derived constraints on the dark
Chandrasekhar mass and showed that $M>0.95\,{\rm GeV}$ with 99.9% C.L., if the
unusual, so called mass-gap, gravitational wave event GW190425 (Abbott et al.,
2020) is interpreted as the merger of an atomic-dark-matter black hole binary.

So far, both Shandera et al. (2018) and Singh et al. (2020) have relied on
estimates of the dark black hole binary parameters obtained by re-scaling the
Population III star binary literature (Abel, 2001; Bromm & Larson, 2004;
Hartwig et al., 2016). To obtain more realistic constraints, we must tackle
directly the complex physics of the formation and evolution of dark-matter
structures in the atomic-dark-matter model. As a stepping stone toward that
complete study, in this work we consider the cosmological evolution of atomic
dark matter. Cyr-Racine & Sigurdson (2013) have previously studied the
recombination of dark hydrogen, which we extend by including the formation of
molecular states. In the Standard Model, molecular hydrogen is the dominant
coolant for pristine (zero metallicity) minihalos with virial temperature
below $10^{4}\,{\rm K}$. So, studying the molecular abundance evolution must
be a crucial step in determining the dark-halo properties within the atomic-
dark-matter paradigm.

Another essential ingredient for these studies is the modification of the
linear matter power spectrum introduced by the atomic dark matter.
Specifically, atomic dark matter produces dark acoustic oscillations (DAO) due
to the propagation of pressure waves in the tightly coupled dark-matter/dark-
photon fluid, as well as dark diffusion damping due to the non-zero mean free
path of the dark photon. Both effects lead to a suppression of power on small
scales, as discussed in detail in Cyr-Racine & Sigurdson (2013). While Cyr-
Racine & Sigurdson (2013) studied a relatively high value for the temperature
ratio, most extensively $\xi=0.37$, we expand the study by varying $\xi$ in
particular to a small value, $\xi<0.1$. We have also added the molecular
processes. We have shown that although molecular processes do change the
evolution of the ionization fraction, as long as dark molecular hydrogen is
rare relative to dark atoms, the effect is too small to alter the atomic dark
matter effects on the linear matter power spectrum.

We find that varying $\xi$ impacts all effects studied. At a given dark photon
temperature, this parameter alters the redshift and hence the dark matter
number density at a given temperature, which controls freeze-out and molecule
formation. Similarly, the diffusion and DAO scales depend on the horizon size
at decoupling, which depends on $\xi$ through the decoupling redshift.

We finally investigate the effects of atomic dark matter on the halo mass
function in the presence of a cutoff in the linear power spectrum. This
approach has been validated by simulations of warm dark matter (Schneider et
al., 2013; Angulo et al., 2013) and is an appropriate ansatz for the formation
of atomic-dark-matter halos at high redshift.

This paper is the second in a series of three. In Paper I, Ryan et al.
(2022a), we defined a re-scaling procedure to compute the molecular chemistry
for atomic dark matter from known Standard Model results. Paper III, Ryan et
al. (2022b), describes DarkKROME, a modification of the chemistry package
KROME (Grassi et al., 2014) to include the chemistry of atomic dark matter. In
Paper III we use DarkKROME  and the cosmological abundances derived in this
paper to explore one-zone collapse scenarios for a range of atomic dark matter
parameters.

This paper is organized as follows. In Sec. 2, we present the numerical
implementation of the cosmological evolution of atomic dark matter states
using Recfast++ modified by including the relevant molecular processes (Sec.
2.1), the computational method for event rates of these reactions (Sec. 2.2),
and the result and an analytic approximation (Sec. 2.5). In Sec. 3, we study
the effects of the addition of atomic dark matter and molecules on the
formation of large-scale structure by calculating the dark diffusion scale,
the DAO scale (Sec. 3.1), and finally estimating the halo mass function at its
formation time (Sec. 3.2). We conclude Sec. 4 by pointing out the path toward
the complete study of the large-scale structure with atomic dark matter.

Throughout, we shall attach the subscript ${\rm D}$ to the Standard Model
symbols to refer to the corresponding dark-matter symbols. For example, ${\rm
e}_{\rm D}+{\rm p}_{\rm D}\to{\rm H}_{\rm D}+\gamma_{\rm D}$ refers to the
dark-recombination reaction, the dark-matter analogue of the Standard Model
reaction: ${\rm e}+{\rm p}\to{\rm H}+\gamma$.

##  2 Dark Recombination and Molecule Formation

###  2.1 The rate equations and implementation

We begin by briefly summarizing the apparatus necessary to solve for the
chemical evolution of atomic-dark-matter states, including the recombination
and molecule formation, in an expanding universe. To do so, we have
implemented the core222That is, we have not implemented extra heating sources
due to dark-matter annihilation (Chen & Kamionkowski, 2004; Padmanabhan &
Finkbeiner, 2005) or magnetic fields (Sethi & Subramanian, 2005; Kunze &
Komatsu, 2014). of Recfast++ (Seager et al., 1999; Chluba et al., 2010) in
Julia (Bezanson et al., 2012), using the DifferentialEquations library
(Rackauckas & Nie, 2017). We present a few details necessary for applying this
code, which computes the Standard Model recombination, to the atomic-dark-
matter scenario in Appendix A.

The Recfast++ code uses Peeble’s effective three-level recombination (Peebles,
1968) approximation, where the recombination is treated entirely by tracking
the non-equilibrium population in the lowest three ($1s$, $2s$ and $2p$)
atomic states and the continuum. High optical depth in the early universe to
Lyman-limit photons prevents the direct recombination to the ground (1s) state
(case-B recombination), so the two $n=2$ states ($2s$ and $2p$) are populated
as meta-stable states. The $n=2$ states are populated with the rate given by
case-B recombination coefficient $\alpha_{H_{\rm D}}$ that includes
recombination to higher levels followed by a cascade down to $n=2$. They leave
the $n=2$ state by one of three processes: (A) the forbidden two photon decay
($2s\to 1s$), (B) redshifting of Lyman-$\alpha$ photons out of the resonance
band, and (C) photoionization. Exploiting the meta-stability (chemical
equilibrium) of the $n=2$ state, $dx_{2}/dt=0$, it is possible to solve for
the population of that state and then for the net change in the ionization
fraction. This leads to the following equation (Seager et al., 1999), whose
right-hand-side terms (in square brackets) represent (i) a change of variables
from $t$ to $z$ (time to redshift), (ii) a branching ratio between the ground-
state recombinations (two photon decay and redshifted Lyman-$\alpha$) and all
exits from the $n=2$ state, including photoionization, and (iii) the net
Case-B recombination rate:

${\frac{dx_{{\rm p}_{\rm
D}}}{dz}}=\left[-\frac{1}{H(z)(1+z)}\right]\left[\frac{\Gamma_{\alpha,{\rm
D}}+A_{2\gamma,{\rm D}}}{\Gamma_{\alpha,{\rm D}}+A_{2\gamma,{\rm
D}}+\beta_{{\rm H}_{\rm D}}}\right]\left[\beta_{{\rm H}_{\rm D}}x_{{\rm
H}_{{\rm D}}\rm I}{e}^{-h\nu_{{\rm Ly\alpha},D}/k_{B}T_{g,{\rm D}}}-x_{{\rm
e}_{\rm D}}x_{{\rm p}_{\rm D}}n_{{\rm H}_{\rm D}}\alpha_{{\rm H}_{\rm
D}}\right],$ (2.1)

with the quantities appearing in the equation defined as follows: $k_{B}$ is
the Boltzmann constant, $h$ is the Planck’s constant, $T_{\rm g,D}$ is the
dark-matter gas temperature and $T_{\gamma,D}$ is the dark-photon temperature,
$x_{{\rm p}_{\rm D}}=n_{{\rm p}_{\rm D}}/n_{{\rm H}_{\rm D}}$ is the free
proton fraction, $x_{{\rm e}_{\rm D}}$ the free electron fraction, $x_{{\rm
H}_{\rm D}\rm I}$ is the neutral hydrogen fraction. Note that $x_{{\rm e}_{\rm
D}}$ and $x_{{\rm p}_{\rm D}}$ differ because we also include the complex
states ${\rm H}_{\rm D}^{-}$, ${\rm H}_{{\rm D},2}^{+}$, ${\rm H}_{{\rm
D},3}^{+}$ in the chemical network. The free electron fraction is calculated
as $x_{{\rm e}_{\rm D}}\equiv x_{{\rm p}_{\rm D}}-x_{{\rm H}^{-}_{\rm
D}}+x_{{\rm H}_{{\rm D},2}^{+}}+x_{{\rm H}_{{\rm D},3}^{+}}$. We also have
$\alpha_{{\rm H}_{\rm D}}$ the case-B recombination rate, $\Gamma_{\alpha,{\rm
D}}\equiv 8\pi H(z)/\left[n_{{\rm H}_{\rm D}\rm I}\lambda_{{\rm Ly\alpha},{\rm
D}}^{3}\right]$ the rate at which redshifted Lyman-$\alpha_{D}$ photons escape
the resonance band, $\nu_{\rm Ly\alpha_{\rm D}}$, the Lyman-$\alpha_{\rm D}$
frequency, $\beta_{{\rm H}_{\rm D}}$ the photoionization coefficient, and
$A_{2\gamma,{\rm D}}$ the 2-photon decay rate.

The recombination coefficient is given by the fitting formula of Pequignot et
al. (1991a). This formula includes an empirical 1.14 “fudge factor” (Seager et
al., 1999) which accounts for the neglected higher ($n>2$) level transitions.
The fudge factor takes into account all electric dipole transitions from
higher $(n,\ell)$ modes. Because dark hydrogen energy and angular momentum
structure is merely re-scaled (Sec. 2.2) from the Standard Model hydrogen, we
expect the same fudge factor to be applied to the atomic dark matter case. We
calculate the photoionization coefficient $\beta_{{\rm H}_{\rm D}}$ by using
detailed balance: $\beta_{{\rm H}_{\rm D}}=\alpha_{{\rm H}_{\rm D}}(2\pi
mk_{B}T_{\rm g,{\rm D}}/h^{2})^{3/2}\exp(-(E_{{\rm H}_{\rm D}}-h\nu_{{\rm
Ly\alpha},{\rm D}})/kT_{\rm g,{\rm D}})$, with $m$ the mass of the light
fermion and $E_{{\rm H}_{\rm D}}$ the dark hydrogen binding energy. As
recombination proceeds, the dark matter and dark photon fall out of thermal
equilibrium due to the inefficiency of the Compton-scattering-type energy
transfer. This process can be described by solving the following equation for
$T_{\rm g,{\rm D}}$ simultaneously with the equation for $x_{{\rm p}_{\rm D}}$
(Seager et al., 1999):

$\frac{dT_{\rm g,D}}{dz}=\frac{8\sigma_{\rm T,D}a_{\rm D}T_{\gamma,{\rm
D}}^{4}}{3H(z)(1+z)mc}\,\frac{n_{{\rm e}_{\rm D}}}{n_{\rm DM}}\,(T_{\rm
g,D}-T_{\gamma,{\rm D}})+\frac{2T_{\rm g,D}}{(1+z)}.$ (2.2)

where $n_{\rm DM}\equiv n_{{\rm H}_{\rm D}}(x_{{\rm e}_{\rm D}}+x_{{\rm
p}_{\rm D}}+x_{\rm HI}+x_{{\rm H}_{{\rm D},2}}+x_{{\rm H}_{{\rm
D},2}^{+}}+x_{{\rm H}_{{\rm D},3}^{+}}+x_{{\rm H}_{{\rm D}}^{-}})$ is the
total dark-matter particle number density, $\sigma_{\rm T,D}$ is the Thomson
cross section and $a_{\rm D}$ the radiation constant. We discuss the
relationship between these constants in the dark-matter sector and the
Standard-Model values in the following subsection (Sec. 2.2).

We solve the ionization fraction and matter temperature equations [Eqs.
(2.1)–(2.2)] simultaneously with the equations for the evolution of ${\rm
H}_{D}^{-}$, ${\rm H}_{{\rm D},2}^{+}$, ${\rm H}_{{\rm D},3}^{+}$, and ${\rm
H}_{{\rm D},2}$. These reactions also require the addition of terms which we
collectively denote by $\frac{d\Delta x_{\rm p_{\rm D}}}{dz}$ to Eq. (2.1) to
account for the effect of the molecular reactions on the free proton fraction.
The reaction rates, $c_{{\rm H}i}$, appearing below are as summarized in Table
1 ($i$ is the reaction number in the Table), and have been obtained by re-
scaling standard model rates according to the procedure describe in Sec. 2.2
and in Ryan et al. (2022a). The number density of neutral dark atoms is
$n_{{\rm H}_{\rm D}\rm I}\equiv n_{{\rm H}_{\rm D}}(1-x_{{\rm p}_{\rm
D}}-2x_{{{\rm H}_{\rm D,2}}}-2x_{{{\rm H}_{{\rm D},2}^{+}}}-3x_{{\rm H}_{{\rm
D},3}^{+}}-x_{{{\rm H}^{-}_{{\rm D}}}})$. The additional evolution equations
are

$\displaystyle\frac{d\Delta x_{\rm p_{\rm D}}}{dz}=-\frac{-c_{{\rm H}7}n_{{\rm
p}_{\rm D}}x_{{{\rm H}^{-}_{{\rm D}}}}-c_{\rm{\rm H}8}n_{{\rm H}_{D}I}x_{{\rm
p}_{D}}+c_{\rm H9}x_{{\rm H}_{{\rm D},2}^{+}}+c_{\rm H10}x_{{\rm H}_{\rm
D,2}^{+}}n_{{\rm H}_{\rm D}\rm I}-c_{\rm H15}x_{{{\rm H}_{\rm D,2}}}n_{{\rm
p}_{\rm D}}}{H(z)(1+z)}$ (2.3) $\displaystyle\frac{dx_{{\rm H}_{\rm
D}^{-}}}{dz}=-\frac{c_{\rm{\rm H}3}n_{{\rm H}_{D}\rm I}x_{{\rm e}_{\rm
D}}-c_{\rm{\rm H}4}x_{{\rm H}^{-}_{\rm D}}-c_{{\rm H}5}n_{{\rm H}_{\rm D}\rm
I}x_{{{\rm H}^{-}_{{\rm D}}}}-c_{{\rm H}7}n_{{\rm p}_{\rm D}}x_{{{\rm
H}^{-}_{{\rm D}}}}}{H(z)(1+z)}$ (2.4) $\displaystyle\frac{dx_{{{\rm H}_{{\rm
D},2}^{+}}}}{dz}=-\frac{c_{\rm{\rm H}8}n_{{\rm H}_{\rm D}I}x_{{\rm p}_{\rm
D}}-c_{\rm H9}x_{{\rm H}_{{\rm D},2}^{+}}-c_{\rm H10}x_{{\rm H}_{\rm
D,2}^{+}}n_{{\rm H}_{\rm D}\rm I}-c_{\rm H13}x_{{{\rm H}_{\rm D,2}}}n_{{{\rm
H}_{{\rm D},2}^{+}}}+c_{\rm H15}n_{{\rm p}_{\rm D}}x_{{{\rm H}_{\rm
D,2}}}+c_{\rm 3B4}n_{{\rm H}_{\rm D}\rm I}^{2}x_{{\rm p}_{\rm D}}}{H(z)(1+z)}$
(2.5) $\displaystyle\frac{dx_{{\rm H}_{{\rm D},3}^{+}}}{dz}=-\frac{c_{\rm
H13}x_{{{\rm H}_{\rm D,2}}}n_{{{\rm H}_{{\rm D},2}^{+}}}-c_{\rm H20}x_{{\rm
H}_{{\rm D},3}^{+}}n_{\rm e_{D}}}{H(z)(1+z)}$ (2.6)
$\displaystyle\frac{dx_{{{\rm H}_{\rm D,2}}}}{dz}=-\frac{c_{\rm H5}n_{{\rm
H}_{\rm D}\rm I}x_{{\rm H}_{\rm D}^{-}}+c_{\rm H10}x_{{{\rm H}_{{\rm
D},2}^{+}}}n_{{\rm H}_{\rm D}\rm I}-c_{\rm H15}x_{{{\rm H}_{\rm D,2}}}n_{{\rm
p}_{\rm D}}+c_{\rm H20}x_{{\rm H}_{{\rm D},3}}^{+}}{H(z)(1+z)}$ (2.7)
$\displaystyle\qquad\qquad-\frac{-c_{{\rm H}*}n_{{\rm H}_{\rm D}\rm I}x_{{{\rm
H}_{\rm D,2}}}+c_{\rm 3B1}n_{{\rm H}_{\rm D}\rm I}^{2}x_{{\rm H}_{\rm D}\rm
I}+c_{\rm 3B2}n_{{\rm H}_{\rm D}\rm I}^{2}x_{{{\rm H}_{\rm D,2}}}+c_{\rm
3B3}n_{{\rm H}_{\rm D}\rm I}^{2}x_{H_{\rm D}}^{+}}{H(z)(1+z)}.$ (2.8)

# | Reaction | Cross section source | $\sigma$ | Re-scaling pre-factor | b | Additional notes
---|---|---|---|---|---|---
$g(r_{\alpha},r_{m},r_{M})$
1 | $p+e\rightarrow{\rm H}+\gamma$ | Mo et al. (2010) | $\frac{\alpha^{5}}{\text{K.E.}(\text{K.E.}+\Delta E)}$ | $r_{\alpha}^{2}\,r_{m}^{-2}$ | $-0.62,-1.15$ |  a , b
2 | ${\rm H}+\gamma\rightarrow p+e$ | Mo et al. (2010) | $\mu\,\alpha^{5}\frac{1}{(\text{K.E.}+\Delta E)^{3}}$ | $r_{\alpha}^{5}\,r_{m}$ | 0.88, 0.35 | c
3 | ${\rm H}+e\rightarrow{\rm H}^{-}+\gamma$ | de Jong (1972) | $\frac{\alpha}{\mu^{2}}\;\frac{\Delta E^{1/2}\,\text{K.E.}^{1/2}}{(\text{K.E.}+\Delta E)}$ | $r_{\alpha}^{2}\,r_{m}^{-2}$ | 0.928 |  a , d
4 | ${\rm H}^{-}+\gamma\rightarrow{\rm H}+e$ | Armstrong (1963) | $\frac{\alpha}{\mu}\frac{\Delta E^{1/2}K^{3/2}}{(\text{K.E.}+\Delta E)^{3}}$ | $r_{\alpha}^{5}\,r_{m}$ | 2.13 |  c , d
5 | ${\rm H}^{-}+{\rm H}\rightarrow{\rm H}_{2}+e$ | Browne & Dalgarno (1969) | $\sqrt{\frac{\alpha\,a_{0}^{3}}{\text{K.E.}}}$ | $r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}$ | 0 |  e , f , g
7 | ${\rm H}^{-}+p\rightarrow 2{\rm H}$ | Bates & Lewis (1955) | $\alpha\,a_{0}^{2}\sqrt{\mu}\,\frac{\sqrt{\text{K.E.}+\Delta E}}{\text{K.E.}\,\Delta E}$ | $r_{\alpha}^{-3}\,r_{m}^{-3}$ | $-\frac{1}{2}$ |  e , h
8 | ${\rm H}+p\rightarrow{\rm H}_{2}^{+}+\gamma$ | Stancil et al. (1993) | $\frac{({\text{K.E.}}+\Delta E)^{3}\alpha^{4}}{E_{\rm H}^{3}{\text{K.E.}}^{3/2}M^{1/2}}$ | $r_{\alpha}^{2}\,r_{m}^{-1}\,r_{M}^{-1}$ | 1.8 | i
9 | ${\rm H}_{2}^{+}+\gamma\rightarrow{\rm H}+p$ | Stancil et al. (1993) | $\left(\frac{\mu\,v}{h\,\nu}\right)^{2}\,\frac{({\text{K.E.}}+\Delta E)^{3}\alpha^{4}}{E_{\rm H}^{3}{\text{K.E.}}^{3/2}M^{1/2}}$ | $r_{\alpha}^{5}\,r_{m}^{1/2}\,r_{M}^{1/2}$ | 1.59 |  c , j
10 | ${\rm H}_{2}^{+}+{\rm H}\rightarrow{\rm H}_{2}+p$ | Galli & Palla (1998) | $\sqrt{\frac{\alpha\,a_{0}^{3}}{\text{K.E.}}}$ | $r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}$ | 0 | g
13 | ${\rm H}_{2}^{+}+{\rm H}_{2}\rightarrow{\rm H}_{3}^{+}+{\rm H}$ | —"— | —"— | —"— | 0 | g
15 | ${\rm H}_{2}+p\rightarrow{\rm H}_{2}^{+}+{\rm H}$ | —"— | —"— | —"— | 0 |  g , j
20 | ${\rm H}_{3}^{+}+e\rightarrow{\rm H}+{\rm H}_{2}$ | Draine (2011) | $\frac{\alpha\,a_{0}}{\text{K.E.}}$ | $r_{\alpha}^{-1}\,r_{m}^{-2}\,r_{M}$ | -0.65 | k
| ${\rm H}_{2}+{\rm H}\rightarrow 3{\rm H}$ | Hard Sphere | $a_{0}^{2}$ | $r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}$ | 0 | l
3B1 | $3{\rm H}\rightarrow{\rm H}_{2}+{\rm H}$ | Hard Sphere/Detailed Balance | $a_{0}^{2}\left[\frac{n_{{\rm H}_{2}}}{n^{2}_{\rm H}}\right]_{\rm LTE}$ | $r_{\alpha}^{-4}\,r_{m}^{-4}\,r_{M}^{-1}$ | -1 |  j , m
3B2 | ${\rm H}_{2}+2{\rm H}\rightarrow 2{\rm H}_{2}$ | —"— | —"— | —"— | -1 |  j , m
3B3 | $2{\rm H}+{\rm H}^{+}\rightarrow{\rm H}_{2}+{\rm H}^{+}$ | —"— | —"— | —"— | -1 |  j , m
3B4 | $2{\rm H}+{\rm H}^{+}\rightarrow{\rm H}_{2}^{+}+{\rm H}$ | —"— | —"— | —"— | -1 |  j , m

Table 1: Table of Reactions. Reactions are numbered according to Galli & Palla
(1998) and the section numbers in parentheses indicate where details can be
found in the text. Only the parametric dependence on key quantities are shown
for the cross sections (no numerical factors). The reactions included in this
table are those that were considered “important" in Galli & Palla (1998) and
additions discussed in the text. To leading order the Standard Model rates can
be approximated as $\gamma\propto g(r_{\alpha},r_{m},r_{M})\;(T/r_{\Delta
E})^{b}$ and values for $b$ are from Galli and Palla unless otherwise noted.
All reactions listed have binding energy proportional to $E_{\rm H}$, so here
$r_{\Delta E}=r_{\alpha}^{2}\,r_{m}$. For reactions 1-4 and 20, the reduced
mass, $\mu$, is proportional to the dark electron mass, $m$. For the remaining
reactions, $\mu$ is proportional to the dark proton mass, $M$.

  * a

Milne Relation

  * b

The $b$ values are an expansion of the Case B recombination coefficient in low
and high temperature regimes (respectively), from Pequignot et al. (1991b).
The re-scaling of the Case A coefficient is identical, although $b$ differs.

  * c

Photoionization re-scaling

  * d

Effective Range approximation

  * e

These rates are uncertain by up to an order of magnitude Glover & Abel (2008)

  * f

Rate is constant for $T\leq 300{\,\rm K}$, and very uncertain for
$T>300{{\,\rm K}}$ Galli & Palla (1998)

  * g

Langevin Reaction

  * h

Landau-Zener method

  * i

$b$ value from Stancil et al. (1998)

  * j

General Detailed Balance

  * k

Coulomb focusing

  * l

The $b$ value is from Lepp & Shull (1983) and does not account for density
dependence.

  * m

The $b$ value is from Yoshida et al. (2006).

Most reactions in Table 1 are those of the minimal model in Galli & Palla
(1998), and the rates are those compiled in the same reference, except that we
replace the more complicated temperature dependence of H7 and H8 in Galli &
Palla (1998) with the simpler power laws taken from Bates & Lewis (1955) and
Shapiro & Kang (1987), respectively.

###  2.2 Re-scaling the reaction rates

We have re-scaled the recombination rate with different values of $m_{e}$,
$m_{p}$, and $\alpha$ as described in Hart & Chluba (2017):

$\displaystyle\begin{split}&\sigma_{\rm T,D}\propto
r_{\alpha}^{2}r_{m}^{-2}\qquad A_{2\gamma,{\rm D}}\propto
r_{\alpha}^{8}r_{m}\qquad\Gamma_{\alpha,{\rm D}}\propto
r_{\alpha}^{6}r_{m}^{3}\\\ &\alpha_{{\rm H}_{\rm D}}\propto
r_{\alpha}^{2}r_{m}^{-2}\qquad\beta_{{\rm H}_{\rm D}}\propto
r_{\alpha}^{5}r_{m}\,\end{split}$ (2.9)

where $r_{x}$ denotes the ratio of a dark quantity to its standard model
value, for example $r_{M}=M/m_{p}$. Here, the re-scalings of $\sigma_{T,D}$,
$A_{2\gamma,D}$, and $\Gamma_{\alpha D}$ are immediate from the scaling of the
hydrogen binding energy, $E_{{\rm H}}\propto r_{m}r_{\alpha}^{2}$. The
scalings of $\alpha_{{\rm H}_{\rm D}}$ and $\beta_{{\rm H}_{\rm D}}$ are
derived in Ryan et al. (2022a) and Hart & Chluba (2017).

Additionally, we take into account the temperature dependence of the
recombination and photoionization coefficients by re-scaling the temperature
to keep the temperature-to-hydrogen-binding-energy ratio invariant:

$T_{\rm D}=\frac{T}{E_{\rm H}}E_{{\rm H},{\rm D}}=r_{\alpha}^{2}r_{m}T\,.$
(2.10)

When the functional form is given in terms of the Standard-Model temperature,
we re-scale the dark-matter temperature to find the corresponding Standard-
Model temperature by $T=r_{\alpha}^{-2}r_{m}^{-1}T_{D}$ for evaluating the
rate coefficient. This re-scaled temperature is denoted $\tilde{T}_{a}$. For
example, we calculate the recombination coefficient by

$\alpha_{{\rm H}_{\rm D}}(T_{\rm D})=r_{\alpha}^{2}r_{m}^{-2}\alpha_{\rm
H}(r_{\alpha}^{-2}r_{m}^{-1}T_{\rm D})=r_{\alpha}^{2}r_{m}^{-2}\alpha_{\rm
H}(\tilde{T}_{a})\,.$ (2.11)

For the reactions involved in the formation of molecular states we adopt re-
scaling procedure described in detail in Ryan et al. (2022a). The reaction
rates in Table 1 are fits to either theoretical calculations or experimental
data. The form emerges from integrating the cross section over a Maxwell-
Boltzmann distribution (or, for radiative reactions, a Planck distribution)
with a lower energy cutoff corresponding to the reaction threshold. As a
result, the form of the re-scaled rate is given as

$c_{\rm Hi,D}={g(r_{\alpha},r_{m},r_{M})}c_{\rm Hi}(\tilde{T}_{a})$ (2.12)

where $g(r_{\alpha},r_{m},r_{M})$ contains the entire parametric dependence of
the cross section.

The reactions 1 to 15 in Table 1 are considered “important" in Galli & Palla
(1998). One might worry that other reactions deemed unimportant in the
Standard Model by Galli & Palla (1998) could also contribute for nonstandard
parameters. This issue is discussed in more detail in Ryan et al. (2022a).
Briefly, the volumetric rates $c_{\rm Hi}$ of the unimportant reactions may be
small either due to a strong parametric suppression of the rate coefficient
$c{\rm Hi}$ or by the small relative abundance of the reactants. Overcoming
these large order-or-magnitude differences requires either an extreme-valued
re-scaling factor $g(r_{\alpha},r_{m},r_{M})$ outside of the range that we
consider here, or a similarly large modification of the relative number
densities. Such a modification of the relative abundances can occur when $\xi$
is small, where recombination occurs in a strongly collisional environment.
For this reason, we have included the reaction $*$ for the collisional
destruction of ${\rm H}_{2}$, as well as the three body reactions relevant to
the formation of ${\rm H}_{2}$ used in Yoshida et al. (2006). We also include
the dominant ${\rm H}_{{\rm D},3}^{+}$ formation and destruction mechanisms,
which are suppressed for the Standard Model only by the low fractional
abundance of ${\rm H}_{2}$. While reasonable, these modifications may be
insufficient in some parameter regimes. In the Standard Model literature, the
minimal reaction list of Galli & Palla (1998) has been largely validated by
later, exhaustive studies (i.e. Glover & Savin (2009)). Pending similar work
for the dark sector, some caution is warranted when the relative abundances of
the species are inverted by orders of magnitude compared to the Standard Model
values.

###  2.3 Limitations of Re-scaling

As we have discussed in Ryan et al. (2022a), the re-scaling procedure has the
following limitations.

First, the theoretical calculations of the molecular cross sections in the
Standard Model rely on a Born-Oppenheimer (stationary nuclei) approximation
that fails when $m\approx M$.

Second, some theoretical cross sections are computed under assumptions of
thermal equilibrium between matter and radiation and/or chemical equilibrium
between excited states. This assumption is not strictly justified either for
the Standard Model or for atomic dark matter, but is adopted regardless for
simplicity. For example, Cyr-Racine & Sigurdson (2013) has pointed out that
weakly coupled dark matter can thermally decouple from the radiation prior to
recombination, clearly breaking the assumption that each reaction is a
function of only one temperature. That work has also derived a recombination
cross-section with stimulated recombination, which explicitly depends on the
radiation temperature, and showed that the stimulated recombination yields
less than 60% change in the final free electron fraction. This level of
difference is small compared to the difference due to the overall uncertainty
in the reaction rates. The suggested rates in the literature for a single
reaction can differ by an order of magnitude or more (Savin et al., 2004).

Finally, the accuracy of the re-scaling method is limited by the accuracy of
the Standard Model interaction rates. In this work we do not attempt to
improve upon the uncertainties in the Standard Model literature. Instead, we
scrutinize the sensitivity of the final free electron fraction and molecular
hydrogen abundance on the interaction rates. In Sec. 2.4, we show that the
uncertainties in the important reaction rates propagate through the
calculation in a simple manner and demonstrate that the uncertainties at worst
linearly affect the final free electron fraction and molecular abundance (Fig.
3).

There are additionally three conditions which limit the range of parameters
within which our results can be trusted. The first is that if recombination
occurs in a sufficiently high density environment, generally quoted around
$10^{8}\,\rm cm^{-3}$ in the Standard Model, three body reactions will become
important (Palla et al., 1983). We can similarly find the density at which
three-body reactions become important by equating the two-body interaction
timescale

$t\sim(n\sigma v)^{-1},$ (2.13)

with the three-body interaction timescale

$t\sim(n^{2}\sigma^{5/2}v)^{-1}.$ (2.14)

Then, scaling the density at recombination as $n\sim
10^{3}\xi^{-3}r_{m}^{3}r_{\alpha}^{6}r_{M}^{-1}$ and $\sigma\propto a_{0}^{2}$
with the Bohr radius $a_{0}$, we find the following condition for three-body
reactions to be unimportant at recombination time:

$\xi^{-3}r_{M}^{-1}r_{\alpha}^{3}<10^{5}.$ (2.15)

We have included a minimal set of three-body reactions, but as discussed above
some caution is still warranted in the parameter regime where Eq. (2.15) is
violated.

A second limiting condition arises from the spectral distortions in the dark
photon field induced by the recombination process. In the Standard Model,
Hirata & Padmanabhan (2006) showed that the excess Lyman-$\alpha$ and two-
photon process contributions to the radiation field suppress the formation of
${\rm H}_{2}$ by a factor of $\sim 5$ compared to naive calculations. This
suppression is primarily due to ionization of ${\rm H}^{-}$ by the redshifted
spectral distortions, which dominate the radiation field at those frequencies.
The size of this effect is thus controlled by the number of spectral
distortion photons compared to the number of ${\rm H}^{-}$ ions. In the
Standard Model, roughly 5 photons are produced per hydrogen atom during the
recombination process (Chluba & Sunyaev, 2006). As long as this number does
not depend sensitively on the dark parameters, we do not expect any
enhancement of this effect and neglect it for simplicity. Since the energy and
angular momentum structure of the dark hydrogen are the same as ordinary
hydrogen, we do not expect this number to change as long as the radiation
field remains approximately black-body.

However, this is not guaranteed for all values of the parameters. The Standard
Model baryon-to-photon ratio $\eta$ is order $10^{-10}$, which means that the
photons produced by the recombination process do not substantially alter the
radiation field, other than through the Lyman-$\alpha$ and two photon
processes which we explicitly include. Manipulating $\xi$ and $M$ can produce
a much larger value of $\eta_{\rm D}$. To handle the recombination and
molecule formation problems in this case would require a full non-equilibrium
treatment of the spectral distortions, since the recombination photons will
destroy the black-body nature of the radiation field. To determine when this
may occur, we start from the fact that for the Standard Model the spectral
distortions due to higher hydrogen transitions are

$\frac{\Delta I_{\nu}}{B_{\nu}}\lesssim 10^{-7},$ (2.16)

with $\Delta I_{\nu}$ the intensity of the spectral distortion and $B_{\nu}$
the black-body intensity. We insist that the distortions remain less than
order unity and insert the scaling of $\eta_{D}$ to obtain

$\xi^{-3}r_{M}^{-1}<\mathcal{O}(10^{7}).$ (2.17)

Finally, in the Standard Model, the recombination temperature is lower than
the binding energy of hydrogen because of the low cosmological number density.
Since the gas temperature is set by the CMB, this density is parameterized by
the baryon-to-photon ratio $\eta$. Even the weak (logarithmic) dependence of
the recombination redshift on $\eta$ eventually pushes the recombination
temperature outside the regime of validity of some of our interaction rates.
We omit the parameters for which the re-scaled temperature (see Eq. (2.11)) at
the recombination redshift is greater than $10^{4}$ K. This is the highest
temperature at which all of the reaction rate laws remain valid. We do not
consider parameters where the re-scaled temperature at recombination is much
less than in the Standard Model.

When $r_{\alpha}\gtrsim 10^{-2/3}\simeq 0.2$, satisfying Eq. (2.15), the three
body limit, automatically satisfies Eq. (2.17), the spectral distortion
condition. For smaller $r_{\alpha}$, the opposite is true. The third
condition, from the temperature range of the reaction rates, seems never to be
a problem when the first two conditions are met. It is worth emphasizing that
these are practical rather than theoretical limitations. In theory, there is
no obstacle to the dark recombination occurring at arbitrarily high $\eta_{\rm
D}$.

###  2.4 Analytic estimation

We now develop some analytic approximations for the final abundance of
molecular states to help us digest the numerical results in Sec. 2.5 and
quantify the propagation of the uncertainties in reaction rates. Our approach
is based on Lepp & Shull (1984) (see also Stiavelli (2009)).

We begin with the two pathways for the formation of molecular hydrogen known
from the Standard Model, the ${{\rm H}_{{\rm D},2}^{+}}$ channel:

$\displaystyle{\rm H}_{\rm D}+{\rm H}_{\rm D}^{+}\rightarrow{\rm H}_{{\rm
D},2}^{+}+\gamma$ (2.18) $\displaystyle{\rm H}_{{\rm D},2}^{+}+{\rm H}_{\rm
D}\rightarrow{\rm H}_{{\rm D},2}+{\rm H}^{+}\,,$ (2.19)

and the ${\rm H}^{-}$ channel:

$\displaystyle{\rm H}_{\rm D}+{\rm e}_{\rm D}\rightarrow{\rm H}^{-}_{\rm
D}+\gamma$ (2.20) $\displaystyle{\rm H}_{\rm D}^{-}+{\rm H}_{\rm
D}\rightarrow{\rm H}_{{\rm D},2}+{\rm e}_{\rm D}\,.$ (2.21)

These reactions occur in chemical equilibrium:

$\frac{dx_{i}}{dt}=0,$ (2.22)

where $x_{i}$ refers to the fraction of ${\rm H}_{{\rm D},2}^{+}$ or ${\rm
H}_{\rm D}^{-}$. The equilibrium condition reduces the differential equation
for the abundance $x_{i}$ to an algebraic equation. We further simplify this
equation by considering only the dominant production and destruction channels
for the species. Initially, the second reaction in each pathway is negligible
compared to the radiative destruction of ${{\rm H}_{{\rm D},2}^{+}}$ and
${{\rm H}^{-}_{{\rm D}}}$. For example, the chemical equilibrium equation for
${\rm H}_{\rm D}^{-}$ is

$\displaystyle\frac{dx_{{\rm H}^{-}_{\rm D}}}{dt}$ $\displaystyle=0$
$\displaystyle=c_{{\rm H}3}x_{{\rm e}}n_{{\rm H}_{\rm D}\rm I}-c_{{\rm
H}4}x_{{\rm H}_{\rm D}^{-}}\,.$ (2.23)

Here, we use the time variable $t$ to emphasize that this chemical equilibrium
holds on time scales much shorter than the background evolution time scale,
the Hubble time.

As the universe cools adiabatically, chemical equilibrium dictates that
$x_{i}$ must increase so that the photo-destruction rate continues to balance
the creation rate. At some critical redshift, however, ionizing photons can no
longer efficiently destroy the intermediate species such as ${\rm H}_{2,{\rm
D}}^{+}$ and ${\rm H}_{\rm D}^{-}$. Then, chemical equilibrium exists not
between creation of the intermediate and its photo-destruction but instead
between creation of the intermediate and the formation of ${{\rm H}_{\rm
D,2}}$. Since the reactions which produce ${{\rm H}_{\rm D,2}}$ are fast
compared to those which produce the intermediate, we assume that all of the
intermediate species present at this critical redshift are converted into
${{\rm H}_{\rm D,2}}$. No additional ${{\rm H}_{\rm D,2}}$ is produced via the
given channel at lower redshifts.

That is, the ${\rm H}_{2}$ produced from the channel $i$ is

$x_{{\rm H}_{{\rm D},2},i}=x_{i,\rm CE}(z_{\rm crit}),$ (2.24)

where $x_{i,CE}$ is the chemical equilibrium solution of Eq. (2.22).

The ${\rm H}_{2}$ production from each channel depends on the free electron
abundance. We approximate the residual electron abundance of recombination as
the free electron abundance available for the molecule formation. The freeze-
out time ($z_{\rm fo}$), when the recombination rate equals the Hubble
expansion rate, sets the residual electron abundance as

$c_{{\rm H}1}(z_{\rm fo})n_{{\rm H}_{\rm D}}x_{{\rm e}}(z_{\rm fo})=H(z_{\rm
fo}).$ (2.25)

Here, we use $x_{{\rm e}}(z)$ from Saha equilibrium, that is, setting the
right hand side of Eq. (2.1) to zero, and neglecting the difference between
$x_{e_{\rm D}}$ and $x_{p_{\rm D}}$:

$\displaystyle\frac{x_{\rm e}^{2}}{1-x_{\rm e}}=$
$\displaystyle\frac{\beta_{{\rm H}_{\rm D}}}{\alpha_{{\rm H}_{\rm
D}}}\frac{1}{n_{{\rm H}_{\rm D}}}e^{-h\nu_{{\rm Ly\alpha},D}/k_{B}T_{\rm
g,D}}$ $\displaystyle=$ $\displaystyle\frac{1}{n_{{\rm H}_{\rm
D}}}\left(\frac{2\pi mk_{B}T}{h^{2}}\right)^{3/2}e^{-E_{{\rm H}_{\rm
D}}/k_{B}T_{\rm g,D}}\,.$ (2.26)

Some molecular hydrogen can be also produced during recombination, but in the
Standard Model the amount is very small and is usually neglected. For some
dark parameters, however, this first burst of formation dominates the final
molecular abundance. In particular, small $\xi$ triggers such a burst because
dark recombination happens earlier, so that the collisional $H_{2}$ formation
processes are enhanced compared to the Standard-Model case. In this case, we
estimate the formation of ${{\rm H}_{\rm D,2}}$ under the assumption of
chemical equilibrium between the dominant creation and destruction reactions:

$\displaystyle\frac{dx_{{{\rm H}_{\rm D,2}}}}{dt}$ $\displaystyle=0$
$\displaystyle=c_{\rm H5}n_{{\rm H}_{\rm D}\rm I}x_{{\rm H}_{\rm
D}^{-}}+c_{\rm H10}n_{{\rm H}_{\rm D}\rm I}x_{{\rm H}_{{\rm D},2}^{+}}$
$\displaystyle\quad-c_{\rm H15}n_{{\rm H}_{\rm D}^{+}}x_{{{\rm H}_{\rm
D,2}}}-c_{{\rm H}*}n_{{\rm H}_{\rm D}\rm I}x_{{{\rm H}_{\rm D,2}}},$ (2.27)

where the abundances of ${\rm H}_{{\rm D},2}^{+}$ and ${\rm H}^{-}_{\rm D}$
are themselves evaluated using chemical equilibrium between their dominant
creation and photodestruction processes and the expression is evaluated at
$z_{\rm eff}=(z_{\rm rec}+z_{\rm dec})/2$. That is, midway between the
redshift $z_{\rm rec}$ at which the Saha ionization fraction [Eq. (2.26)] is
$x_{\rm e}=0.5$ and the redshift of last scattering $z_{\rm dec}$, the peak of
the visibility function $g(z)=(d\tau/dz)e^{-\tau(z)}$ with the optical depth
$\tau(z)=\int{\rm d}\ell n_{e}\sigma_{T}$. This approximation relies on the
fact that the ${{\rm H}_{\rm D,2}}$ production is still slow relative to the
production and other (besides ${{\rm H}_{\rm D,2}}$ formation) destruction
reactions for the intermediates. That is, the equilibrium point between
creation and destruction of ${{\rm H}_{{\rm D},2}^{+}}$ and ${{\rm
H}^{-}_{{\rm D}}}$ is not strongly affected by this burst of ${{\rm H}_{\rm
D,2}}$ formation.

For the efficiency of molecule formation during recombination, the particle
number density $n_{{\rm H}_{\rm D}}$ at recombination turns out to be the
single strongest controlling parameter. This is because, unlike atomic
hydrogen, molecular hydrogen is mainly created and destroyed collisionally. In
the Standard Model, by the time of recombination ($T\sim 0.25$ eV), it is
already energetically favorable to form molecules, and in fact we do see a
burst of molecule formation in all cases we study in Sec. 2.5 [Fig. 1], as
soon as neutral atoms are available. Nonetheless, this formation channel is
completely negligible for the parameters close to the Standard Model. In the
Standard Model, this channel is inhibited by two factors: the radiative
destruction of ${\rm H}_{{\rm D},2}^{+}$ and ${\rm H}_{\rm D}^{-}$ (they are
fragile enough to be destroyed even at low temperatures), and the simple
unlikeliness of the necessary encounters between species with only a small
fractional abundance. Both problems are remedied by increasing the particle
number density. For a fixed $\rho_{\rm DM}$, the dark-particle-to-photon ratio
scales as $\xi^{3}/r_{M}$. That is, for low $\xi$ recombination happens at
very high redshift, and since $n\propto(1+z)^{3}$ the number density can be
much higher compared to the Standard Model. The dependence on the other
parameters is similar: the higher $m$ and $\alpha$ the higher the atomic
binding energy and hence the redshift of recombination. As the density at
recombination increases, we find that $x_{{\rm H}_{{\rm D},2}}\to 0.5$; that
is, the dark-matter gas approaches a completely molecular state.

###  2.5 Numerical results

In Fig. 1, we show the cosmological evolution from the modified Recfast++ of
the atomic dark matter states (ionized, atomic, and molecular) as a function
of redshift, renormalized in the unit of $z_{\rm rec}$, the recombination
redshift. For each panel, we show the abundances of ${{\rm H}}_{\rm D}^{+}$
(black, solid line), ${\rm H}_{{\rm D},2}^{+}$ (cyan, dot-dashed line), ${\rm
H}_{{\rm D}}^{-}$ (red, dashed line), and ${\rm H}_{{\rm D},2}$ (gold, dotted
line). The redshift range we show here spans from the beginning of
recombination to the end of molecule formation. The four panels in Fig. 1 show
a range of parameters, which are representative of the range of possible
evolution histories.

The dots correspond to our analytic estimates of the electron freeze out [Eq.
(2.25)] and each of the three channels for ${\rm H}_{{\rm D},2}$ production
[Eqs. (2.24)–(2.27)]. Note that each dot indicates only the additional ${{\rm
H}_{\rm D,2}}$ generated from a particular channel, so if this amount is much
smaller than the preexisting ${{\rm H}_{\rm D,2}}$ abundance the dot will lie
far from the gold curve. In that case, we have verified only that the estimate
correctly predicts little enough ${{\rm H}_{\rm D,2}}$ production to
negligibly alter the final abundance—not that this tiny fractional increase in
$x_{{\rm H}_{\rm D,2}}$ agrees exactly with the numerical result. As shown in
Fig. 1, the analytic approximation is generally accurate to within an order of
magnitude for the final ${\rm H}_{{\rm D},2}$ fraction.

The two top panels with $\xi=0.02$ are qualitatively quite similar to Standard
Model (e.g. Galli & Palla (1998); Lepp & Shull (1984)) tracks. Note however
that the recombination occurs at $z_{\rm rec}\sim 10^{5}$ in the top left
panel, due to the reduced temperature. In the right panel, this is largely
offset ($z_{\rm rec}\sim 2000$) by the low mass $m$ reducing the hydrogen
binding energy. By lowering $m$ or $\alpha$ (or increasing $\xi)$,
recombination can be delayed to arbitrarily low redshift or pushed into the
future. These parameter choices are less observationally relevant, since,
unless the coupling $\alpha$ is extremely small, the ionized dark matter with
the dark-electromagnetic interaction would fail to behave like CDM on large
scales.

The bottom two panels show more extreme parameters. The bottom left shows
“heavy” dark matter with $\xi=0.005$. Here, the recombination happens at
$z\sim 10^{9}$ in a very high density environment. This precedes even Big-Bang
Nucleosynthesis (BBN), but for completely secluded DM, BBN is not affected
because the contribution to the total relativistic species from the dark
matter is extremely small with $\xi=0.005$. In this case we find the
asymptotic value of $x_{{{\rm H}_{\rm D,2}}}=0.13$. As we have discussed at
the end of Sec. 2.4, this reflects a general trend toward forming more
molecular hydrogen when the formation occurs at high density. Here, the very
large ${{\rm H}_{\rm D,2}}$ abundance allows the formation of ${\rm H}_{{\rm
D},3}^{+}$ to compete with the formation of ${{\rm H}_{\rm D,2}}$ once the
${{\rm H}_{{\rm D},2}^{+}}$ photo-dissociation becomes ineffective. The ${\rm
H}_{{\rm D},3}^{+}$ is in turn efficiently converted into neutral hydrogen and
${{\rm H}_{\rm D,2}}$ by dissociative recombination. This leads to a net
reduction in $x_{{\rm H}_{\rm D}^{+}}$.

Finally, the bottom right panel shows $r_{\alpha}=0.274$ and $\xi=0.005$, and
$r_{M}=10$. Here we see that even for small values of $\xi$ we can end up with
the same qualitative story as in the Standard Model provided that the other
parameters are adjusted to control the number density at recombination. Note
however that the weak coupling leads to a relatively high freeze-out
ionization fraction and relatively low ${{\rm H}_{\rm D,2}}$ fraction.

Figure 1: The evolution of the dark matter gas for several choices of
parameters. The circles colored to match each line show our analytic
approximations for the freeze out free proton fraction and ${\rm H}_{{\rm
D},2}$ produced from each of the three colors (gold for the recombination
channel, cyan for the ${\rm H}_{{\rm D},2}^{+}$ channel and red for the ${\rm
H}_{\rm D}^{-}$ channel. If a previous channel has already created a large
amount of molecular hydrogen these dots can be far from the ${\rm H}_{{\rm
D},2}$ curve and still accurately depict the amount of ${\rm H}_{{\rm D},2}$
formed through that channel.

In Fig. 2, we show explicitly the evolution of the free ion fraction in our
full treatment (black, solid), the 3-level treatment neglecting the formation
of molecules (green, dashed), and in Saha equilibrium (blue, dot-dashed) for
the same four choices of parameters. Except in the bottom left, the
recombination proceeds initially nearly in equilibrium, until the expansion of
the universe drives the freeze-out. In the bottom left, the recombination rate
is suppressed by the large electron mass ($\alpha_{{\rm H}_{\rm D}}\propto
r_{m}^{-2}r_{\alpha}^{2}$), so the recombination rate drops out of equilibrium
(i.e. falls below the Hubble rate) early in the recombination. This leads to
early departure from the Saha ionization fraction and a large freeze-out
ionization fraction.

We also plot $T_{\rm g,D}/T_{\gamma,D}$, illustrating the thermal decoupling
of the gas from the dark photons. In the Standard Model, the thermal
decoupling is delayed compared to recombination by the low baryon-to-photon
ratio, $\eta$. The decoupling in the bottom left is qualitatively most similar
to the Standard Model, but the delay is due instead to the large freeze-out
free electron fraction and high temperature around recombination. In the
bottom right, the weak coupling means that the decoupling precedes
recombination, while in the top panels the decoupling occurs during
recombination.

Note that except in the bottom left the molecules have a minimal effect on the
ionization fraction at all times. The ${\rm H}_{{\rm D},2}$ production
pathways do not lead to a net depletion of electrons and (for the parameters
considered) they proceed rapidly compared to all other reactions. However,
both the formation and destruction of ${\rm H}_{{\rm D},3}^{+}$ lead to a net
recombination. A difference is also discernible in the bottom right, due to
the mutual neutralization reaction $c_{\rm H7}$ converting some of the ${\rm
H}_{\rm D}^{-}$ and ${\rm p}_{\rm D}$ into neutral hydrogen rather than
molecular hydrogen. That reaction rate is enhanced by
$r_{\alpha}^{3}r_{m}^{-3}\approx 10^{5}$, and still leads to only a factor
$\lesssim 2$ difference in the final ionization fraction. We conclude that
except when ${{\rm H}_{\rm D,2}}$ dominates over $x_{{\rm H}_{\rm D}^{+}}$,
the cosmological scales which depend on opacity to Thomson scattering
(discussed below) are weakly affected by the inclusion of molecules

Figure 2: The evolution of the ionization fraction for several choices of
parameters: in our full treatment, neglecting molecular states, and in Saha
equilibrium. Note that the inclusion of molecules has a quite small effect on
the free proton fraction, except in the bottom left.

We conclude this section by analyzing the propagation of uncertainties in the
rates through our calculation, using our analytic approximation. In Fig. 3,
for each reaction in turn we have artificially first increased and then
decreased the rate by a factor of 10. We use the reaction variations which
produce the highest and lowest final $x_{{\rm H}_{\rm D}^{+}}$, $x_{e_{\rm
D}}$, and $x_{{\rm H}_{{\rm D},2}}$ to plot an uncertainty envelope around
each track. In the same way, by feeding these artificial rates to our analytic
model, we can understand how sensitively our results depend on the accuracy of
the reaction rates.

Except in the bottom left panel, the analytic error bars closely match the
numerical error bars. In the bottom left, the estimate for the ${{\rm H}_{\rm
D,2}}$ abundance fails because the high density, $n_{{\rm H}_{\rm D}}$, allows
significant molecule formation for a large window around recombination, even
in chemical equilibrium. This is not captured by our estimate assuming a short
burst of formation with the electron fraction inferred from Saha equillibrium.
The estimate for $x_{{\rm H}_{\rm D}^{+}}$ fails because it does not take into
account the ${\rm H}_{{\rm D},3}^{+}$ recombination channel: the reaction
which produces the largest effect numerically (${{\rm H}_{{\rm D},2}^{+}}$
formation, $c_{{\rm H}8}$) has no effect in the analytic treatment. We also
see a uniquely sensitive dependence of $x_{{\rm H}_{\rm D}^{+}}$ on this
reaction rate. The linear dependence of the freeze-out abundance on the
reaction rate revealed in Eq. (2.23) holds only when the abundances of the
reactants can be treated as constant. Here we see instead the exponential rate
dependence of the solution $x_{{\rm H}_{\rm D}^{+}}\propto\exp[{-t}/{(c_{{\rm
H}i}n_{{\rm H}_{{\rm D}}\rm I}})]$. The electron abundance is not so radically
affected, because only the lesser of $x_{{\rm H}_{\rm D}^{+}}$ and $x_{{\rm
H}_{{\rm D},3}^{+}}$ is subject to exponential depletion, and once $x_{{\rm
H}_{{\rm D},3}^{+}}$ exceeds $x_{{\rm p}_{\rm D}}$, we have $x_{\rm e_{\rm
D}}\approx x_{{\rm H}_{{\rm D},3}^{+}}$. Fig. 3 tells us that we can be
confident in the accuracy of our results to almost exactly the extent we are
confident in the rates of these four key reactions. There are no non-linear
dependencies where small errors in certain rates leading to large changes in
the final results. This mitigates the impact of uncertainties in the Standard-
Model reaction coefficients.

Figure 3: Effects of increasing and decreasing individually the reactions with
the largest effect on the final ${{\rm H}_{\rm D,2}}$ fraction (purple,
shaded), $x_{{\rm p}_{\rm D}}$ (gray, shaded), and $x_{\rm e_{{\rm D}}}$
(gray, shaded, hatched), along with the effect inferred from varying those
rates by the same amount in our analytic approximation (dots with bars). The
results do not sensitively depend on most of the rates involved. See text for
discussion of the large difference between $x_{{\rm p}_{\rm D}}$ and $x_{\rm
e_{{\rm D}}}$ in the bottom left.

Now we are ready to show the dependence of the molecular-hydrogen fraction at
redshift 30 on the model parameters in Fig. 4, by first fixing $\xi$ and
$\alpha$ (left panel) while varying the masses, and then fixing the masses
while varying $\xi$ and $\alpha$ (right panel). The effect of varying each
parameter can be understood through its effect on the number density $n_{{\rm
H}_{\rm D}}$ at fixed photon temperature: a higher number density leads to
more molecular hydrogen.

Figure 4: The dependence of the molecular hydrogen fraction $x_{{\rm
H}_{2,{\rm D}}}$ at redshift 30 on the model parameters.

##  3 Cosmological Consequences

A full exploration of the cosmological consequences of dark matter which can
cool is an ambitious program. Generating predictions for the merger rate of
dark compact objects will require further analytic work and simulations at a
wide range of scales. In this section, however, armed with a detailed picture
of the homogeneous evolution of atomic dark matter in the previous sections,
we use simple analytic tools to sketch out the large-scale structure in this
model. We shall demonstrate that the model can easily satisfy existing
observational constraints and is consistent with related work while estimating
the results of future, more detailed investigations. We begin by examining the
alterations to the linear matter power spectrum in our model and then proceed
to estimate the halo mass function in the Press-Schecter formalism (Press &
Schechter, 1974).

###  3.1 Linear Power Spectrum: Acoustic Oscillations and Diffusion Damping

On scales much larger than the horizon size at the dark matter/dark photon
decoupling, dissipative dark matter behaves like ordinary cold dark matter. On
small scales this is not the case. An exact treatment of the alterations to
the linear matter power spectrum in a dissipative dark matter model requires
solving the Einstein-Boltzmann equations for the given dark matter model, as
in Cyr-Racine & Sigurdson (2013). Alternatively, one can adopt a
phenomenological parameterization of the matter power spectrum and then map
the parameters of this model back to the dark matter microphysics, as in Bohr
et al. (2020). Our approach hews closer to that of Bohr et al. (2020), but
simpler still.

Atomic dark matter can alter the linear matter power spectrum on small scales
primarily through two effects: dark acoustic oscillations (DAO) and dark
diffusion (Silk) damping. The following two length scales define the pivot
length scales below which the two effects alter the matter power spectrum.

The DAO scale is given as

$d_{\rm DAO}=\int_{0}^{a_{\rm dec}}c_{s}(a)\frac{da}{a^{2}H(a)},$ (3.1)

where we calculate the physical sound speed from

$c_{s}(a)=\frac{c}{\sqrt{3(1+R(a))}}.$ (3.2)

as a function of the scale factor $a=1/(1+z)$. The quantity $R(a)$ is the
ratio of the relativistic inertias, $\bar{\rho}+\bar{P}$, with the background
density $\bar{\rho}$ and pressure $\bar{P}$, of the dark matter and dark
photon, given as

$R(a)=\frac{3\bar{\rho}_{DM}}{4\bar{\rho}_{\gamma,D}}=2.32\times
10^{10}a\left(\frac{\Omega_{DM}}{0.25}\right)\left(\frac{\xi}{0.02}\right)^{-4}.\,$
(3.3)

Clearly, for a low dark photon temperature ($\xi$), $R(a)$ can be very large.
The DAO scale depends only on $R(a)$ (which depends most strongly on $\xi$)
and on the scale factor at decoupling, which is controlled dominantly by
$E_{{\rm H}_{\rm D}}$ compared to $\xi$.

The diffusion damping scale is given by (Zaldarriaga & Harari, 1995)

$\frac{1}{k_{D}^{2}}=\int_{0}^{a_{\rm
dec}}\frac{da}{a^{2}H(a)}\frac{1}{6(1+R)n_{{\rm e}_{\rm D}}\sigma_{\rm
T,D}a}\left[\frac{16}{15}+\frac{R^{2}}{1+R}\right]\,.$ (3.4)

Note that even for very small or large values of $R(a)$ the diffusion scale is
the geometric mean of the photon mean free path and the horizon size, modified
by factors of order unity. A weak dependence on $M$ enters through the number
density.

Since the DAO scale is suppressed by $R(a)$ while the diffusion scale is not,
for low $\xi$ the diffusion scale can far exceed the DAO scale. In this case,
the DAOs are completely washed out by diffusion damping. We show the
dependence of these two scales on the model parameters in Fig. 5.

Figure 5: The dependence of the diffusion damping scale (black, dashed) and
the acoustic scale (color bars) on the model parameters.

Since we find that for the parameters studied in Ryan et al. (2022b) the
diffusion scale is larger than the acoustic scale, we calculate the matter
power spectrum assuming CDM using CAMB (Lewis et al., 2000) and then apply a
Gaussian cutoff (as in e.g. Seljak (1994)):

$P(k,z)=P(k,z)_{\rm CAMB}e^{-(k/k_{D})^{2}}.$ (3.5)

On scales smaller than the diffusion scale, any other modifications to the
power spectrum are buried under the exponential diffusion damping. This
treatment is sufficient to determine whether the formation of the halos
studied in Ryan et al. (2022b) is suppressed by the diffusion scale, and
whether this suppression will appear at observable scales in the matter power
spectrum.

###  3.2 Implications for Dissipative Dark Halos

Equipped with the modified matter power spectrum, we can estimate the halo
mass function at various redshifts. In addition to the aforementioned
modification to the linear power spectrum, atomic dark matter does not
virialize exactly like CDM on small scales. In particular, the dynamics of the
collapse will be altered by scattering and subsequent radiative transitions if
some of the dark matter is ionized.

However, early in the nonlinear collapse (which takes approximately a Hubble
time), the number density does not differ much from background. Therefore, the
ionization fraction will not differ much from the freeze-out value. Since the
gas is predominantly neutral, it will behave as CDM through this phase. Just
prior to virialization, the number density and temperature rise rapidly, but
at this point dissipative effects cannot prevent the formation of the halo or
alter its mass.

However, the eventual shape and structure of these halos may be quite
different if they can cool efficiently subsequent to their formation.
Hydrodynamical simulations with cooling will be necessary to resolve these
effects.

Motivated by this discussion, we adopt the approach which has been developed
for warm dark matter (WDM), and has recently been applied to gravity-only
simulations of various self interacting dark matter models (Bohr et al.,
2021). Atomic dark matter is similar to WDM in that both models generate a
suppression of power at small scales while behaving like CDM on large scales.

In principle, thermal velocities in WDM models can not only lead to a
suppression of power on small scales but can also interfere with the
virialization of halos. This has been treated in Benson et al. (2012) by
employing a “moving barrier”, raising the critical density for halos smaller
than the suppression scale. Later work (Schneider et al., 2013; Angulo et al.,
2013), however, shows that for simulations initialized in the nonlinear regime
the thermal velocity effect can be safely neglected compared to the
gravitational velocities. That is, the only modification necessary is to the
linear matter power spectrum, exactly as in atomic dark matter.

We begin with the analytic prediction for the halo mass function in Press-
Schechter theory (Press & Schechter, 1974). There, the cumulative mass
function of halos of mass greater than $M_{\rm h}$ is given by the fraction of
the smoothed linear density field (with smoothing radius $R=[3M_{\rm
h}/(4\pi\bar{\rho}_{M})]^{1/3}$) exceeding the critical value
$\delta_{c}\simeq 1.69$ calculated in the spherical collapse model. Here,
$\bar{\rho}_{M}$ is the comoving background matter density. Then, the
differential number of halos of mass $M_{h}$ per volume is given as

$\frac{dn}{dM_{\rm h}}=\frac{\bar{\rho}}{M_{\rm
h}^{2}}\left\rvert\frac{d\operatorname{ln}\sigma_{M_{\rm
h}}}{d\operatorname{ln}M_{\rm h}}\right\rvert f(\sigma_{M_{\rm h}}),$ (3.6)

with $\sigma_{M_{\rm h}}$ is the root-mean-square of the overdensity smoothed
over a region with length scale $R$:

$\sigma^{2}_{M_{h}}=\frac{1}{2\pi^{2}}\int dkk^{2}P(k)W_{R}^{2}(k).$ (3.7)

The window function $W_{R}(k)$ is the Fourier transform of the smoothing
filter $W_{R}(r)$, and the Press-Schechter multiplicity function is

$f(\sigma_{M_{\rm h}})=\sqrt{\frac{2}{\pi}}\frac{\delta_{c}}{\sigma_{M_{\rm
h}}}{\rm exp}\left(-\frac{\delta_{c}^{2}}{2\sigma_{M_{\rm h}}^{2}}\right)\,.$
(3.8)

Strictly speaking (Bond et al., 1991), the formalism is only valid for the
sharp-$k$ filter:

$W_{R}(k)=\begin{cases}1&k\leq 1/R\\\ 0&k>1/R\end{cases},$ (3.9)

which has the property that the contributions to $\sigma_{M_{h}}$ from each
successive decrease in $M_{h}$ are uncorrelated. However, for CDM the results
are quite insensitive to the choice of filter, so a real space tophat filter
is often employed. Finally, Sheth & Tormen (1999) derived a multiplicity
function by fitting to simulation results,

$f_{ST}(\sigma)=A\sqrt{\frac{2a}{\pi}}\left[1+\left(\frac{\sigma^{2}}{a\delta_{c}^{2}}\right)^{p}\right]\frac{\delta_{c}}{\sigma}\exp{\left(-\frac{a\delta_{c}^{2}}{2\sigma^{2}}\right)},$
(3.10)

with $p=0.3$, $A=0.322$, $a=0.707$ (Stiavelli, 2009), which has been widely
used in the literature.

The problem with applying the Press-Schechter formalism to atomic dark matter
is that, for models with a cut-off in the power spectrum, the halo mass
function becomes sensitive to the choice of the smoothing filter. Schneider et
al. (2013) found that the sharp-$k$ filter (which is anyway the only
theoretically justified choice) produces the best agreement with simulations
when modifying the mass-to-smoothing-radius relation as

$M_{sk}=\frac{4\pi}{3}\bar{\rho}[cR_{sk}]^{3},$ (3.11)

where $c=2.7$ from calibration to their simulations, and setting $a=1$ for the
$a$ parameter in the Sheth-Tormen multiplicity function. Note that warm dark
matter simulations including Schneider et al. (2013) have been plagued by
force resolution issues which have led to spurious small-scale structures.
This was largely resolved in Angulo et al. (2013), which justifies using the
Press-Schechter formalism with modifications described by Schneider et al.
(2013) to obtain a qualitatively correct halo mass function. However, that
work cautions that further effort is needed to obtain a halo mass function
which is quantitatively accurate for scales near and below the cutoff in the
matter power spectrum. These caveats apply equally to atomic dark matter.

In Fig. 6, we show the halo mass function $z=0$ (top panel) and $z=5$ (bottom
panel) for the usual CDM case and for the atomic dark matter case with
$k_{d}=10{h/{\rm Mpc}}$ (blue, dot-dashed line) and $k_{d}=1{h/{\rm Mpc}}$
(green, dashed line). The diffusion scale is the only relevant variable
affecting the halo mass function.

Figure 6: The halo mass function at $z=0$ (upper) and $z=5$ (lower). The
Gaussian cutoff in the power spectrum leads to a corresponding suppression in
the halo mass function at small scales. Clearly, $k_{d}=1$ is ruled out by the
existence of $\sim 10^{11}\,\rm M_{\odot}$ galaxies. We include it here to
illustrate the dependence of the halo mass function on $k_{d}$.

##  4 Conclusion

In this paper, we have calculated the background evolution of the abundance of
dark molecules and ions in the atomic-dark-matter model. To do so, we use the
Recfast++ software and the interaction-rate coefficients re-scaled by the
methods derived in Ryan et al. (2022a). Previous work in the mirror matter
model (equivalent to atomic dark matter with $r_{\alpha}=r_{m}=r_{M}=1$ for
our purposes) (Latif et al., 2019; D’Amico et al., 2018) assumed the ratio
$x_{e}/x_{H_{2}}\approx 10^{2}$ carried over from the Standard Model,
independent of the dark matter parameters. Here, we calculate the primordial
molecular abundance directly and find that this assumption often fails,
because the primordial molecular abundance depends not only on the primordial
ionization fraction but also on the dark-particle number density at the time
of molecule formation.

We have also presented an analytic calculation for the final dark-molecular
abundance estimate, which agrees well with the full numerical result. The
analytic estimate also indicates that the final abundance is at worst linearly
sensitive to the uncertainties in the reaction rates.

We find that when dark recombination happens at high redshift, large dark
particle number density promotes the molecular processes so that most of the
atomic dark matter particles ends up in molecules. However, as long as the
molecule formation is fast compared to the recombination (which is true for
all cases studied here) the molecular processes hardly change the
recombination history, thereby, the DAO scale and the dark diffusion (Silk)
damping scale.

We have also estimated the linear power spectrum of atomic dark matter by
calculating the DAO scale and dark-diffusion (Silk) damping scale. When the
dark photon temperature is much lower than the CMB temperature, the dark
diffusion scale is much larger than the DAO scale, so that we neglect the DAO
effect and modify the linear power spectrum by a Gaussian damping factor due
to diffusion.

Finally, using the modified linear matter power spectrum, we have estimated
the halo mass function using the Press-Schechter formalism modified to fit the
warm dark matter simulation results. This work provides the initial conditions
and some useful estimates for the results of full atomic dark matter
simulations with chemistry, beginning with the companion paper, Ryan et al.
(2022b).

In the future, hydrodynamical simulations with radiative cooling will provide
insight into the implications of atomic dark matter at scales ranging from the
largest structures in the universe to the smallest compact objects. Although
precise results for such complex nonlinear physics will require significant
investment, this program will begin to promote a broad class of alternative
dark matter models—all those in which the dark matter forms hydrogen-like
bound states—to testable theories whose parameters can be constrained by
observation.

The authors thank Jens Chluba for providing a special version of the Recfast++
code and for useful discussions about the calculation of cosmic recombination
history and spectral distortions. We thank the anonymous referee for providing
a truly excellent and professional report, which contributed significantly to
the depth of analysis presented in the revised version. Funding for this work
was provided by the Charles E. Kaufman Foundation of the Pittsburgh
Foundation.

##  Appendix A Atomic dark matter implementation of Recfast++

The Recfast++ code calculates the temperature and ionization history of
Standard Model hydrogen and helium. In order to implement the evolution of
atomic dark matter states using the Recfast++, besides changing atomic rate
coefficients as described in Sec. 2.1, we have re-defined the input variables
as follows. With these redefinitions, the same Recfast++ code can be applied
to both the standard recombination and the atomic dark matter scenario.

  * •

We have preserved the basic Recfast++ functionality, including helium
recombination, but set the helium fraction to zero.

  * •

The Recfast++ variable representing the CMB temperature $T_{0}$ is
reinterpreted as the dark photon temperature today, parameterized by
$T_{\gamma,D}=\xi T_{0}$.

  * •

The number density of the dark matter is calculated as

$n_{H_{D}}(z)=\frac{3H_{0}^{2}\Omega_{DM}}{8\pi GM}(1+z)^{3}\,,$ (A1)

where $M$ is the mass of dark matter proton. In Recfast++, the hydrogen number
density is parameterized by $\Omega_{b}$. To match Eq. (A1), we set

$\Omega_{b}\rightarrow\Omega_{DM}r_{M}^{-1}\,,$ (A2)

where $r_{M}=M/m_{p}$. This is equivalent to setting the input dark-matter
density parameter

$\Omega_{DM}\to\Omega_{b}+(1-r_{M}^{-1})\Omega_{DM}\,,$ (A3)

because we need to keep the total matter density
$\Omega_{m}=\Omega_{DM}+\Omega_{b}$ equal to the usual $\Lambda$CDM cosmology
value. Since Recfast++ takes $\Omega_{m}$ as an input, the modification to
$\Omega_{DM}$ is only implicit.

  * •

Finally, we have included all the relativistic degrees of freedom other than
the dark photons into $N_{\rm eff}$:

$N_{\rm eff}\rightarrow\frac{1}{\xi^{4}}\left[N_{\rm
eff}+\frac{8}{7}\left(\frac{T_{\rm CMB}}{T_{\nu}}\right)^{4}\right]\,,$ (A4)

so that the total radiation energy density becomes

$\displaystyle\rho_{\rm rad}=\rho_{\gamma}+\rho_{\nu}+\rho_{{\rm D},\gamma}=$
$\displaystyle\frac{\pi^{2}}{15}\left[1+\frac{7}{8}N_{\rm
eff}\left(\frac{T_{\nu}}{T_{\rm CMB}}\right)^{4}+\left(\frac{T_{\gamma,{\rm
D}}}{T_{\rm CMB}}\right)^{4}\right]T_{\rm CMB}^{4}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{15}\left[1+\frac{7}{8}\frac{1}{\xi^{4}}\left(N_{\rm
eff}+\frac{8}{7}\left(\frac{T_{\rm
CMB}}{T_{\nu}}\right)^{4}\right)\left(\frac{T_{\nu}}{T_{\rm
CMB}}\right)^{4}\right]T_{\gamma,{\rm D}}^{4}\,.$ (A5)

## References

  * Abbott et al. (2020) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, The Astrophysical Journal, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abel (2001) Abel, T. 2001, Science, 295, 93–98, doi: 10.1126/science.1063991
  * Ackerman et al. (2009) Ackerman, L., Buckley, M. R., Carroll, S. M., & Kamionkowski, M. 2009, Phys. Rev. D, 79, doi: 10.1103/physrevd.79.023519
  * Ade et al. (2016) Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, Astronomy & Astrophysics, 594, A13, doi: 10.1051/0004-6361/201525830
  * Aghanim et al. (2020) Aghanim, N., Akrami, Y., Arroja, F., et al. 2020, Astronomy & Astrophysics, 641, A1, doi: 10.1051/0004-6361/201833880
  * Agnese et al. (2018) Agnese, R., Aramaki, T., Arnquist, I., et al. 2018, Physical Review Letters, 120, doi: 10.1103/physrevlett.120.061802
  * Agrawal et al. (2017a) Agrawal, P., Cyr-Racine, F.-Y., Randall, L., & Scholtz, J. 2017a, JCAP , 2017, 021, doi: 10.1088/1475-7516/2017/08/021
  * Agrawal et al. (2017b) —. 2017b, Journal of Cosmology and Astroparticle Physics, 2017, 022, doi: 10.1088/1475-7516/2017/05/022
  * Akerib et al. (2017) Akerib, D., Alsum, S., Araújo, H., et al. 2017, Physical Review Letters, 118, doi: 10.1103/physrevlett.118.021303
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Physical Review D, 103, doi: 10.1103/physrevd.103.083533
  * Albert et al. (2017) Albert, A., Anderson, B., Bechtol, K., et al. 2017, The Astrophysical Journal, 834, 110, doi: 10.3847/1538-4357/834/2/110
  * Angulo et al. (2013) Angulo, R. E., Hahn, O., & Abel, T. 2013, Monthly Notices of the Royal Astronomical Society, 434, 3337–3347, doi: 10.1093/mnras/stt1246
  * Armstrong (1963) Armstrong, B. H. 1963, Physical Review, 131, 1132, doi: 10.1103/physrev.131.1132
  * Bates & Lewis (1955) Bates, D. R., & Lewis, J. T. 1955, Proceedings of the Physical Society. Section A, 68, 173, doi: 10.1088/0370-1298/68/3/306
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, The Astrophysical Journal Supplement Series, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Benson et al. (2012) Benson, A. J., Farahi, A., Cole, S., et al. 2012, Monthly Notices of the Royal Astronomical Society, 428, 1774–1789, doi: 10.1093/mnras/sts159
  * Bergström et al. (2013) Bergström, L., Bringmann, T., Cholis, I., Hooper, D., & Weniger, C. 2013, Physical Review Letters, 111, doi: 10.1103/physrevlett.111.171101
  * Bezanson et al. (2012) Bezanson, J., Karpinski, S., Shah, V. B., & Edelman, A. 2012, Julia: A Fast Dynamic Language for Technical Computing.  https://arxiv.org/abs/1209.5145
  * Boddy et al. (2016) Boddy, K. K., Kaplinghat, M., Kwa, A., & Peter, A. H. G. 2016, Physical Review D, 94, 123017, doi: 10.1103/PhysRevD.94.123017
  * Bode et al. (2001) Bode, P., Ostriker, J. P., & Turok, N. 2001, The Astrophysical Journal, 556, 93–107, doi: 10.1086/321541
  * Bohr et al. (2021) Bohr, S., Zavala, J., Cyr-Racine, F.-Y., & Vogelsberger, M. 2021, The halo mass function and inner structure of ETHOS haloes at high redshift.  https://arxiv.org/abs/2101.08790
  * Bohr et al. (2020) Bohr, S., Zavala, J., Cyr-Racine, F.-Y., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 3403, doi: 10.1093/mnras/staa2579
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, Astrophys. J., 379, 440, doi: 10.1086/170520
  * Boveia & Doglioni (2018) Boveia, A., & Doglioni, C. 2018, Annual Review of Nuclear and Particle Science, 68, 429–459, doi: 10.1146/annurev-nucl-101917-021008
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, Monthly Notices of the Royal Astronomical Society: Letters, 415, L40–L44, doi: 10.1111/j.1745-3933.2011.01074.x
  * Boylan-Kolchin et al. (2012) —. 2012, Monthly Notices of the Royal Astronomical Society, 422, 1203–1218, doi: 10.1111/j.1365-2966.2012.20695.x
  * Bromm & Larson (2004) Bromm, V., & Larson, R. B. 2004, Annual Review of Astronomy and Astrophysics, 42, 79–118, doi: 10.1146/annurev.astro.42.053102.134034
  * Brooks (2019) Brooks, A. M. 2019, Illuminating Dark Matter, 19–28, doi: 10.1007/978-3-030-31593-1_3
  * Browne & Dalgarno (1969) Browne, J. C., & Dalgarno, A. 1969, Journal of Physics B: Atomic and Molecular Physics, 2, 885, doi: 10.1088/0022-3700/2/8/309
  * Buckley & Peter (2018) Buckley, M. R., & Peter, A. H. 2018, Physics Reports, 761, 1–60, doi: 10.1016/j.physrep.2018.07.003
  * Bullock (2010) Bullock, J. S. 2010, Notes on the Missing Satellites Problem.  https://arxiv.org/abs/1009.4505
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, Annual Review of Astronomy and Astrophysics, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Chen & Kamionkowski (2004) Chen, X., & Kamionkowski, M. 2004, Phys. Rev. D, 70, 043502, doi: 10.1103/PhysRevD.70.043502
  * Chluba & Sunyaev (2006) Chluba, J., & Sunyaev, R. A. 2006, Astronomy & Astrophysics, 458, L29–L32, doi: 10.1051/0004-6361:20066191
  * Chluba et al. (2010) Chluba, J., Vasil, G. M., & Dursi, L. J. 2010, Monthly Notices of the Royal Astronomical Society, 407, 599, doi: 10.1111/j.1365-2966.2010.16940.x
  * Cline et al. (2014) Cline, J. M., Liu, Z., Moore, G. D., & Xue, W. 2014, Phys. Rev. D, 90, 015023
  * Collaboration et al. (2013) Collaboration, C., Agnese, R., Ahmed, Z., et al. 2013, Silicon Detector Dark Matter Results from the Final Exposure of CDMS II.  https://arxiv.org/abs/1304.4279
  * Cyr-Racine et al. (2014) Cyr-Racine, F.-Y., de Putter, R., Raccanelli, A., & Sigurdson, K. 2014, Phys. Rev. D, 89, 063517, doi: 10.1103/PhysRevD.89.063517
  * Cyr-Racine & Sigurdson (2013) Cyr-Racine, F.-Y., & Sigurdson, K. 2013, Physical Review D, 87, 103515, doi: 10.1103/PhysRevD.87.103515
  * D’Amico et al. (2018) D’Amico, G., Panci, P., Lupi, A., Bovino, S., & Silk, J. 2018, Monthly Notices of the Royal Astronomical Society, 473, 328, doi: 10.1093/mnras/stx2419
  * de Blok (2010) de Blok, W. J. G. 2010, Advances in Astronomy, 2010, 1–14, doi: 10.1155/2010/789293
  * de Jong (1972) de Jong, T. 1972, Astron. Astrophys., 20, 263.  https://ui.adsabs.harvard.edu/abs/1972A&A....20..263D
  * Diemand et al. (2008) Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Nature, 454, 735–738, doi: 10.1038/nature07153
  * Draine (2011) Draine, B. T. 2011, Physics of the interstellar and intergalactic medium (Princeton University Press)
  * Fan et al. (2013a) Fan, J., Katz, A., Randall, L., & Reece, M. 2013a, PDU, 2, 139, doi: 10.1016/j.dark.2013.07.001
  * Fan et al. (2013b) —. 2013b, Physics of the Dark Universe, 2, 139, doi: 10.1016/j.dark.2013.07.001
  * Feng et al. (2009) Feng, J. L., Kaplinghat, M., Tu, H., & Yu, H.-B. 2009, JCAP , 2009, 004, doi: 10.1088/1475-7516/2009/07/004
  * Foot & Vagnozzi (2015) Foot, R., & Vagnozzi, S. 2015, Phys. Rev. D, 91, 023512, doi: 10.1103/physrevd.91.023512
  * Foot & Vagnozzi (2016) —. 2016, JCAP , 2016, 013, doi: 10.1088/1475-7516/2016/07/013
  * Galli & Palla (1998) Galli, D., & Palla, F. 1998, arXiv:astro-ph/9803315.  http://arxiv.org/abs/astro-ph/9803315
  * Ghalsasi & McQuinn (2018) Ghalsasi, A., & McQuinn, M. 2018, Phys. Rev. D, 97, 123018, doi: 10.1103/PhysRevD.97.123018
  * Glover & Abel (2008) Glover, S. C. O., & Abel, T. 2008, Mon. Not. R. Astron. Soc., 388, 1627, doi: 10.1111/j.1365-2966.2008.13224.x
  * Glover & Savin (2009) Glover, S. C. O., & Savin, D. W. 2009, Monthly Notices of the Royal Astronomical Society, 393, 911, doi: 10.1111/j.1365-2966.2008.14156.x
  * Goldberg & Hall (1986) Goldberg, H., & Hall, L. J. 1986, Physics Letters B, Volume 174, Issue 2, p. 151-155., 174, 151, doi: 10.1016/0370-2693(86)90731-8
  * Grassi et al. (2014) Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, Mon. Not. R. Astron. Soc., 439, 2386, doi: 10.1093/mnras/stu114
  * Hart & Chluba (2017) Hart, L., & Chluba, J. 2017, Monthly Notices of the Royal Astronomical Society, 474, 1850–1861, doi: 10.1093/mnras/stx2783
  * Hartwig et al. (2016) Hartwig, T., Volonteri, M., Bromm, V., et al. 2016, Monthly Notices of the Royal Astronomical Society: Letters, 460, L74–L78, doi: 10.1093/mnrasl/slw074
  * Hirata & Padmanabhan (2006) Hirata, C. M., & Padmanabhan, N. 2006, Monthly Notices of the Royal Astronomical Society, 372, 1175–1186, doi: 10.1111/j.1365-2966.2006.10924.x
  * Hogan & Dalcanton (2000) Hogan, C. J., & Dalcanton, J. J. 2000, Phys. Rev. D, 62, 063511, doi: 10.1103/PhysRevD.62.063511
  * Kaplan et al. (2010) Kaplan, D. E., Krnjaic, G. Z., Rehermann, K. R., & Wells, C. M. 2010, Journal of Cosmology and Astroparticle Physics, 2010, 021, doi: 10.1088/1475-7516/2010/05/021
  * Kaplan et al. (2011) —. 2011, Journal of Cosmology and Astroparticle Physics, 2011, 011, doi: 10.1088/1475-7516/2011/10/011
  * Kunze & Komatsu (2014) Kunze, K. E., & Komatsu, E. 2014, JCAP , 2014, 009, doi: 10.1088/1475-7516/2014/01/009
  * Latif et al. (2019) Latif, M. A., Lupi, A., Schleicher, D. R. G., et al. 2019, Monthly Notices of the Royal Astronomical Society, 485, 3352, doi: 10.1093/mnras/stz608
  * Lepp & Shull (1983) Lepp, S., & Shull, J. M. 1983, Astrophys. J., 270, 578, doi: 10.1086/161149
  * Lepp & Shull (1984) —. 1984, The Astrophysical Journal, 280, 465, doi: 10.1086/162013
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, The Astrophysical Journal, 538, 473–476, doi: 10.1086/309179
  * Markevitch et al. (2004) Markevitch, M., Gonzalez, A. H., Clowe, D., et al. 2004, The Astrophysical Journal, 606, 819, doi: 10.1086/383178
  * Mitsou (2015) Mitsou, V. A. 2015, Journal of Physics: Conference Series, 651, 012023, doi: 10.1088/1742-6596/651/1/012023
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge New York: Cambridge University Press).  https://ui.adsabs.harvard.edu/abs/2010gfe..book.....M
  * Ostriker (2003) Ostriker, J. P. 2003, Science, 300, 1909–1913, doi: 10.1126/science.1085976
  * Padmanabhan & Finkbeiner (2005) Padmanabhan, N., & Finkbeiner, D. P. 2005, Phys. Rev. D, 72, 023508, doi: 10.1103/PhysRevD.72.023508
  * Palla et al. (1983) Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, Astrophys. J., 271, 632, doi: 10.1086/161231
  * Peebles (1968) Peebles, P. J. E. 1968, Astrophys. J., 153, 1, doi: 10.1086/149628
  * Pequignot et al. (1991a) Pequignot, D., Petitjean, P., & Boisson, C. 1991a, Astron. Astrophys., 251, 680
  * Pequignot et al. (1991b) —. 1991b, Astron. Astrophys., 251, 680
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, Astrophys. J., 187, 425, doi: 10.1086/152650
  * Rackauckas & Nie (2017) Rackauckas, C., & Nie, Q. 2017, The Journal of Open Research Software, 5, doi: 10.5334/jors.151
  * Randall & Scholtz (2015) Randall, L., & Scholtz, J. 2015, JCAP , 2015, 057, doi: 10.1088/1475-7516/2015/09/057
  * Ryan et al. (2022a) Ryan, M., Gurian, J., Shandera, S., & Jeong, D. 2022a, Astrophys. J., 934, 120, doi: 10.3847/1538-4357/ac75ef
  * Ryan et al. (2022b) Ryan, M., Shandera, S., Gurian, J., & Jeong, D. 2022b, Astrophys. J., 934, 122, doi: 10.3847/1538-4357/ac75e5
  * Savin et al. (2004) Savin, D. W., Krstic, P. S., Haiman, Z., & Stancil, P. C. 2004, arXiv:astro-ph/0404288, doi: 10.1086/421108
  * Schneider et al. (2013) Schneider, A., Smith, R. E., & Reed, D. 2013, Monthly Notices of the Royal Astronomical Society, 433, 1573–1587, doi: 10.1093/mnras/stt829
  * Seager et al. (1999) Seager, S., Sasselov, D. D., & Scott, D. 1999, The Astrophysical Journal, 523, L1–L5, doi: 10.1086/312250
  * Seljak (1994) Seljak, U. 1994, The Astrophysical Journal, 435, L87, doi: 10.1086/187601
  * Sethi & Subramanian (2005) Sethi, S. K., & Subramanian, K. 2005, Mon. Not. R. Astron. Soc., 356, 778, doi: 10.1111/j.1365-2966.2004.08520.x
  * Shandera et al. (2018) Shandera, S., Jeong, D., & Gebhardt, H. S. G. 2018, Physical Review Letters, 120, 241102, doi: 10.1103/PhysRevLett.120.241102
  * Shapiro & Kang (1987) Shapiro, P. R., & Kang, H. 1987, Astrophys. J., 318, 32, doi: 10.1086/165350
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, Monthly Notices of the Royal Astronomical Society, 308, 119–126, doi: 10.1046/j.1365-8711.1999.02692.x
  * Singh et al. (2020) Singh, D., Ryan, M., Magee, R., et al. 2020, arXiv:2009.05209 [astro-ph].  http://arxiv.org/abs/2009.05209
  * Sommer-Larsen & Dolgov (2001) Sommer-Larsen, J., & Dolgov, A. 2001, The Astrophysical Journal, 551, 608–623, doi: 10.1086/320211
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Physical Review Letters, 84, 3760–3763, doi: 10.1103/physrevlett.84.3760
  * Springel et al. (2008) Springel, V., Wang, J., Vogelsberger, M., et al. 2008, Mon. Not. R. Astron. Soc., 391, 1685, doi: 10.1111/j.1365-2966.2008.14066.x
  * Stancil et al. (1993) Stancil, P. C., Babb, J. F., & Dalgarno, A. 1993, Astrophys. J., 414, 672, doi: 10.1086/173113
  * Stancil et al. (1998) Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, Astrophys. J., 509, 1, doi: 10.1086/306473
  * Stiavelli (2009) Stiavelli, M. 2009, The First Stars (John Wiley & Sons, Ltd), 5–38, doi: https://doi.org/10.1002/9783527627363.ch2
  * Yoshida et al. (2006) Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, The Astrophysical Journal, 652, 6–25, doi: 10.1086/507978
  * Zaldarriaga & Harari (1995) Zaldarriaga, M., & Harari, D. D. 1995, Phys. Rev. D, 52, 3276, doi: 10.1103/PhysRevD.52.3276