Molecular Chemistry for Dark Matter

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.


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 mzr55@psu.edu jhg5248@psu.edu ses47@psu.edu djeong@psu.edu 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. 2018Abbott et al. , 2019Singh 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, α. We take m 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 "mir-ror" 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. 2010Kaplan et al. , 2011Fan et al. 2013;Cyr-Racine & Sigurdson 2013;Cyr-Racine et al. 2014;Cline et al. 2014;Foot & Vagnozzi 2015Randall & 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 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 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(Bromm et al. , 2002Hirano 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 abun-dance 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: 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 semianalytic results used in the literature. 1 Sections 3 and 4 present the cooling functions and reaction rates obtained using these re-scaling techniques. Note that we assume c = = 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 longrange 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 . 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 α 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).

BASIC TOOLS FOR DARK MOLECULAR PHYSICS
The hydrogenic molecule is a system of two heavy particles of mass M and charge +1, at positions X A , X B , and two light particles of mass m and charge −1, at positions x 1 , x 2 , all interacting electromagnetically, with the fine-structure constant α. 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 Here, Ψ(X A , X B ; x 1 , x 2 ) is the wave function for the stationary state with energy E, and we use = 1. We use the shorthanded notation for the distance X AB ≡ |X A − X B |, and so on. As long as M 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: (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 = (α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: where X 0 is the nuclear separation at the minimum of the potential. The typical scale for the parameters are 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: 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 (2.6) At this level of approximation, the molecular rotational energies and low-lying vibrational energies ( ν K0b 1) for generic masses m M and coupling (α 1) can be found from the Standard Model (SM) values via 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 internuclear distances. In this basis, the ansatz for the electronic wave functions u m,X0 is written in terms of coordinates y ij ≡ 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., where P AB (P 12 ) indicates the permutation of A and B (1 and 2), and the coefficients c {n} and β are variational parameters. The individual indices n i are integers ≥ 0, and the summation is over all sets of indices that satisfy 5 i=1 n i ≤ Ω for Ω = 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.7) should remain accurate. Heuristically, we can obtain the re-scaling in Eq. (2.7) by using the spring constant for the vibrational mode which yields the vibrational energy level as (2.10) and the moment of inertia, I = M X 2 0 , which yields the rotational energy level as (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.7). 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 ∆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 /λ u ) 2 where λ u is the wavelength corresponding to the energy gap.

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 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): (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 H , and the constant does not depend on mass (m, M ) nor coupling constant α. This method calculates a binding energy of 3.2 eV (Heitler & London 1927) compared to the measured value of 4.74 eV, and the same approach applied to H + 2 yields an energy of 1.77 eV compared to the true value of 2.64 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 where u 1s is the hydrogenic wavefunction for an effective (screened) nuclear charge Z, and u is the the wave function of the 2p state in a hydrogenic atom with charge 2Z. Both Z and σ 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 eV, which differs by 15.1% from the actual value, and again clearly scales as E 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 H is the "dimensionless temperature" relevant to these reactions. Experimental data suggests that this conclusion is correct: the dissociation energy of H 2 differs from that of D 2 by less than 2% (Herzberg & Monfils 1961), consistent with predominantly E H dependence, insensitive to the mass M .

Spontaneous emission rate
Next, we consider the spontaneous emission rates, or Einstein A coefficients.
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 Ψ ν,J |Q|Ψ ν ,J+2 . 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 ∝ ω 6 Q 2 (Jackson 1998), so equating that with A quad. ω we find the transition rate between states separated by energy ∆E scales as 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 ∝ (length) 2 ∝ a 2 0 regardless. In the equation above we have explicitly written the factors of so that the units on the right hand side are obviously s −1 .
Then, the dark matter Einstein A coefficient for transitions between rotational levels can be estimated by Similarly, for transitions involving only the vibrational level, we estimate the rate by (2.18)

Collisional excitation scattering rates
In the simplest scattering between hydrogen molecules (H 2 -H 2 ), or between a hydrogen molecule and a hydrogen atom (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.7), 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 H 2 -H scattering. The scattering rate coefficients, γ ≡ σv , are determined by the velocity-averaged cross section σ 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 where E is the center of mass collision energy, µ is the reduced mass of the H 2 − H pair, R c ∝ a 0 is the distance of closest approach in a classical head-on collision. Assuming the Morse-type potential, F 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 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 γ comes from the two factors of length in the cross section (R c ∝ a 0 ) and the thermal/kinematic factor from the velocity average, k B T /µ. 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 = ∆E/(k B T ). While F 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 (∆E/y 2 ) = f ∆E (y), and using the dependence of R C ∝ a 0 on fundamental parameters and µ ∝ M , the complete parameteric dependence is (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 deexcitation of rotational levels in the vibrational ground state (ν = 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 Then, the rotational transition rate may be written as ( 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 2 rot as the temperature in ∆E rot units, the relevant energy scale of the problem. For a given dark-matter temperature T DM , we have to evaluate the function f rot at the corresponding Standard Model temperature T SM yielding the same y 2 rot : (2.23) Using the rotational-energy re-scaling in Eq. (2.7), then, we can define a re-scaled temperature, and obtain the final dark matter rate coefficient as 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 rescaling 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 ∝ 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 H 2 -H 2 and 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 assuming the collisional excitation rate similarly contains a dimensionless function that now depends on y vib , f vib (y vib ), and including the remaining terms from Eq. (2.20), we find that for which the re-scaling for the atomic dark-matter becomes Again, we apply this formula for vibrational transition rates from both H 2 -H 2 and H 2 -H collisions.

Summary
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 H 2 has no dipole moment (d = 0), 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 ∝ a 3 0 ∝ 1/(m 3 α 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 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.

APPLICATION: ROVIBRATIONAL COOLING FOR DARK MOLECULAR HYDROGEN
Molecular hydrogen cooling becomes important at low temperatures and high H 2 concentrations in the Standard Model (Glover 2012). In these conditions, molecular line cooling is the dominant process, with certain 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), Table 1. Table of primary quantities important for reaction rates, their dependence on the parameters m, M , and α, and the powers of ratios of parameters (see Eq. (1.1)) needed to re-scale the quantities.
who considered only on H 2 − H and H 2 − 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.

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 erg s −1 cm −3 ) can be written as the sum over all radiative transitions: where n u is the number density of molecules in the upper-energy state u, the are the lower-energy states available for the transition, A u is the Einstein A coefficient, and ∆E u = E u − E 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 rescaled energy of each eigenstate. We have already derived the energies in Section 2, so the energy differences are given approximately as where the vibrational and rotational energies are approximated by Eq. (2.7). 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 gas , then the ratio of the population of the two states is given as the Boltzmann factor: where g u, are the number of internal degrees of freedom for each state. Considering those same levels as a twostate 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 where n X is the number density of the colliding species and γ u (γ u ) is the rate coefficient for the collisional excitation (de-exitation). The steady state solution (dn u /dt = 0) is There are two limiting cases. First, when the radiative decay dominates over the collisional de-excitation, (A u n X γ u ), the population density in the excited states is given as n u = n n X γ u /A u , or the collisionally excited molecules stay in the excited states for the mean lifetime of A −1 u . On the other hand, when the collisional de-excitation dominates over the radiative decay, (n X γ u A u ), the population at the excited state returns to the thermal solution in Eq. (3.3) because n u = n γ u /γ u = n g u /g exp (−∆E u /k B T gas ) from the principal of detailed balance. These limits are commonly denoted as the n → 0 (since n X A u /γ u ) and the 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 (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 crit,u (X), as a benchmark for when collisional deexcitation 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 .
( 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: where L has units of erg cm 3 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 crit,u (X). For example, for a simple two-level system, C = n n X γ u ∆E u at the low-density limit, n X n crit , and C = n u A u ∆E u at the high-density limit, n X n crit .

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 ( 3.9) with which the cooling rate can be written as 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 → 0 limit) so that the de-excitation happens only radiatively, the cooling rate becomes, (3.11) or every collision leads to an emission of photon with energy ∆E u .
In the opposite (n → ∞) limit, we recover the LTE cooling rate as follows. Because the 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: That means when collisional de-excitation dominates over the radiative decay, we recover the LTE cooling rate 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 (3.15) so that C approaches to the correct limit at both lowdensity 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 ( 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 ( 3.17) In reality, this approximation may work better for T < ∆E so that only a few lowest-level states are populated. This is particularly true for the rotational state where ∆E ∝ J, but may be little worse for the vibrational modes where the energy gap stays constant.

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 ν = 0. Let γ H JJ (T ), γ H2 JJ (T ) be the rate coefficient for the (J → J ) transition to occur due to a collision with, respectively, a hydrogen atom and molecule.
At the low-density (n → 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.8), we have where we use the low-density limit of Eq. (3.5), n u = (n n(X)γ u )/A u , to calculate the upper-level population, and set n J=0 /n(H 2 ) = 1/4, n J=1 /n(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 C = X=H,H2 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 where T 3 = T /1000 K. Because the J-dependent part in the second parenthesis is near unity for J = 3, we can see that γ H2 γ 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 for which we evaluate the Standard Model rates at rescaled temperatureT r = r M /(r 2 α r 2 m T ). Finally, substituting into Eq. (3.18) along with the energy re-scaling in Eq. (2.7), gives the re-scaling for the cooling rate (3.21) At the high-density limit, collisional processes dominate so that all rotational levels are in thermal equilib-rium, and the cooling rate is given by Eq. (3.14): The re-scaling results from Section 2 for the Einstein A coefficient [Eq. (2.16)] and the energy levels [Eq. (2.7)] can be used to determine the re-scaled cooling rate in the LTE regime to account for changes in the microphysical parameters as following: We then calculate the cooling rate in the intermediate density regime by interpolation following Eq. (3.15).

Vibrational cooling rates
For cooling by vibrational transitions, we follow Hollenbach & McKee (1979) who approximate the states as a three-level system (ν = 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 ν (dropping the terms higher order in ν): ∆E 10 = ∆E 21 = 0.5∆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: which is nearly the expression given by Hollenbach & McKee (1979) except that we have included the γ 21 term dropped there. The expression including the γ 21 term provides a better fit to the later results of Glover & Abel (2008). As for the collisional de-excitation rates, we use where T is in units of K, then re-scale these rates by using Eq. (2.28): Finally, substituting into Eq. (3.18) along with the energy re-scaling in Eq. (2.7), gives the re-scaling for the cooling rate (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): From the re-scaling results from Section 2 for the Einstein A coefficient [Eq. (2.18)] and the energy levels [Eq. (2.7)], we find the re-scaled cooling rate in the LTE regime:

Rovibrational cooling
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 with X = H, 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) C tot = C rot,LTE + C vib,LTE 1 + (C rot,LTE + C vib,LTE )/C tot,n→0 , 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 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 C H,H2 rot,n→0 and C H,H2 vib,n→0 are usually given in the form of an analytic fit to the sum C H,H2 rot,n→0 + C H,H2 vib,n→0 , 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 H 2 . For literature comparison, we have plotted Λ defined as (n(H) + n(p) + 2 n(H 2 )) 2 . (3.38) For this plot, we choose the values n(H 2 ) = 10 −4 cm −3 , n(H) = 1 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 K ≤ T ≤ 10 4 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 eV; hence, we do not expect to have significant amounts of hydrogen molecules at that temperature.

Collisions with other species
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 rescaling the cooling rate involving collisions with protons and electrons.
The 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 γ DM,v→v (T ) = r −1 α r −3/2 m r −1/2 M γ SM,v→v (T v ). The cooling rates are then derived from the appropriate substitutions into Eqs. (3.18) and (3.24): We also use the Langevin-type cross section for the H 2 − e collision, which leads to a rate rescaling of γ DM,j→j (T ) = r −1 α r −2 m γ SM,j→j (T r ) and γ DM,v→v (T ) = r −1 α r −2 m γ SM,v→v (T v ). Again, substituting appropriately into Eqs. (3.42) Like the H 2 − {H, H 2 } collisions, the rotational and vibrational cooling rates for both the H 2 − p and 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 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 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 H 2 − {e, p} collisional cooling does not contribute meaningfully, regardless of the re-scaling used, until ∼1300 K, where much stronger atomic processes become relevant.

Summary of cooling results
Finally, in Figure 2, we show the entire atomic-darkmatter cooling rate Λ, including both the molecular rovibrational cooling defined above and the combination of cooling processes provided in Rosenberg & Fan (2017). . The atomic cooling rates were derived in Rosenberg & Fan (2017) and the H2 fraction, xH 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.
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 H2 to match the cosmological abundance obtained in Gurian et al. (2022).

APPLICATION: REACTION RATES
The molecular line cooling in the previous section is proportional to the abundance of H and 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, where the first sum represents the source term including all chemical reactions that contain H 2 as one of the products (Formation reactions), and the second sum represents the sink term including all chemical reactions that contain 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 (γ) 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.

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 standardmodel rate γ SM , the corresponding rate for the darkmatter model can be estimated by an overall re-scaling factor multiplying γ SM evaluated at a re-scaled temper-atureT . That is, where g(r α , r m , r M ) is the overall re-scaling, andT ∆ = T /r ∆E with r ∆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 Here, we reuse the two dimensionless parameters defined in Section 2: where ∆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 α.
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"R (with dimensions of length) squared, and a dimensionless integral, analogously to Eq. (2.19). The reactiondependent part of the pre-factor, g(r α , 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, α, 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 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 ofR.
For photochemistry reactions, we have different factors and integrals, but the overall approach is identical, with the substitution of the radiation temperature, T γ , instead of the particle temperature, T , and a slightly different universal thermal factor. Thus, we have with g(r α , r m , r M ) = r 3 ∆E r Rα α r Rm m r R M M r R T ∆E .
(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 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 powerlaw 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 γ SM ∝ 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 subsections 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: While p b and p f are reaction specific, recombination/photoionization-type reactions all have p b = (h ν)/c and p f = m ej v and we obtain the Milne Relation (Mo et al. 2010), where σ rec is the recombination cross section, σ pi is the photonionization cross section, h ν is the photon energy, and m 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 ν = K.E. + ∆E = k B T x 2 + y 2 through conservation of energy (with x and y defined in Eqs. (4.4) -(4.5)) and usually m ej ≈ µ, such that σ rec ∝ T µ x 2 + y 2 2 x 2 σ pi (4.11)     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  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. The b value is from Lepp & Shull (1983) and does not account for density dependence. The b value is from Yoshida et al. (2006).

Reaction 1: Hydrogen recombination
Let us start with one of the simplest reactions, hydrogen recombination: p + e → H + γ. 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) (4.12) where E H = α/(2a 0 ) is the binding energy of hydrogen, g 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 ν = K.E. + ∆E, and we have (4.13) Now that we have the photoionization cross section, application of the Milne relation [Eq. (4.10)] with µ = m gives . (4.14) Substituting this into Eq. (4.6) and definingT a = T /(r 2 α r m ) as the atomic energy analog ofT r andT v , we obtain our final re-scaling which matches the results in Rosenberg & Fan (2017) and Hart & Chluba (2017).
A simplified Standard Model recombination rate is given by Cen (1992) 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, 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.

Reaction 3 and 4: Hydrogen electron attachment and H − photodetachment
The hydrogen electron attachment reaction, H+e − → H − + γ, 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 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 ν = ∆E + K.E. = ∆E + p 2 e /(2m), where ∆E is simply the H − ground state energy, proportional to E H ) such that with γ 2 /2 = 0.027 75 a.u. 2 the electron affinity of hydrogen, and = 2.646 a.u. the effective range (de Jong 1972) in the Standard Model. The effective range , having dimensions of length, must be proportional to a 0 . The electron affinity is defined in Bethe (1949) as γ 2 /2 = (µ ∆E)/(2 ) ∝ m 2 α 2 , so γ ∝ 1/a 0 . Thus, the re-scaling in the two terms effectively cancels, and we can ignore the parametric dependence of the 1/(1 − γ ) term.
We then find the cross-section for hydrogen electron attachment (reaction 3 on our Table) from the Milne relation [Eq. (4.10)] and µ = m, , (4.20) which re-scales as [Eq. (4.6)] (4.21) The full Standard Model rate is given as (Galli & Palla 1998)  Comparing this result with Eq. (4.17), we note that, while both rates use the same overall pre-factor g(r α , r m , r M ), the final re-scalings are quite different. This divergence highlights the necessity of including the r ∆E = r 2 α r m temperature re-scaling factor. Conveniently, Eq. (4.8), (4.19) and µ = m are sufficient to find the re-scaling for the H − photodetachment rate, giving γ 4,DM ≈ r 5 α r m γ 4,SM (T a ) . (4.24)

Reaction 7: mutual neutralization
The cross section for the mutual neutralization reaction, H − + p → 2H, 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, H(1s)H(2s or p) and H(1s)H(3s, p, or d), for our purposes the resulting cross section is effectively where R ∝ 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 2 i ∝ µ K.E. and k 2 f ∝ k 2 i + µ∆E = µ (K.E. + ∆E), F 3 (ζ) = ζ n−1 Γ(n − 1, ζ) − (2ζ) n−1 Γ(n − 1, 2ζ) and Γ(n − 1, ζ) is the incomplete gamma function. The term ζ is a dimensionless ratio of momenta, which, after some simplification, can be written as ζ ∝ k −1 f (αµ∆U )/∆E 2 , and ∆U and ∆E are the zeroth order and true potential energy curve differences. In our parameter regime of interest, ζ 1 and a series expansion gives F 3 (ζ) ≈ ζ. Then we have where, in the second line, we used the fact that ∆U and ∆E are proportional. Then, since ∆E is proportional to the atomic binding energy and µ here is the proton mass, we obtain a net re-scaling of γ 7,DM ≈ r −3 α r −3 m γ 7,SM (T a ) . An interaction without a short-range potential barrier will always occur as long as the incoming particle 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 α = 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. % Difference(Ros, Cen) = 100% |γ Ros −γ Cen | 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: (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 ef f (R) = V (R) + V centrifugal (R), where we include a centrifugal term, V 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: where we have used the property that p ∝ a 3 0 . Then using our tools from Section 4.1 and the definition µ ≈ M , we obtain an approximate rate re-scaling of With the standard assumption that the internal degrees of freedom g 10 and g 15 stay constant, we have σ 15 ∝ σ 10 , and we only need the temperature re-scaling, which, from reaction 10, is proportional to the hydrogen binding energy E H .

H 2 Dissociation Scaling
We approximate the H 2 dissociation reaction, H 2 + H → 3H, at the lowest level as a hard-sphere collision. Then, since the hard-sphere impact parameter, b, scales as a 0 ∝ (αm) −1 , we have (4.31) The dominant binding energy is simply the dissociation energy, D e which depends on α 2 m, so r De = r 2 α r m . The overall re-scaled rate is then given by where the first term in brackets comes from the velocity average (with the approximation µ = m H2 = 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 H 2 + H → 3H 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.

Three-Body-Reactions
We re-scale the reaction rate for the following threebody reactions, taken from Yoshida et al. (2006): 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 γ a n 2 H = γ d n H2 , (4.38) where γ a is the rate coefficient for three-body association reaction and γ d is the rate coefficient for the inverse process. Then, using the Saha equation: where Z = i g i e −Ei/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) Flower (2007) gave Z H2 ≈ .028T , and at the low temperatures where these reactions are significant, this scales as the rotational energy, T →T r . We substitute z H2 (T r ) = r m r −1 M z H2 (T a ) and apply the scaling of γ diss,DM from Section 4.6 to substitute γ diss,DM (T ) = r −1 α r (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 H + 2 and H 2 as well as between H and H + .

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.

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 H − + p → 2H was only known to within an order of magnitude until recently (Glover 2015b), while higher vibrational levels in the H 2 + e → H + 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 H + p → H + 2 + γ has a reaction rate defined on the interval 1 K ≤ T ≤ 32 000 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 γ . 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 M . Trivially, in many cases we approximate the reduced mass µ = (m M )/(m + M ) ≈ m for the e, p system, or µ ≈ m H ≈ M for the H, H (or p) system, which can be easily adjusted as m → 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 → 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 (H + e → H − + γ) assumed a massive stationary, screened proton captures a free electron. In the extreme case, m ≈ M , the physics should become much closer to positronium production, e − + e + → P s or radiative association of H and p, H + p → H + 2 + γ, 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).

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. D + + e D + γ), one is based directly on experimental data. The reaction D 2 + D + → D + 2 + D has a rate derived from the cross section listed in Wang & Stancil (2002) and given by While this appears significantly different to our rescaled result from Galli & Palla (1998), γ 15,DM ≈ 2.12 × 10 −10 cm 3 s −1 exp (−21050/T ), 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 H 2 + H + → H + 2 + 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 K to 10 4.5 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.
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,known → a 2,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 α , r m , r M ) = r −1/2−a/2 µ and ∆E at least approximately independent of µ. Since µ ∝ M for these cases, this holds reasonably well for those reactions with simple cross sections and atomic binding energy, for example, mutual neutralization, H − + p → 2H. Our more complete re-scaling procedure (γ 7,DM ∝ r −0.5 M γ 7,SM ) gives equivalent results to the mass re-scaling (γ 7,DM ∝ r −0.487 M γ 7,SM ). Lastly, we consider the cooling rates and demonstrate why we do not validate using HD chemistry. Unfortunately, we were unable to find a D 2 cooling function based on experimental data or theoretical calculations and can only compare our re-scaled rates to HD cooling  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.
rates, setting M = 1.5m p . However, HD has significantly different behavior, as compared to H 2 , due to its mass asymmetry-induced dipole moment and allowance of J ±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) H 2 cooling rates versus the HD cooling rate from Glover & Abel (2008). Note that while the differences between the HD curve and the rescaled curves are small in the high temperature region, at low temperatures the HD curve is orders of magnitude higher.

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. It may be out-competed by a different channel involving identical reactants, for example γ 6 compared to γ 7 .
2. It may be suppressed below some threshold energy which is large compared to competing reactions. For example reaction 18 requires a photon 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 ±1 transitions enable HD to have an orders of magnitude higher cooling rate.
with enough energy to excite the electronic degrees of freedom of the 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. 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 H 2 .
For the first case, since the basic energy and angular momentum structure of the molecule is preserved under rescaling, 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 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) 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.

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 rescaling 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)  and to model the low-virial temperature halo formation in analogous fashion to the minihalo formation in Standard Model . 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 darkmatter density may reach high enough values so that the re-scaling of other process rates, especially many-body processes such as 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.
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 with temperature T and constants a i , then in the mass-scaled rate the coefficient a 1 takes the form where µ known (µ 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, D − + D + → 2D can be approximately obtained by re-scaling the rate for the regular hydrogen reaction, H − + p → 2H. The hydrogen reaction is known to be of the form of Eq.(A1) with a 1 = 1.4 × 10 −7 cm 3 s −1 and a 2 = −0.487 (Stancil et al. 1998).
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 where B is a mass independent constant with units [length] 2−a [time] a . The rate is given by the thermal average, The integral can be performed exactly, giving Then re-scaling µ gives back the expression found in Eq. (A2), where a 1 = A(a)µ −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: 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 α and m to vary from their Standard Model values. Second, even in the case where r α = r m = 1, the B term in Eq. (A3) is dimensionful, and often contains factors of µ that are not accounted for. Likewise, the dimensionful a 2 and a 3 terms in Eq. (A1) can depend on m, M , and α, significantly changing the rate's temperature dependence. The re-scaling procedure we have presented in this paper overcomes these limitations.

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π 0 = ≡ 1, masses are measured in m e ≡ 1, and energy is measured in units of the Hartree energy, E H = α 2 m e c 2 = 1. With these definitions, c = 1/α, and the unit of time is /E 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 (ν , J + 2) to the state with (ν, J) is (Turner et al. 1977;Flower 2007) A(ν, J ← ν , J + 2) = 1.43 × 10 4 s −1 (E ν ,J+2 − E ν,J ) 5 3(J + 2)(J + 1) 2(2J + 5)(2J + 3) Q 2 ν ,J+2;ν,J = 1 60c 5 (E ν ,J+2 − E ν,J ) 5 3(J + 2)(J + 1) 2(2J + 5)(2J + 3) Q 2 ν ,J+2;ν,J The rates for other transitions, with ∆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, /E H , ( /E H = 2.418 884 × 10 −17 s in S.I. units), c = 137, and the factor of 60 (note the typo in Flower (2007)). From Eq.(B1), 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 H /c 5 = E H α 5 in the pre-factor in Eq. (B1), 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 2 0 , there is no additional parametric dependence from Q 2 . That is, suppose a Standard Model quadrupole in S.I. units is Q ≡Qa 2 0 , whereQ is value of the quadrupole in Standard Model atomic units. Then the dark matter quadrupole can be obtained by Q DM = Q (a 0,DM /a 0 ) 2 . Using dark atomic units on the left hand side, Finally, then, the dark matter Einstein A coefficient for transitions between rotational levels can be estimated by the re-scaling These expressions agree with those obtained in the main text.

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 H 2 + p and H 2 + 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 Γ ≈ 10 −21 s −1 (Pachucki & Komasa 2008) is several orders of magnitude less than the equivalent proton rate, n p k p,O→P (50 K) ≈ 10 −10 s −1 , assuming n p = 1 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.
This condition holds for the example cases given in the main article, but can be violated relatively easily, for example (M, m, α) = (1 GeV, 5.11 MeV, 0.02),T r ≤ 10 4 K and n p = n H /10 4 = 1 cm −3 .

D. RE-SCALING COMBINED COOLING RATE
The literature on molecular hydrogen cooling rates tends to present the low-density limit as a single analytic fit of the full experimental or numerical rovibrational results, C rovib (see e.g. Galli & Palla (1998); Glover & Abel (2008)). This presents a difficulty in applying the re-scaling procedure presented here because the cooling process is composed of different physical processes, apparently complicating our assumption of a single energy scale by which the process can be non-dimensionalized and the fit re-scaled (a problem not unique to molecular cooling). Here, we show how these complications can be resolved by applying piecewise re-scalings in the temperature regime relevant to each physical process.
The re-scaling of the temperature regimes themselves may lead to a gap between the rescaled curves. In this case, the higher temperature process is strongly suppressed below threshold and dominates above threshold, so we extrapolate each re-scaling out to the re-scaled threshold temperature and allow these curves to meet with a discontinuous derivative.
For the case of n → 0 molecular cooling, there are only two processes to consider: rotational and vibrational level transitions. The net rate is the simple sum, C tot = C rot + C vib . Further, rotational transitions dominate at low temperatures (insufficient energy to induce a vibrational transition) and vibrational transitions dominate at high temperatures (larger energy transition), with a clear changeover between the regimes. Defining the pivot temperature T 0 as the temperature where the rotational cooling rate equals the vibrational cooling rate, we can divide the total rovibrational rate into C rovib (T ) = C rovib,T <T0 (T ) + C rovib,T ≥T0 (T ) = R(T ) + V (T ).
Note that R(T ) is only defined on the interval 0 < T ≤ T 0 and V (T ) is only defined on the interval T 0 ≤ T < ∞. From Section 3, we know R and V can be re-scaled as , R DM (T ) = g(r α , r m , r M )R SM (T r ) (D2) V DM (T ) = g(r α , r m , r M )V SM (T v ).
But, since R and V are only defined for T ≤ T 0 and T ≥ T 0 , R DM is only defined for T ≤ (r 2 α r 2 m )/r M T 0 = T 0,r and V DM is only defined for T ≥ (r 2 α r 3/2 m )/r 1/2 M T 0 = T 0,v . Thus, the sum of R DM and V DM only gives the re-scaled rovibrational cooling rate across all T if T 0,r ≥ T 0,v (effectively r M /r m ≤ 1), i.e.
The other case, r M > r m , requires additional information to fill the two gaps, T 0,r < T ≤ T 0,DM and T 0,DM ≤ T < T 0,v , where T 0,DM is the re-scaled pivot point. One option exploits the high (low) temperature behavior of the rotational (vibrational) analytic curves from Hollenbach & McKee (1979). Essentially, the rotational and vibrational curves both behave logarithmically at the temperature extremes, so a linear extrapolation of R DM from T 0,r to T 0,DM and of V DM from T 0,v to T 0,DM provides an acceptable fit. 4 The final full rovibrational cooling curve is given by where linExtrap(F (T )) is the linear extrapolation of the function F (T ). This is the approach taken in Figures 1 and 2 collisions, we find T 0 = 5.40 × 10 3 K and T 0,DM ≈ r 2 α r 1.58 m r −0.586 M T 0,SM . We have been unable to determine an appropriate pivot temperature and corresponding re-scaling for the H 2 − {e, p} collisions as of publication. As such, we will use rotational scaling for the entire temperature range, equivalent to T 0 → ∞. See Section 3.6 for more information.