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

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


INTRODUCTION
About 84% of the matter in the universe is dark, revealing itself only through its gravitational charge, which is its mass. Various large-scale structure observation campaigns (Bennett et al. 2013;Aghanim et al. 2020;Alam et al. 2021) show that the observed spectrum of large-scale structure is consistent with the cold dark matter (CDM) model, and detailed study of the Bullet Cluster collision shows that the dark matter selfinteraction is much weaker than gas particles (Marke- . There are small-scale phenomena such as cored density-profiles of dwarf galaxies (Ostriker 2003;de Blok 2010), the missing satellite problem (Diemand et al. 2008;Springel et al. 2008;Bullock 2010) and the too-big-to-fail problem (Boylan-Kolchin et al. 2011) that might suggest a deviation from prevailing CDM model, such as self-interacting dark matter (Spergel & Steinhardt 2000) or warm dark matter with non-negligible velocity dispersion (Bode et al. 2001;Sommer-Larsen & Dolgov 2001;Hogan & Dalcanton 2000). None of them, however, have conclusive evidence, and more accurate accounting for baryonic physics can also yield some of these phenomena (see Brooks (2019) for a comprehensive review).
On the other hand, despite decades of searching through direct detection (Agnese et al. 2018;Collaboration et al. 2013;Akerib et al. 2017), indirect detection (Albert et al. 2017;Bergström et al. 2013) and pair cre-ation (Boveia & Doglioni 2018;Mitsou 2015), there is no definitive evidence that dark matter interacts with Standard Model particles. It may be that the dark matter is for all practical purposes completely secluded from the Standard Model. Although this prospect would doom many current detection programs, it still leaves open the possibility of a rich dark sector consisting of one or more new particles interacting via new forces. Because such a model predicts distinctive dark-matter halo structure, it is possible to constrain the dark matter self interactions using only the gravitational interaction which the dark matter must possess (Buckley & Peter 2018).
Here, we consider the particular possibility often referred to as "atomic dark matter" (Goldberg & Hall 1986;Ackerman et al. 2009;Feng et al. 2009;Kaplan et al. 2010Kaplan et al. , 2011Cyr-Racine & Sigurdson 2013;Cyr-Racine et al. 2014;Fan et al. 2013a;Cline et al. 2014;Boddy et al. 2016;Foot & Vagnozzi 2015Randall & Scholtz 2015;Agrawal et al. 2017a;Ghalsasi & Mc-Quinn 2018). In this model, the dark matter consists of a heavy and a light fermion charged under a "dark electromagnetism". This dark matter can form hydrogenlike bound states including molecules. From one perspective, this can be seen as a toy model which captures some of the essential physics of a more complexly interacting dark sector. This model can satisfy constraints on the dark matter self-interaction cross section from, for example, colliding galaxy clusters and dwarf galaxy density profiles (Agrawal et al. 2017b). In fact, dark matter models which introduce a mechanism to suppress the linear matter power spectrum on small scales (including warm dark matter and atomic dark matter) may alleviate challenges facing the ΛCDM model on small scales (Bullock & Boylan-Kolchin 2017).
The parameters of the atomic-dark-matter model are the mass of the heavy particle M , that of the light particle m, the coupling constant α and the ratio of the dark-photon temperature to the cosmic microwave background (CMB) temperature, ξ = T γ,D /T γ . Observations of Big-Bang Nucleosynthesis (BBN) and the CMB temperature anisotropies constrain the energy density in relativistic species other than the three known neutrino generations, parameterized as the change in the effective number of neutrino species ∆N eff . The most recent such constraint is ∆N eff 0.2 (Ade et al. 2016), which restricts ξ 0.4. This upper limit is easy to accommodate since the dark-photon temperature need not be the same as the CMB temperature. In some regions of the parameter space, the atomic dark matter model allows formation of a double disk (Fan et al. 2013b), but this high dark-photon temperature region is strongly constrained by the lack of dark acoustic oscil-lations (DAO) in the observed galaxy power spectrum (Cyr-Racine et al. 2014).
One particularly interesting aspect of atomic dark matter is that the dark matter is dissipative: atomic dark matter halos can radiate away kinetic energy. This means that these halos can undergo collapse, fragmentation, and can eventually form dark compact objects. This raises the prospect of probing an atomic dark matter model through gravitational wave events from their compact-binary coalescence (Shandera et al. 2018). Crucially, the Chandrasekhar mass M Ch ∝ 1/M 2 only depends on the mass of the dark 'proton', and it can be much smaller than the nominal value M Ch 1.4M in the Standard Model. That is, the mass spectrum of merger events, particularly for the sub-solar mass region, can carry information about the physics of even a fully secluded dark sector. Shandera et al. (2018) has showed that a range of model parameters can produce merger rates accessible to the current generation of advanced-LIGO detectors. More recently, Singh et al. (2020) has derived constraints on the dark Chandrasekhar mass and showed that M > 0.95 GeV with 99.9% C.L., if the unusual, so called mass-gap, gravitational wave event GW190425 (Abbott et al. 2020) is interpreted as the merger of an atomic-dark-matter black hole binary.
So far, both Shandera et al. (2018) and Singh et al. (2020) have relied on estimates of the dark black hole binary parameters obtained by re-scaling the Population III star binary literature (Abel 2001;Bromm & Larson 2004;Hartwig et al. 2016). To obtain more realistic constraints, we must tackle directly the complex physics of the formation and evolution of dark-matter structures in the atomic-dark-matter model. As a stepping stone toward that complete study, in this work we consider the cosmological evolution of atomic dark matter. Cyr-Racine & Sigurdson (2013) have previously studied the recombination of dark hydrogen, which we extend by including the formation of molecular states. In the Standard Model, molecular hydrogen is the dominant coolant for pristine (zero metallicity) minihalos with virial temperature below 10 4 K. So, studying the molecular abundance evolution must be a crucial step in determining the dark-halo properties within the atomic-dark-matter paradigm.
Another essential ingredient for these studies is the modification of the linear matter power spectrum introduced by the atomic dark matter. Specifically, atomic dark matter produces dark acoustic oscillations (DAO) due to the propagation of pressure waves in the tightly coupled dark-matter/dark-photon fluid, as well as dark diffusion damping due to the non-zero mean free path of the dark photon. Both effects lead to a suppression of power on small scales, as discussed in detail in Cyr-Racine & Sigurdson (2013). While Cyr-Racine & Sigurdson (2013) studied a relatively high value for the temperature ratio, most extensively ξ = 0.37, we expand the study by varying ξ in particular to a small value, ξ < 0.1. We have also added the molecular processes. We have shown that although molecular processes do change the evolution of the ionization fraction, as long as dark molecular hydrogen is rare relative to dark atoms, the effect is too small to alter the atomic dark matter effects on the linear matter power spectrum.
We find that varying ξ impacts all effects studied. At a given dark photon temperature, this parameter alters the redshift and hence the dark matter number density at a given temperature, which controls freeze-out and molecule formation. Similarly, the diffusion and DAO scales depend on the horizon size at decoupling, which depends on ξ through the decoupling redshift.
We finally investigate the effects of atomic dark matter on the halo mass function in the presence of a cutoff in the linear power spectrum. This approach has been validated by simulations of warm dark matter (Schneider et al. 2013;Angulo et al. 2013) and is an appropriate ansatz for the formation of atomic-dark-matter halos at high redshift.
This paper is the second in a series of three. In Paper I, Ryan et al. (2022a), we defined a re-scaling procedure to compute the molecular chemistry for atomic dark matter from known Standard Model results. Paper III, Ryan et al. (2022b), describes DarkKROME, a modification of the chemistry package KROME (Grassi et al. 2014) to include the chemistry of atomic dark matter. In Paper III we use DarkKROME and the cosmological abundances derived in this paper to explore one-zone collapse scenarios for a range of atomic dark matter parameters. This paper is organized as follows. In Sec. 2, we present the numerical implementation of the cosmological evolution of atomic dark matter states using Rec-fast++ modified by including the relevant molecular processes (Sec. 2.1), the computational method for event rates of these reactions (Sec. 2.2), and the result and an analytic approximation (Sec. 2.5). In Sec. 3, we study the effects of the addition of atomic dark matter and molecules on the formation of large-scale structure by calculating the dark diffusion scale, the DAO scale (Sec. 3.1), and finally estimating the halo mass function at its formation time (Sec. 3.2). We conclude Sec. 4 by pointing out the path toward the complete study of the large-scale structure with atomic dark matter.
Throughout, we shall attach the subscript D to the Standard Model symbols to refer to the corresponding dark-matter symbols. For example, e D + p D → H D + γ D refers to the dark-recombination reaction, the darkmatter analogue of the Standard Model reaction: e + p → H + γ.

The rate equations and implementation
We begin by briefly summarizing the apparatus necessary to solve for the chemical evolution of atomic-darkmatter states, including the recombination and molecule formation, in an expanding universe. To do so, we have implemented the core 1 of Recfast++ (Seager et al. 1999;Chluba et al. 2010) in Julia (Bezanson et al. 2012), using the DifferentialEquations library (Rackauckas & Nie 2017). We present a few details necessary for applying this code, which computes the Standard Model recombination, to the atomic-dark-matter scenario in Appendix A.
The Recfast++ code uses Peeble's effective three-level recombination (Peebles 1968) approximation, where the recombination is treated entirely by tracking the nonequilibrium population in the lowest three (1s, 2s and 2p) atomic states and the continuum. High optical depth in the early universe to Lyman-limit photons prevents the direct recombination to the ground (1s) state (case-B recombination), so the two n = 2 states (2s and 2p) are populated as meta-stable states. The n = 2 states are populated with the rate given by case-B recombination coefficient α HD that includes recombination to higher levels followed by a cascade down to n = 2. They leave the n = 2 state by one of three processes: (A) the forbidden two photon decay (2s → 1s), (B) redshifting of Lyman-α photons out of the resonance band, and (C) photoionization. Exploiting the meta-stability (chemical equilibrium) of the n = 2 state, dx 2 /dt = 0, it is possible to solve for the population of that state and then for the net change in the ionization fraction. This leads to the following equation (Seager et al. 1999), whose right-hand-side terms (in square brackets) represent (i) a change of variables from t to z (time to redshift), (ii) a branching ratio between the ground-state recombinations (two photon decay and redshifted Lyman-α) and all exits from the n = 2 state, including photoionization, and (iii) the net Case-B recombination rate: with the quantities appearing in the equation defined as follows: k B is the Boltzmann constant, h is the Planck's constant, T g,D is the dark-matter gas temperature and T γ,D is the dark-photon temperature, x pD = n pD /n HD is the free proton fraction, x eD the free electron fraction, x HDI is the neutral hydrogen fraction. Note that x eD and x pD differ because we also include the complex states H − D , H + D,2 , H + D,3 in the chemical network. The free electron fraction is calculated as . We also have α HD the case-B recombination rate, Γ α,D ≡ 8πH(z)/ n HDI λ 3 Lyα,D the rate at which redshifted Lyman-α D photons escape the resonance band, ν LyαD , the Lyman-α D frequency, β HD the photoionization coefficient, and A 2γ,D the 2-photon decay rate.
The recombination coefficient is given by the fitting formula of Pequignot et al. (1991a). This formula includes an empirical 1.14 "fudge factor" (Seager et al. 1999) which accounts for the neglected higher (n > 2) level transitions. The fudge factor takes into account all electric dipole transitions from higher (n, ) modes. Because dark hydrogen energy and angular momentum structure is merely re-scaled (Sec. 2.2) from the Standard Model hydrogen, we expect the same fudge factor to be applied to the atomic dark matter case. We calculate the photoionization coefficient β HD by using detailed balance: β HD = α HD (2πmk B T g,D /h 2 ) 3/2 exp(−(E HD − hν Lyα,D )/kT g,D ), with m the mass of the light fermion and E HD the dark hydrogen binding energy. As recombination proceeds, the dark matter and dark photon fall out of thermal equilibrium due to the inefficiency of the Compton-scattering-type energy transfer. This process can be described by solving the following equation for T g,D simultaneously with the equation for x pD (Seager et al. 1999): is the total dark-matter particle number density, σ T,D is the Thomson cross section and a D the radiation constant. We discuss the relationship between these constants in the dark-matter sector and the Standard-Model values in the following subsection (Sec. 2.2).
We solve the ionization fraction and matter temperature equations [Eqs. (2.1)-(2.2)] simultaneously with the equations for the evolution of H − D , H + D,2 , H + D,3 , and H D,2 . These reactions also require the addition of terms which we collectively denote by d∆xp D dz to Eq. (2.1) to account for the effect of the molecular reactions on the free proton fraction. The reaction rates, c Hi , appearing below are as summarized in Table 1 (i is the reaction number in the Table), and have been obtained by re-scaling standard model rates according to the procedure describe in Sec. 2.2 and in Ryan et al. (2022a). The number density of neutral dark atoms is n HDI ≡ . The additional evolution equations are (2.8) Most reactions in Table 1 are those of the minimal model in Galli & Palla (1998), and the rates are those compiled in the same reference, except that we replace the more complicated temperature dependence of H7 and H8 in Galli & Palla (1998) with the simpler power laws taken from Bates & Lewis (1955) and Shapiro & Kang (1987), respectively.

Re-scaling the reaction rates
We have re-scaled the recombination rate with different values of m e , m p , and α as described in Hart & Chluba (2017): where r x denotes the ratio of a dark quantity to its standard model value, for example r M = M/m p . Here, the re-scalings of σ T,D , A 2γ,D , and Γ αD are immediate from the scaling of the hydrogen binding energy, E H ∝ r m r 2 α . The scalings of α HD and β HD are derived in Ryan et al. (2022a) and Hart & Chluba (2017).
Additionally, we take into account the temperature dependence of the recombination and photoionization coefficients by re-scaling the temperature to keep the temperature-to-hydrogen-binding-energy ratio invariant: When the functional form is given in terms of the Standard-Model temperature, we re-scale the darkmatter temperature to find the corresponding Standard-Model temperature by T = r −2 α r −1 m T D for evaluating the rate coefficient. This re-scaled temperature is denoted T a . For example, we calculate the recombination coefficient by (2.11) For the reactions involved in the formation of molecular states we adopt re-scaling procedure described in detail in Ryan et al. (2022a). The reaction rates in Table  1 are fits to either theoretical calculations or experimental data. The form emerges from integrating the cross section over a Maxwell-Boltzmann distribution (or, for radiative reactions, a Planck distribution) with a lower energy cutoff corresponding to the reaction threshold. As a result, the form of the re-scaled rate is given as where g(r α , r m , r M ) contains the entire parametric dependence of the cross section.
The reactions 1 to 15 in Table 1 are considered "important" in Galli & Palla (1998). One might worry that other reactions deemed unimportant in the Standard Model by Galli & Palla (1998) could also contribute for nonstandard parameters. This issue is discussed in more detail in Ryan et al. (2022a). Briefly, the volumetric rates c Hi of the unimportant reactions may be small either due to a strong parametric suppression of the rate coefficient cHi or by the small relative abundance of the reactants. Overcoming these large order-ormagnitude differences requires either an extreme-valued re-scaling factor g(r α , r m , r M ) outside of the range that we consider here, or a similarly large modification of the relative number densities. Such a modification of the relative abundances can occur when ξ is small, where Hard Sphere  Galli & Palla (1998) and the section numbers in parentheses indicate where details can be found in the text. Only the parametric dependence on key quantities are shown for the cross sections (no numerical factors). The reactions included in this table are those that were considered "important" in Galli & Palla (1998) and additions discussed in the text. To leading order the Standard Model rates can be approximated as The b value is from Lepp & Shull (1983) and does not account for density dependence. The b value is from Yoshida et al. (2006).
recombination occurs in a strongly collisional environment. For this reason, we have included the reaction * for the collisional destruction of H 2 , as well as the three body reactions relevant to the formation of H 2 used in Yoshida et al. (2006). We also include the dominant H + D,3 formation and destruction mechanisms, which are suppressed for the Standard Model only by the low fractional abundance of H 2 . While reasonable, these modifications may be insufficient in some parameter regimes. In the Standard Model literature, the minimal reaction list of Galli & Palla (1998) has been largely validated by later, exhaustive studies (i.e. Glover & Savin (2009)). Pending similar work for the dark sector, some caution is warranted when the relative abundances of the species are inverted by orders of magnitude compared to the Standard Model values.

Limitations of Re-scaling
As we have discussed in Ryan et al. (2022a), the rescaling procedure has the following limitations.
First, the theoretical calculations of the molecular cross sections in the Standard Model rely on a Born-Oppenheimer (stationary nuclei) approximation that fails when m ≈ M .
Second, some theoretical cross sections are computed under assumptions of thermal equilibrium between matter and radiation and/or chemical equilibrium between excited states. This assumption is not strictly justified either for the Standard Model or for atomic dark matter, but is adopted regardless for simplicity. For example, Cyr-Racine & Sigurdson (2013) has pointed out that weakly coupled dark matter can thermally decouple from the radiation prior to recombination, clearly breaking the assumption that each reaction is a function of only one temperature. That work has also derived a recombination cross-section with stimulated recombination, which explicitly depends on the radiation temperature, and showed that the stimulated recombination yields less than 60% change in the final free electron fraction. This level of difference is small compared to the difference due to the overall uncertainty in the reaction rates. The suggested rates in the literature for a single reaction can differ by an order of magnitude or more (Savin et al. 2004).
Finally, the accuracy of the re-scaling method is limited by the accuracy of the Standard Model interaction rates. In this work we do not attempt to improve upon the uncertainties in the Standard Model literature. Instead, we scrutinize the sensitivity of the final free electron fraction and molecular hydrogen abundance on the interaction rates. In Sec. 2.4, we show that the uncertainties in the important reaction rates prop-agate through the calculation in a simple manner and demonstrate that the uncertainties at worst linearly affect the final free electron fraction and molecular abundance (Fig. 3).
There are additionally three conditions which limit the range of parameters within which our results can be trusted. The first is that if recombination occurs in a sufficiently high density environment, generally quoted around 10 8 cm −3 in the Standard Model, three body reactions will become important (Palla et al. 1983). We can similarly find the density at which three-body reactions become important by equating the two-body interaction timescale with the three-body interaction timescale t ∼ (n 2 σ 5/2 v) −1 . (2.14) Then, scaling the density at recombination as n ∼ 10 3 ξ −3 r 3 m r 6 α r −1 M and σ ∝ a 2 0 with the Bohr radius a 0 , we find the following condition for three-body reactions to be unimportant at recombination time: We have included a minimal set of three-body reactions, but as discussed above some caution is still warranted in the parameter regime where Eq. (2.15) is violated. A second limiting condition arises from the spectral distortions in the dark photon field induced by the recombination process. In the Standard Model, Hirata & Padmanabhan (2006) showed that the excess Lyman-α and two-photon process contributions to the radiation field suppress the formation of H 2 by a factor of ∼ 5 compared to naive calculations. This suppression is primarily due to ionization of H − by the redshifted spectral distortions, which dominate the radiation field at those frequencies. The size of this effect is thus controlled by the number of spectral distortion photons compared to the number of H − ions. In the Standard Model, roughly 5 photons are produced per hydrogen atom during the recombination process (Chluba & Sunyaev 2006). As long as this number does not depend sensitively on the dark parameters, we do not expect any enhancement of this effect and neglect it for simplicity. Since the energy and angular momentum structure of the dark hydrogen are the same as ordinary hydrogen, we do not expect this number to change as long as the radiation field remains approximately black-body.
However, this is not guaranteed for all values of the parameters. The Standard Model baryon-to-photon ratio η is order 10 −10 , which means that the photons produced by the recombination process do not substantially alter the radiation field, other than through the Lymanα and two photon processes which we explicitly include. Manipulating ξ and M can produce a much larger value of η D . To handle the recombination and molecule formation problems in this case would require a full nonequilibrium treatment of the spectral distortions, since the recombination photons will destroy the black-body nature of the radiation field. To determine when this may occur, we start from the fact that for the Standard Model the spectral distortions due to higher hydrogen transitions are ∆I ν B ν 10 −7 , (2.16) with ∆I ν the intensity of the spectral distortion and B ν the black-body intensity. We insist that the distortions remain less than order unity and insert the scaling of η D to obtain ξ −3 r −1 M < O(10 7 ).
( 2.17) Finally, in the Standard Model, the recombination temperature is lower than the binding energy of hydrogen because of the low cosmological number density. Since the gas temperature is set by the CMB, this density is parameterized by the baryon-to-photon ratio η. Even the weak (logarithmic) dependence of the recombination redshift on η eventually pushes the recombination temperature outside the regime of validity of some of our interaction rates. We omit the parameters for which the re-scaled temperature (see Eq. (2.11)) at the recombination redshift is greater than 10 4 K. This is the highest temperature at which all of the reaction rate laws remain valid. We do not consider parameters where the re-scaled temperature at recombination is much less than in the Standard Model.
When r α 10 −2/3 0.2, satisfying Eq. (2.15), the three body limit, automatically satisfies Eq. (2.17), the spectral distortion condition. For smaller r α , the opposite is true. The third condition, from the temperature range of the reaction rates, seems never to be a problem when the first two conditions are met. It is worth emphasizing that these are practical rather than theoretical limitations. In theory, there is no obstacle to the dark recombination occurring at arbitrarily high η D .

Analytic estimation
We now develop some analytic approximations for the final abundance of molecular states to help us digest the numerical results in Sec. 2.5 and quantify the propagation of the uncertainties in reaction rates. Our approach is based on Lepp & Shull (1984) (see also Stiavelli (2009)).
We begin with the two pathways for the formation of molecular hydrogen known from the Standard Model, the H + D,2 channel: (2.19) and the H − channel: (2.21) These reactions occur in chemical equilibrium: where x i refers to the fraction of H + D,2 or H − D . The equilibrium condition reduces the differential equation for the abundance x i to an algebraic equation. We further simplify this equation by considering only the dominant production and destruction channels for the species. Initially, the second reaction in each pathway is negligible compared to the radiative destruction of H + D,2 and H − D .

For example, the chemical equilibrium equation for H
(2.23) Here, we use the time variable t to emphasize that this chemical equilibrium holds on time scales much shorter than the background evolution time scale, the Hubble time.
As the universe cools adiabatically, chemical equilibrium dictates that x i must increase so that the photodestruction rate continues to balance the creation rate. That is, the H 2 produced from the channel i is where x i,CE is the chemical equilibrium solution of Eq. (2.22). The H 2 production from each channel depends on the free electron abundance. We approximate the residual electron abundance of recombination as the free electron abundance available for the molecule formation. The freeze-out time (z fo ), when the recombination rate equals the Hubble expansion rate, sets the residual electron abundance as c H1 (z fo )n HD x e (z fo ) = H(z fo ). (2.25) Here, we use x e (z) from Saha equilibrium, that is, setting the right hand side of Eq. (2.1) to zero, and neglecting the difference between x eD and x pD : Some molecular hydrogen can be also produced during recombination, but in the Standard Model the amount is very small and is usually neglected. For some dark parameters, however, this first burst of formation dominates the final molecular abundance. In particular, small ξ triggers such a burst because dark recombination happens earlier, so that the collisional H 2 formation processes are enhanced compared to the Standard-Model case. In this case, we estimate the formation of H D,2 under the assumption of chemical equilibrium between the dominant creation and destruction reactions: where the abundances of H + D,2 and H − D are themselves evaluated using chemical equilibrium between their dominant creation and photodestruction processes and the expression is evaluated at z eff = (z rec + z dec )/2. That is, midway between the redshift z rec at which the Saha ionization fraction [Eq. (2.26)] is x e = 0.5 and the redshift of last scattering z dec , the peak of the visibility function g(z) = (dτ /dz)e −τ (z) with the optical depth τ (z) = d n e σ T . This approximation relies on the fact that the H D,2 production is still slow relative to the production and other (besides H D,2 formation) destruction reactions for the intermediates. That is, the equilibrium point between creation and destruction of H + D,2 and H − D is not strongly affected by this burst of H D,2 formation. For the efficiency of molecule formation during recombination, the particle number density n HD at recombination turns out to be the single strongest controlling parameter. This is because, unlike atomic hydrogen, molecular hydrogen is mainly created and destroyed collisionally. In the Standard Model, by the time of recombination (T ∼ 0.25 eV), it is already energetically favorable to form molecules, and in fact we do see a burst of molecule formation in all cases we study in Sec. 2.5 [ Fig. 1], as soon as neutral atoms are available. Nonetheless, this formation channel is completely negligible for the parameters close to the Standard Model. In the Standard Model, this channel is inhibited by two factors: the radiative destruction of H + D,2 and H − D (they are fragile enough to be destroyed even at low temperatures), and the simple unlikeliness of the necessary encounters between species with only a small fractional abundance. Both problems are remedied by increasing the particle number density. For a fixed ρ DM , the darkparticle-to-photon ratio scales as ξ 3 /r M . That is, for low ξ recombination happens at very high redshift, and since n ∝ (1 + z) 3 the number density can be much higher compared to the Standard Model. The dependence on the other parameters is similar: the higher m and α the higher the atomic binding energy and hence the redshift of recombination. As the density at recombination increases, we find that x HD,2 → 0.5; that is, the dark-matter gas approaches a completely molecular state.

Numerical results
In Fig. 1, we show the cosmological evolution from the modified Recfast++ of the atomic dark matter states (ionized, atomic, and molecular) as a function of redshift, renormalized in the unit of z rec , the recombination redshift. For each panel, we show the abundances of H + D (black, solid line), H + D,2 (cyan, dot-dashed line), H − D (red, dashed line), and H D,2 (gold, dotted line). The redshift range we show here spans from the beginning of recombination to the end of molecule formation. The four panels in Fig. 1 show a range of parameters, which are representative of the range of possible evolution histories.
The dots correspond to our analytic estimates of the electron freeze out [Eq. (2.25)] and each of the three channels for H D,2 production [Eqs. (2.24)-(2.27)]. Note that each dot indicates only the additional H D,2 generated from a particular channel, so if this amount is much smaller than the preexisting H D,2 abundance the dot will lie far from the gold curve. In that case, we have verified only that the estimate correctly predicts little enough H D,2 production to negligibly alter the final abundance-not that this tiny fractional increase in x HD,2 agrees exactly with the numerical result. As shown in Fig. 1, the analytic approximation is generally accurate to within an order of magnitude for the final H D,2 fraction.
The two top panels with ξ = 0.02 are qualitatively quite similar to Standard Model (e.g. Galli & Palla (1998); Lepp & Shull (1984)) tracks. Note however that the recombination occurs at z rec ∼ 10 5 in the top left panel, due to the reduced temperature. In the right panel, this is largely offset (z rec ∼ 2000) by the low mass m reducing the hydrogen binding energy. By lowering m or α (or increasing ξ), recombination can be delayed to arbitrarily low redshift or pushed into the future. These parameter choices are less observationally relevant, since, unless the coupling α is extremely small, the ionized dark matter with the dark-electromagnetic interaction would fail to behave like CDM on large scales.
The bottom two panels show more extreme parameters. The bottom left shows "heavy" dark matter with ξ = 0.005. Here, the recombination happens at z ∼ 10 9 in a very high density environment. This precedes even Big-Bang Nucleosynthesis (BBN), but for completely secluded DM, BBN is not affected because the contribution to the total relativistic species from the dark matter is extremely small with ξ = 0.005. In this case we find the asymptotic value of x HD,2 = 0.13. As we have discussed at the end of Sec. 2.4, this reflects a general trend toward forming more molecular hydrogen when the formation occurs at high density. Here, the very large H D,2 abundance allows the formation of H + D,3 to compete with the formation of H D,2 once the H + D,2 photo-dissociation becomes ineffective. The H + D,3 is in turn efficiently converted into neutral hydrogen and H D,2 by dissociative recombination. This leads to a net reduction in x H + D . Finally, the bottom right panel shows r α = 0.274 and ξ = 0.005, and r M = 10. Here we see that even for small values of ξ we can end up with the same qualitative story as in the Standard Model provided that the other parameters are adjusted to control the number density at recombination. Note however that the weak coupling leads to a relatively high freeze-out ionization fraction and relatively low H D,2 fraction.
In Fig. 2, we show explicitly the evolution of the free ion fraction in our full treatment (black, solid), the 3level treatment neglecting the formation of molecules (green, dashed), and in Saha equilibrium (blue, dotdashed) for the same four choices of parameters. Except in the bottom left, the recombination proceeds initially nearly in equilibrium, until the expansion of the universe drives the freeze-out. In the bottom left, the recombination rate is suppressed by the large electron mass (α HD ∝ r −2 m r 2 α ), so the recombination rate drops out of equilibrium (i.e. falls below the Hubble rate) early in the recombination. This leads to early departure from the Saha ionization fraction and a large freeze-out ionization fraction.
We also plot T g,D /T γ,D , illustrating the thermal decoupling of the gas from the dark photons. In the Standard Model, the thermal decoupling is delayed compared to recombination by the low baryon-to-photon ratio, η. The decoupling in the bottom left is qualitatively most similar to the Standard Model, but the delay is due instead to the large freeze-out free electron fraction and high temperature around recombination. In the bottom right, the weak coupling means that the decoupling precedes recombination, while in the top panels the decoupling occurs during recombination.
Note that except in the bottom left the molecules have a minimal effect on the ionization fraction at all times. The H D,2 production pathways do not lead to a net depletion of electrons and (for the parameters considered) they proceed rapidly compared to all other reactions. However, both the formation and destruction of H + D,3 lead to a net recombination. A difference is also discernible in the bottom right, due to the mutual neutralization reaction c H7 converting some of the H − D and p D into neutral hydrogen rather than molecular hydrogen. That reaction rate is enhanced by r 3 α r −3 m ≈ 10 5 , and still leads to only a factor 2 difference in the final ionization fraction. We conclude that except when H D,2 dominates over x H + D , the cosmological scales which depend on opacity to Thomson scattering (discussed below) are weakly affected by the inclusion of molecules We conclude this section by analyzing the propagation of uncertainties in the rates through our calculation, using our analytic approximation. In Fig. 3, for each reaction in turn we have artificially first increased and then decreased the rate by a factor of 10. We use the reaction variations which produce the highest and lowest final x H + D , x eD , and x HD,2 to plot an uncertainty envelope around each track. In the same way, by feeding these artificial rates to our analytic model, we can understand how sensitively our results depend on the accuracy of the reaction rates.
Except in the bottom left panel, the analytic error bars closely match the numerical error bars. In the bottom left, the estimate for the H D,2 abundance fails because the high density, n HD , allows significant molecule formation for a large window around recombination, even in chemical equilibrium. This is not captured by our estimate assuming a short burst of formation with the electron fraction inferred from Saha equillibrium. The estimate for x H + D fails because it does not take into account the H + D,3 recombination channel: the reaction which produces the largest effect numerically (H + D,2 formation, c H8 ) has no effect in the analytic treatment. The circles colored to match each line show our analytic approximations for the freeze out free proton fraction and HD,2 produced from each of the three colors (gold for the recombination channel, cyan for the H + D,2 channel and red for the H − D channel. If a previous channel has already created a large amount of molecular hydrogen these dots can be far from the HD,2 curve and still accurately depict the amount of HD,2 formed through that channel. r m = 1000, r M = 1000, r α = 1.37, ξ = 0.005 We also see a uniquely sensitive dependence of x H + D on this reaction rate. The linear dependence of the freeze-out abundance on the reaction rate revealed in Eq. (2.23) holds only when the abundances of the reactants can be treated as constant. Here we see instead the exponential rate dependence of the solution The electron abundance is not so radically affected, because only the lesser of x H + D and x H + D,3 is subject to exponential depletion, and once Fig. 3 tells us that we can be confident in the accuracy of our results to almost exactly the extent we are confident in the rates of these four key reactions. There are no non-linear dependencies where small errors in certain rates leading to large changes in the final results. This mitigates the impact of uncertainties in the Standard-Model reaction coefficients.
Now we are ready to show the dependence of the molecular-hydrogen fraction at redshift 30 on the model parameters in Fig. 4, by first fixing ξ and α (left panel) while varying the masses, and then fixing the masses while varying ξ and α (right panel). The effect of varying each parameter can be understood through its effect on the number density n HD at fixed photon temperature: a higher number density leads to more molecular hydrogen.

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

Linear Power Spectrum: Acoustic Oscillations and Diffusion Damping
On scales much larger than the horizon size at the dark matter/dark photon decoupling, dissipative dark matter behaves like ordinary cold dark matter. On small scales this is not the case. An exact treatment of the alterations to the linear matter power spectrum in a dissi-pative dark matter model requires solving the Einstein-Boltzmann equations for the given dark matter model, as in Cyr-Racine & Sigurdson (2013). Alternatively, one can adopt a phenomenological parameterization of the matter power spectrum and then map the parameters of this model back to the dark matter microphysics, as in Bohr et al. (2020). Our approach hews closer to that of Bohr et al. (2020), but simpler still.
Atomic dark matter can alter the linear matter power spectrum on small scales primarily through two effects: dark acoustic oscillations (DAO) and dark diffusion (Silk) damping. The following two length scales define the pivot length scales below which the two effects alter the matter power spectrum.
The DAO scale is given as where we calculate the physical sound speed from .
( 3.2) as a function of the scale factor a = 1/(1 + z). The quantity R(a) is the ratio of the relativistic inertias,ρ + P , with the background densityρ and pressureP , of the dark matter and dark photon, given as (3.3) Clearly, for a low dark photon temperature (ξ), R(a) can be very large. The DAO scale depends only on R(a) (which depends most strongly on ξ) and on the scale factor at decoupling, which is controlled dominantly by E HD compared to ξ.
The diffusion damping scale is given by (Zaldarriaga & Harari 1995) (3.4) Note that even for very small or large values of R(a) the diffusion scale is the geometric mean of the photon mean free path and the horizon size, modified by factors of order unity. A weak dependence on M enters through the number density.
Since the DAO scale is suppressed by R(a) while the diffusion scale is not, for low ξ the diffusion scale can far exceed the DAO scale. In this case, the DAOs are completely washed out by diffusion damping. We show the dependence of these two scales on the model parameters in Fig. 5.   Since we find that for the parameters studied in Ryan et al. (2022b) the diffusion scale is larger than the acoustic scale, we calculate the matter power spectrum assuming CDM using CAMB (Lewis et al. 2000) and then apply a Gaussian cutoff (as in e.g. Seljak (1994)): (3.5) On scales smaller than the diffusion scale, any other modifications to the power spectrum are buried under the exponential diffusion damping. This treatment is sufficient to determine whether the formation of the halos studied in Ryan et al. (2022b) is suppressed by the diffusion scale, and whether this suppression will appear at observable scales in the matter power spectrum.

Implications for Dissipative Dark Halos
Equipped with the modified matter power spectrum, we can estimate the halo mass function at various redshifts. In addition to the aforementioned modification to the linear power spectrum, atomic dark matter does not virialize exactly like CDM on small scales. In particular, the dynamics of the collapse will be altered by scattering and subsequent radiative transitions if some of the dark matter is ionized.
However, early in the nonlinear collapse (which takes approximately a Hubble time), the number density does not differ much from background. Therefore, the ionization fraction will not differ much from the freeze-out value. Since the gas is predominantly neutral, it will behave as CDM through this phase. Just prior to virialization, the number density and temperature rise rapidly, but at this point dissipative effects cannot prevent the formation of the halo or alter its mass.
However, the eventual shape and structure of these halos may be quite different if they can cool efficiently subsequent to their formation. Hydrodynamical simulations with cooling will be necessary to resolve these effects.
Motivated by this discussion, we adopt the approach which has been developed for warm dark matter (WDM), and has recently been applied to gravity-only simulations of various self interacting dark matter models (Bohr et al. 2021). Atomic dark matter is similar to WDM in that both models generate a suppression of power at small scales while behaving like CDM on large scales.
In principle, thermal velocities in WDM models can not only lead to a suppression of power on small scales but can also interfere with the virialization of halos. This has been treated in Benson et al. (2012) by employing a "moving barrier", raising the critical density for halos smaller than the suppression scale. Later work (Schneider et al. 2013;Angulo et al. 2013), however, shows that for simulations initialized in the nonlinear regime the thermal velocity effect can be safely neglected compared to the gravitational velocities. That is, the only modification necessary is to the linear matter power spectrum, exactly as in atomic dark matter.
We begin with the analytic prediction for the halo mass function in Press-Schechter theory (Press & Schechter 1974). There, the cumulative mass function of halos of mass greater than M h is given by the fraction of the smoothed linear density field (with smoothing radius R = [3M h /(4πρ M )] 1/3 ) exceeding the critical value δ c 1.69 calculated in the spherical collapse model. Here,ρ M is the comoving background matter density. Then, the differential number of halos of mass M h per volume is given as with σ M h is the root-mean-square of the overdensity smoothed over a region with length scale R: The window function W R (k) is the Fourier transform of the smoothing filter W R (r), and the Press-Schechter multiplicity function is (3.8) Strictly speaking (Bond et al. 1991), the formalism is only valid for the sharp-k filter: (3.10) with p = 0.3, A = 0.322, a = 0.707 (Stiavelli 2009), which has been widely used in the literature.
The problem with applying the Press-Schechter formalism to atomic dark matter is that, for models with a cut-off in the power spectrum, the halo mass function becomes sensitive to the choice of the smoothing filter. Schneider et al. (2013) found that the sharp-k filter (which is anyway the only theoretically justified choice) produces the best agreement with simulations when modifying the mass-to-smoothing-radius relation as where c = 2.7 from calibration to their simulations, and setting a = 1 for the a parameter in the Sheth-Tormen multiplicity function. Note that warm dark matter simulations including Schneider et al. (2013) have been plagued by force resolution issues which have led to spurious small-scale structures. This was largely resolved in Angulo et al. (2013), which justifies using the Press-Schechter formalism with modifications described by Schneider et al. (2013) to obtain a qualitatively correct halo mass function. However, that work cautions that further effort is needed to obtain a halo mass function which is quantitatively accurate for scales near and below the cutoff in the matter power spectrum. These caveats apply equally to atomic dark matter. In Fig. 6, we show the halo mass function z = 0 (top panel) and z = 5 (bottom panel) for the usual CDM case and for the atomic dark matter case with k d = 10h/Mpc (blue, dot-dashed line) and k d = 1h/Mpc (green, dashed line). The diffusion scale is the only relevant variable affecting the halo mass function.

CONCLUSION
In this paper, we have calculated the background evolution of the abundance of dark molecules and ions in the atomic-dark-matter model. To do so, we use the Recfast++ software and the interaction-rate coefficients re-scaled by the methods derived in Ryan et al. (2022a). Previous work in the mirror matter model (equivalent to atomic dark matter with r α = r m = r M = 1 for our purposes) (Latif et al. 2019;D'Amico et al. 2018) assumed the ratio x e /x H2 ≈ 10 2 carried over from the Standard Model, independent of the dark matter parameters. Here, we calculate the primordial molecular abundance directly and find that this assumption often  Figure 6. The halo mass function at z = 0 (upper) and z = 5 (lower). The Gaussian cutoff in the power spectrum leads to a corresponding suppression in the halo mass function at small scales. Clearly, k d = 1 is ruled out by the existence of ∼ 10 11 M galaxies. We include it here to illustrate the dependence of the halo mass function on k d .
fails, because the primordial molecular abundance depends not only on the primordial ionization fraction but also on the dark-particle number density at the time of molecule formation. We have also presented an analytic calculation for the final dark-molecular abundance estimate, which agrees well with the full numerical result. The analytic estimate also indicates that the final abundance is at worst linearly sensitive to the uncertainties in the reaction rates.
We find that when dark recombination happens at high redshift, large dark particle number density promotes the molecular processes so that most of the atomic dark matter particles ends up in molecules. However, as long as the molecule formation is fast compared to the recombination (which is true for all cases studied here) the molecular processes hardly change the recombination history, thereby, the DAO scale and the dark diffusion (Silk) damping scale.
We have also estimated the linear power spectrum of atomic dark matter by calculating the DAO scale and dark-diffusion (Silk) damping scale. When the dark photon temperature is much lower than the CMB temperature, the dark diffusion scale is much larger than the DAO scale, so that we neglect the DAO effect and modify the linear power spectrum by a Gaussian damping factor due to diffusion.
Finally, using the modified linear matter power spectrum, we have estimated the halo mass function using the Press-Schechter formalism modified to fit the warm dark matter simulation results. This work provides the initial conditions and some useful estimates for the results of full atomic dark matter simulations with chemistry, beginning with the companion paper, Ryan et al. (2022b).
In the future, hydrodynamical simulations with radiative cooling will provide insight into the implications of atomic dark matter at scales ranging from the largest structures in the universe to the smallest compact objects. Although precise results for such complex nonlinear physics will require significant investment, this program will begin to promote a broad class of alternative dark matter models-all those in which the dark matter forms hydrogen-like bound states-to testable theories whose parameters can be constrained by observation.

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