Molecular Chemistry for Dark Matter III: DarkKROME

Dark matter that is dissipative may cool sufficiently to form compact objects, including black holes. Determining the abundance and mass spectrum of those objects requires an accurate model of the chemistry relevant for the cooling of the dark matter gas. Here we introduce a chemistry tool for dark matter, DarkKROME, an extension of the KROME software package. DarkKROME is designed to include all atomic and molecular processes relevant for dark matter with two unequal-mass fundamental fermions, interacting via a massless-photon mediated $U(1)$ force. We use DarkKROME to perform one-zone collapse simulations and study the evolution of temperature-density phase diagrams for various dark-sector parameters.

The traditional dark-matter model of the WIMP (weakly interacting massive particle), where dark matter has a very simple particle content, has an attractive simplicity. But, dark matter need not be so minimal. Some observations (Bullock & Boylan-Kolchin 2017;de Martino et al. 2020;Cyr-Racine et al. 2021) and some theoretical considerations (Zurek 2014;Petraki & Volkas 2013;Arkani-Hamed et al. 2016;Chacko et al. 2018) suggest that more complex physics may be at work. If dark matter has a richer particle content, it likely has internal chemistry. And, if dark matter has chemistry, it may be able to dissipate sufficient kinetic energy through scattering and atomic-or molecular-like transitions to cool and form compact objects, including black holes (Cline et al. 2014;D'Amico et al. 2017;Shandera et al. 2018;Latif et al. 2019). In that case, fully modeling the evolution of dark-matter structures (and the baryonic structures that trace, albeit biased, the dark structure) requires new numerical tools to evolve the dark matter gas, including all relevant scattering and chemical processes. The modeling in turn will enable data on the abundance and mass spectrum of black holes from gravitational wave observatories to be used to constrain the particle properties of dark matter.
A particular model of dissipative dark matter that is complex yet calculable is the "atomic" dark matter scenario (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. 2013;Cline et al. 2014;Foot & Vagnozzi 2015Randall & Scholtz 2015;Boddy et al. 2016;Agrawal et al. 2017;Ghalsasi & McQuinn 2018). Here dark matter consists of a heavy fermion with mass M and a light fermion with mass m. These particles are oppositely charged by a U (1) force of strength α D , mediated by a dark photon, γ D , that allows the particles to form atoms and molecules, H D and H D,2 , nearly analogous to atomic and molecular hydrogen. When M m, dark-molecular processes can be obtained by a simple re-scaling of Standard Model processes ). There are no weak or strong force analogs in this model, so no quarks, neutrons, muons, etc. While this simplicity may seem ad hoc, it allows the model to be treated in precise numerical detail, providing an important benchmark scenario to calibrate phenomenological treatments applicable to dissipative scenarios more generally.
Cosmologically, an important additional parameter is the ratio of the dark-photon background temperature to the standard-photon background temperature, ξ = T γ,D /T CM B . We allow data considerations to drive the choice of temperature (that is, no input assumption that the dark matter and Standard Model were thermalized at any point) and use ξ 0.4 to be consistent with constraints on additional light degrees of freedom and the lack of observed dark acoustic oscillations (Cyr-Racine et al. 2014). We also assume the existence of a dark matter/anti-matter asymmetry, such that matter dominates, and the net dark U (1) charge of the universe is ≈ 0 (Kaplan et al. 2010).
In two companion papers we have derived the molecular physics of atomic dark matter , and computed the cosmological abundances of the atoms and molecules for a wide range of parameter values ). Here, we present an application of those works: an extension of the software package KROME (Grassi et al. 2014), used to implement chemical and thermal evolution in astrophysical and cosmological simulations, to evolve the chemical network of atomic dark matter. This extension, DarkKROME, is a tool that can help bridge the gap between the semianalytic models of dissipative dark matter, simulations of structure on cluster and galaxy scales, and future gravitational-wave observations. This article is organized as follows. In Section 2, we introduce the DarkKROME extension, explaining how it implements the dissipative dark-matter model. In Section 3 we show how the results of a simple one-zone collapse simulation can vary depending on the masses and coupling strength of the dark matter. We conclude in Section 4.
2. INTRODUCTION TO DarkKROME KROME 1 (Grassi et al. 2014) is a tool to evaluate the thermal and chemical evolution of astrophysical gas clouds, assuming Standard Model particle content. It is used to generate a library of function calls which implement the user-specified reactions at the current gas composition, temperature, and timestep. This library can then be included in simulations of the gravitational dynamics of the gas. KROME provides several simple example simulations, including one-zone spherical cloud collapse and a one dimensional spherical shock, or it can be embedded into more complex simulations (Suazo et al. 2019;Latif et al. 2019;Capelo et al. 2018;Prieto et al. 2015). KROME leaves all dark-matter related physics, including any possible chemistry, to the exterior simulation. DarkKROME extends KROME by enabling it to include dissipative dark reactions and thermal processes in the atomic dark matter model (i.e. the dark chemistry), without changing the overall library structure. DarkKROME is publicly available under the GNU GPLv3 license and can be found at https://bitbucket.org/mtryan83/darkkrome.
The KROME software package provides a Python pre-processor that takes as input a list of chemical reactions and heating and cooling processes, and produces FORTRAN subroutines. These are used by the krome function call to solve the set of (usually) sparse, stiff ordinary differential equations for the time evolution of particle number densities n i and the gas temperature T , In the first equation, the change in number density of species i is determined by all the formation reactions F i , with rates k j and reactants R j , and all the destruction reactions D i , with reactants R j . In the second equation, the change in temperature depends on the adiabatic index, γ, the collective heating, Γ, and cooling, Λ, (both in erg cm −3 s −1 ) in addition to the particle number densitiesn = {n i } and Boltzmann Constant, k B (Grassi et al. 2014). DarkKROME extends the KROME framework to support a new class of chemical species, designated with a Q (i.e. H D would be QH), representing a second sector with a chemistry entirely decoupled from that of the Standard Model, i.e. the dark sector. 2 The sub-atomic particles included are a dark photon, QG, and the two fundamental fermions in the atomic dark matter model, with masses qp mass (M earlier), and qe mass (m). These masses, along with the parameters Dalpha (α D ) and xi (ξ), are additional inputs to DarkKROME with default values equal to their Standard Model analogs. As an example usage, the masses and Dalpha are used in the reaction rates found in the chemical network, react dark, and all four parameters are part of several new thermal process rates. We have not built any dark antimatter parameters into Dark-KROME, although they could be included by the user in the normal KROME fashion. The list of dark matter model parameters can be found in the first section of Table 1.
As appropriate for the atomic-dark-matter model, we have duplicated or extended the KROME subroutines to account for these new species, including reaction number and charge balancing, computing mean molecular weight and adiabatic index, etc. KROME provides significant additional machinery to compute grain physics, cosmic ray chemistry, advanced photochemistry, and other optional features, but we have not implemented corresponding dark versions beyond including a Q-photon, dark CMB flux (darkCMB), and the reaction and heating rates for H D photoionization. We have included three sets of thermal processes: rescaled dark atomic cooling (with -cooling=DARKATOM), analytic dark atomic cooling and heating (with -cooling=ADARKATOM and -heating=ADARKATOM), and dark molecular cooling and heating (with -cooling=DARKMOL and -heating=DARKMOL). The rescaled rates are obtained by extracting the dominant parametric dependence on m, M , and α for each process. Then, given a temperature-dependent rate Λ SM (T ) for a Standard Model process, the re-scaled rate appropriate for the corresponding process in the dark matter is given by where r Λ is a dimensionless product of ratios of dark-tostandard-model parameters, andT is the temperature re-scaled by ratio of dark to Standard Model (predominantly atomic) energy scales. Detailed expressions for the atomic processes, with Standard-Model processes from Cen (1992) as used in KROME, can be found in Appendix A. The analytic rates are those derived in Rosenberg & Fan (2017), with some implementation details in Appendix A. The analytic rates suffer from increased computational complexity and run time. Note that since the -cooling=DARKATOM and -cooling=ADARKATOM options include the same atomic processes, they are mutually exclusive.
The dark atomic cooling rates, both re-scaled and analytic, include contributions from inverse Compton scattering, bremsstrahlung, collisional ionization, collisional excitation, and recombination, and are the dark analogs of the -cooling=ATOMIC plus -cooling=COMPTON and -cooling=FF rates in KROME . The analytic dark atomic heating rate only includes photoionization and is analogous to -heating=PHOTO. Figure 1 shows the slight difference in the net cooling rate, Λ, between the re-scaled and analytic approaches. The difference is predominantly due to the expressions for the collisional excitation cooling rate, and is of a similar level to the difference in Standard Model rates from Cen (1992) and those found in other literature (e.g., Abel et al. (1997)). More detailed comparison is provided in Appendix A.
The dark molecular cooling rate includes contributions from dark molecular hydrogen (H D,2 ) rovibrational cooling and endoergic reactions, analogous to -cooling=H2 plus -cooling=CHEM, while the heating rate includes several exoergic reactions, similar to -heating=CHEM. These re-scaled dark-matter rates were derived in Ryan et al. (2022). We give additional implementation details in Appendix B. Since the rates are computed by re-scaling the pre-factors and temperature dependence of the built-in KROME H 2 rovibrational cooling rates from Glover (2015), Glover & Abel (2008) and chemical thermal rates from Omukai (2000), the rates are identical to the KROME rates when Standard Model parameter values are used. Model parameter values, ξ = 0.01, and assuming chemical equilibrium at a particle density of 1 cm −3 . The difference between the re-scaled and analytic collisional excitation rates dominates the variation between the two rates, with more detail in Appendix A.
As with the original KROME, chemical reaction networks are included as external files, specified during the call to the darkkrome Python pre-processor. We have included an example chemical reaction network, react dark, that contains a minimal chemical network for primordial cloud collapse, as well as some ancillary variable definitions. More details of the react dark file can be found in Appendix C, but there are two important limitations as compared to common primordial chemical networks like the examples provided by KROME (e.g. react primordial*), or used in early universe literature (e.g. Galli & Palla (1998); Glover (2015)). First, the network does not include several subdominant H D,2 and H + D,2 destruction reactions, such as H D,2 + e D → H D + H − D , which reduces the accuracy at temperatures approaching the H D,2 dissociation temperature and above. Second, the network only includes the minimal set of 3-body reactions from Ryan et al. (2022) and does not include any reactions that involve H D,3 , which become important at high (n tot > 10 8 cm −3 (m/511 keV) 3 (α/137 −1 ) 3 ) densities (Glover 2012;Gurian et al. 2022). We leave the addition of these reactions to future work.

ONE-ZONE COLLAPSE
To verify DarkKROME , we explored a simple density evolution model: a one-zone, uniform density cloudcollapse. First, we demonstrate that DarkKROME can reproduce the results from running the one-zone collapse model built into KROME (-test=CollapseZ) with zero metallicity (i.e. only hydrogen and helium in the Standard Model), where the density, ρ, follows freefall or adiabatic evolution, where the free-fall time t f f = (3π)/(32Gρ), and G is the gravitational constant. The thermal processes include DARKATOM/ADARKATOM and DARKMOL heating and cooling, and compressional (adiabatic) heating, defined as We end the simulation before the optically thick regime and so ignore the continuum cooling included in the Col-lapseZ test. For the CollapseZ comparison, we use the initial chemistry parameters including total particle density and species abundances specified in Grassi et al. (2014). Second, we demonstrate that DarkKROME can reproduce results found in dissipative dark matter literature (D'Amico et al. 2017), which previously used the one-zone collapse model provided with KROME to model a "mirror" dark sector, where particle content and parameter values are the same as in the Standard Model, except for ξ. In that example, the density evolution follows Equation 4 unless the sound-crossing time t s = R/c s is shorter than the free-fall time, wherein it follows isobaric evolution In these equations, c s is the speed of sound in the gas, R is the cloud radius, and T 0 is the initial temperature. The compressional heating term is turned off during isobaric evolution. Essentially, this allows for a more accurate evolution of the cloud wherein collapse only occurs if sound waves cannot traverse the cloud faster than a free-fall time.
The initial conditions of these simulations are determined by the cosmological parameters and the primary dark parameters, which set the initial chemistry variables including total particle density and abundance. We use the values for the cosmological parameters from Planck Collaboration et al. (2016), and set ξ = 0.01. We further define the fraction of dissipative dark matter (out of total dark matter), . Given ξ and , the primordial species abundances at the time of structure formation provide the initial species abundances in the cloud, where the primordial species abundance is either taken directly from D'Amico et al. (2017) for comparison with that work or computed using the results from Gurian et al. (2022). Lastly, we need to specify the initial temperature, provided as an input to the simulation. The required cosmological parameters are listed in part 1 of Table 2, with the initial chemistry parameters in part 2.
For further comparison, we show the results of varying the dark parameters, including ξ, and . Doing so introduces additional categories of behaviors and indicates the potential for a wide range in dark, collapsed, halo mass scales.

Verification Results
To check the reproduction, we consider the temperature evolution as a function of total particle density, or n tot vs T . This is common practice in the literature (see e.g. Glover & Abel (2008); Grassi et al. (2014); Yoshida et al. (2006); D'Amico et al. (2017)). The Standard Model one-zone behavior has a direct analog in the radial temperature profile in full 3D hydrodynamical simulations, as seen in Yoshida et al. (2006) and Latif et al. (2019), and consists of initial virialization heating, followed by efficient H 2 rovibrational cooling until density saturation, followed by rapid molecularization due to three-body reactions. Figure 2 demonstrates that DarkKROME can reproduce the results of running the primordial onezone collapse simulation, -test=CollapseZ, provided in KROME, using the initial fractional abundances of x eD = x H + D = 10 −4 , x HD,2 = 10 −6 , and x HD ≈ 1 and starting gas temperature T gas = 450K. DarkKROME very nearly matches the results of KROME, running with zero metallicity and the amount of helium set to zero. If helium is included, there are additional collisional channels that cool the gas at low densities. Dark-KROME does not currently contain the exothermic H D,3 reactions that significantly heat the gas at densities above 10 8 cm −3 (m/511 keV) 3 (α/137 −1 ) 3 . Figure 2. Comparison of KROME and DarkKROME primordial one-zone collapse. Note that the DarkKROME curves exhibit less cooling than the CollapseZ test, as they lack collisional interactions with elements heavier than hydrogen, but otherwise demonstrate the expected qualitative behavior. The figure also demonstrates that there is minimal discernible difference between the different dark atomic cooling options.
In Figure 3, we reproduce the ξ = 0.01 results of D'Amico et al. (2017), with the built-in Standard Model KROME chemical and thermal processes on the left and our added atomic dark matter model chemistry and thermal processes on the right. For both simulations we have used a starting total particle number density of n tot = 2.6 cm −3 , corresponding to z = 40 and Ω ADM = Ω b ( = 0.18) and initial abundances of x eD = x H + D = 10 −8 , x HD,2 = 10 −10 , and x HD ≈ 1, following D' Amico et al. (2017). While the atomic dark matter model does not contain dark helium, the amount of mirror helium is negligible at ξ = 0.01 and does not factor into the comparison (Berezhiani et al. 2001). Each line corresponds to a different initial virial temperature, in the range T 0 =300 K to 15 000 K. Both simulations show three categories of trajectories: efficient molecular cooling (red lines), low-temperature quasiisothermal collapse due to low free ionization (yellow lines), and high-temperature quasi-isothermal collapse due to delayed H D,2 formation (blue lines). The dark chemical network used here contains the relevant reactions to reproduce the main behavior categories, and in general matches the trajectories. There are, however, some differences shown in Figure 3, because we did not include some of the reactions used in the more complete network from D' Amico et al. (2017). In particular, there are fewer H + D,2 destruction channels, leading to slightly higher H D,2 formation and cooling rates in some cases. This difference explains the discrepancy where the highest (lowest) high-temperature quasi-isothermal Standard Model trajectories in the left column of the figure are converted into efficient molecular cooling (lowtemperature quasi-isothermal) dark trajectories, in the right column.

Parameter Exploration Results
The initial H D,2 abundance has a critical role in the evolution of these halos. As mentioned, D'Amico et al.
(2017) used a low initial abundance of x HD,2 = 10 −10 , assuming for simplicity that x eD /x HD,2 = 100, independent of ξ as in Standard Model cosmological recombination, and expecting that this was an overestimate of the true H D,2 fraction (Latif et al. 2019). By solving the background evolution equations, however, Gurian et al. (2022) have shown that, for Standard Model values of m, M , and α but ξ = 0.01, the primordial H D,2 abundance is comparable to the Standard Model value, x HD,2 ≈ 10 −6 at the time of structure formation, even though the free dark ion fraction is much lower. Thus, when using the primordial abundances from Gurian et al. (2022) we obtain the results in Figure 4, where we have also varied . The O(10 4 ) increase in H D,2 computed in Gurian et al. (2022), over the value used in D' Amico et al. (2017) and Figure 3, ensures that halos with sufficient free electrons will undergo efficient molec-ular cooling before reaching the cooling-behavior density transition at n tot ≈ 10 5 cm −3 . That is, the trajectories that exhibited high-temperature quasi-isothermal behavior (blue) in Figure 3 are converted to trajectories with efficient cooling (red) in the top right panel of Figure 4.
The effect of varying is a bit more subtle. Primarily, the initial value of n tot scales linearly with , as does the total particle number density at recombination, n rec . The freeze-out free ion abundance is inversely related to n rec . With the other parameters (and hence z rec ) fixed, this implies x e ( ) ≈ x e ( = 1)/ . Meanwhile, x HD,2 ( ) ∝ x eD n form , where n form is the number density at molecule formation. Since z rec /z form is a constant, so is n rec /n form . Thus, x HD,2 ( ) ≈ x HD,2 ( = 1): decreasing increases the free ion abundance without changing the H D,2 abundance. Varying may also cause the halo to enter isobaric evolution, resulting in cooling. Essentially, the sound above which isobaric evolution occurs. Thus, if the halo is heating adiabatically and crosses the temperature threshold, it generally cools until it can evolve adiabatically again. As the threshold temperature has an and n-dependent minimum, the threshold-crossing behavior produces a characteristic dip in the temperature at low densities, highly noticeable in the low temperature trajectories in the small-panels of Figure 4. The = 0.356 panel demonstrates the maximum value of epsilon where this threshold-crossing behavior occurs, for the given initial conditions. Above this value, the trajectories never heat enough to cross the temperature threshold and so experience purely adiabatic evolution.
Lastly, in Figures 5a-5b we consider dark matter that is fully dissipative ( = 1) and demonstrate how varying the values of m, M , α, and ξ can drastically change the evolution of the halos. In Figure 5a, with m = 250 keV, M = 20 GeV, α = 2/137, and ξ = 0.02 four different behaviors emerge. For these parameters, atomic collisional excitation cooling becomes efficient at approximately 3000 K, so for halos with initial temperatures  2017)) and atomic dark matter for the ξ = 0.01 case. We have assumed xH D,2 = 10 −10 , which is substantially lower than the actual cosmological abundance, for purposes of comparison with D' Amico et al. (2017). In the left column we run the one-zone collapse simulation for a "mirror" dark sector, using the standard KROME processes for a dissipative dark matter fraction of = 0.18. In the right column we use the reduced chemical network of atomic dark matter, described in Section C, and the ADARKATOMIC and DARKMOL cooling options ( Table 2). The first row displays the temperature evolution as a function of total particle density and the middle and bottom row display the dark molecule HD,2(QH2) and dark electron eD (QE) abundances. Each line shows the evolution with different initial temperature T0 and line colors follow the definitions in D' Amico et al. (2017), with red denoting efficient gas cooling, yellow denoting quasi-isothermal collapse in the range (500 − 900)K, and blue denoting quasi-isothermal collapse at 9000K. Also plotted are various temperature thresholds: the lowest HD,2 rotational and vibrational energy transitions, and the low-temperature peak of atomic collisional excitation cooling. Not shown is the HD,2 dissociation temperature, at approximately 5 × 10 4 K. Figure 4. The temperature and density evolution of the dissipative dark matter component in a one-zone collapse simulations, varying the dissipative dark matter fraction , and using the primordial abundances from Gurian et al. (2022). The simulations have otherwise identical initializations to the right column of Figure 3. As in the previous figure, red lines denote efficient gas cooling, and yellow is quasi-isothermal collapse in the range (500 − 900)K. Significantly, no trajectories exhibit the high-temperature quasi-isothermal behaviour seen there, due to the greater initial HD,2 abundance. Varying most directly affects the initial total number density, with low temperature trajectories at lower values also demonstrating the alternating adiabatic/isobaric evolution that results in low-density cooling. The threshold value of = 0.356, below which the isobaric phase of cooling is significant, is derived from Eq.(7).
below that threshold (in cyan), the halo heats adiabatically until atomic cooling and compressional heating balance. Likewise, for the majority of temperatures above the threshold (in blue), the halo immediately cools until the balance is achieved. In both cases, insufficient H D,2 production prevents efficient molecular cooling. In a (relatively) small temperature range (in magenta), however, enough H D,2 is produced to cool the cloud, at least down to below the lowest vibrational transition. The halo is unable to cool down to the rotational regime before transitioning from low-density rovibrational cooling (which scales as n 2 ) to less-densityefficient high-density rovibrational cooling (which scales as n), around n tot ≈ 1 cm −3 . Since compressional heating scales as n 3/2 , as the density increases, it begins to dominate the thermal evolution. With minimal 3-body processes, and without tracking photons and H D,3 reactions, the behavior above n tot ≈ 10 8 cm −3 , i.e. in the high density and opacity regime, is uncertain. Lastly, in the lowest temperature halos, the trace H D,2 production is sufficient for some initial cooling, seen when the trajectory switches from n 2/3 compressional heating to the adiabatic/isobaric oscillatory behavior (which follows n 1/3 behavior, as seen in Equation 7). Unlike the low-temperature trajectories of Figure 4 and the cyan trajectories however, here the temperature is too low for atomic cooling to be relevant, and the molecular cooling has already entered the inefficient, high-density regime. The dissipative-dark-matter component of the halo thus enters a pseudo-equilibrium, where it is cooling back below the isobaric temperature threshold very inefficiently, taking longer and longer to re-enter the adiabatic phase. In this case, the halo remains in said state between 1 and 10 Gyr(black with diamond) or more than 10 Gyr(black with star), at time of simulation termination.
In Figure 5b, with m = 1 MeV, M = 0.1 GeV, α = 137 −1 , and ξ = 0.05, we observe five categories of behavior, with only some of the behaviors in common with the previous parameter set. For starting temperatures below approximately 6000 K, the halo collapses and heats adiabatically until n tot ≈ 600 cm −3 , at which point halos begin alternating between adiabatic and isobaric evolution and the amount of generated H D,2 becomes critical. At the lowest temperatures, the halos enter pseudo-equilibrium, as before. However, for the starting temperatures in the approximately 150 K to 800 K range, the halo does not form sufficient H D,2 early enough to prevent further heating and undergoes the low-temperature quasi-isothermal evolution seen in the Standard Model halos instead (mustard). Note that the later cooling from the long tail of the thermal population is occurring much further below the lowest ro-tational transition than in the Standard Model, as the rovibrational cooling channels (which re-scale inversely to the dark proton mass, see Equations B22-B24) have significantly increased magnitude. At higher starting temperatures, halos achieve the high-temperature quasiisothermal evolution seen in Figure 5a. At the highest starting temperatures, however, we observe the effect of coexisting atomic and molecular processes (lavender). Essentially, both the H − D and H + D,2 paths contribute to H D,2 formation, rapidly forming large amounts of H D,2 , which, combined with the increased atomic cooling, drops the halo once again to the point of lowtemperature quasi-isothermal evolution.

CONCLUSION AND OUTLOOK
We have created an extension of the KROME (Grassi et al. 2014) software package that enables the inclusion of dark sector chemistry in simulations. DarkKROME aims to provide a drop-in replacement for KROME with expanded functionality that can flexibly add dark reactions and thermal processes while still solving the rate equations and providing additional KROME features like charge balancing and reaction checking.
We demonstrated that we can reproduce results found in the literature on dissipative dark matter simulations and can use DarkKROME to explore the parameter space. We showed that the thermal evolution populations of a one-zone cloud collapse model are dependent on the initial abundances as well as the overall dissipative dark matter fraction, and differ from prior literature. Finally, we provided examples of other possible behavior populations that may be encountered at low particle densities in the dark parameter landscape.
As the high-density dark chemistry has not been fully determined at time of publication, we do not speculate on the end points of these clouds and leave more thorough cloud collapse simulations to future work. We anticipate DarkKROME will find significant utilization in those simulations, and, due to the high extensibility of both it and KROME, will further assist other simulations involving dark chemistry, either in the atomic dark matter model or in other dissipative dark matter models.  Figure 5. Temperature versus total particle number density for the specified parameters. In the first panel, four behaviors are observed: pseudo-equilibrium (black with diamonds and stars), initial heating/cooling until atomic cooling balances compressional heating (cyan/blue), and efficient molecular cooling at low densities (magenta). In the second panel we see additional behaviors: low-temperature quasi-isothermal evolution due to insufficient HD,2 production (mustard) and low-temperature quasi-isothermal cooling from rapid HD,2 production (lavender). Also plotted are various temperature thresholds as described in Figure 3, re-scaled as described in Ryan et al. (2022). In the second panel, since T vib /T C.E.peak ≈ 0.4 (1836 m)/M ≈ 1.7, the lowest vibrational energy transition occurs at a higher temperature than the atomic collisional-excitation-cooling peak.

ACKNOWLEDGMENTS
Funding for this work was provided by the Charles E. Kaufman Foundation of the Pittsburgh Foundation. We thank Guido D'Amico for his input on how he used KROME in his simulations. 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.

A. ATOMIC COOLING PROCESS RATES
The KROME software package uses the chemical cooling and heating rates primarily found in Cen (1992) for standard, baryonic matter. DarkKROME provides two alternate sets of cooling rates for dissipative dark matter: re-scaled versions of the Cen rates either directly from or based on the procedure in Ryan et al. (2022) (the -cooling=DARKATOM option), or the analytical expressions from Rosenberg & Fan (2017) (the -cooling=ADARKATOM option). When m = m e = 511 keV, M = m p = 0.938 GeV, and α D = α = 137 −1 , the rates very nearly agree.
In the following subsections, we discuss how all of the rates are re-scaled along with some implementation details for the recombination, collisional ionization, collisional excitation, and bremsstrahlung analytic rates. In the equations below, the cooling rate Λ has units of erg cm −3 s −1 and we assume chemical equilibrium in the figures. We use In addition, we define a temperature re-scaled by the atomic energy scale as In Figure 6, we demonstrate how the individual components contribute to the overall cooling rate, Λ, for Standard Model values and the two sets of dark parameters used in Section 3.1, (m, M, α, ξ) = {(250 keV, 20 GeV, 2/137, 0.02), (1 MeV, 0.1 GeV, 137 −1 , 0.05)}.  Figures {2,5a,5b}. In all cases, the dark gas is assumed to be in chemical equilibrium, the total particle density is ntot = 1 cm −3 , and redshift is z = 40. The slight kink in the recombination rate results from the low-to-high temperature-limit transition. Note the increased cooling in panel (b), the temperature shifts from the change in the dark hydrogen binding energy, and the strong ξ dependence of inverse Compton scattering.

A.1. Recombination
The dark atomic recombination (H + D + e D → H D + γ D ) cooling rate is simply given as the thermal average of the collision kinetic energy summed over all energy levels (Rosenberg & Fan 2017), or To compute the re-scaled rate then, we need the σ rec,n re-scaling, which, from Ryan et al. (2022), is σ rec,n,DM = r 5 α r −2 ∆E σ rec,n .
The kinetic energy term introduces an additional factor of T , so the final re-scaling is where the first term in the first line comes from the T /m term in the thermal average and the second from the collision kinetic energy, and in the third line we used ∆E ∝ mα 2 . For the analytic rate, we avoid the full integral calculation from Rosenberg & Fan (2017)  (5 + y 2 (2.860 + 14/3 log y 2 )) y 1 .
Transitioning at y 2 = 1/4 provides a good fit to the full integral, varying at most by approximately 40 percent at the transition point.

A.2. Collisional Ionization
Like the recombination cooling rate, the collisional ionization cooling rate is simply the reaction rate multiplied by the energy lost, ∆E ∝ mα 2 . To compute the re-scaled reaction rate, we use the simple binary encounter approximation (Peterkops 1977) to obtain the overall parametric dependence of the cross section, followed by the re-scaling procedure. This gives a re-scaled cross section of σ ci,D = r −2 α r −2 m σ ci and the final re-scaling is then While Rosenberg & Fan (2017) also provides the more accurate binary-encounter-Bethe model for the analytic rate, we continue to use the binary encounter approximation, which has an analytic solution involving exponential integrals (Ei(z) = − ∞ −z e −t /t dt), Λ ci n eD n HD = 3.9 × 10 −18 r 2 α r −1/2 As a well-known special function, Ei(z) can be rapidly calculated from library functions, avoiding the computational slow-down of the more complicated binary-encounter-Bethe model, with only an O(1) difference in magnitude (Rosenberg & Fan 2017), comparable to the difference between different Standard Model rate sources. In Figure 7 we compare the re-scaled Cen (1992) rate used in the -cooling=DARKATOM option, the analytic rate based on the binary encounter approximation used in the -cooling=ADARKATOM option, the full binary-encounter-Bethe model, and a more recent rate from Abel et al. (1997) for Standard Model values of m and α. Figure 7. Comparison of the DARKATOM, ADARKATOM, and binary-encounter-Bethe collisional ionization rates for m = 511 keV and α = 137 −1 . The re-scaled KROME rate in DARKATOM (solid blue) is based on Cen (1992), while the analytic rates (red dashed and yellow dot-dashed) are from Rosenberg & Fan (2017). We also plot a more recent rate from Abel et al. (1997) (dotted purple), demonstrating that the differences between the various dark rates are comparable to the difference between the Cen (1992) and Abel et al. (1997) rates.

A.3. Collisional Excitation
To compute the re-scaled collisional excitation cooling rate, we use the same overall cross-sectional re-scaling as collisional ionization, σ ce,D = r −2 α r −2 m σ ce . This matches the dependence found using the Born approximation combined with empirical scaling (Schiff 1968;Kim 2001;Rosenberg & Fan 2017), and leads to the same overall rate re-scaling, Analytically, the 1s → 2p transition dominates over all other ground state transitions and so we ignore other transitions and keep only the leading order term when computing the final rate (with y 2 = (mα 2 )/(2k B T )), Λ ce n H + D n HD = 3.9 × 10 −18 r 2 α r −1/2 m 10 5 T g(y 2 ) (A14) Note that this rate has a higher magnitude than the Cen (1992) rate when calculated using Standard Model values, but as demonstrated by Rosenberg & Fan (2017), the analytic formula has the same level of agreement with the Cen rate as a newer Standard Model rate from Callaway (1994).

A.4. Bremsstrahlung
The nonrelativistic, thermal bremsstrahlung cooling rate, up to a Gaunt factor, has an analytic expression of the form (Rybicki & Lightman 1985), That is, the analytic expression from Rosenberg & Fan (2017) and the re-scaled rate are identical. The default value of g f f,D = 1.5 (darkGauntFF) matches the value used in KROME .

A.5. Inverse Compton Scattering
Like bremsstrahlung, in the nonrelativistic, low-energy limit, the cooling rate due to inverse Compton scattering has an analytic form (Mo et al. 2010 so the re-scaled rate is identical to the rate given in Rosenberg & Fan (2017). We assume the dark photon temperature is set to background, i.e. T γ D = (1 + z)T γ,D = (1 + z) ξ T CMB for redshift z.
whereT R = r 2 α r 2 m r −1 M T andT V = r 2 α r the Λ HD,2,eD cooling rate as Λ HD,2,eD = r 1 α r −1 M Λ H2,e (T r ), Here we have defined T 0 as 855.833 K for H 2 −H and 5402.44 K for H 2 −H 2 , T 0,r = r 2 α r 2 m r −1 M T 0 , T 0,v = r 2 α r 3/2 m r −1/2 M T 0 , and lerp(f (x), x 0 , x) as the linear extrapolation of the function f (x) from the point x 0 . This is the dark equivalent to the cooling found in -cooling=H2.
The net low-density rovibrational cooling rate consists of the sum of the H D,2 − X i collisional rates with various species {X i }. Of note, the H 2 rovibrational cooling from Glover (2015) used in KROME also includes the cooling terms from collisions with helium, assuming it is present in the simulation. Since our dark matter model does not contain neutrons, and thus no helium, this cooling channel has been omitted.

B.2. Endo-and Exoergic Processes
As described in Omukai (2000); Grassi et al. (2014), certain endo-and exoergic reactions contribute significantly to the thermal evolution of the gas. Some of these reactions are included in other heating and cooling options, like atomic collisional ionization and recombination, but the remainder are considered the "chemical" cooling and heating processes and listed in Table 3.

Reaction
Energy/eV Notes Reaction Energy/eV Notes  Table 3. Reactions considered part of the "chemical" heating and cooling processes. Reactions with negative energies cool the gas, positive energies heat. From Ryan et al. (2022), all reactions have energy re-scalings given by eD = r 2 α r m ESM . 1 Forward reaction provides cooling, inverse reaction provides heating.
2 Reaction not included in minimal reaction network. 3 Forward reaction not included in minimal reaction network.
The basic rate equation is quite simple, with where E j the energy consumed/produced and k j is the rate for reaction j, and the net rate is just Λ chem = j Λ j . The dark analogs of E j and k j have already been computed in Ryan et al. (2022). The rates in KROME are weighted by a critical density factor, f = (1 + n cr /n tot ), however, following Hollenbach & McKee (1979), whose parametric dependence must still be determined. The critical density is approximated in the Standard Model rate as with A vib the Einstein A coefficient, x H , x 2 the relative abundances of H and H 2 , and γ H,H2 ∆v the collisional de-excitation rate coefficients for the ∆v vibrational transition. From Ryan et al. (2022), these quantities all re-scale with the dark