# Molecular Chemistry for Dark Matter

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
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  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  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

###### Abstract

Molecular cooling is essential for studying the formation of sub-structure of
dissipative dark-matter halos that may host compact objects such as black
holes. Here, we analyze the reaction rates relevant for the formation,
dissociation, and transition of hydrogenic molecules while allowing for
different values of the physical parameters: the coupling constant, the proton
mass, and the electron mass. For all cases, we re-scale the reaction rates for
the standard molecular hydrogen, so our results are valid as long as the dark
matter is weakly coupled and one of the fermions is much heavier than the
other. These results will allow a robust numerical treatment of cosmic
structure, in particular for mini-halos for which molecular cooling is
important, in a dissipative dark matter scenario.

cosmology: theory – dark matter – molecular processes

##  1 Introduction

If the particle content of dark matter has enough complexity to support its
own chemistry, the universe may contain a much richer array of structures than
surveys of stars and galaxies have already revealed. Dark matter with
chemistry can be dissipative, so that in sufficiently dense environments it
will radiate kinetic energy and cool, allowing compact structures to form.
While gravitational evidence on galactic scales indicates that any
interactions among dark-matter particles must be somewhat weaker than those
observed in the visible matter (Markevitch et al., 2004), some small-scale
data may be better explained by dark matter with self interactions (Bullock &
Boylan-Kolchin, 2017). Future observations of smaller-scale structure,
including the possible compact-object detections by gravitational wave
observatories, will further probe for dissipative dark-matter interactions
(Kouvaris & Tinyakov, 2011; de Lavallaz & Fairbairn, 2010; Bramante & Elahi,
2015; Bramante & Linden, 2014; Kouvaris et al., 2018; Shandera et al., 2018;
Latif et al., 2019; Abbott et al., 2018, 2019; Singh et al., 2021; Hippert et
al., 2021) and may ultimately resolve the dark matter puzzle.

As an illustrative and calculable example, we consider a dissipative dark
matter model consisting of two oppositely charged fundamental fermions with
masses $m$ and $M$, interacting via a $U(1)$ force mediated by a dark photon.
The interaction strength is determined by a coupling constant, $\alpha$. We
take $m\ll M$ and the dark photon to be massless. For simplicity, we assume
that any non-gravitational interactions between dark matter and the Standard
Model particles are too weak to be relevant. Here, Standard Model refers to
the Standard Model of particle physics, but, hereafter, we shall use the term
Standard Model for referring to atomic and molecular physics of just electrons
and protons with Standard Model masses and electromagnetic interaction
strength. Other dissipative scenarios, especially “mirror" dark matter where
the dark sector has a copy of the particle content of the Standard Model, have
also been studied in the literature (Berezhiani et al., 2001; Chang et al.,
2019; Dessert et al., 2019; D’Amico et al., 2017; Choquette et al., 2019).

In the scenario we consider, pairs of oppositely charged particles ($M$ and
$m$) can form bound states analogous to hydrogen atoms, we call them dark
atoms, and a gas of these atoms can undergo all of the familiar
electromagnetic processes, including energy-level transition, ionization, and
recombination. In addition to the free-free (bremsstrahlung) emission from the
ionized states, dark atoms can dissipate energy by radiating dark photons
through recombination or collisional excitation followed by radiative decay
(Rosenberg & Fan, 2017; Buckley & DiFranzo, 2018). Because the neutral, bound
state atoms play an important role, this dark matter model is often called
“atomic" dark matter (Goldberg & Hall, 1986; Ackerman et al., 2009; Feng et
al., 2009; Kaplan et al., 2010, 2011; Fan et al., 2013; Cyr-Racine &
Sigurdson, 2013; Cyr-Racine et al., 2014; Cline et al., 2014; Foot & Vagnozzi,
2015, 2016; Randall & Scholtz, 2015; Agrawal et al., 2017a; Boddy et al.,
2016; Ghalsasi & McQuinn, 2018).

However, the energy dissipation allowed by atomic processes requires excited-
state atoms and free electrons, which are rare at low gas temperature where
essentially all atoms are neutral and in the ground state. So, atomic cooling
is inefficient at low temperatures. In the Standard Model, for example,
cooling below ${10}^{4}\text{\,}\mathrm{K}$ proceeds mostly through line
cooling of molecular hydrogen (Glover, 2012; Mo et al., 2010). The lower-lying
and more closely spaced molecular levels (see Sec. 2) allow the gas to cool to
approximately $100\text{\,}\mathrm{K}$. The Jeans mass at this minimum
temperature sets the mass scale of the collapsing gas cloud, which in the case
of Pop. III stars is correlated with the eventual mass of the proto-stars
(Bromm et al., 1999, 2002; Hirano et al., 2014).

The molecular and atomic cooling processes for hydrogen are well-studied in
the Standard Model. Of particular importance to us are the atomic and
molecular processes for Population III star formation that take place within a
pristine gas that is 92% hydrogen. In this paper we derive the chemical
reaction rates and cooling functions relevant for the analogous processes in
an atomic-dark-matter gas containing all species (free particles, atoms, ions,
and molecules). Many of the atomic processes were derived in Rosenberg & Fan
(2017), so our focus is on molecular physics.

The following three ingredients are necessary to understand the cooling
processes of a gas of atoms and molecules (Flower, 2007; Draine, 2011): (1)
the abundance of the gas comprising each species, (2) the complete set of
quantum states available to each component of the gas and how those states are
populated, and (3) the interaction potentials between each species. To obtain
the information from first principles requires solving multi-particle
Schrödinger equations. Even within the Standard Model the solutions to these
equations, and the molecular states and scattering rates, can only be found
approximately and numerically. The cooling functions used in simulations of
structure formation are typically semi-analytic, and significant uncertainties
in various rates remain. The state of the art is reviewed, for example, in
Galli & Palla (2013); Glover & Abel (2008); Glover (2015a).

Instead of obtaining the molecular wave functions, energies, and interactions
from first principles, in this paper we model dark molecular processes based
on the known results for Standard Model hydrogen. Fortunately, a little
dimensional analysis goes a long way. For example, in the next section
(Section 2), we review the Born-Oppenheimer approximation for the quantum
mechanical calculation of hydrogenic molecules and show how those results can
be re-scaled by ratios of the dark matter parameter values to the Standard
Model values:

$\displaystyle
r_{m}=\frac{m}{$511\text{\,}\mathrm{keV}$}\,,r_{M}=\frac{M}{$0.938\text{\,}\mathrm{GeV}$}\,,r_{\alpha}=\frac{\alpha}{137^{-1}}\,.$
(1.1)

We then argue that the re-scaling that we derive for the analytic but
approximate Born-Oppenheimer treatment largely carries over to the more
accurate but semi-analytic results used in the literature.111Appendix A
compares our re-scaling procedure to a simpler mass re-scaling technique used
in literature on Standard Model reaction rates. Sections 3 and 4 present the
cooling functions and reaction rates obtained using these re-scaling
techniques. Note that we assume $c=\hbar=1$ throughout except where noted.

The results we present in this work have a broad range of potential
applications including a more accurate treatment of Standard Model deuterated
chemistry. However, we developed the re-scaling scheme to compute two specific
aspects of dissipative dark matter scenarios: (A) the cosmological evolution
of the dark molecular hydrogen fraction and (B) formation of compact objects
from dark molecular cooling. We apply the outcome of this paper to these cases
in two companion papers, determining the early-time makeup of the cosmological
gas in Gurian et al. (2022) and modeling the dark-matter halo cooling from
both dark hydrogen atoms and molecules in Ryan et al. (2022).

While many of the results in this paper are independent of cosmology, our
eventual goal is to work out the observational consequences if some or all of
the dark matter is atomic dark matter. For that purpose, we make the following
four assumptions: (A) dark charge neutrality, (B) sufficiently fast
recombination for “atomic" processes to be relevant, (C) thermalized particle
distribution, and (D) collisional-ionization equilibrium. That is, we assume
an equal number of $M$ and $m$ particles, as doing otherwise would result in
long-range $U(1)$ interactions overwhelming the known bulk gravitational
behavior to alter, for example, the Friedmann equation. We assume a
recombination rate large enough that there is a cosmologically significant
fraction of dark atoms. We undertake a detailed treatment of the cosmological
recombination of atomic dark matter in a companion paper (Gurian et al.,
2022). Briefly, Saha equilibrium implies that the dark atoms will recombine
when the dark radiation temperature drops sufficiently far below the atomic
binding energy, as long as the recombination rate exceeds the Hubble rate.
This is the case for the entire range of parameters studied in our companion
papers. Third, we assume the dark particles have approximately-Maxwellian
velocity distributions at a common temperature. Agrawal et al. (2017b)
generically validate this assumption in the cosmological context. In the case
of cooling atomic dark matter halos, Fan et al. (2013); Rosenberg & Fan (2017)
delineate a regime where equi-partition is maintained by elastic collisions
with the dominant dark hydrogen species. We expect equi-partition must be
maintained well beyond this regime by other mechanisms (i.e. Landau damping
(Shu, 1992)). Fundamentally, as long as $\alpha\gg GM^{2}$, the U(1) force
will prevent the large scale charge-separations implied by differential
cooling of heavy and light particles in a gravitational potential. Lastly, we
calculate the cooling rates assuming collisional-ionization equilibrium and
ignoring the effects of photons produced through cooling processes, like
potential reabsorption or reionization. In the absence of strong radiation
sources (Wiersma et al., 2009), this treatment is equivalent to assuming the
optically thin gas where the photons escape the medium. This turns out to be a
good approximation in the Standard Model (Lykins et al., 2013).

##  2 Basic tools for dark molecular physics

The hydrogenic molecule is a system of two heavy particles of mass $M$ and
charge $+1$, at positions ${\bf X}_{A}$, ${\bf X}_{B}$, and two light
particles of mass $m$ and charge $-1$, at positions ${\bf x}_{1}$, ${\bf
x}_{2}$, all interacting electromagnetically, with the fine-structure constant
$\alpha$. The Schrödinger equation for stationary states of this system
depends on the distance between each pair of particles (e.g., $x_{1A}$ between
light particle 1 and heavy particle $A$) as

$\displaystyle H\Psi({\bf X}_{A},{\bf X}_{B};{\bf x}_{1},{\bf
x}_{2})=E\Psi({\bf X}_{A},{\bf X}_{B};{\bf x}_{1},{\bf x}_{2})\,,$ (2.1)

with Hamiltonian

$\displaystyle H$
$\displaystyle=-\frac{1}{2M}(\nabla^{2}_{A}+\nabla^{2}_{B})-\frac{1}{2m}(\nabla^{2}_{1}+\nabla^{2}_{2})$
$\displaystyle+\alpha\left(\frac{1}{X_{AB}}+\frac{1}{x_{12}}-\frac{1}{x_{1A}}-\frac{1}{x_{2A}}-\frac{1}{x_{1B}}-\frac{1}{x_{2B}}\right)\,.$
(2.2)

Here, $\Psi({\bf X}_{A},{\bf X}_{B};{\bf x}_{1},{\bf x}_{2})$ is the wave
function for the stationary state with energy $E$, and we use $\hbar=1$. We
use the shorthanded notation for the distance $X_{AB}\equiv|{\bf X}_{A}-{\bf
X}_{B}|$, and so on. As long as $M\gg m$, we can separate the Schrödinger
equation into electronic part and nucleic parts by using the Born-Oppenheimer
approximation as follows. First, an ansatz for the wave functions of the light
particles (the electrons) is introduced assuming the heavy particles (the
protons) remain stationary (Bethe & Jackiw, 1997). The electronic wave
functions are only parametrically dependent on the nuclear separation $X_{AB}$
and can be used to find an effective potential between the two nuclei, which
is in turn solved for the nucleic part of the wave functions.

The simplest ansatz for the electronic ground state is the symmetric
combination of products of ground state wave functions $u_{is}$ for the
individual hydrogen atoms:

$\psi=\frac{1}{\sqrt{2}}\left(u_{1s}(x_{1A})u_{1s}(x_{2B})+u_{1s}(x_{2A})u_{1s}(x_{1B})\right)\,.$
(2.3)

This is the Heitler-London (Heitler & London, 1927) wave function for the
molecule, which depends on fundamental parameters only in the usual Bohr
radius combination $a_{0}=(\alpha m)^{-1}$. This result is exact when the two
hydrogen atoms are infinitely far apart, and provides the asymptotic form for
more accurate electronic wave functions at large nuclear separation. Assuming
this solution for the electronic wave functions, the effective potential for
the nuclei (called the Heitler-London potential) can be written analytically.
Its shape is very well approximated by the simpler Morse potential, a function
of the nuclear separation:

$\displaystyle V_{\rm
Morse}(X_{AB})=V_{0}(e^{-2(X_{AB}-X_{0})/b}-2e^{-(X_{AB}-X_{0})/b})\,.$ (2.4)

where $X_{0}$ is the nuclear separation at the minimum of the potential. The
typical scale for the parameters are $X_{0},b\propto a_{0}$,
$V_{0}\propto\alpha/a_{0}\propto m\alpha^{2}$.

The advantage of the Morse potential is that its energy levels can be found
exactly, and analytic expressions for the corresponding nuclear wave functions
can be obtained (Dong, 2007). The energies of the vibrational modes of the
molecule are given approximately by the bound state energies in the Morse
potential:

$E_{{\rm
vib},\nu}=-V_{0}\left(1-\frac{\nu+\frac{1}{2}}{K_{0}b}\right)^{2}\,,\;\;\;0\leq\nu\leq
K_{0}b-\frac{1}{2}\;,$ (2.5)

where $K_{0}=\sqrt{2MV_{0}}$. Below the energy required to excite a
vibrational mode, the separation of the two nuclei can be treated as fixed,
$X_{AB}=X_{0}$, and then the Hamiltonian becomes that of a rigid rotor. By
changing coordinates to those for the center of mass and relative motion, the
usual rotational states are found with energy levels

$E_{{\rm rot},J}=\frac{J(J+1)}{MX_{0}^{2}}\,.$ (2.6)

At this level of approximation, the molecular rotational energies and low-
lying vibrational energies ($\frac{\nu}{K_{0}b}\ll 1$) for generic masses
$m\ll M$ and coupling ($\alpha\ll 1$) can be found from the Standard Model
(SM) values via

$\displaystyle E_{\bf mol.\,vib}$
$\displaystyle=\left[\frac{r_{\alpha}^{2}r_{m}^{3/2}}{r_{M}^{1/2}}\right]E_{\bf
mol.\,vib,SM}\,,$ $\displaystyle E_{\bf mol.\,rot}$
$\displaystyle=\left[\frac{r_{\alpha}^{2}r_{m}^{2}}{r_{M}}\right]E_{\bf
mol.\,rot,SM}\,.$ (2.7)

The literature on molecular cooling and scattering uses more accurate
solutions for the molecular hydrogen energy eigenstates, found by introducing
more sophisticated ansatze for the electronic wave functions. These are
typically constructed from a basis set of trial wave functions, and then the
task is to solve for the coefficients of the basis states that contribute to
each stationary state. For example, the James-Coolidge basis (James &
Coolidge, 1933) for solving the Schrödinger equation is relatively simple and
is valid for small inter-nuclear distances. In this basis, the ansatz for the
electronic wave functions $u_{m,X_{0}}$ is written in terms of coordinates
$y_{ij}\equiv x_{ij}/X_{0}$ in units of the fixed nuclear separation, which up
to a constant is the same as units of the Bohr radius (e.g.,
$y_{1A}=x_{1A}/X_{0}\propto x_{1A}/a_{0}$):

$\displaystyle u_{m,X_{0}}$
$\displaystyle=\sum_{\\{n\\}}c_{\\{n\\}}(1+P_{AB})(1+P_{12})y_{12}^{n_{1}}(y_{1A}-y_{1B})^{n_{2}}$
$\displaystyle\times(y_{2A}-y_{2B})^{n_{3}}(y_{1A}+y_{1B})^{n_{4}}(y_{2A}+y_{2B})^{n_{5}}$
$\displaystyle\times e^{-\beta(y_{1A}+y_{1B})}e^{-\beta(y_{2A}+y_{2B})}$ (2.8)

where $P_{AB}$ ($P_{12}$) indicates the permutation of $A$ and $B$ (1 and 2),
and the coefficients $c_{\\{n\\}}$ and $\beta$ are variational parameters. The
individual indices $n_{i}$ are integers $\geq 0$, and the summation is over
all sets of indices that satisfy $\sum_{i=1}^{5}n_{i}\leq\Omega$ for
$\Omega=3,4,5...$, up to the desired precision. This basis was used in Kołos &
Roothaan (1960) in an early determination of electronic wave functions, and
more recently (e.g, Sims & Hagstrom (2006)) to make high-accuracy calculations
for the hydrogen molecule. Pachucki (2010) gives analytic expressions for many
integrals needed in molecular calculations using these wave functions. Other
classic literature, beginning with a paper by Kołos & Wolniewicz (1965)
generalizes the James-Coolidge basis to be more accurate at large inter-
nuclear distances, where it should asymptote to the Heitler-London solution.
But, none of these prescriptions introduce a new parametric dependence on the
fundamental parameters into the molecular potential, and so the re-scaling
shown in Eq. (2) should remain accurate.

Heuristically, we can obtain the re-scaling in Eq. (2) by using the spring
constant for the vibrational mode

$k=\left(\frac{{\rm d}^{2}V_{\rm Morse}}{{\rm
d}X^{2}}\right)_{X=X_{0}}=\frac{V_{0}}{b^{2}}\propto m^{3}\alpha^{4},$ (2.9)

which yields the vibrational energy level as

$E_{\bf
mol.\,vib}=\nu\hbar\omega=\nu\hbar\sqrt{\frac{k}{M}}\propto\sqrt{\frac{m}{M}}m\alpha^{2}\,,$
(2.10)

and the moment of inertia, $I=MX_{0}^{2}$, which yields the rotational energy
level as

$E_{\bf
mol.\,rot}=\frac{L^{2}}{2I}\propto\frac{1}{Ma_{0}^{2}}=\frac{m^{2}\alpha^{2}}{M}\,.$
(2.11)

Note that these heuristic results rely on the facts that the typical length
scale and energy scale for the hydrogen molecule are determined by the
electronic states in the framework of Born-Oppenheimer approximation. That is,
it is the Bohr radius and the hydrogen binding energy that set the relevant
scales, and the heuristic re-scaling must hold beyond the specific electronic
wave function being used. That argument further supports our extension of the
re-scaling in Eq. (2).

The spin part of the nuclear wave functions is also important for
understanding molecular properties. There is a singlet state, parahydrogen,
with total nuclear spin $S_{N}=0$, and a triplet state, orthohydrogen, with
total nuclear spin $S_{N}=1$. Since the total wave function must be
antisymmetric under the exchange of the two heavy fermions (protons) the anti-
symmetric singlet state (parahydrogen) can only have $J$-even rotational
states, while the symmetric triplet state (orthohydrogen) must have $J$-odd.
Because the majority of radiative transitions do not change the spin state,
and a photon carries angular momentum $\Delta J=1$, the leading order
radiative transitions of hydrogenic molecules are electric quadrupole
transitions, which are slower compared to the dipole transition by a factor of
$(a_{0}/\lambda_{u\ell})^{2}$ where $\lambda_{u\ell}$ is the wavelength
corresponding to the energy gap.

###  2.1 Molecular binding energy

Molecular transitions and reaction rates generally depend on the ratio of the
temperature to the energy scale associated with the reaction. For rotational
and vibrational transitions, we presented the parametric dependence of the
relevant energy scales in the previous section. On the other hand, formation
and dissociation processes of molecular states involve a sum of partial cross
sections over rovibrational states and so there is at least a formal
dependence on multiple energy scales. Nonetheless, we argue that the most
relevant energy scale for the formation and dissociation of hydrogenic
molecules is the binding energy $E_{\rm H}$ of hydrogen.

This is most obvious in the Heitler-London treatment of the problem. Writing
down the symmetric combination of ground state wave functions for the
individual hydrogen atoms, Eq. (2.3) leads to an expression for the
expectation value of the total energy of the molecule of the form (Heitler &
London, 1927; Dushman, 1936):

$E_{\text{H}_{2}}=E_{\rm H}[1+f(X_{AB}/a_{0})],$ (2.12)

from which the binding energy is determined by minimizing the energy over
$X_{AB}$ and comparing to the total energy of the separated hydrogen energies.
There are two key observations. First, although molecular dissociation
involves separating the nuclei, the binding energy is determined principally
by the electronic configuration. Second, this binding energy is (at lowest
order in perturbation theory) just a constant fraction of $E_{\rm H}$, and the
constant does not depend on mass ($m$, $M$) nor coupling constant $\alpha$.
This method calculates a binding energy of $3.2\text{\,}\mathrm{eV}$ (Heitler
& London, 1927) compared to the measured value of $4.74\text{\,}\mathrm{eV}$,
and the same approach applied to $\text{H}_{2}^{+}$ yields an energy of
$1.77\text{\,}\mathrm{eV}$ compared to the true value of
$2.64\text{\,}\mathrm{eV}$ (Li et al., 2007).

Unfortunately, the accuracy of this method is clearly less than excellent.
This is because the molecule is too tightly coupled for the Heitler-London
wave function to be a good approximation to the molecular electronic wave
function. Historically, this problem was first addressed by attaching some
ansatz together with a variational parameter to the individual electron wave
functions. One early such attempts is Rosen (1931). There, the individual
electronic wave functions are given as for example

$\psi_{1A}=u_{1s}(x_{1A})+\sigma u^{\prime}({\bf x}_{1A}),$ (2.13)

where $u_{1s}$ is the hydrogenic wavefunction for an effective (screened)
nuclear charge $Z$, and $u^{\prime}$ is the the wave function of the $2p$
state in a hydrogenic atom with charge $2Z$. Both $Z$ and $\sigma$ are
determined variationally. The full electronic wavefunction is taken as the
Heitler-London combination of the individual electronic wavefunctions. This
approach ultimately yields a ground state energy of
$4.02\text{\,}\mathrm{eV}$, which differs by 15.1% from the actual value, and
again clearly scales as $E_{\rm H}$.

To achieve percent-level accuracy requires a fully variational calculation,
while also including the previously neglected electron-electron separation, in
the James-Coolidge basis (James & Coolidge, 1933) (a calculation which James &
Coolidge (1933) note “can be evaluated and checked by an experienced computer
in a little over two hours.”). This variational calculation may well contain
some subdominant higher order energy parameter dependence which cannot be
easily extracted. Still, we can be encouraged by two facts. First, the
calculation is still concerned entirely with determining the electronic wave
function as a function of fixed internuclear separation (and then minimizing
the energy of that state with respect to the separation). Second, the
variational ansatz carries dimensions only through the Bohr radius. All this
makes clear that the molecular dissociation energy scales dominantly with the
atomic binding energy, and that $T/E_{\rm H}$ is the “dimensionless
temperature” relevant to these reactions. Experimental data suggests that this
conclusion is correct: the dissociation energy of $\text{H}_{2}$ differs from
that of $\text{D}_{2}$ by less than 2% (Herzberg & Monfils, 1961), consistent
with predominantly $E_{\rm H}$ dependence, insensitive to the mass $M$.

###  2.2 Spontaneous emission rate

Next, we consider the spontaneous emission rates, or Einstein $A$
coefficients. $\text{H}_{2}$ is a homonuclear molecule with no permanent
electric dipole moment, so transitions between states are quadrupolar, and the
hard work is to compute the quadrupole transition
$\langle\Psi_{\nu,J}|\hat{Q}|\Psi_{\nu^{\prime},J+2}\rangle$. Here, we bypass
this calculation by re-scaling known results from the Standard Model to
estimate the $A$ coefficient for the dark-atomic model.

The electric quadrupole radiation power is proportional to
$P\propto\omega^{6}Q^{2}$ (Jackson, 1998), so equating that with $A^{\rm
quad.}\hbar\omega$ we find the transition rate between states separated by
energy $\Delta E$ scales as222In the literature, this rate is often given in
atomic units. See Appendix B for a derivation of the re-scaling in those
units.

$\displaystyle A^{\rm quad.}$ $\displaystyle\propto\frac{\alpha(\Delta
E)^{5}}{\hbar^{5}c^{4}}\langle Q\rangle^{2}\,.$ (2.14)

The literature contains at least two definitions for the quadrupole moment
that differ by a factor of 2 (Poll & Wolniewicz, 1978; Aannestad, 1973), but
$Q\propto{\rm(length)}^{2}\propto a_{0}^{2}$ regardless. In the equation above
we have explicitly written the factors of $\hbar$ so that the units on the
right hand side are obviously $\rm{s}^{-1}$.

Then, the dark matter Einstein $A$ coefficient for transitions between
rotational levels can be estimated by

$\displaystyle A_{\rm rot,\,DM}\propto$ $\displaystyle\,\alpha E_{\rm\bf
mol.~{}rot}^{5}(a_{0}^{2})^{2}$ (2.15) $\displaystyle=$ $\displaystyle
r_{\alpha}\left[\frac{r_{\alpha}^{2}r_{m}^{2}}{r_{M}}\right]^{5}\left(\frac{1}{(r_{\alpha}r_{m})^{2}}\right)^{2}\;A_{\rm
rot,\;SM}$ $\displaystyle=$
$\displaystyle\left[\frac{r_{\alpha}^{7}r_{m}^{6}}{r_{M}^{5}}\right]\;A_{\rm
rot,\;SM}\,.$ (2.16)

Similarly, for transitions involving only the vibrational level, we estimate
the rate by

$\displaystyle A_{\rm vib,\,DM}\propto$ $\displaystyle\,\alpha E_{\rm\bf
mol.~{}vib}^{5}(a_{0}^{2})^{2}$ (2.17) $\displaystyle=$ $\displaystyle
r_{\alpha}\left[\frac{r_{\alpha}^{2}r_{m}^{3/2}}{r_{M}^{1/2}}\right]^{5}\left(\frac{1}{(r_{\alpha}r_{m})^{2}}\right)^{2}\;A_{\rm
vib,\;SM}$ $\displaystyle=$
$\displaystyle\left[\frac{r_{\alpha}^{7}r_{m}^{7/2}}{r_{M}^{5/2}}\right]\;A_{\rm
vib,\;SM}\,.$ (2.18)

###  2.3 Collisional excitation scattering rates

In the simplest scattering between hydrogen molecules
($\text{H}_{2}$-$\text{H}_{2}$), or between a hydrogen molecule and a hydrogen
atom ($\text{H}_{2}$-H), there is a change in the state of the molecule but no
exchange of particles, generally called collisional excitation. Here, we
discuss how to re-scale the rates for molecular scattering. We then generalize
the tools developed here to general scattering processes involved in molecular
chemical reactions in Section 4.

With molecular states and energies in hand from Eq. (2), the additional
ingredient needed for scattering rates is the interaction potentials between
two molecules, or between a molecule and a hydrogen atom. Takayanagi &
Nishimura (1960) used a Morse-type potential in a classic analytic treatment
of $\text{H}_{2}$-H scattering. The scattering rate coefficients,
$\gamma\equiv\left<\sigma v\right>$, are determined by the velocity-averaged
cross section $\sigma$ for changing the rotational states, where velocity $v$
is assumed to be drawn from a Maxwell-Boltzmann distribution at temperature
$T$. Takayanagi & Nishimura (1960) have denoted

$\displaystyle\gamma_{J,J+2}=$ $\displaystyle\langle\sigma
v\rangle_{J\rightarrow J+2}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{8}{\pi\mu
k_{B}^{3}T^{3}}}\int_{0}^{\infty}\sigma_{J\rightarrow
J+2}(E^{\prime})e^{-E^{\prime}/k_{B}T}E^{\prime}dE^{\prime}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{8k_{B}T}{\pi\mu}}\pi R_{c}^{2}\langle F\rangle$
(2.19)

where $E^{\prime}$ is the center of mass collision energy, $\mu$ is the
reduced mass of the $\text{H}_{2}-\text{H}$ pair, $R_{c}\propto a_{0}$ is the
distance of closest approach in a classical head-on collision. Assuming the
Morse-type potential, $\langle F\rangle$ is the dimensionless result of an
integral and is known analytically. Although this is a simple approach, it
remains reasonably accurate even as calculations for the interaction potential
have continued to improve. For example, Wrathmall & Flower (2007) compared a
Morse-type potential, although not identical to the one used in Takayanagi &
Nishimura (1960), to a very accurate, recent $\text{H}_{2}$-H potential and
found reasonable agreement. So, we will use the Eq.(2.19) to re-scale the
scattering rates.

The dimensionality of $\gamma$ comes from the two factors of length in the
cross section ($R_{c}\propto a_{0}$) and the thermal/kinematic factor from the
velocity average, $\sqrt{k_{B}T/\mu}$. Since the energy difference between the
initial and final levels remains parametrically important after the
integration, it is useful to define the dimensionless combination
$y^{2}=\Delta E/(k_{B}T)$. While $\langle F\rangle$ is dimensionless, it does
become a function of $y^{2}$, due, at minimum, to the presence of $T$ in the
exponential term. Then, if a scattering rate is known to scale with
temperature as $f(T)=f(\Delta E/y^{2})=f_{\Delta E}(y)$, and using the
dependence of $R_{C}\propto a_{0}$ on fundamental parameters and $\mu\propto
M$, the complete parameteric dependence is

$\gamma\propto\sqrt{\frac{\Delta E}{M}}\frac{1}{\alpha^{2}m^{2}}f_{\Delta
E}(y).$ (2.20)

Collisions in general cause transition between both vibrational and rotational
states. At each gas temperature, however, we can approximate the collisional
transitions by using the dominant transition: Cooling via de-excitation of
rotational levels in the vibrational ground state ($\nu=0$) dominates at lower
temperatures, and cooling via transitions between the lowest few vibrational
levels (averaging over the more closely spaced rotational levels) dominates at
higher temperatures.

In this approximation, we first consider transitions between rotational levels
within a fixed vibrational level. Let us define the dimensionless parameter

$y^{2}_{\rm rot}=\frac{\Delta E_{\rm
rot}}{k_{B}\,T}=\frac{m^{2}\alpha^{2}}{Mk_{B}\,T}\,.$ (2.21)

Then, the rotational transition rate may be written as

$\displaystyle\gamma\propto$
$\displaystyle\frac{1}{\alpha^{2}m^{2}}\sqrt{\frac{\Delta E_{\rm
rot}}{M}}f_{\rm rot}(y_{\rm rot})$ $\displaystyle\propto$
$\displaystyle\frac{1}{\alpha^{2}m^{2}}\sqrt{\frac{\alpha^{2}m^{2}}{M^{2}}}f_{\rm
rot}(y_{\rm rot})$ $\displaystyle=$ $\displaystyle\frac{f_{\rm rot}(y_{\rm
rot})}{\alpha mM}\,.$ (2.22)

The re-scaling for the dimensionful prefactor in the dark-matter model rate
coefficient is then obvious from Eq. (2.22). We can think of $y_{\rm rot}^{2}$
as the temperature in $\Delta E_{\rm rot}$ units, the relevant energy scale of
the problem. For a given dark-matter temperature $T_{\rm DM}$, we have to
evaluate the function $f_{\rm rot}$ at the corresponding Standard Model
temperature $T_{\rm SM}$ yielding the same $y^{2}_{\rm rot}$:

$T_{\rm SM}=T_{\rm DM}\;\left(\frac{\Delta E_{\rm rot,SM}}{\Delta E_{\rm
rot,DM}}\right).$ (2.23)

Using the rotational-energy re-scaling in Eq. (2), then, we can define a re-
scaled temperature,

$\displaystyle\tilde{T}_{r}\equiv$
$\displaystyle\left(\frac{r_{M}}{r_{m}^{2}r_{\alpha}^{2}}\right)\;T,$ (2.24)

and obtain the final dark matter rate coefficient as

$\displaystyle\gamma_{J,{\rm DM}}(T)\approx$
$\displaystyle\left[r_{\alpha}^{-1}\,r_{m}^{-1}\,r_{M}^{-1}\right]\;\gamma_{J,SM}(\tilde{T}_{r})\,.$
(2.25)

The considerations for molecule-molecule scattering are more complicated than
the molecular-atomic case, as the para-para, ortho-ortho, and ortho-para cases
must all be considered (McMahan et al., 1974; Silvera, 1980). In addition, the
relative orientations of the molecules is important. Still, it seems
reasonable that the basic re-scaling for molecule-molecule scattering should
be similar to the molecule-atom case: the cross section should increase as the
size of the molecule ($R_{c}\propto a_{0}$) increases, and the reduced mass in
the denominator be nearly $M$. We therefore use Eq. (2.25) to re-scale rates
involving a change in rotational level for both $\text{H}_{2}$-$\text{H}_{2}$
and $\text{H}_{2}$-H processes.

For collisions involving a change in the vibrational level, the calculation
proceeds along the same lines: the wave functions for the initial and final
states are different from above, but the interaction potential is the same.
So, defining

$y^{2}_{\rm vib}=\frac{m^{3/2}\alpha^{2}}{M^{1/2}k_{B}T}\,,$ (2.26)

assuming the collisional excitation rate similarly contains a dimensionless
function that now depends on $y_{\rm vib}$, $f_{\rm vib}(y_{\rm vib})$, and
including the remaining terms from Eq. (2.20), we find that

$\displaystyle\gamma_{\rm vib}\propto$
$\displaystyle\frac{1}{\alpha^{2}m^{2}}\sqrt{\frac{\Delta E_{\rm
vib}}{M}}f_{\rm vib}\left(y_{\rm vib}\right)$ $\displaystyle\propto$
$\displaystyle\frac{1}{\alpha^{2}m^{2}}\sqrt{\frac{m^{3/2}\alpha^{2}}{M^{3/2}}}f_{\rm
vib}\left(y_{\rm vib}\right)$ $\displaystyle=$ $\displaystyle\frac{f_{\rm
vib}(y_{\rm vib})}{\alpha\,m^{5/4}\,M^{3/4}}\,,$ (2.27)

for which the re-scaling for the atomic dark-matter becomes

$\displaystyle\gamma_{\rm vib,DM}(T)\approx$
$\displaystyle\left[r_{\alpha}^{-1}\,r_{m}^{-5/4}\,r_{M}^{-3/4}\right]\;\gamma_{\rm
vib,SM}(\tilde{T}_{v})\,,$ (2.28)

with

$\tilde{T}_{v}\equiv\left(\frac{r_{M}^{1/2}}{r_{m}^{3/2}r_{\alpha}^{2}}\right)\,T.$
(2.29)

Again, we apply this formula for vibrational transition rates from both
$\text{H}_{2}$-$\text{H}_{2}$ and $\text{H}_{2}$-H collisions.

###  2.4 Summary

Quantity | Dependence | Re-scaling
---|---|---
$a_{0}$ (Bohr radius) | $\frac{1}{\alpha\,m}$ | $r_{\alpha}^{-1}\,r_{m}^{-1}$
$E_{\rm H}$ (atomic energy level spacing) | $m\,\alpha^{2}$ | $r_{\alpha}^{2}\,r_{m}$
$E_{\rm rot}$ (Molecular rotational energy) | $\frac{\alpha^{2}\,m^{2}}{M}$ | $r_{\alpha}^{2}\,r_{m}^{2}\,r_{M}^{-1}$
$E_{\rm vib}$ (Molecular vibrational energy) | $\frac{\alpha^{2}\,m^{3/2}}{M^{1/2}}$ | $r_{\alpha}^{2}\,r_{m}^{3/2}\,r_{M}^{-1/2}$
${\bf d}$ (dipole moment) | $0$ | $0$
$A_{\rm rot}$ (quadrupolar rotational Einstein coefficient) | $\frac{\alpha^{7}\,m^{6}}{M^{5}}$ | $r_{\alpha}^{7}\,r_{m}^{6}\,r_{M}^{-5}$
$A_{\rm vib}$ (quadrupolar vibrational Einstein coefficient) | $\frac{\alpha^{7}\,m^{7/2}}{M^{5/2}}$ | $r_{\alpha}^{7}\,r_{m}^{7/2}\,r_{M}^{-5/2}$
$p_{ij}$ (polarizability) | $a_{0}^{3}=\frac{1}{m^{3}\,\alpha^{3}}$ | $r_{\alpha}^{-3}\,r_{m}^{-3}$
Table 1: Table of primary quantities important for reaction rates, their
dependence on the parameters $m$, $M$, and $\alpha$, and the powers of ratios
of parameters (see Eq.(1.1)) needed to re-scale the quantities.

We can extract the parametric dependence of many relevant processes in the
atomic-dark-matter model by re-scaling. As we have shown, reaction and
transition rates can at lowest order be treated as dependent only on some
length scale (proportional to $a_{0}$) and some energy scale (atomic,
rotational, or vibrational). Since $\text{H}_{2}$ has no dipole moment (${\bf
d}=0$), $\text{H}_{2}$ state transitions are quadrupolar, with the vibrational
and rotational Einstein coefficients derived above. We neglect higher moments,
which are less important for scattering rates and cooling processes. Lastly,
several reaction rates will require the ground state polarizability, a tensor
mediating the external electric field and the induced electric dipole moment,
$p_{ij}\propto a_{0}^{3}\propto 1/(m^{3}\,\alpha^{3})$. This is
straightforward to show for the H atom (see e.g. Sakurai & Napolitano (2017)
for a discussion using perturbation theory), but is more complicated for the
$\text{H}_{2}$ molecule (Kołos & Wolniewicz, 1967). We have listed the
parametric dependence and re-scaling equation for all of the aforementioned
quantities in Table 1.

##  3 Application: Rovibrational cooling for dark molecular hydrogen

Molecular hydrogen cooling becomes important at low temperatures and high
$\text{H}_{2}$ concentrations in the Standard Model (Glover, 2012). In these
conditions, molecular line cooling is the dominant process, with certain
$\text{H}_{2}$ production and destruction channels providing auxiliary
cooling. In this section, we calculate the molecular line cooling rate for the
atomic-dark-matter model by applying the re-scaling that we have derived in
Section 2.2-2.3. To start, we briefly derive the molecular line cooling
equations in Sections 3.1 and 3.2, connecting our notation to that in the
classic paper of Hollenbach & McKee (1979), whose analytic model forms the
basis of comparison with the more modern computational and empirical results
(Glover & Abel, 2008; Glover, 2015b; Galli & Palla, 1998) that we would prefer
to use. In Sections 3.3 and 3.4 we re-scale the equations based on the
dominant physical process, and in Section 3.5 we combine the equations to
obtain the full, re-scaled, rovibrational cooling rates. We present full
results based on Hollenbach & McKee (1979) as well as Glover & Abel (2008).
Since it is useful to gain an understanding of the rates using the analytic
results of Hollenbach & McKee (1979), who considered only on
$\text{H}_{2}-\text{H}$ and $\text{H}_{2}-\text{H}_{2}$ collisions, Sections
3.3 -3.5 concentrate on those as collision partners. However, at low particle
densities molecular line cooling can be driven by collisions with various
other species (Glover & Abel, 2008), which we discuss in Section 3.6.

###  3.1 Molecular line cooling

Molecular line cooling occurs when the coolant molecules in excited states
decay to lower-level states by spontaneously emitting photons. The cooling
rate per volume (in units of
$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{cm}}^{-3}$) can be
written as the sum over all radiative transitions:

${\cal C}=\sum_{u}n_{u}\sum_{\ell<{u}}A_{u\ell}\Delta E_{u\ell}\,,$ (3.1)

where $n_{u}$ is the number density of molecules in the upper-energy state
$u$, the $\ell$ are the lower-energy states available for the transition,
$A_{u\ell}$ is the Einstein $A$ coefficient, and $\Delta
E_{u\ell}=E_{u}-E_{\ell}$ is the energy carried away by the photon.

Evaluating the re-scaled cooling rate from Eq. (3.1) requires knowing three
pieces of information: the occupation numbers of the various
rotational/vibrational levels, the re-scaled Einstein $A$ coefficient, and the
re-scaled energy of each eigenstate. We have already derived the energies in
Section 2, so the energy differences are given approximately as

$\displaystyle\Delta E(\nu^{\prime},J\pm 2,\nu,J)$ $\displaystyle=\,$
$\displaystyle(E_{\rm vib}(\nu^{\prime})-E_{\rm vib}(\nu))+(E_{\rm rot}(J\pm
2)-E_{\rm rot}(J)),$ (3.2)

where the vibrational and rotational energies are approximated by Eq. (2). We
have also derived the re-scaling of the Einstein $A$ coefficient for the
electric quadrupole transition in Eqs. (2.16)-(2.18). The remaining quantity
to compute is the number density of molecules in the excited states, which
depends on the temperature and the collision rate with the colliding particle
$X$.

To simplify the discussion, let us first consider the simple two-level system.
If a gas is in (local) thermal equilibrium (LTE) at temperature $T_{\rm gas}$,
then the ratio of the population of the two states is given as the Boltzmann
factor:

$\frac{n_{u}^{\rm LTE}}{n_{\ell}^{\rm LTE}}=\frac{g_{u}}{g_{\ell}}\exp(-\Delta
E_{u\ell}/k_{B}T_{\rm gas})\,,$ (3.3)

where $g_{u,\ell}$ are the number of internal degrees of freedom for each
state. Considering those same levels as a two-state system subject to
collisional excitation, collisional de-excitation, and radiative decay
(without ambient radiation, or ignoring the stimulated emission or
absorption), then the population of the excited state changes according to

$\frac{dn_{u}}{dt}=n_{X}n_{\ell}\gamma_{\ell
u}-n_{X}n_{u}\gamma_{u\ell}-n_{u}A_{u\ell}$ (3.4)

where $n_{X}$ is the number density of the colliding species and $\gamma_{\ell
u}$ ($\gamma_{u\ell}$) is the rate coefficient for the collisional excitation
(de-exitation). The steady state solution ($dn_{u}/dt=0)$ is

$\frac{n_{u}}{n_{\ell}}=\frac{n_{X}\gamma_{\ell
u}}{n_{X}\gamma_{u\ell}+A_{u\ell}}\,.$ (3.5)

There are two limiting cases. First, when the radiative decay dominates over
the collisional de-excitation, ($A_{u\ell}\gg n_{X}\gamma_{u\ell}$), the
population density in the excited states is given as
$n_{u}=n_{\ell}n_{X}\gamma_{\ell u}/A_{u\ell}$, or the collisionally excited
molecules stay in the excited states for the mean lifetime of
$A_{u\ell}^{-1}$. On the other hand, when the collisional de-excitation
dominates over the radiative decay, ($n_{X}\gamma_{u\ell}\gg A_{u\ell}$), the
population at the excited state returns to the thermal solution in Eq. (3.3)
because $n_{u}=n_{\ell}\gamma_{\ell
u}/\gamma_{u\ell}=n_{\ell}g_{u}/g_{\ell}\exp\left(-\Delta
E_{u\ell}/k_{B}T_{\rm gas}\right)$ from the principal of detailed balance.
These limits are commonly denoted as the $n\rightarrow 0$ (since $n_{X}\ll
A_{u\ell}/\gamma_{u\ell}$) and the $\rm LTE$ cases, respectively.

For a general system with more than two levels, collisional excitation, de-
excitation and radiative decay can happen between any pairs of levels, so the
equilibrium condition for the $i$-state becomes

$\sum_{j\neq i}n_{j}n_{X}\gamma_{ji}+\sum_{j>i}n_{j}A_{ji}=\sum_{j\neq
i}n_{i}n_{X}\gamma_{ij}+\sum_{j<i}n_{i}A_{ij}\,.$ (3.6)

The solutions for the three-level case are given in Section 17.5 of Draine
(2011). For systems with more than three levels, the equations are typically
solved numerically.

As the number density of the colliding particle determines the collisional de-
excitation rate, it is useful to define a critical density for level $u$ of a
collision partner $X$, $n_{\rm crit,u}(X)$, as a benchmark for when
collisional de-excitation is important for determining the excited-level
population. Although different definitions exist, the one used in Draine
(2011), for example, (in the limit of no ambient radiation) is

$n_{\rm
crit,u}(X)\equiv\frac{\sum_{\ell<u}A_{u\ell}}{\sum_{\ell<u}\gamma_{u\ell}(X)}\,.$
(3.7)

That is, this is the density of a collision particle for which collisional de-
excitation equals radiative de-excitation.

The relationship between the cooling rate in Eq. (3.1) and the cooling
functions $L$ used in Hollenbach & McKee (1979), specific to a particular
colliding particle $X$, is:

$\displaystyle L_{u\ell}^{X}=$
$\displaystyle\frac{1}{n(X)}\frac{1}{n(\text{H}_{2})}\sum_{u}n_{u}\sum_{\ell<u}A_{u\ell}\Delta
E_{u\ell}$ $\displaystyle=$ $\displaystyle\frac{{\cal
C}}{n(X)n(\text{H}_{2})}\,,$ (3.8)

where $L$ has units of
$\mathrm{erg}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$. Note that
this definition is motivated by the fact that the upper-level population is
proportional to the density of the colliding particle $n_{X}$ in the low-
density limit, $n_{X}<n_{\rm crit,u}(X)$. For example, for a simple two-level
system, ${\cal C}=n_{\ell}n_{X}\gamma_{\ell u}\Delta E_{u\ell}$ at the low-
density limit, $n_{X}\ll n_{\rm crit}$, and ${\cal C}=n_{u}A_{u\ell}\Delta
E_{u\ell}$ at the high-density limit, $n_{X}\gg n_{\rm crit}$.

###  3.2 Molecular line cooling with multiple colliding particles

In the Standard Model, atomic hydrogen is the main collisional partner of
molecular hydrogen, simply due to the high abundances, but several other
species are also relevant for cooling, including other hydrogen molecules and
helium.. When there is more than one colliding species then, analogous to Eq.
(3.5), the steady-state population of the upper state (in the two-level
approximation where we assume every pair of states separately satisfy the
stationary condition) may be approximated as

$\frac{n_{u}}{n_{\ell}}=\frac{\sum_{X}n_{X}\gamma_{\ell
u}(X)}{\sum_{X}n_{X}\gamma_{u\ell}(X)+A_{u\ell}}\,,$ (3.9)

with which the cooling rate can be written as

${\cal C}=\sum_{u}\sum_{\ell<u}\frac{\sum_{X}n_{X}\gamma_{\ell
u}(X)}{\sum_{X}n_{X}\gamma_{u\ell}(X)+A_{u\ell}}n_{\ell}A_{u\ell}\Delta
E_{u\ell}\,.$ (3.10)

Let us begin by considering the two limiting cases for the density of
colliding partners: the low-density and high-density limits.

When the density of colliding particles is small ($n\to 0$ limit) so that the
de-excitation happens only radiatively, the cooling rate becomes,

${\cal C}_{n\to 0}\equiv\sum_{X}{\cal C}_{n\to
0}^{X}=\sum_{X}\sum_{u}\sum_{\ell<u}n_{X}n_{\ell}\gamma_{\ell u}(X)\Delta
E_{u\ell}\,,$ (3.11)

or every collision leads to an emission of photon with energy $\Delta
E_{u\ell}$.

In the opposite ($n\to\infty$) limit, we recover the LTE cooling rate as
follows. Because the $\text{H}_{2}$ molecules must be in the LTE state when
the collision is frequent enough, and this must be true even with a single
colliding particle species, the collisonal excitation rate and the collisional
de-excitation rate must satisfy individually the detailed balance relation:

$\frac{\gamma_{\ell u}(X)}{\gamma_{u\ell}(X)}=\frac{n_{u}^{\rm
LTE}}{n_{\ell}^{\rm LTE}}=\frac{g_{u}}{g_{\ell}}e^{-E_{u\ell}/k_{B}T}\,.$
(3.12)

That means when collisional de-excitation dominates over the radiative decay,

$A_{u\ell}\ll\sum_{X}n_{X}\gamma_{u\ell}(X)\,,$ (3.13)

we recover the LTE cooling rate

$\displaystyle{\cal C}_{\rm LTE}=$
$\displaystyle\,\sum_{u}n_{u}\sum_{\ell<u}A_{u\ell}\Delta E_{u\ell}$
$\displaystyle=$ $\displaystyle\,\frac{n(\text{H}_{2})}{\cal
Z}\sum_{u}g_{u}e^{-E_{u}/k_{B}T}\sum_{\ell<u}A_{u\ell}\Delta E_{u\ell}\,,$
(3.14)

where $n_{u}=(n(\text{H}_{2})/{\cal Z})g_{u}e^{-E_{u}/k_{B}T}$ with the
partition function ${\cal Z}=\sum_{i}g_{i}e^{-E_{i}/k_{B}T}$.

In the intermediate density regime, rather than fully solving for the level
population as a function of density of all colliding particles [Eq. (3.6) or
its two-level approximation in Eq. (3.9)], it is a common practice, for
example in Hollenbach & McKee (1979); Galli & Palla (1998); Glover & Abel
(2008), to interpolate between the low-density and high-density (LTE) cases by

$\displaystyle{\cal C}\simeq\frac{{\cal C}_{\rm LTE}}{1+{\cal C}_{\rm
LTE}/\sum_{X}{\cal C}_{n\to 0}^{X}}\,,$ (3.15)

so that ${\cal C}$ approaches to the correct limit at both low-density and
high-density limits.

To understand this approximation better, let us explore the case when the
approximation becomes exact. Starting from the two-level approximation in Eq.
(3.10), and using Eq. (3.12), we find the cooling rate as

$\displaystyle{\cal C}=$
$\displaystyle\sum_{u}\sum_{\ell<u}\frac{\sum_{X}n_{X}\gamma_{\ell
u}(X)}{\sum_{X}n_{X}\gamma_{u\ell}(X)+A_{u\ell}}n_{\ell}A_{u\ell}\Delta
E_{u\ell}$ $\displaystyle=$
$\displaystyle\sum_{u}\sum_{\ell<u}\frac{[n_{u}^{\rm LTE}/n_{\ell}^{\rm
LTE}]n_{\ell}/n_{u}}{1+A_{u\ell}/\sum_{X}[n_{X}\gamma_{u\ell}(X)]}n_{u}A_{u\ell}\Delta
E_{u\ell}\,.$ (3.16)

That is, in two-level approximation, the interpolation in Eq. (3.15) is exact
if the level population follows the LTE value, for which

$\displaystyle\frac{A_{u\ell}}{\sum_{X}[n_{X}\gamma_{u\ell}(X)]}=$
$\displaystyle\,\frac{A_{u\ell}n_{u}\Delta
E_{u\ell}}{\sum_{X}[n_{X}\gamma_{u\ell}(X)n_{u}\Delta E_{u\ell}]}$
$\displaystyle=$ $\displaystyle\,\frac{{\cal C}_{\rm LTE}^{\rm
2-level,\,LTE}}{{\cal C}_{n\to 0}^{\rm 2-level,\,LTE}}\,.$ (3.17)

In reality, this approximation may work better for $T<\Delta E$ so that only a
few lowest-level states are populated. This is particularly true for the
rotational state where $\Delta E\propto J$, but may be little worse for the
vibrational modes where the energy gap stays constant.

###  3.3 Rotational cooling rates

We must also account for the type of energy transition (rotational versus
vibrational), separate from colliding particle density, since the results of
Section 2 have demonstrated they have different re-scalings. We first consider
rotational transitions between levels $J$ and $J-2$, in the vibrational ground
state $\nu=0$. Let $\gamma_{JJ^{\prime}}^{\text{H}}(T)$,
$\gamma_{JJ^{\prime}}^{\text{H}_{2}}(T)$ be the rate coefficient for the
($J\to J^{\prime}$) transition to occur due to a collision with, respectively,
a hydrogen atom and molecule.

At the low-density ($n\to 0$) limit, with much suppressed collisional
excitation, only the ground states (of both the para, $J=0$, and ortho, $J=1$,
type molecules) are populated, and the dominant excitations are to the $J=2$
and $J=3$ states. Beginning with Eq. (3.1) and Eq. (3.1), we have

$\displaystyle{\cal C}_{\rm{rot},n\rightarrow 0}^{X}$ $\displaystyle=\,$
$\displaystyle n_{J=2}A_{20}\Delta E_{20}+n_{J=3}A_{31}\Delta E_{31}$
$\displaystyle=\,$ $\displaystyle n_{J=0}n(X)\gamma^{X}_{02}\Delta
E_{20}+n_{J=1}n(X)\gamma^{X}_{13}\Delta E_{31}$ $\displaystyle=\,$
$\displaystyle\left(\frac{1}{4}\gamma^{X}_{02}\Delta
E_{20}+\frac{3}{4}\gamma^{X}_{13}\Delta E_{31}\right)n(\text{H}_{2})n(X)$
$\displaystyle=\,$ $\displaystyle\left(\frac{5}{4}\gamma^{X}_{20}e^{-\Delta
E_{20}/k_{B}T}\Delta E_{20}+\frac{7}{4}\gamma^{X}_{31}e^{-\Delta
E_{31}/k_{B}T}\Delta E_{31}\right)$ $\displaystyle\times
n(\text{H}_{2})n(X)\,,$ (3.18)

where we use the low-density limit of Eq. (3.5),
$n_{u}=(n_{\ell}n(X)\gamma_{\ell u})/A_{u\ell}$, to calculate the upper-level
population, and set $n_{J=0}/n(\text{H}_{2})=1/4$,
$n_{J=1}/n(\text{H}_{2})=3/4$ from the ortho-para ratio. The final equality
follows from the detailed-balance relation in Eq. (3.12). The total cooling
rate is given as the sum ${\cal C}=\sum_{X=\text{H},\text{H}_{2}}{\cal
C}^{X}$. Here we have assumed the commonly used ortho-para ratio of 3/1 and
use that for the comparison of cooling rates in Fig. (1). In general, however,
the balance between the two spin states may be different. In Appendix C we
derive the re-scaling for several processes that affect the ortho-para ratio
and show how they depend on the dark parameters.

Hollenbach & McKee (1979) provide an analytic form for a fit to the low-
density limit collisional de-excitation rates

$\displaystyle\gamma_{J,J-2}^{\text{H}}(T)$ $\displaystyle=$
$\displaystyle${10}^{-12}$\left(\frac{10T_{3}^{1/2}}{1+60T_{3}^{-4}}+T_{3}\right)$
$\displaystyle\times\left(0.33+0.9\exp\left[-\left(\frac{J-3.5}{0.9}\right)^{2}\right]\right)${\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$
$\displaystyle\gamma_{J,J-2}^{\text{H}_{2}}(T)$ $\displaystyle=$
$\displaystyle$3.3\text{\times}{10}^{-12}$\left(1+2T_{3}\right)$ (3.19)
$\displaystyle\times\left(0.276J^{2}\exp\left[-\left(\frac{J}{3.18}\right)^{1.7}\right]\right)${\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$

where $T_{3}=T/1000~{}{\rm K}$. Because the $J$-dependent part in the second
parenthesis is near unity for $J=3$, we can see that
$\gamma^{\text{H}_{2}}\gg\gamma^{\text{H}}$ for $T_{3}<1$. Hollenbach & McKee
(1979) reported at the time that Eq. (3.19) was $50~{}\%$ accurate for
$0.1<T_{3}<5$, above which we start to see the vibrational line cooling,
though subsequent results (e.g. (Glover & Abel, 2008; Lee et al., 2008; Lique,
2015)) have shown the rate actually differs significantly at low temperatures.
However, these rates, and the vibrational rates below, are still useful as an
analytic comparison to demonstrate the basic re-scaling behavior.

Using our result in Eq. (2.25), the dependence on microphysical parameters can
be captured by the re-scaling

$\displaystyle\gamma_{J,J-2,{\rm
DM}}^{\text{H}(\text{H}_{2})}(T)=\left[r_{\alpha}^{-1}\,r_{m}^{-1}\,r_{M}^{-1}\right]\;\gamma^{\text{H}(\text{H}_{2})}_{J,J-2,\rm{SM}}(\tilde{T}_{r}),$
(3.20)

for which we evaluate the Standard Model rates at re-scaled temperature
$\tilde{T}_{r}=r_{M}/(r_{\alpha}^{2}r_{m}^{2}T)$. Finally, substituting into
Eq. (3.3) along with the energy re-scaling in Eq. (2), gives the re-scaling
for the cooling rate

$\displaystyle\frac{{\cal C}_{{\rm rot},n\to 0,\rm
DM}^{X}(T)}{n(\text{H}_{2})n(X)}=\left[r_{\alpha}\,r_{m}\,r_{M}^{-2}\right]\;\frac{{\cal
C}_{{\rm rot},n\to 0,\rm DM}^{X}(\tilde{T}_{r})}{n(\text{H}_{2})n(X)}.$ (3.21)

At the high-density limit, collisional processes dominate so that all
rotational levels are in thermal equilibrium, and the cooling rate is given by
Eq. (3.14):

$\displaystyle{\cal C}_{\rm rot,LTE}(T,n(\text{H}_{2}))$ $\displaystyle=$
$\displaystyle\,n(\text{H}_{2})\sum_{J\geq
2}g_{J}\frac{\exp[-E_{J}/k_{B}T]}{{\cal Z}(T)}A_{J,J-2}\Delta E_{J,J-2}\,.$
(3.22)

The re-scaling results from Section 2 for the Einstein $A$ coefficient [Eq.
(2.16)] and the energy levels [Eq. (2)] can be used to determine the re-scaled
cooling rate in the LTE regime to account for changes in the microphysical
parameters as following:

${\cal
C}_{\rm{rot},\,LTE,DM}(T)=\left[r_{\alpha}^{9}\,r_{m}^{8}\,r_{M}^{-6}\right]\;{\cal
C}_{\rm rot,\,LTE,SM}(\tilde{T}_{r}).$ (3.23)

We then calculate the cooling rate in the intermediate density regime by
interpolation following Eq. (3.15).

###  3.4 Vibrational cooling rates

For cooling by vibrational transitions, we follow Hollenbach & McKee (1979)
who approximate the states as a three-level system ($\nu=2,1,0$), assuming all
levels have the same number of internal degrees of freedom
($g_{0}=g_{1}=g_{2}=1$). We also follow them in approximating the vibrational
motion as the simple harmonic oscillator, so that the energy levels are
proportional to $\nu$ (dropping the terms higher order in $\nu$): $\Delta
E_{10}=\Delta E_{21}=0.5\Delta E_{20}$. Because the rotational energy gap is
much smaller than the vibrational gap, we ignore any differences in energy if
the $J$ value changes in the transition.

The low-density limit for this system comes from a similar calculation as was
used above for the rotational case, except that the factor from degrees of
freedom is one (compared to $(2J+1)/(2J-3)$) and there is no distinction
between para and ortho states:

$\displaystyle\frac{{\cal C}_{\rm{vib},n\rightarrow
0}^{X}}{n(\text{H}_{2})n(X)\Delta E_{10}}$ (3.24) $\displaystyle=$
$\displaystyle\,\gamma_{10}^{X}e^{-\Delta
E_{10}/k_{B}T}+2\gamma^{X}_{20}e^{-\Delta
E_{20}/k_{B}T}+\gamma^{X}_{21}e^{-\Delta E_{21}/k_{B}T}\,,$

which is nearly the expression given by Hollenbach & McKee (1979) except that
we have included the $\gamma_{21}$ term dropped there. The expression
including the $\gamma_{21}$ term provides a better fit to the later results of
Glover & Abel (2008). As for the collisional de-excitation rates, we use
Hollenbach & McKee (1979)333Although the
$\gamma_{21}^{\text{H}_{2}},\,\gamma_{10}^{\text{H}_{2}}$ rate coefficients
have been updated in Hollenbach & McKee (1989), we use the original here and
in Figures 1 and 2 for simplicity.(rates are in
${\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$):

$\displaystyle\gamma_{10}^{\text{H}}$ $\displaystyle=$
$\displaystyle$1.0\text{\times}{10}^{-12}$T^{1/2}\exp\left[-\frac{1000}{T}\right]$
(3.25) $\displaystyle\gamma_{20}^{\text{H}}$ $\displaystyle=$
$\displaystyle$1.6\text{\times}{10}^{-12}$T^{1/2}\exp\left[-\left(\frac{400}{T}\right)^{2}\right]$
(3.26) $\displaystyle\gamma_{21}^{\text{H}}$ $\displaystyle=$
$\displaystyle$4.5\text{\times}{10}^{-12}$T^{1/2}\exp\left[-\left(\frac{500}{T}\right)^{2}\right]$
(3.27) $\displaystyle\gamma_{10}^{\text{H}_{2}}$ $\displaystyle=$
$\displaystyle$1.4\text{\times}{10}^{-12}$T^{1/2}\exp\left[-\frac{12000}{T+1200}\right]$
(3.28) $\displaystyle\gamma_{21}^{\text{H}_{2}}$ $\displaystyle=$
$\displaystyle\gamma_{10}^{\text{H}_{2}}$ (3.29)
$\displaystyle\gamma_{20}^{\text{H}_{2}}$ $\displaystyle=$ $\displaystyle 0$
(3.30)

where $T$ is in units of $\mathrm{K}$, then re-scale these rates by using Eq.
(2.28):

$\displaystyle\gamma_{\rm vib,DM}^{X}(T)\approx$
$\displaystyle\left[r_{\alpha}^{-1}\,r_{m}^{-5/4}\,r_{M}^{-3/4}\right]\;\gamma_{\rm
vib,SM}^{X}(\tilde{T}_{v})\,,$ (3.31)

with
$\tilde{T}_{v}\equiv\left(r_{M}^{1/2}/(r_{m}^{3/2}r_{\alpha}^{2})\right)T$.
Finally, substituting into Eq. (3.3) along with the energy re-scaling in Eq.
(2), gives the re-scaling for the cooling rate

$\displaystyle\frac{{\cal C}_{{\rm vib},n\to 0,\rm
DM}^{X}(T)}{n(\text{H}_{2})n(X)}=\left[r_{\alpha}\,r_{m}^{1/4}\,r_{M}^{-5/4}\right]\;\frac{{\cal
C}_{{\rm vib},n\to 0,\rm DM}^{X}(\tilde{T}_{v})}{n(\text{H}_{2})n(X)}.$ (3.32)

At the high-density limit, collisional processes dominate so that all
vibrational levels are in thermal equilibrium, and the cooling rate is given
by Eq. (3.14):

$\displaystyle{\cal C}_{\rm vib,LTE}(T,n(\text{H}_{2}))$ (3.33)
$\displaystyle=$
$\displaystyle\,n(\text{H}_{2})\sum_{\nu,\nu^{\prime}<\nu}g_{J}\frac{\exp[-E_{\nu}/k_{B}T]}{{\cal
Z}(T)}A_{\nu\nu^{\prime}}\Delta E_{\nu\nu^{\prime}}\,.$

From the re-scaling results from Section 2 for the Einstein $A$ coefficient
[Eq. (2.18)] and the energy levels [Eq. (2)], we find the re-scaled cooling
rate in the LTE regime:

${\cal
C}_{\rm{vib},\,LTE,DM}(T)=\left[r_{\alpha}^{9}\,r_{m}^{5}\,r_{M}^{-3}\right]\;{\cal
C}_{\rm vib,\,LTE,SM}(\tilde{T}_{v}).$ (3.34)

###  3.5 Rovibrational cooling

|
---|---
|
Figure 1: $\text{H}_{2}-\\{\text{H},\text{H}_{2}\\}$ collisional cooling
curves obtained from Hollenbach & McKee (1979) and Glover & Abel (2008) in the
low-density limit ($\Lambda_{n\rightarrow 0}^{\text{H},\text{H}_{2}}$). All
plots have $n(\text{H}_{2})=${10}^{-4}\text{\,}{\mathrm{cm}}^{-3}$$ and
$n(\text{H})=$1\text{\,}{\mathrm{cm}}^{-3}$$. The top left panel shows a
comparison of two cooling curves with Standard Model values. The top right
panel demonstrates the necessity of including both rotational and vibrational
re-scaling for atomic dark matter parameters
$(M,\,m,\,\alpha)=($40\text{\,}\mathrm{GeV}$,\,$40\text{\,}\mathrm{keV}$,\,$0.01$)$:
the Glover curves are shown using rotational only (green x), vibrational only
(cyan +), or both rotational and vibrational re-scaling (red dash-dot).
Clearly, using only one re-scaling fails in temperature regions with a
different dominant process. The lower Hollenbach and McKee curves in the upper
left plot, and all curves in the lower two plots, use the combined re-scaling
detailed in Appendix D. The “kink" in the $\text{H}_{2}$ curves arises from a
combination of approximating the $\text{H}_{2}$ transition temperature re-
scaling and the simple gap-filling method we use. The other parameter values
shown are
$(M,\,m,\,\alpha)=($100\text{\,}\mathrm{GeV}$,\,$10\text{\,}\mathrm{keV}$,\,$0.01$)$
(bottom left), and
$(M,\,m,\,\alpha)=($40\text{\,}\mathrm{GeV}$,\,$40\text{\,}\mathrm{keV}$,\,$0.001$)$
(bottom right).

To calculate the total cooling rate, we have to combine the rotational cooling
rate and vibrational cooling rate. First of all, using the approximation in
Eq. (3.15), we can combine the results of the previous sections by

$\displaystyle{\cal C}_{\rm tot}=$ $\displaystyle{\cal C}_{\rm rot}+{\cal
C}_{\rm vib}$ $\displaystyle=$ $\displaystyle\frac{{\cal C}_{\rm
rot,LTE}}{1+{\cal C}_{\rm rot,LTE}/\sum_{X}{\cal C}_{{\rm rot},n\to 0}^{X}}$
$\displaystyle+\frac{{\cal C}_{\rm vib,LTE}}{1+{\cal C}_{\rm
vib,LTE}/\sum_{X}{\cal C}_{{\rm vib},n\to 0}^{X}}\,,$ (3.35)

with $X={\text{H}},~{}{\text{H}}_{2}$ being the colliding particles. As
described in the previous sections, we use the collisional excitation rates
from Hollenbach & McKee (1979).

On the other hand, most recent literature combines the rotational and
vibrational cooling rates in the manner specified in Galli & Palla (1998)

${\cal C}_{\rm tot}=\frac{{\cal C}_{\rm{rot},LTE}+{\cal
C}_{\rm{vib},LTE}}{1+({\cal C}_{\rm{rot},LTE}+{\cal C}_{\rm{vib},LTE})/{\cal
C}_{{\rm tot},n\to 0}}\,,$ (3.36)

with

${\cal C}_{{\rm tot},n\to 0}=\sum_{X}\left({\cal C}_{\rm{rot},n\rightarrow
0}^{X}+{\cal C}_{\rm{vib},n\rightarrow 0}^{X}\right)\,.$ (3.37)

Since more recent rates have improved collisional coefficients (Galli & Palla,
1998; Lee et al., 2008; Lique, 2015), and include ortho-para conversions and
other improvements (Glover & Abel, 2008), we would prefer to use the rate
formulation given in Eq. (3.36). Our companion works (Gurian et al., 2022;
Ryan et al., 2022), for example, use the rates from Glover & Abel (2008);
Glover (2015b). This introduces complications, however, as the low-density
limit terms ${\cal C}_{\rm{rot},n\rightarrow 0}^{\text{H},\text{H}_{2}}$ and
${\cal C}_{\rm{vib},n\rightarrow 0}^{\text{H},\text{H}_{2}}$ are usually given
in the form of an analytic fit to the sum ${\cal C}_{\rm{rot},n\rightarrow
0}^{\text{H},\text{H}_{2}}+{\cal C}_{\rm{vib},n\rightarrow
0}^{\text{H},\text{H}_{2}}$, while the LTE terms are generally given
individually. This poses a non-trivial problem for our purpose, because, as
demonstrated in Sections 3.3 and 3.4, rotational and vibrational processes are
subject to different re-scalings. In Appendix D, we present a method for re-
scaling the total cooling rate with the assumption that rotational cooling
rate and vibrational cooling rate dominate at opposite ends of the temperature
range.

In Figure 1, we show the resulting low-density-limit total cooling rate of the
hydrogenic molecule colliding with H or $\text{H}_{2}$. For literature
comparison, we have plotted $\Lambda$ defined as

$\Lambda^{X}=\frac{{\cal C}^{X}_{{\rm rot},n\to 0}+{\cal C}^{X}_{{\rm
vib},n\to 0}}{\left(n(\text{H})+n(p)+2\,n(\text{H}_{2})\right)^{2}}\,.$ (3.38)

For this plot, we choose the values
$n(\text{H}_{2})=${10}^{-4}\text{\,}{\mathrm{cm}}^{-3}$$,
$n(\text{H})=$1\text{\,}{\mathrm{cm}}^{-3}$$, and neglect the ionized state
$n(p)$. While simply using just the rotational or just the vibrational re-
scaling is insufficient (top-right panel), the method that we outlined in
Appendix D provides results close to what we obtain by re-scaling the
individual rotational and vibrational formulas in Hollenbach & McKee (1979)
(and outlined in the previous section). This method, therefore, allows us to
re-scale the modern and more accurate cooling function from Glover & Abel
(2008). Note that the fit from Glover & Abel (2008) is only well-defined in
the range $$10\text{\,}\mathrm{K}$\leq T\leq${10}^{4}\text{\,}\mathrm{K}$$ and
that this range shifts proportionally to the dark parameters. The high
temperature behaviors are not a significant issue, however, as they occur
above the disassociation temperature, $4.48\text{\,}\mathrm{eV}$; hence, we do
not expect to have significant amounts of hydrogen molecules at that
temperature.

###  3.6 Collisions with other species

Figure 2: Combined atomic and molecular cooling rate for atomic dark matter
for the dark parameters listed, in the format
$(M[$\mathrm{GeV}$],m[$\mathrm{keV}$],\alpha,x_{\text{H}_{2}})$. The atomic
cooling rates were derived in Rosenberg & Fan (2017) and the $\text{H}_{2}$
fraction, $x_{\text{H}_{2}}$, represents the cosmological fraction found in
Gurian et al. (2022). For the Standard Model case, we have highlighted the
atomic (green circles) and molecular (cyan stars) contributions. Note that the
molecular cooling can be significantly reduced compared to atomic cooling as
the parameters are varied from the Standard Model values.

Atomic and molecular hydrogen are not the only collisional partners in the
Standard Model primordial gas. Glover & Abel (2008) showed that collisions
with He, protons, and electrons also contribute to cooling, although with
slower rates than the hydrogenic cooling due to their lower abundances. Even
at low abundances, however, the proton collisions can be important at low
temperatures, both as a direct additional source of cooling and through an
effect on the ortho-para ratio. The atomic dark matter model has no dark
neutrons and no dark helium. This section, therefore, focuses on re-scaling
the cooling rate involving collisions with protons and electrons.

The $\text{H}_{2}-p$ scattering rate re-scaling for a rotational transition
and the effect of proton collisions on the ortho-para ratio are covered in
more detail in Appendix C, where we assume the cross-section can be
approximated by a Langevin rate (Section 4.5 for more details). The re-scaling
for a vibrational transition follows from using the vibrational energy re-
scaling in Eq. C1, giving $\gamma_{{\rm DM},v\rightarrow
v^{\prime}}(T)=\left[r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\right]\gamma_{{\rm
SM},v\rightarrow v^{\prime}}(\tilde{T}_{v})$. The cooling rates are then
derived from the appropriate substitutions into Eqs. (3.3) and (3.24):

$\displaystyle\frac{{\cal C}_{{\rm rot},n\to 0,\rm
DM}^{p}(T)}{n(\text{H}_{2})n(p)}$
$\displaystyle=\left[r_{\alpha}\,r_{m}^{1/2}\,r_{M}^{-3/2}\right]\;\frac{{\cal
C}_{{\rm rot},n\to 0,\rm SM}^{p}(\tilde{T}_{r})}{n(\text{H}_{2})n(p)}$ (3.39)
$\displaystyle\frac{{\cal C}_{{\rm vib},n\to 0,\rm
DM}^{p}(T)}{n(\text{H}_{2})n(p)}$
$\displaystyle=\left[r_{\alpha}\,r_{M}^{-1}\right]\;\frac{{\cal C}_{{\rm
vib},n\to 0,\rm SM}^{p}(\tilde{T}_{v})}{n(\text{H}_{2})n(p)}.$ (3.40)

We also use the Langevin-type cross section for the $\text{H}_{2}-e$
collision, which leads to a rate re-scaling of $\gamma_{{\rm DM},j\rightarrow
j^{\prime}}(T)=\left[r_{\alpha}^{-1}\,r_{m}^{-2}\right]\gamma_{{\rm
SM},j\rightarrow j^{\prime}}(\tilde{T}_{r})$ and $\gamma_{{\rm
DM},v\rightarrow
v^{\prime}}(T)=\left[r_{\alpha}^{-1}\,r_{m}^{-2}\right]\gamma_{{\rm
SM},v\rightarrow v^{\prime}}(\tilde{T}_{v})$. Again, substituting
appropriately into Eqs. (3.3) and (3.24) we have

$\displaystyle\frac{{\cal C}_{{\rm rot},n\to 0,\rm
DM}^{e}(T)}{n(\text{H}_{2})n(e)}$
$\displaystyle=\left[r_{\alpha}\,r_{M}^{-1}\right]\;\frac{{\cal C}_{{\rm
rot},n\to 0,\rm SM}^{e}(\tilde{T}_{r})}{n(\text{H}_{2})n(e)}$ (3.41)
$\displaystyle\frac{{\cal C}_{{\rm vib},n\to 0,\rm
DM}^{e}(T)}{n(\text{H}_{2})n(p)}$
$\displaystyle=\left[r_{\alpha}\,r_{m}^{-1/2}\,r_{M}^{-1/2}\right]\;\frac{{\cal
C}_{{\rm vib},n\to 0,\rm SM}^{e}(\tilde{T}_{v})}{n(\text{H}_{2})n(e)}.$ (3.42)

Like the $\text{H}_{2}-\\{\text{H},\text{H}_{2}\\}$ collisions, the rotational
and vibrational cooling rates for both the $\text{H}_{2}-p$ and
$\text{H}_{2}-e$ collisions have different re-scalings, yet the Standard Model
rates are fits to the sum. We would prefer to account for this using the
approach in Appendix D, however, we have been unable to obtain a clear
separation between the two temperature regimes in the Standard Model rate. As
such, we have used the rotational scaling [Eqs. 3.39 and 3.41] for the
entirety of the temperature range. This leads to cooling that is stronger by a
factor of $\sqrt{r_{m}/r_{M}}$ in the vibrational temperature regime. We
believe this an acceptable approximation for two reasons: first, the rates
have the largest relative contribution in the low temperature regime (Glover &
Abel, 2008) where the vibrational portion is unlikely to be relevant, and more
importantly, we expect the overall contribution from $\text{H}_{2}-\\{e,p\\}$
cooling to be small due to relative abundances for the majority of cases. For
example, for the parameter values used in the top right panel of Figure 1, the
$\text{H}_{2}-\\{e,p\\}$ collisional cooling does not contribute meaningfully,
regardless of the re-scaling used, until $\sim$$1300\text{\,}\mathrm{K}$,
where much stronger atomic processes become relevant.

###  3.7 Summary of cooling results

Finally, in Figure 2, we show the entire atomic-dark-matter cooling rate
$\Lambda$, including both the molecular rovibrational cooling defined above
and the combination of cooling processes provided in Rosenberg & Fan (2017).
We have computed the rate for the same dark parameter sets used in Figure 1 to
demonstrate the significant behavioral changes induced by varying the
parameters. We have also varied the ionization fraction, assuming chemical
equilibrium, and the dark-atomic hydrogen fractions $x_{\text{H}_{2}}$ to
match the cosmological abundance obtained in Gurian et al. (2022).

##  4 Application: Reaction rates

The molecular line cooling in the previous section is proportional to the
abundance of H and $\text{H}_{2}$. In order to complete the molecular line
cooling calculation for the dark-atomic model, therefore, we have to determine
their abundances. These abundances are generally computed numerically from a
set of coupled differential rate equations similar to Eq. (3.4), for example,

$\frac{dn_{\text{H}_{2}}}{dt}=\sum_{j\in F_{i}}\left(\gamma_{j}\prod_{k\in
R_{j}}n_{k}\right)-\sum_{j\in D_{i}}\left(\gamma_{j}\prod_{k\in
R_{j}}n_{k}\right),$ (4.1)

where the first sum represents the source term including all chemical
reactions that contain $\text{H}_{2}$ as one of the products (Formation
reactions), and the second sum represents the sink term including all chemical
reactions that contain $\text{H}_{2}$ as one of the reactants (Destruction
reactions). Thus, to determine the species’ abundance evolution in the atomic-
dark-matter model, we need to know how to re-scale the reaction rates
($\gamma$) to give the leading dependence on the dark-matter parameters.

We shall present the re-scaling of these rates in this section. In the first
sub-section, we give an overview of how the re-scaling procedure works and
present a summary table of our results. In the remaining sub-sections we
provide the details.

###  4.1 Overview of method and summary of result

Just like the case for the collisional excitation rates that we studied in
Section 2.3, from a given standard-model rate $\gamma_{\rm SM}$, the
corresponding rate for the dark-matter model can be estimated by an overall
re-scaling factor multiplying $\gamma_{\rm SM}$ evaluated at a re-scaled
temperature $\tilde{T}$. That is,

$\displaystyle\gamma_{\rm DM}(T)=g(r_{\alpha},r_{m},r_{M})\;\gamma_{\rm
SM}\left(\tilde{T}_{\Delta}\right)\,,$ (4.2)

where $g(r_{\alpha},r_{m},r_{M})$ is the overall re-scaling, and
$\tilde{T}_{\Delta}=T/r_{\Delta E}$ with $r_{\Delta E}$ the re-scaling of the
primary binding energy relevant for the reaction.

For non-photochemistry reactions, the calculations of interaction rates are
similar to Eq. (2.19) from Section 2, but with the added possibility of
chemical reactions instead of just pure scattering. Assuming that all
particles in the gas are at the same temperature, the reaction rate is

$\gamma_{\rm DM,non-\gamma}=\langle\sigma
v\rangle\propto\sqrt{\frac{T}{\mu}}\int_{0}^{\infty}\sigma(x;y)x^{3}e^{-x^{2}}\text{d}x.$
(4.3)

Here, we reuse the two dimensionless parameters defined in Section 2:

$\displaystyle x^{2}$ $\displaystyle=$
$\displaystyle\frac{\text{K.E.}}{k_{B}\,T}=\frac{\mu\,v^{2}}{2\,k_{B}\,T}\,,$
(4.4) $\displaystyle y^{2}$ $\displaystyle=$ $\displaystyle\frac{\Delta
E}{k_{B}\,T}\,,$ (4.5)

where $\Delta E$ is the relevant binding energy of the chemical reactions.
Note that, for chemical reactions, we cannot assume that the cross-section
scales only as some effective size of the particles. Instead, we need to
determine in more detail how the cross section depends on kinetic energy and
binding energies, and how strongly the reaction depends on $\alpha$.

# | Reaction | Cross section source | $\sigma$ | Re-scaling pre-factor | b | Additional notes
---|---|---|---|---|---|---
$g(r_{\alpha},r_{m},r_{M})$
1 |  $p+e\rightarrow\text{H}+\gamma$(4.2) | 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 |  $\text{H}+\gamma\rightarrow p+e$(4.2) | 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 |  $\text{H}+e\rightarrow\text{H}^{-}+\gamma$ (4.3) | 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 | $\text{H}^{-}+\gamma\rightarrow\text{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 | $\text{H}^{-}+\text{H}\rightarrow\text{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 |  $\text{H}^{-}+p\rightarrow 2\text{H}$ (4.4) | 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 | $\text{H}+p\rightarrow\text{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 | $\text{H}_{2}^{+}+\gamma\rightarrow\text{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 |  $\text{H}_{2}^{+}+\text{H}\rightarrow\text{H}_{2}+p$(4.5) | 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 | $\text{H}_{2}^{+}+\text{H}_{2}\rightarrow\text{H}_{3}^{+}+\text{H}$ | —"— | —"— | —"— | 0 | g
15 |  $\text{H}_{2}+p\rightarrow\text{H}_{2}^{+}+\text{H}$(4.5) | —"— | —"— | —"— | 0 |  g , j
20 | $\text{H}_{3}^{+}+e\rightarrow\text{H}+\text{H}_{2}$ | Draine (2011) | $\frac{\alpha\,a_{0}}{\text{K.E.}}$ | $r_{\alpha}^{-1}\,r_{m}^{-2}\,r_{M}$ | -0.65 | k
|  $\text{H}_{2}+\text{H}\rightarrow 3\text{H}$ (4.6) | Hard Sphere | $a_{0}^{2}$ | $r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}$ | 0 | l
3B1 |  $3\text{H}\rightarrow\text{H}_{2}+\text{H}$ (4.7) | Hard Sphere/Detailed Balance | $a_{0}^{2}\left[\frac{n_{\text{H}_{2}}}{n^{2}_{\text{H}}}\right]_{\rm LTE}$ | $r_{\alpha}^{-4}\,r_{m}^{-4}\,r_{M}^{-1}$ | -1 |  j , m
3B2 |  $\text{H}_{2}+2\text{H}\rightarrow 2\text{H}_{2}$ (4.7) | —"— | —"— | —"— | -1 |  j , m
3B3 |  $2\text{H}+\text{H}^{+}\rightarrow\text{H}_{2}+\text{H}^{+}$ (4.7) | —"— | —"— | —"— | -1 |  j , m
3B4 |  $2\text{H}+\text{H}^{+}\rightarrow\text{H}_{2}^{+}+\text{H}$ (4.7) | —"— | —"— | —"— | -1 |  j , m

Table 2: 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 only those that were considered “important" in Galli & Palla (1998),
are required for our companion paper (Gurian et al., 2022), or are three-body
reactions. 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. (1991). 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\text{\,}\mathrm{K}$$, and very uncertain for
$T>$300\text{\,}\mathrm{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).

Even for chemical reactions, however, we can organize the scattering rate, Eq.
(4.3), as a product of the universal thermal factor, an effective “impact
parameter" $\tilde{R}$ (with dimensions of length) squared, and a
dimensionless integral, analogously to Eq. (2.19). The reaction-dependent part
of the pre-factor, $g(r_{\alpha},r_{m},r_{M})$, then depends on the re-scaling
of the “impact parameter". In general, this will be an algebraic combination
of the masses, $\alpha$, and $T$ (we do not need to track whether the $T$
comes from $x$ or $y$ substitutions). Considering only the dominant parametric
dependence as power laws, the pre-factor can be written

$g(r_{\alpha},r_{m},r_{M})=\left(r_{\mu}^{-1/2}r_{\Delta
E}^{1/2}\right)\left(r_{\alpha}^{R_{\alpha}}\,r_{m}^{R_{m}}\,r_{M}^{R_{M}}r_{\Delta
E}^{R_{T}}\right),$ (4.6)

where the terms in the first parenthesis come from the universal thermal
factor, and the terms in the second are from the parametric dependencies of
$\tilde{R}$.

For photochemistry reactions, we have different factors and integrals, but the
overall approach is identical, with the substitution of the radiation
temperature, $T_{\gamma}$, instead of the particle temperature, $T$, and a
slightly different universal thermal factor. Thus, we have

$\displaystyle\gamma_{\rm DM,\gamma}(T_{\gamma})$
$\displaystyle=2c\int_{\Delta E/h}^{\infty}d\nu\frac{\nu^{2}\sigma_{\rm
photo}(\nu)}{e^{h\nu/k_{B}T_{\gamma}}-1}$
$\displaystyle=g(r_{\alpha},r_{m},r_{M})\;\gamma_{\rm
SM,\gamma}\left(\frac{T_{\gamma}}{r_{\Delta E}}\right),$ (4.7)

with

$g(r_{\alpha},r_{m},r_{M})=\left(r_{\Delta
E}^{3}\right)\left(r_{\alpha}^{R_{\alpha}}\,r_{m}^{R_{m}}\,r_{M}^{R_{M}}r_{\Delta
E}^{R_{T}}\right).$ (4.8)

Table 2 summarizes the main results of this section. This table contains “the
most important” molecular hydrogen reactions as designated by Galli & Palla
(1998) (see also Abel et al. (1997)), as well as several additional reactions
needed for accurate results in our companion papers (Gurian et al., 2022; Ryan
et al., 2022). The table shows the parametric dependence of the cross
sections, the re-scaling pre-factor (capturing only the dominant parametric
dependence), and the dominant power-law temperature dependence of the rate as
$b$. Standard Model rates frequently have a temperature dependence that can be
captured to a first approximation by a power law $\gamma_{\rm SM}\propto
T^{b}$, and it is often useful to see this power law in order to understand
the primary parametric dependence of a reaction. However, the full temperature
dependence known for the Standard Model can of course be used (as long as the
re-scaled temperature is used throughout) and is recommended. In the following
sub-sections we go through detailed derivations of a few representative cases.

#### Principle of detailed balance

Several reactions in Table 2 are inverses of each other, and we can exploit
that property to compute some re-scalings using the principal of detailed
balance (Mo et al., 2010). In thermal equilibrium, forward and backward
reaction rates must match, leading to the following relationship between their
cross sections:

$\sigma_{f}\propto\left(\frac{g_{b}}{g_{f}}\right)\,\left(\frac{p_{b}}{p_{f}}\right)^{2}\,\sigma_{b},$
(4.9)

While $p_{b}$ and $p_{f}$ are reaction specific,
recombination/photoionization-type reactions all have $p_{b}=(h\,\nu)/c$ and
$p_{f}=m_{\rm ej}\,v$ and we obtain the Milne Relation (Mo et al., 2010),

$\sigma_{\rm rec}(v,n)=\frac{g_{n}}{g_{n+1}}\left(\frac{h\,\nu}{m_{\rm
ej}\,v\,c}\right)^{2}\sigma_{\rm pi}(\nu,n),$ (4.10)

where $\sigma_{\rm rec}$ is the recombination cross section, $\sigma_{\rm pi}$
is the photonionization cross section, $h\,\nu$ is the photon energy, and
$m_{\rm ej}\,v$ is the momentum of the ejected particle (e.g. an electron in
standard hydrogen recombination; see Section 4.2). Generally, we have
$h\,\nu=\text{K.E.}+\Delta E=k_{B}\,T\left(x^{2}+y^{2}\right)$ through
conservation of energy (with $x$ and $y$ defined in Eqs. (4.4) - (4.5)) and
usually $m_{\rm ej}\approx\mu$, such that

$\sigma_{\rm
rec}\propto\frac{T}{\mu}\frac{\left(x^{2}+y^{2}\right)^{2}}{x^{2}}\sigma_{\rm
pi}$ (4.11)

###  4.2 Reaction 1: Hydrogen recombination

Let us start with one of the simplest reactions, hydrogen recombination:
$p+e\rightarrow{\text{H}}+\gamma$. This reaction is often divided into Case A,
where the hydrogen atoms are sparse enough to recombine directly to the ground
state, and Case B, where the hydrogen atoms are too dense for the Lyman
continuum photons to escape, restricting recombination to excited states. The
dependence on dark-sector parameters has been worked out already in Rosenberg
& Fan (2017) (Case A) and in Hart & Chluba (2017) (Case B), and so this is a
good example to demonstrate the validity of the re-scaling procedure in this
work. The cross section scales in the same way for Case A and Case B, but the
reaction rates of course differ. In Table 2, we have used $b$ values derived
from the Case B rate of Pequignot et al. (1991) in order to have a relevant
comparison to the photoionization rate (in Case A, the appropriate comparison
would be with the collisional ionization rate). In this section, we obtain the
expression for the recombination cross-section using the principal of detailed
balance and the photoionization cross section, and apply this expression to
the Case A rate of Cen (1992).

The photoionization cross section from the ground state of a hydrogenic atom
is (Mo et al., 2010)

$\sigma_{\rm pi}(\nu)=\frac{g_{\rm
bf}}{Z^{2}}\frac{64\pi}{3\sqrt{3}}\frac{h^{2}}{\alpha
m^{2}c^{2}}\left[\frac{E_{\text{H}}}{h\,\nu}\right]^{3}\propto\alpha\,a_{0}^{2}\,\left(\frac{E_{\text{H}}}{h\,\nu}\right)^{3},$
(4.12)

where $E_{\text{H}}=\alpha/(2a_{0})$ is the binding energy of hydrogen,
$g_{\rm bf}$ is the bound-free Gaunt factor, and $Z$ is the nuclear charge
($Z=1$ in the atomic dark matter case). From conservation of energy,
$h\,\nu=\text{K.E.}+\Delta E$, and we have

$\sigma_{\rm pi}\propto\alpha
a_{0}^{2}\frac{\left(\frac{\alpha}{a_{0}}\right)^{3}}{(\text{K.E.}+\Delta
E)^{3}}=m\,\alpha^{5}\frac{1}{(k_{B}\,T(x^{2}+y^{2}))^{3}}.$ (4.13)

Now that we have the photoionization cross section, application of the Milne
relation [Eq. (4.10)] with $\mu=m$ gives

$\sigma_{\rm
rec}\propto\left(\frac{\alpha^{5}}{T^{2}}\right)\,\left(\frac{1}{x^{2}\,(x^{2}+y^{2})}\right).$
(4.14)

Substituting this into Eq. (4.6) and defining
$\tilde{T}_{a}=T/(r_{\alpha}^{2}\,r_{m})$ as the atomic energy analog of
$\tilde{T}_{r}$ and $\tilde{T}_{v}$, we obtain our final re-scaling

$\gamma_{\rm rec,DM}(T)=\left[r_{\alpha}^{2}\,r_{m}^{-2}\right]\;\gamma_{\rm
rec,SM}(\tilde{T}_{a}).$ (4.15)

which matches the results in Rosenberg & Fan (2017) and Hart & Chluba (2017).

A simplified Standard Model recombination rate is given by Cen (1992); Mo et
al. (2010):

$\displaystyle\gamma_{\rm rec,SM}(T)$
$\displaystyle=$8.40\text{\times}{10}^{-11}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$T^{-1/2}\left(\frac{T}{${10}^{3}\text{\,}\mathrm{K}$}\right)^{-0.2}\left(1+\left(\frac{T}{${10}^{6}\text{\,}\mathrm{K}$}\right)^{0.7}\right)^{-1},$
(4.16)

and we find the corresponding dark recombination rate following the procedure
in Section 4.1 as

$\displaystyle\gamma_{\rm rec,DM}(T)$
$\displaystyle=$2.66\text{\times}{10}^{-12}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$\left[r_{\alpha}^{2}\,r_{m}^{-2}\right]\left(\frac{T/(r_{\alpha}^{2}\,r_{m})}{${10}^{3}\text{\,}\mathrm{K}$}\right)^{-0.7}\left(1+\left(\frac{T/(r_{\alpha}^{2}\,r_{m})}{${10}^{6}\text{\,}\mathrm{K}$}\right)^{0.7}\right)^{-1}$
$\displaystyle\approx\begin{dcases}\begin{rcases}$2.66\text{\times}{10}^{-12}$\left[r_{\alpha}^{3.4}\,r_{m}^{-1.3}\right]\left(\frac{T}{${10}^{3}\text{\,}\mathrm{K}$}\right)^{-0.7}&T\ll
r_{\alpha}^{2}\,r_{m}${10}^{6}\text{\,}\mathrm{K}$\\\
$2.11\text{\times}{10}^{-14}$\left[r_{\alpha}^{4.8}\,r_{m}^{-0.6}\right]\left(\frac{T}{${10}^{6}\text{\,}\mathrm{K}$}\right)^{-1.4}&T\gg
r_{\alpha}^{2}\,r_{m}${10}^{6}\text{\,}\mathrm{K}$\end{rcases}\end{dcases}${\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$
(4.17)

In Figure 3, we show that this re-scaling provides a good approximation of the
true dark-recombination rate and that we can improve on these results simply
by using the same re-scaling but a more accurate Standard Model rate with the
form,

$\displaystyle\gamma_{\rm
rec,SM}(T)=\exp\left(\sum_{n=0}^{9}A_{n}\log(T)^{n}\right),$ (4.18)

where the coefficients are given in Abel et al. (1997). We compare the re-
scaled recombination rates, Eq. (4.17) and the re-scaled version of Eq. (4.18)
with the first-principle calculation from Rosenberg & Fan (2017) in Fig. 3.

|
---|---
Figure 3:  Comparisons between our re-scaled recombination rates from Eq.
(4.17) and (4.18) based on Cen (1992) and Abel et al. (1997) and the full rate
given in Rosenberg & Fan (2017). On the left we plot the Standard Model rates,
equivalent to setting $r_{\alpha}=r_{m}=r_{M}=1$. On the right we compute the
percent difference between the Rosenberg and Cen and Rosenberg and Abel
curves, i.e. $\text{\%
Difference}(\rm{Ros},\rm{Cen})=100\%\frac{\left|\gamma_{\rm Ros}-\gamma_{\rm
Cen}\right|}{\gamma_{\rm Ros}}$. We have plotted the rates using
$m=\\{511,40,10\\}$\mathrm{keV}$$, $M=\\{0.938,40,100\\}$\mathrm{GeV}$$, and
$\alpha=\\{$0.007\,29$,$0.01$,$0.01$\\}$, labeled {SM,DM1,DM2}, as a function
of the re-scaled temperature, $\tilde{T}_{a}=r_{\alpha}^{-2}\,r_{m}^{-1}\,T$.
Note that the increased difference between the Rosenberg and Cen curves in
both plots arises from Cen being a simplified analytic fit.

###  4.3 Reaction 3 and 4: Hydrogen electron attachment and $\text{H}^{-}$
photodetachment

The hydrogen electron attachment reaction,
$\text{H}+e^{-}\rightarrow{\text{H}}^{-}+\gamma$, can be effectively treated
as a screened proton capturing a free electron, in a similar process to
standard hydrogen recombination (Ohmura & Ohmura, 1960). The reaction rate can
be computed in analogous fashion, but where the $\text{H}^{-}$ photodetachment
cross section is used instead of the H photoionization cross section. Thus, we
start with the photodetachment cross section, where we use results from
effective range theory (Armstrong, 1963; Ohmura & Ohmura, 1960) and
conservation of energy ($h\,\nu=\Delta E+\text{K.E.}=\Delta E+p_{e}^{2}/(2m)$,
where $\Delta E$ is simply the $\text{H}^{-}$ ground state energy,
proportional to $E_{\rm H}$) such that

$\displaystyle\sigma_{\rm pd}(\nu)$
$\displaystyle=\frac{1}{1-\gamma\varrho}\frac{16\pi}{3}\frac{\alpha\hbar^{2}}{\mu}\frac{\Delta
E^{1/2}(h\nu-\Delta E)^{3/2}}{\left(p_{e}^{2}/2m+\Delta E\right)^{3}}$
$\displaystyle\propto\frac{1}{1-\gamma\varrho}\frac{\alpha}{\mu}\frac{y\,x^{3}}{T\left(x^{2}+y^{2}\right)^{3}}$
(4.19)

with $\gamma^{2}/2=$0.027\,75\text{\,}{\mathrm{a.u.}}^{2}$$ the electron
affinity of hydrogen, and $\varrho=$2.646\text{\,}\mathrm{a.u.}$$ the
effective range (de Jong, 1972) in the Standard Model. The effective range
$\varrho$, having dimensions of length, must be proportional to $a_{0}$. The
electron affinity is defined in Bethe (1949) as $\gamma^{2}/2=(\mu\,\Delta
E)/(2\hbar)\propto m^{2}\alpha^{2}$, so $\gamma\propto 1/a_{0}$. Thus, the re-
scaling in the two terms effectively cancels, and we can ignore the parametric
dependence of the $1/(1-\gamma\varrho)$ term.

We then find the cross-section for hydrogen electron attachment (reaction 3 on
our Table) from the Milne relation [Eq. (4.10)] and $\mu=m$,

$\displaystyle\sigma_{3}(x,y)$
$\displaystyle\propto\frac{\alpha}{\mu^{2}}\frac{x\,y}{\left(x^{2}+y^{2}\right)},$
(4.20)

which re-scales as [Eq. (4.6)]

$\gamma_{3,\rm DM}(T)\approx r_{\alpha}^{2}\,r_{m}^{-2}\,\gamma_{3,\rm
SM}(\tilde{T}_{a}).$ (4.21)

The full Standard Model rate is given as (Galli & Palla, 1998)

$\gamma_{3,\rm
SM}(T)=$1.4\text{\times}{10}^{-18}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-3}\text{\,}{\mathrm{s}}^{-1}$\;T^{0.928}\exp\left(\frac{T}{16200}\right),$
(4.22)

and so the final, re-scaled dark rate is

$\displaystyle\gamma_{3,\rm DM}(T)=$
$\displaystyle$1.4\text{\times}{10}^{-18}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-3}\text{\,}{\mathrm{s}}^{-1}$\,\left[r_{\alpha}^{0.1444}\,r_{m}^{-2.928}\right]$
$\displaystyle\,\times
T^{0.928}\exp\left(\frac{T}{16200\,r_{\alpha}^{2}\,r_{m}}\right).$ (4.23)

Comparing this result with Eq. (4.17), we note that, while both rates use the
same overall pre-factor $g(r_{\alpha},r_{m},r_{M})$, the final re-scalings are
quite different. This divergence highlights the necessity of including the
$r_{\Delta E}=r_{\alpha}^{2}r_{m}$ temperature re-scaling factor.

Conveniently, Eq. (4.8), (4.19) and $\mu=m$ are sufficient to find the re-
scaling for the ${\text{H}}^{-}$ photodetachment rate, giving

$\gamma_{4,{\rm DM}}\approx r_{\alpha}^{5}\,r_{m}\gamma_{4,{\rm
SM}}(\tilde{T}_{a})\,.$ (4.24)

###  4.4 Reaction 7: mutual neutralization

The cross section for the mutual neutralization reaction,
$\text{H}^{-}+p\rightarrow 2\text{H}$, was computed from the Landau-Zener
transition probability at the initial and final potential energy curve psuedo-
crossing (Bates & Lewis, 1955). Although there are two possible final states,
$\text{H}(1s)\text{H}(2s\text{ or }p)$ and $\text{H}(1s)\text{H}(3s,p,\text{or
}d)$, for our purposes the resulting cross section is effectively

$\sigma_{7}\propto
R^{2}\left(\frac{k_{f}^{2}}{k_{i}^{2}}\right)\,F_{3}(\zeta),$ (4.25)

where $R\propto a_{0}$ is the internuclear distance at the potential crossing
point, $k_{i}$ and $k_{f}$ are the momentum of the initial state and final
states, with $k_{i}^{2}\propto\mu\,\text{K.E.}$ and $k_{f}^{2}\propto
k_{i}^{2}+\mu\Delta E=\mu\left(\text{K.E.}+\Delta E\right)$,
$F_{3}(\zeta)=\zeta^{n-1}\Gamma(n-1,\zeta)-(2\zeta)^{n-1}\Gamma(n-1,2\zeta)$
and $\Gamma(n-1,\zeta)$ is the incomplete gamma function. The term $\zeta$ is
a dimensionless ratio of momenta, which, after some simplification, can be
written as $\zeta\propto k_{f}^{-1}(\alpha\mu\Delta U)/\Delta E^{2}$, and
$\Delta U$ and $\Delta E$ are the zeroth order and true potential energy curve
differences. In our parameter regime of interest, $\zeta\ll 1$ and a series
expansion gives $F_{3}(\zeta)\approx\zeta$. Then we have

$\displaystyle\sigma_{7}$ $\displaystyle\propto
a_{0}^{2}\left(\frac{\text{K.E.}+\Delta
E}{\text{K.E.}}\right)\left(\frac{\alpha\,\sqrt{\mu}\,\Delta U}{\Delta
E^{2}\sqrt{\text{K.E.}+\Delta E}}\right)$ $\displaystyle\propto\alpha
a_{0}^{2}\sqrt{\mu}\frac{\sqrt{\text{K.E.}+\Delta E}}{\text{K.E.}\,\Delta E}$
$\displaystyle\propto\frac{\alpha
a_{0}^{2}\sqrt{\mu}}{T^{3/2}}\frac{\sqrt{x^{2}+y^{2}}}{x^{2}y^{2}}$ (4.26)

where, in the second line, we used the fact that $\Delta U$ and $\Delta E$ are
proportional. Then, since $\Delta E$ is proportional to the atomic binding
energy and $\mu$ here is the proton mass, we obtain a net re-scaling of

$\gamma_{7,\rm DM}\approx r_{\alpha}^{-3}\,r_{m}^{-3}\gamma_{7,\rm
SM}(\tilde{T}_{a})\,.$ (4.27)

###  4.5 Reactions 10 and 15: Langevin and detailed balance

An interaction without a short-range potential barrier will always occur as
long as the incoming particle passes the centrifugal potential barrier and
spirals into the target. These reactions can be well-treated by classical
capture models, the theory of which is presented in e.g. Hirasawa (1969);
Zhang & Willitsch (2018). In this approach, the long-range potential which
governs the reaction can be written as a series in powers of the separation
$R$:

$V(R)=-\sum_{n}\frac{C_{n}}{R^{n}}.$ (4.28)

The leading order term in this series dominates in general, and in the
Langevin (charge-induced dipole) case, $n=4$, $C_{4}=pe^{2}/2$ where $p$ is
the isotropic polarizability of the molecule, proportional to the trace of
$p_{ij}$. The reaction occurs if the collision energy exceeds the maximum of
$V_{eff}(R)=V(R)+V_{\rm centrifugal}(R)$, where we include a centrifugal term,
$V_{\rm centrifugal}(R)$, to account for a non-zero impact parameter. Then,
since the impact parameter is proportional to the angular momentum, we obtain
the cross section in the Langevin case:

$\displaystyle\sigma_{\rm
Lang}(x)\propto\sqrt{\frac{pe^{2}}{x^{2}k_{B}\,T}}\propto\sqrt{\frac{a_{0}^{3}\,\alpha}{T}}\frac{1}{x},$
(4.29)

where we have used the property that $p\propto a_{0}^{3}$. Then using our
tools from Section 4.1 and the definition $\mu\approx M$, we obtain an
approximate rate re-scaling of

$\gamma_{\rm Lang,DM}(T)\approx
r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\gamma_{\rm
Lang,SM}(\tilde{T}_{a}),$ (4.30)

which we use for both of the Langevin-type reactions: associative detachment
of H and $\text{H}^{-}$ ($\text{H}^{-}+\text{H}\rightarrow\text{H}_{2}+e$) and
$\text{H}_{2}$ formation due to $\text{H}_{2}^{+}$
($\text{H}_{2}^{+}+\text{H}\rightarrow\text{H}_{2}+p$).

The computation for the inverse reaction, $\text{H}_{2}^{+}$ charge exchange
($\text{H}_{2}+p\rightarrow\text{H}_{2}^{+}+\text{H}$), cross section uses the
more general detailed balance results combined with the Langevin cross
section. In this case, we have
$p_{10}=\mu_{10}v_{10}\approx\mu_{15}v_{15}=p_{15}$ (where the labels $10$ and
$15$ correspond to those reactions), since
$\mu_{10}=(m_{\text{H}_{2}^{+}}\,m_{\text{H}})/(m_{\text{H}_{2}^{+}}+m_{\text{H}})\approx(2\,M^{2})/(3\,M)\approx(m_{\text{H}_{2}}M)/(m_{\text{H}_{2}}+M)\approx\mu_{15}$.
With the standard assumption that the internal degrees of freedom $g_{10}$ and
$g_{15}$ stay constant, we have $\sigma_{15}\propto\sigma_{10}$, and we only
need the temperature re-scaling, which, from reaction 10, is proportional to
the hydrogen binding energy $E_{\rm H}$.

###  4.6 $\text{H}_{2}$ Dissociation Scaling

We approximate the $\text{H}_{2}$ dissociation reaction,
$\text{H}_{2}+\text{H}\rightarrow 3\text{H}$, at the lowest level as a hard-
sphere collision. Then, since the hard-sphere impact parameter, $b$, scales as
$a_{0}\propto(\alpha m)^{-1}$, we have

$\sigma_{\rm diss}\propto b^{2}\propto\left(\frac{1}{\alpha m}\right)^{2}.$
(4.31)

The dominant binding energy is simply the dissociation energy, $D_{e}$ which
depends on $\alpha^{2}m$, so $r_{D_{e}}=r_{\alpha}^{2}\,r_{m}$. The overall
re-scaled rate is then given by

$\displaystyle\gamma_{\rm diss,DM}(T)$
$\displaystyle\approx\left(\sqrt{\frac{r_{D_{e}}}{r_{M}}}\right)\left(r_{\alpha}^{-2}\,r_{m}^{-2}\right)\;\gamma_{\rm
diss}(T/r_{D_{e}})$ (4.32)
$\displaystyle=r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\;\gamma_{\rm
diss,SM}(\tilde{T}_{a}),$ (4.33)

where the first term in brackets comes from the velocity average (with the
approximation $\mu=m_{\text{H}_{2}}=M$) and the second from the cross section.

Note that there is a major limitation to this re-scaled rate: the hard-sphere
approximation introduces significant simplifications into the cross section
parametric dependence that are not necessarily justified. As Martin et al.
(1996) points out, the $\text{H}_{2}+\text{H}\rightarrow 3\text{H}$ reaction
rate requires accounting for all of the possible energy level transitions and
populations, including not just collisions but also spontaneous dissociation
via tunneling from quasi-bound states. This introduces myriad locations for
the various parameters to creep in, including, at a minimum, the resulting
temperature dependence in the cross section.

###  4.7 Three-Body-Reactions

We re-scale the reaction rate for the following three-body reactions, taken
from Yoshida et al. (2006):

$\displaystyle 3\text{H}$ $\displaystyle\rightarrow\text{H}_{2}+\text{H}$
(4.34) $\displaystyle\text{H}_{2}+2\text{H}$ $\displaystyle\rightarrow
2\text{H}_{2}$ (4.35) $\displaystyle 2\text{H}+\text{H}^{+}$
$\displaystyle\rightarrow\text{H}_{2}+\text{H}^{+}$ (4.36) $\displaystyle
2\text{H}+\text{H}^{+}$ $\displaystyle\rightarrow\text{H}_{2}^{+}+\text{H}.$
(4.37)

The re-scaling is common to all four reactions, and can be derived from the
principle of detailed balance (see, e.g. Flower (2007)). The chemical
equilibrium relation gives

$\gamma_{a}n_{\text{H}}^{2}=\gamma_{d}n_{\text{H}_{2}},$ (4.38)

where $\gamma_{a}$ is the rate coefficient for three-body association reaction
and $\gamma_{d}$ is the rate coefficient for the inverse process. Then, using
the Saha equation:

$n(\text{H}_{2})=n^{2}(\text{H})\frac{Z_{\text{H}_{2}}}{Z_{\text{H}}^{2}}\left(\frac{1}{\pi
MT}\right)^{3/2}e^{D_{e,D}/T},$ (4.39)

where $Z=\sum_{i}g_{i}e^{-E_{i}/T}$ is the partition function with $g_{i}$ the
statistical weight.

Substituting, and noting that the partition function for atomic hydrogen
(assumed predominantly in the ground state) is a constant (whose value is
given as $1$ in Flower (2007) and $2$ in Forrey (2013)):

$\gamma_{a,\text{D}}(T)=\gamma_{d,\text{D}}(T)\frac{2}{Z_{\text{H}_{2}}}\left(\frac{1}{\pi
MT}\right)^{3/2}e^{D_{\rm e}/T}.$ (4.40)

Flower (2007) gave $Z_{\text{H}_{2}}\approx.028T$, and at the low temperatures
where these reactions are significant, this scales as the rotational energy,
$T\rightarrow\tilde{T}_{r}$. We substitute
$z_{\text{H}_{2}}(\tilde{T}_{r})=r_{m}\,r_{M}^{-1}z_{H_{2}}(T_{a})$ and apply
the scaling of $\gamma_{\rm diss,DM}$ from Section 4.6 to substitute
$\gamma_{\rm
diss,DM}(T)=r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\gamma_{\rm
diss,SM}(\tilde{T}_{a})$, yielding

$\displaystyle\gamma_{a,\text{D}}(T)$
$\displaystyle=r_{\alpha}^{-4}\,r_{m}^{-4}\,r_{M}^{-1}$
$\displaystyle\times\frac{2}{Z_{\text{H}_{2}}(\tilde{T}_{a})}\left(\frac{1}{\pi
m_{\text{H}}\tilde{T}_{a}}\right)^{3/2}e^{D_{e}/\tilde{T}_{a}}\gamma_{\rm
diss,SM}(\tilde{T}_{a})$ (4.41)
$\displaystyle=r_{\alpha}^{-4}\,r_{m}^{-4}\,r_{M}^{-1}\gamma_{a,\rm
SM}(\tilde{T}_{a}).$ (4.42)

The reaction rates are taken from Yoshida et al. (2006) and Krstic et al.
(2003). This derivation can be applied to the final reaction of the four
(which involves a charge exchange) by considering the Saha factors between
$\text{H}_{2}^{+}$ and $\text{H}_{2}$ as well as between H and $\text{H}^{+}$.

##  5 Validity of results

There are several limitations to the results presented above, some of which
are based on assumptions we have made, whereas others follow from the
dependence on Standard Model chemistry. Therefore, many of the results so far
presented are only valid in certain regions of the parameter space and we
discuss the reasons and regions here.

###  5.1 Limitations inherited from Standard Model results

The largest source of inaccuracy of our results follows from our dependence on
the Standard Model rates. Many of the reactions, especially those considered
“unimportant" (in Table 3) are poorly known, while others are only known for
specific cases. For example, the reaction $\text{H}^{-}+p\rightarrow
2\text{H}$ was only known to within an order of magnitude until recently
(Glover, 2015b), while higher vibrational levels in the
$\text{H}_{2}+e\rightarrow\text{H}+\text{H}^{-}$ reaction were not thought to
have a significant contribution to the net rate until Capitelli et al. (2007).

Most of the reactions listed in Table 2 are obtained from older cross sections
that may not account for these corrections and special cases found more
recently. However, we expect that as long as the approximations used in the
older rates remain reflective of the underlying physics, the rate re-scalings
obtained will remain valid, as demonstrated in Figures 1 and 3.

Additionally, the majority of reaction cross sections listed in section 4 are
analytic fits to numerical computations or even experimental data. As such,
they are only defined for certain temperature regimes. For example, the
reaction $\text{H}+p\rightarrow\text{H}_{2}^{+}+\gamma$ has a reaction rate
defined on the interval $$1\text{\,}\mathrm{K}$\leq
T\leq$32\,000\text{\,}\mathrm{K}$$, but is undefined elsewhere. Extrapolation
past those temperature limits can result in significant divergences,
especially with some of the more accurate fits that have strong exponential
terms, like the rovibrational cooling rates mentioned in Section 3. These
temperature limits will scale with the dark parameters but extra care should
still be taken when exploring more extreme parameter space ranges.

A more subtle issue exists for radiative reactions because the gas and photons
need not share a temperature. In principle, the reaction rates should then
depend on both the gas and photon temperature, but all of our rates depend on
only one of these temperatures. Cyr-Racine & Sigurdson (2013) points out that
in particular the standard calculation of the recombination coefficient
neglects stimulated recombination, which explicitly depends on the photon
temperature. That work found a change of at most a factor of a few in the
recombination coefficient due to this effect when $T_{g}=0.01T_{\gamma}$. In
the particular case of hydrogen recombination, this issue is relevant only for
weakly coupled dark matter models, since otherwise the gas and radiation are
thermally coupled at recombination. For the problem of halo formation (Grassi
et al., 2014), the matter and radiation do not in general share a temperature,
but even so standard calculations consider each rate to be a function of only
one temperature.

A final limitation derives from several of the rates listed here depending on
the assumption from the Standard Model that $m\ll M$. Trivially, in many cases
we approximate the reduced mass $\mu=(m\,M)/(m+M)\approx m$ for the $e$, $p$
system, or $\mu\approx m_{\text{H}}\approx M$ for the $H$, $H$ (or $p$)
system, which can be easily adjusted as $m\rightarrow M$. On the other hand,
in section 2 we assumed stationary nuclei, with electronic motion resulting in
minor corrections at most. This no longer holds as $m\rightarrow M$, as we can
no longer use the Born-Oppenheimer approximation. This affects both the
fundamental wave functions and cross section calculations. For example, the
derivation for the hydrogen electron attachment reaction
($\text{H}+e\rightarrow\text{H}^{-}+\gamma$) assumed a massive stationary,
screened proton captures a free electron. In the extreme case, $m\approx M$,
the physics should become much closer to positronium production,
$e^{-}+e^{+}\rightarrow Ps$ or radiative association of H and $p$,
$\text{H}+p\rightarrow\text{H}_{2}^{+}+\gamma$, with the corresponding re-
scaling, but where and how the transition would occur is beyond the scope of
this work. Even the basic hydrogen atomic physics are modified; with a
significantly larger electron mass, the energy eigenvalues would have a
significant dependence on the nuclear (proton) mass (Poszwa et al., 2016).

###  5.2 Comparison with deuterated chemistry

As mentioned in the introduction, we can compare our rates with those obtained
for deuterium reactions. Gay et al. (2011) has compiled a significant list of
deuterium reaction rates; while the majority are computed using the mass re-
scaling described in Appendix A, or assumed identical to the hydrogen rate
(e.g. $\text{D}^{+}+e\rightleftharpoons\text{D}+\gamma$), one is based
directly on experimental data. The reaction
$\text{D}_{2}+\text{D}^{+}\rightarrow\text{D}_{2}^{+}+\text{D}$ has a rate
derived from the cross section listed in Wang & Stancil (2002) and given by

$\displaystyle\gamma_{D_{2}+D^{+}\rightarrow D_{2}^{+}+D}$ $\displaystyle=$
$\displaystyle\,$2.46\text{\times}{10}^{-9}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$\,T_{3}^{-0.439}\exp\left(\frac{-{30500}}{T}\right).$
(5.1)

While this appears significantly different to our re-scaled result from Galli
& Palla (1998), $\gamma_{15,\rm
DM}\approx$2.12\text{\times}{10}^{-10}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$\,\exp\left(-21050/T\right)$,
the discrepancy is not fatal to our argument. Savin et al. (2004) points out
that there is significant variation between the reported rates for the
$\text{H}_{2}+\text{H}^{+}\rightarrow\text{H}_{2}^{+}+\text{H}$ reaction, with
sometimes drastically different values and temperature behavior even between
rates based on the same source data. While the rates may differ significantly,
they share common features, including an exponential drop off dependent on the
threshold energy at low temperatures, and minimal temperature dependence in
the ${10}^{4}\text{\,}\mathrm{K}e4.5\text{\,}\mathrm{K}$ range. In Figure 4,
we plot several re-scaled hydrogen rates (from Galli & Palla (1998); Abel et
al. (1997); Savin et al. (2004)), as well as the rate from Gay, to demonstrate
that the rates are sufficiently similar, especially at higher temperatures.

Figure 4: $\text{D}_{2}+\text{D}^{+}\rightarrow\text{D}_{2}^{+}+\text{D}$
reaction rate and re-scaled H rates for comparison. Re-scaled
$\text{H}_{2}+\text{H}^{+}\rightarrow\text{H}_{2}^{+}+\text{H}$ reaction rates
are from Galli & Palla (1998); Abel et al. (1997); Savin et al. (2004) and the
deuterium rate is from Gay et al. (2011). The deuterium rate is clearly
comparable to the re-scaled hydrogen rates to the extent that they are
comparable to each other.

As stated, the remaining pure deuterium reactions in Gay et al. (2011) not
identical to hydrogen rates have been computed using the mass re-scaling
method described in section A, with a minor change: in Eq. (A2), $a_{2,\rm
known}\rightarrow a_{2,\rm known}+0.5$. To the extent that those rates have
been used in many astrochemistry calculations, this further validates our
approach of re-scaling known Standard Model rates, and comparison of Eq. (A3)
with Eq. (4.2) and Eq. (4.6) reveals that the simplified version is equivalent
to $g(r_{\alpha},r_{m},r_{M})=r_{\mu}^{-1/2-a/2}$ and $\Delta E$ at least
approximately independent of $\mu$. Since $\mu\propto M$ for these cases, this
holds reasonably well for those reactions with simple cross sections and
atomic binding energy, for example, mutual neutralization,
$\text{H}^{-}+p\rightarrow 2\text{H}$. Our more complete re-scaling procedure
($\gamma_{7,\rm DM}\propto r_{M}^{-0.5}\gamma_{7,\rm SM}$) gives equivalent
results to the mass re-scaling ($\gamma_{7,\rm DM}\propto
r_{M}^{-0.487}\gamma_{7,\rm SM}$).

Lastly, we consider the cooling rates and demonstrate why we do not validate
using HD chemistry. Unfortunately, we were unable to find a $\text{D}_{2}$
cooling function based on experimental data or theoretical calculations and
can only compare our re-scaled rates to HD cooling rates, setting
$M=1.5m_{p}$. However, HD has significantly different behavior, as compared to
$\text{H}_{2}$, due to its mass asymmetry-induced dipole moment and allowance
of $J\pm 1$ transitions. Thus, the cooling comparison breaks down in the low
temperature regime where those aspects become relevant. We demonstrate this
issue in Figure 5, where we plot the re-scaled Hollenbach & McKee (1979) and
Glover & Abel (2008) $\text{H}_{2}$ cooling rates versus the HD cooling rate
from Glover & Abel (2008). Note that while the differences between the HD
curve and the re-scaled curves are small in the high temperature region, at
low temperatures the HD curve is orders of magnitude higher.

Figure 5: Demonstration of the inability to utilize the HD molecule in
validation. While the re-scaled Hollenbach & McKee (1979) and Glover & Abel
(2008) cooling curves are similar to the HD cooling curve in the high
temperature region, at low temperatures the allowed $J\pm 1$ transitions
enable HD to have an orders of magnitude higher cooling rate.

###  5.3 Additional reactions

The minimal reaction list of Galli & Palla (1998) informs our reaction list.
However, we neglect several of these reactions (briefly listed in Table 3).
This list includes three-body reactions and reactions deemed irrelevant by
Galli & Palla (1998). Since different reactions are re-scaled by different
amounts, one might worry that reactions which are insignificant in the
Standard Model could contribute. There are three distinct reasons a reaction
might fail to contribute:

  1. 1.

It may be out-competed by a different channel involving identical reactants,
for example $\gamma_{6}$ compared to $\gamma_{7}$.

  2. 2.

It may be suppressed below some threshold energy which is large compared to
competing reactions. For example reaction 18 requires a photon with enough
energy to excite the electronic degrees of freedom of the $\rm H_{2}$
molecule, but except in the case of a cold gas exposed to an external Lyman-
Werner flux (which we never encounter) collisional dissociation mechanisms
will efficiently destroy the molecules at a much lower temperature.

  3. 3.

It may be suppressed by the relative abundances of the species involved. For
example, reaction 10 dominates over reaction 13 because in general there is
more neutral hydrogen than $\text{H}_{2}$.

For the first case, since the basic energy and angular momentum structure of
the molecule is preserved under re-scaling, whatever physical principle (i.e.
selection rules, charge-dipole interactions) prefers one channel over another
continues to operate. In the second case, the sharp cutoff in the reaction
rate below the threshold energy cannot be overcome by parametric re-scaling of
the rate coefficient. In the third case, serious care is required because
changes to the abundances or the reaction rates can alter the relative
importance of the various reactions. Note also that these categories are not
completely disjoint. In the example mentioned for (2), if one encounters a
case where $\text{H}_{2}$ molecules and energetic photons coexist, then the
considerations of case (3) might still prevent the reaction from contributing.

We have re-scaled and included all reactions omitted due to (3) except for
reaction 11, $\text{H}_{2}^{+}+e\rightarrow 2\text{H}$. This reaction is
outcompeted by reaction 10,
$\text{H}_{2}^{+}+\text{H}\rightarrow\text{H}_{2}+p$, in all cases considered
because we always find $x_{\text{H}}\gg x_{e}$ by the time $\text{H}_{2}^{+}$
forms.

# | Reaction | Justification
---|---|---
6 | $\text{H}^{-}+p\rightarrow\text{H}_{2}^{+}+e$ | 1
11 | $\text{H}_{2}^{+}+e\rightarrow 2\text{H}$ | 3
12 | $\text{H}_{2}^{+}+\gamma\rightarrow 2p+e$ | 1
14 | $\text{H}_{2}^{+}+\text{H}\rightarrow\text{H}_{3}^{+}+\gamma$ | 1
16 | $\text{H}_{2}+e\rightarrow\text{H}+\text{H}^{-}$ | 2
17 | $\text{H}_{2}+e\rightarrow 2\text{H}+e$ | 1,2
18 | $\text{H}_{2}+\gamma\rightarrow\text{H}_{2}^{+}+e$ | 2
19 | $\text{H}_{3}^{+}+\text{H}\rightarrow\text{H}_{2}^{+}+\text{H}_{2}$ | 2
21 | $\text{H}_{2}+\text{H}^{+}\rightarrow\text{H}_{3}^{+}+\gamma$ | 1
22 | $\text{H}_{3}^{+}+\gamma\rightarrow\text{H}_{2}^{+}+\text{H}$ | 2

Table 3: Reactions from Galli and Palla that are not part of the minimal model
and were not included in the re-scaling procedure with the justification for
omission, as enumerated in the text.

##  6 Conclusion

In order to extend the study of atomic dark matter models and determine their
impact on the formation and evolution of large-scale structure, galactic
halos, and compact object, we have extended previous studies (Rosenberg & Fan,
2017; Buckley & DiFranzo, 2018; Boddy et al., 2016) by calculating the
chemistry of dark molecular hydrogen. As dark molecular hydrogen processes
dominate the chemistry at low temperature, this work provides an essential
tool for studying the formation and evolution of low-mass halos with much
lower virial temperatures, and the end point of cooling processes within
halos.

We have provided some basic tools for working with dark molecules, including
the re-scaling for several key atomic and molecular quantities (Table 1), and
the procedure for obtaining approximately re-scaled rates when the physics of
the analogous Standard Model interaction is known (e.g. Eq. (2.25), Eq. (4.2),
and Eq. (4.6)). We have applied these tools to two important areas in
molecular chemistry, thermal processes and chemical reactions, and determined
how the rates in a minimal subset of these processes can be re-scaled to
obtain expressions valid for a wide range of dark matter parameters (Eq.
(3.23), Eq. (3.34), Eq. (3.36) and Table 2). Importantly, these re-scaled
rates are only as accurate as the original Standard Model rates and neither
correct nor transcend any limitations or assumptions therein, as demonstrated
with our comparison with deuterated reactions (Eq. (5.1) and Figures 4 & 5).

Additionally, we note that many astrochemical databases and chemical models
either neglect deuterium reactions entirely, or when present either use the
hydrogen rate directly or use the simplified mass re-scaling method from
Appendix A (e.g. the deuspin chemical network in the KIDA database, from
Majumdar et al. (2016)). Therefore, at a minimum we expect the re-scaling
method presented here to provide improved rates for deuterium reactions that
can easily be included in said models.

We have provided the first calculations of dark molecular processes with
sufficient accuracy to be used in simulations. These are used in our companion
papers to derive cosmic abundances of various states (ionized, atomic, and
molecular hydrogen) (Gurian et al., 2022) and to model the low-virial
temperature halo formation in analogous fashion to the minihalo formation in
Standard Model (Ryan et al., 2022). Eventually, these rates and the re-scaling
procedure can be used to validate phenomenological models used in large-scale
structure simulations for studying the formation and evolution of dark-matter
halos, which increasingly depend on full numerical calculations that depend on
chemistry rates that we study here.

At the core of the atomic-dark-matter halo, the dark-matter density may reach
high enough values so that the re-scaling of other process rates, especially
many-body processes such as $\text{H}_{3}$ reactions, might also be necessary.
Of course, to get the dark-chemistry rates, one can also take the parallel
approach of computing the interaction rates from the first-principle quantum
calculations. We leave these computations for future work.

Funding for this work was provided by the Charles E. Kaufman Foundation of the
Pittsburgh Foundation. 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.

##  Appendix A Brief Review of Mass Re-scaling

Applications of chemistry in astrophysics frequently require chemical reaction
rates that are difficult to compute from first-principles, and are required in
temperature and density regimes inaccessible in laboratory experiments on
Earth. Sometimes, rates may be known for other isotopes or related species and
in these cases a common approach in the literature (see e.g. Stancil et al.
(1998); Walker et al. (2014)) is similar to what we present in this work: re-
scaling the known rate by a ratio of the masses.

Assuming the known (endoergic) reaction rate can be approximated by the form

$\gamma_{\rm known}=a_{1}\;T^{a_{2}}\;\exp\left(\frac{a_{3}}{T}\right),$ (A1)

with temperature $T$ and constants $a_{i}$, then in the mass-scaled rate the
coefficient $a_{1}$ takes the form

$a_{1,\rm new}=a_{1,\rm known}\left(\frac{\mu_{\rm known}}{\mu_{\rm
new}}\right)^{a_{2,\rm known}},$ (A2)

where $\mu_{\rm known}$ ($\mu_{\rm new}$) is the reduced mass of the known
(new) reaction. All other terms remain the same and the same formula holds for
exoergic rates. For example, the rate coefficient for deuterium mutual
neutralization, $\text{D}^{-}+\text{D}^{+}\rightarrow 2\text{D}$ can be
approximately obtained by re-scaling the rate for the regular hydrogen
reaction, $\text{H}^{-}+p\rightarrow 2\text{H}$. The hydrogen reaction is
known to be of the form of Eq.(A1) with
$a_{1}=$1.4\text{\times}{10}^{-7}\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$
and $a_{2}=$-0.487$$ (Stancil et al., 1998). Using
$\mu_{\text{D}}\approx$1\text{\,}\mathrm{GeV}$$ and
$\mu_{\text{H}}\approx$0.5\text{\,}\mathrm{GeV}$$, the mass re-scaling
procedure suggests the corresponding deuterium process is proportional to
$a_{1,\text{D}}=a_{1,\text{H}}\left(\mu_{\text{H}}/\mu_{\text{D}}\right)^{a_{2,\text{H}}}\approx$1.4$\,a_{1,\text{H}}=$1.96\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$.

The derivation of this re-scaling relationship is straightforward (Walker et
al., 2014). Assume the reaction cross section depends on the collision
velocity $v$ as

$\sigma(v)=Bv^{a},$ (A3)

where $B$ is a mass independent constant with units $[{\rm length}]^{2-a}[{\rm
time}]^{a}$. The rate is given by the thermal average,

$\gamma=\langle\sigma\,v\rangle=\left(\frac{2}{\pi}\right)^{1/2}\left(\frac{\mu}{k_{B}T}\right)^{3/2}\int_{0}^{\infty}\sigma(v)\exp(-\mu
v^{2}/2k_{B}T)v^{3}dv\,.$ (A4)

The integral can be performed exactly, giving

$\gamma=\frac{B}{\sqrt{\pi}}2^{\frac{a+3}{2}}\,\left(\frac{k_{B}T}{\mu}\right)^{\frac{a+1}{2}}\Gamma\left(\frac{a}{2}+2\right)=A(a)\left(\frac{T}{\mu}\right)^{b}.$
(A5)

Then re-scaling $\mu$ gives back the expression found in Eq. (A2), where
$a_{1}=A(a)\mu^{-b}$ and $a_{2}=b=(a+1)/2$. Note that the exponential cutoff
in Eq. (A1) not derived in Walker et al. (2014) can arise naturally within
this formalism in at least two ways. First, by imposing a lower energy cutoff
in the integration:

$\displaystyle\gamma$
$\displaystyle=\left(\frac{2}{\pi}\right)^{1/2}\left(\frac{\mu}{k_{B}T}\right)^{3/2}\int_{v_{min}}^{\infty}\sigma(v)\exp(-\mu
v^{2}/2k_{B}T)v^{3}dv$ (A6)
$\displaystyle=\frac{B}{\sqrt{\pi}}2^{\frac{a+3}{2}}\left(\frac{k_{B}T}{\mu}\right)^{\frac{1+a}{2}}\Gamma\left(2+\frac{a}{2},\frac{\mu
v^{2}_{min}}{2k_{b}T}\right),$ (A7)

and approximating the incomplete gamma function in $1/T$ by an exponent in
$-1/T$ (which is exact for $a=-2$). Second, if there is no minimum energy
cutoff in the forward reaction (i.e. Langevin reactions), an exponential term
will appear in the reverse reaction due to the detailed balance factor. This
can give us some confidence that the form Eq. (A2) is physically well
motivated and could be reasonable to re-scale empirical fits of this form over
a wide temperature range.

Unfortunately, Eq. (A2) is insufficient for computing dark reaction rates for
two reasons. First, it clearly cannot capture the effect of allowing $\alpha$
and $m$ to vary from their Standard Model values. Second, even in the case
where $r_{\alpha}=r_{m}=1$, the $B$ term in Eq. (A3) is dimensionful, and
often contains factors of $\mu$ that are not accounted for. Likewise, the
dimensionful $a_{2}$ and $a_{3}$ terms in Eq. (A1) can depend on $m$, $M$, and
$\alpha$, significantly changing the rate’s temperature dependence. The re-
scaling procedure we have presented in this paper overcomes these limitations.

##  Appendix B Quadrupolar spontaneous emission in Atomic Units

The quadrupolar spontaneous emission rate is usually given in atomic units in
the literature. In these units, $e^{2}=4\pi\epsilon_{0}=\hbar\equiv 1$, masses
are measured in $m_{e}\equiv 1$, and energy is measured in units of the
Hartree energy, $E_{\rm H}=\alpha^{2}m_{e}c^{2}=1$. With these definitions,
$c=1/\alpha$, and the unit of time is $\hbar/E_{\rm H}$.

When the energies and quadrupole moments are given in atomic units, the rate
for a standard model hydrogen atom to transition from a state with
(vibrational, rotational) quantum numbers ($\nu^{\prime},J+2$) to the state
with ($\nu,J$) is (Turner et al., 1977; Flower, 2007)

$\displaystyle A(\nu,J\leftarrow\nu^{\prime},J+2)=$
$\displaystyle\left($1.43\text{\times}{10}^{4}\text{\,}{\mathrm{s}}^{-1}$\right)(E_{\nu^{\prime},J+2}-E_{\nu,J})^{5}\frac{3(J+2)(J+1)}{2(2J+5)(2J+3)}Q^{2}_{\nu^{\prime},J+2;\nu,J}$
$\displaystyle=$
$\displaystyle\frac{1}{60c^{5}}(E_{\nu^{\prime},J+2}-E_{\nu,J})^{5}\frac{3(J+2)(J+1)}{2(2J+5)(2J+3)}Q^{2}_{\nu^{\prime},J+2;\nu,J}$
(B1)

The rates for other transitions, with $\Delta J=0,-2$, have the same form. The
numerical pre-factor in the top line can be obtained from the second line
using the conversion to S.I. units from the unit of time in atomic units,
$\hbar/E_{\rm H}$, ($\hbar/E_{\rm
H}=$2.418\,884\text{\times}{10}^{-17}\text{\,}\mathrm{s}$$ in S.I. units),
$c=137$, and the factor of 60 (note the typo in Flower (2007)).

From Eq.(B), the re-scaling relationships needed to approximate dark matter
transition rates can be determined as a change in units from Standard Model
atomic units to dark matter atomic units. Given a numerical value of some
transition rate in the Standard Model, the parametric dependence comes from
changing units in the factor of $E_{\rm H}/c^{5}=E_{\rm H}\alpha^{5}$ in the
pre-factor in Eq. (B), and from the fundamental constants appearing in the
energy difference of vibrational or rotational modes when those energies are
expressed in atomic units. Since the quadrupole moment is proportional to
$a_{0}^{2}$, there is no additional parametric dependence from $Q^{2}$. That
is, suppose a Standard Model quadrupole in S.I. units is
$Q\equiv\tilde{Q}a_{0}^{2}$, where $\tilde{Q}$ is value of the quadrupole in
Standard Model atomic units. Then the dark matter quadrupole can be obtained
by $Q^{\rm DM}=Q\left(a_{0,{\rm DM}}/a_{0}\right)^{2}$. Using dark atomic
units on the left hand side, $\tilde{Q}^{\rm DM}=Q\left(1/a_{0}\right)^{2}$,
or $\tilde{Q}^{\rm DM}=\tilde{Q}$.

Finally, then, the dark matter Einstein $A$ coefficient for transitions
between rotational levels can be estimated by the re-scaling

$\displaystyle A_{\rm quad,rot\,DM.}=$
$\displaystyle\left[\left(\frac{\alpha_{\rm DM}}{\alpha_{\rm
SM}}\right)^{5}\frac{E_{\rm h,\,DM}}{E_{\rm
h\,SM}}\right]\left[\left(\frac{\Delta E_{\rm rot,DM}}{\alpha_{\rm
DM}^{2}mc^{2}}\right)\left(\frac{\alpha_{\rm SM}^{2}m_{e}c^{2}}{\Delta E_{\rm
rot,SM}}\right)\right]^{5}\;A_{\rm quad,rot\;SM}$ (B2) $\displaystyle=$
$\displaystyle\left[\frac{r_{\alpha}^{7}r_{m}^{6}}{r_{M}^{5}}\right]\;A_{\rm
quad,rot,\;SM}$

and, for vibrational transitions,

$\displaystyle A_{\rm quad,vib\,DM.}=$
$\displaystyle\left[\left(\frac{\alpha_{\rm DM}}{\alpha_{\rm
SM}}\right)^{5}\frac{E_{\rm h,\,DM}}{E_{\rm h\,SM}}\right]\left[\frac{\Delta
E_{\rm vib,DM}}{\alpha_{\rm DM}^{2}mc^{2}}\frac{\alpha_{\rm
SM}^{2}m_{e}c^{2}}{\Delta E_{\rm vib,SM}}\right]^{5}\;A_{\rm quad,vib\;SM}$
$\displaystyle=$
$\displaystyle\left[\frac{r_{\alpha}^{7}r_{m}^{7/2}}{r_{M}^{5/2}}\right]\;A_{\rm
quad,vib\;SM}\,.$ (B3)

These expressions agree with those obtained in the main text.

##  Appendix C The Dark Ortho- to Para-Hydrogen Ratio

In low-temperature, low-density hydrogen gas, the ratio of ortho- (nuclear
spin $I=1$, corresponding to rotational level $J=1$) to para- ($I,J=0$) forms
of molecular hydrogen plays a key role in rovibrational cooling, due to the
differing energies of the lowest transition (Glover & Abel, 2008). In the
Standard model, this ratio is set by the collisional processes
$\text{H}_{2}+p$ and $\text{H}_{2}+\text{H}$, with the proton collision
dominating in primordial halos, and is generally assumed to be 3 to 1 (Glover
& Abel, 2008; Lique et al., 2012). Radiative transitions do not contribute, as
their rate of $\Gamma\approx${10}^{-21}\text{\,}{\mathrm{s}}^{-1}$$(Pachucki &
Komasa, 2008) is several orders of magnitude less than the equivalent proton
rate, $n_{p}k_{p,O\rightarrow
P}($50\text{\,}\mathrm{K}$)\approx${10}^{-10}\text{\,}{\mathrm{s}}^{-1}$$,
assuming $n_{p}=$1\text{\,}{\mathrm{cm}}^{-3}$$(Gerlich, 1990). In the dark
sector, however, all three interaction rates depend on the dark parameters,
and their relative dominance need not match the Standard-Model case. In this
section, we study the re-scaling of these rates relevant for setting the
ortho-para ratio.

###  C.1 Collisional and radiative scaling

The collisional reactions are straightforward to re-scale by using the results
in Gerlich (1990). First, we approximate the $\text{H}_{2}+p$ cross section
using the Langevin model and reuse our results from Section 4.5. The overall
rate re-scales as $r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}$. The energy
scale associated with the transition is no longer $E_{\rm H}$, however, and is
instead $E_{\rm rot}$, giving the total re-scaling of

$\gamma_{{\rm DM},j\rightarrow
j^{\prime}}(T)=r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\gamma_{{\rm
SM},j\rightarrow j^{\prime}}(\tilde{T}_{r}).$ (C1)

The $\text{H}_{2}-\text{H}$ collision rate for a rotational transition has
already been re-scaled in Section Eq. 2.25,

$\gamma_{{\rm DM},j\rightarrow
j^{\prime}}(T)=r_{\alpha}^{-1}\,r_{m}^{-1}\,r_{M}^{-1}\gamma_{{\rm
SM},j\rightarrow j^{\prime}}(\tilde{T}_{r}).$ (C2)

The radiative transition is more complicated, involving several corrections
and effects beyond the spin-orbit interaction, but Pachucki & Komasa (2008)
provides an excellent analytic guide. The process is composed of two
transitions, the $J+I=2\rightarrow 0$ transition with rate

$\Gamma_{2}\approx\frac{1}{120}\alpha E_{\rm
rot}^{5}\left(\frac{g_{p}R_{0}}{M}\right)^{2}$ (C3)

and the $J+I=1\rightarrow 0$ transition, with rate

$\Gamma_{1}=\frac{2}{9}\alpha E_{\rm
rot}|\langle\Sigma|Q^{k}|\Pi^{k}\rangle|^{2}.$ (C4)

Here, $g_{p}$ is the proton g-factor (independent of the dark parameters),
$R_{0}$ is the average inter-proton distance (scales as $a_{0}$), and the
matrix element $\langle\Sigma|Q^{k}|\Pi^{k}\rangle$ can be written in terms of
dimensionless functions $F_{i}=F_{i}(m\alpha R_{0})$ as

$\displaystyle\langle\Sigma|Q^{k}|\Pi^{k}\rangle$ $\displaystyle=\frac{E_{\rm
rot}^{2}R_{0}}{2M}\left[\left(\frac{g_{p}}{2}-1\right)-\frac{9}{2}\frac{g_{p}F_{1}}{(m\alpha
R_{0})^{3}}\right]-\frac{E_{\rm
rot}}{M^{2}R_{0}}\left[\left(g_{p}-1\right)-g_{p}F_{3}+\frac{9}{4}\frac{g_{p}F_{2}}{(m\alpha
R_{0})^{3}}\right].$ (C5)

Using the definitions from Section 2, and noting that the re-scaling of the
quantity $m\alpha R_{0}$ is $1$, a little algebra clearly shows that
$\Gamma_{2}$ and $\Gamma_{1}$ both re-scale the same way, with the angular-
momentum-averaged rate, $\Gamma_{\rm avg}=(5\Gamma_{2}+3\Gamma_{1})/9$ then
simply re-scaling as

$\Gamma_{\rm avg,DM}=r_{\alpha}^{9}\,r_{m}^{8}\,r_{M}^{-7}\Gamma_{\rm
avg,SM}.$ (C6)

###  C.2 The ortho-para ratio

As shown in Appendix C.1, the $\text{H}_{2}+\text{H}$ and $\text{H}_{2}+p$
ortho-para conversion rates [Eqs. (C1)-(C2)] have a substantially weaker re-
scaling than the radiative conversion rate [Eq. (C6)]. As such, all three need
to be considered when calculating the dark ortho-para ratio.

Assuming radiative transitions remain irrelevant, the behavior of the ortho-
para ratio will depend on the free proton fraction, as is true in the Standard
Model. Effectively, well below the $\text{H}_{2}-\text{H}$ activation
temperature
($\tilde{T}_{r}=r_{\alpha}^{-2}\,r_{m}^{-2}\,r_{M}T\sim$5000\text{\,}\mathrm{K}$$
(Lique et al., 2012)), the ortho-para ratio can be approximated using the
expression from Gerlich (1990) at re-scaled temperature $\tilde{T}_{r}$,

$K_{\rm
DM}(T)\approx\begin{cases}$9.35$\exp(\frac{-169.4}{\tilde{T}_{r}})&\tilde{T}_{r}<$180\text{\,}\mathrm{K}$\\\
3&\tilde{T}_{r}\geq$180\text{\,}\mathrm{K}$.\end{cases}$ (C7)

Near the activation temperature, and if
$n_{p}/n_{\text{H}}\gtrsim${10}^{-3}$\sqrt{r_{m}/r_{M}}$, the $\text{H}_{2}+p$
reaction dominates to keep the ortho-para ratio 3 to 1. If not, the ratio
approaches the $\text{H}_{2}-\text{H}$ value, approximately $7$, as the
temperature increases. Note, however, that the discrepancy between the fully-
ortho cooling rate and the fully-para cooling rate approaches $\mathcal{O}(1)$
at such high temperatures (Glover & Abel, 2008), so the total cooling rate is
insensitive to the actual ortho-para ratio.