# Molecular Chemistry for Dark Matter III: DarkKROME

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

###### Abstract

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. DarkKROME is publicly available at
https://bitbucket.org/mtryan83/darkkrome.

cosmology: theory – dark matter – molecular processes

††software: KROME (Grassi et al., 2014)

##  1 Introduction

Around 84 percent of all matter in the universe appears to be non-baryonic
(Planck Collaboration et al., 2016); its composition is still unknown. Current
constraints on the particle nature of this dark matter are inferred from
cosmological data (Aghanim et al., 2020), from data on the gravitationally
bound structures in the universe on scales ranging from about 10 Mpc - 100 pc
(Sofue & Rubin, 2001; McConnachie, 2012; Strigari, 2013; Collins et al.,
2021), and via dark-matter particle detection experiments on Earth (Schumann,
2019; Lin, 2019; Agrawal et al., 2021; Billard et al., 2021). Gravitational
wave observations from compact object mergers can also constrain dark matter
models, an approach that is increasingly constructive given the recent success
of current ground-based observatories (Goldman & Nussinov, 1989; Abbott et
al., 2005, 2008; Kouvaris & Tinyakov, 2011; de Lavallaz & Fairbairn, 2010;
Bramante & Linden, 2014; Bramante & Elahi, 2015; Bird et al., 2016; Sasaki et
al., 2016; Bramante et al., 2018; Kouvaris et al., 2018; Abbott et al., 2018;
Gresham & Zurek, 2019; Diego, 2020; Abbott et al., 2019; Gow et al., 2020; De
Luca et al., 2020; Singh et al., 2021; Nitz & Wang, 2021a, b, c). To date,
only gravitational interactions between dark and visible matter have been
observed.

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., 2010, 2011; Cyr-Racine & Sigurdson,
2013; Cyr-Racine et al., 2014; Fan et al., 2013; Cline et al., 2014; Foot &
Vagnozzi, 2015, 2016; Randall & 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 $\alpha_{D}$, mediated by a
dark photon, $\gamma_{D}$, that allows the particles to form atoms and
molecules, $\text{H}_{\rm D}$ and $\text{H}_{\rm D,2}$, nearly analogous to
atomic and molecular hydrogen. When $M\gg m$, dark-molecular processes can be
obtained by a simple re-scaling of Standard Model processes (Ryan et al.,
2022). 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,
$\xi=T_{\gamma,D}/T_{CMB}$. 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 $\xi\ll 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 $\approx 0$
(Kaplan et al., 2010).

In two companion papers we have derived the molecular physics of atomic dark
matter (Ryan et al., 2022), and computed the cosmological abundances of the
atoms and molecules for a wide range of parameter values (Gurian et al.,
2022). 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 semi-analytic 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

KROME111The KROME software package and documentation are available at
kromepackage.org.(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$,

$\displaystyle\frac{dn_{i}}{dt}=$ $\displaystyle\sum_{j\in
F_{i}}\left(k_{j}\prod_{r\in R_{j}}n_{r(j)}\right)-\sum_{j\in
D_{i}}\left(k_{j}\prod_{r\in R_{j}^{\prime}}n_{r(j)}\right)$ (1)
$\displaystyle\frac{dT}{dt}=$
$\displaystyle(\gamma-1)\frac{\Gamma(T,\bar{n})-\Lambda(T,\bar{n})}{k_{B}\sum_{i}n_{i}}\,.$
(2)

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}^{\prime}$. In the second equation, the change in temperature depends on
the adiabatic index, $\gamma$, the collective heating, $\Gamma$, and cooling,
$\Lambda$, (both in
$\mathrm{erg}\text{\,}{\mathrm{cm}}^{-3}\text{\,}{\mathrm{s}}^{-1}$) in
addition to the particle number densities $\bar{n}=\\{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. $\text{H}_{\rm D}$ would be QH),
representing a second sector with a chemistry entirely decoupled from that of
the Standard Model, i.e. the dark sector.222We have chosen $Q$ such that there
is no possible overlap with a Standard Model element name. 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 ($\alpha_{D}$) and xi
($\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 DarkKROME, 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.

Parameter |  | Description
---|---|---
qe_mass |  | dark electron mass: $m$
qp_mass |  | dark proton mass: $M$
Dalpha |  | dark fine structure constant: $\alpha_{D}$
xi |  | ratio of dark photon background
| temperature to CMB: $T_{\gamma,D}/T_{CMB}$
DARKATOM |  | Re-scaled atomic cooling
ADARKATOM |  | Analytic atomic cooling/heating
DARKMOL |  | Re-scaled molecular cooling/heating
Table 1: Table of primary constants and rates added in DarkKROME . The first
section lists the main parameters of the dissipative dark matter model used,
and the second lists the sets of cooling/heating processes included. The main
constants are defined as possible additional parameters in the chemistry
network files.

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 $\text{H}_{\rm D}$ photoionization.

We have included three sets of thermal processes: re-scaled 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 re-scaled rates
are obtained by extracting the dominant parametric dependence on $m$, $M$, and
$\alpha$ for each process. Then, given a temperature-dependent rate
$\Lambda_{\rm SM}(T)$ for a Standard Model process, the re-scaled rate
appropriate for the corresponding process in the dark matter is given by

$\Lambda_{\rm DM}(T)=r_{\Lambda}\Lambda_{\rm SM}(\tilde{T})\,,$ (3)

where $r_{\Lambda}$ is a dimensionless product of ratios of dark-to-standard-
model parameters, and $\tilde{T}$ 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, $\Lambda$, 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 ($\text{H}_{\rm 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
$\text{H}_{\rm 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.

Figure 1: Comparison of the DARKATOM (re-scaled) and ADARKATOM(analytic)
atomic cooling functions with Standard Model parameter values, $\xi=0.01$, and
assuming chemical equilibrium at a particle density of
$1\text{\,}{\mathrm{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 $\text{H}_{\rm D,2}$ and $\text{H}_{\rm D,2}^{+}$
destruction reactions, such as $\text{H}_{\rm D,2}+\text{e}_{\rm
D}\rightarrow\text{H}_{\rm D}+\text{H}_{\rm D}^{-}$, which reduces the
accuracy at temperatures approaching the $\text{H}_{\rm 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 $\text{H}_{D,3}$, which become important at high ($n_{\rm
tot}>${10}^{8}\text{\,}{\mathrm{cm}}^{-3}$(m/$511\text{\,}\mathrm{keV}$)^{3}(\alpha/137^{-1})^{3}$)
densities (Glover, 2012; Gurian et al., 2022). We leave the addition of these
reactions to future work.

##  3 One-Zone Collapse

To verify DarkKROME , we explored a simple density evolution model: a one-
zone, uniform density cloud-collapse. 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, $\rho$, follows free-fall or
adiabatic evolution,

$\frac{d\rho}{dt}=\frac{\rho}{t_{ff}},$ (4)

where the free-fall time $t_{ff}=\sqrt{(3\pi)/(32G\rho)}$, and $G$ is the
gravitational constant. The thermal processes include DARKATOM/ADARKATOM and
DARKMOL heating and cooling, and compressional (adiabatic) heating, defined as

$\Gamma_{\rm compress}=\frac{\sum_{i}n_{i}k_{B}T}{t_{ff}}.$ (5)

We end the simulation before the optically thick regime and so ignore the
continuum cooling included in the CollapseZ 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 $\xi$. 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

$\frac{d\rho}{dt}=-\frac{T_{0}\rho_{0}}{T^{2}}\frac{dT}{dt}.$ (6)

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 $\xi=0.01$. We further define the fraction of dissipative dark matter (out
of total dark matter), $\epsilon$. Given $\xi$ and $\epsilon$, 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 $\xi$, and $\epsilon$. Doing so introduces additional categories of
behaviors and indicates the potential for a wide range in dark, collapsed,
halo mass scales.

Parameter |  | Description
---|---|---
$h$ |  | reduced Hubble constant
$\delta_{V}$ |  | virialization overdensity
$\Omega_{M}$ |  | cosmological matter density
$\Omega_{DM}$ |  | cosmological dark matter density
$\epsilon$ |  | dissipative dark matter to all
| dark matter fraction: $\Omega_{DDM}/\Omega_{DM}$
$z_{s}$ |  | redshift at structure formation
$n_{\rm tot}$ |  | total particle number density of DDM
$n_{0}(Q)=x_{0}(Q)\;n_{\rm tot}$ |  | initial number density of
| dark species $Q$
$T_{0}$ |  | initial temperatue
Table 2: Table of required cosmological parameters and initial chemical
parameters required for the one-zone collapse simulations. Cosmological values
combined with the dark parameters are used to compute the initial chemical
parameters if not specified directly. Cosmological values are from Planck
Collaboration et al. (2016).

###  3.1 Verification Results

To check the reproduction, we consider the temperature evolution as a function
of total particle density, or $n_{\rm 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
$\text{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 one-zone collapse simulation, -test=CollapseZ, provided in KROME,
using the initial fractional abundances of $x_{\text{e}_{\rm
D}}=x_{\text{H}_{\rm D}^{+}}=${10}^{-4}$$, $x_{\text{H}_{\rm
D,2}}=${10}^{-6}$$, and $x_{\text{H}_{\rm D}}\approx 1$ and starting gas
temperature $T_{\rm gas}=450\text{K}$. 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. DarkKROME does not currently contain the
exothermic $H_{D,3}$ reactions that significantly heat the gas at densities
above
$${10}^{8}\text{\,}{\mathrm{cm}}^{-3}$(m/$511\text{\,}\mathrm{keV}$)^{3}(\alpha/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 $\xi=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_{\rm tot}=$2.6\text{\,}{\mathrm{cm}}^{-3}$$, corresponding to
$z=40$ and $\Omega_{\rm ADM}=\Omega_{\rm b}$ ($\epsilon=$0.18$$) and initial
abundances of $x_{\text{e}_{\rm D}}=x_{\text{H}_{\rm D}^{+}}=${10}^{-8}$$,
$x_{\text{H}_{\rm D,2}}=${10}^{-10}$$, and $x_{\text{H}_{\rm D}}\approx 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 $\xi=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\text{\,}\mathrm{K}15\,000\text{\,}\mathrm{K}$. Both simulations
show three categories of trajectories: efficient molecular cooling (red
lines), low-temperature quasi-isothermal collapse due to low free ionization
(yellow lines), and high-temperature quasi-isothermal collapse due to delayed
$\text{H}_{\rm 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 $\text{H}_{\rm D,2}^{+}$ destruction channels,
leading to slightly higher $\text{H}_{\rm 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 (low-
temperature quasi-isothermal) dark trajectories, in the right column.

|
---|---
|
|
Figure 3:  Comparison between mirror dark matter results (equivalent to
D’Amico et al. (2017)) and atomic dark matter for the $\xi=0.01$ case. We have
assumed $x_{\text{H}_{\rm 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 $\epsilon=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 $\text{H}_{\rm D,2}$(QH2)
and dark electron $\text{e}_{\rm D}$ (QE) abundances. Each line shows the
evolution with different initial temperature $T_{0}$ 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)$\mathrm{K}$$, and blue denoting quasi-isothermal collapse at
$9000$\mathrm{K}$$. Also plotted are various temperature thresholds: the
lowest $\text{H}_{\rm D,2}$ rotational and vibrational energy transitions, and
the low-temperature peak of atomic collisional excitation cooling. Not shown
is the $\text{H}_{\rm D,2}$ dissociation temperature, at approximately
$5\text{\times}{10}^{4}\text{\,}\mathrm{K}$.

###  3.2 Parameter Exploration Results

The initial $\text{H}_{\rm 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_{\text{H}_{\rm D,2}}=${10}^{-10}$$, assuming for
simplicity that $x_{\text{e}_{\rm D}}/x_{\text{H}_{\rm D,2}}=100$, independent
of $\xi$ as in Standard Model cosmological recombination, and expecting that
this was an overestimate of the true $\text{H}_{\rm 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
$\alpha$ but $\xi=0.01$, the primordial $\text{H}_{\rm D,2}$ abundance is
comparable to the Standard Model value, $x_{\text{H}_{\rm
D,2}}\approx${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 $\epsilon$. The $\mathcal{O}(10^{4})$ increase in
$\text{H}_{\rm 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 molecular cooling before reaching the
cooling-behavior density transition at $n_{\rm
tot}\approx${10}^{5}\text{\,}{\mathrm{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 $\epsilon$ is a bit more subtle. Primarily, the initial
value of $n_{\rm tot}$ scales linearly with $\epsilon$, as does the total
particle number density at recombination, $n_{\rm rec}$. The freeze-out free
ion abundance is inversely related to $n_{\rm rec}$. With the other parameters
(and hence $z_{\rm rec}$) fixed, this implies $x_{e}(\epsilon)\approx
x_{e}(\epsilon=1)/\epsilon$. Meanwhile, $x_{\text{H}_{\rm
D,2}}(\epsilon)\propto x_{\text{e}_{\rm D}}n_{\rm form}$, where $n_{\rm form}$
is the number density at molecule formation. Since $z_{\rm rec}/z_{\rm form}$
is a constant, so is $n_{\rm rec}/n_{\rm form}$. Thus, $x_{\text{H}_{\rm
D,2}}(\epsilon)\approx x_{\text{H}_{\rm D,2}}(\epsilon=1)$: decreasing
$\epsilon$ increases the free ion abundance without changing the
$\text{H}_{\rm D,2}$ abundance.

Varying $\epsilon$ may also cause the halo to enter isobaric evolution,
resulting in cooling. Essentially, the sound-crossing time behaves as
$t_{s}\propto\left(\epsilon_{M}/\rho_{\rm DDM}\right)^{1/3}T^{-1/2}$, where
$\rho_{\rm DDM}$ is the local dissipative dark matter density,
$\epsilon_{M}=\epsilon\,\Omega_{DM}/\Omega_{M}$, while the free-fall time
behaves as $t_{ff}\propto\left(\rho_{M}(1-\epsilon_{M})+\rho_{\rm
DDM}\right)^{-1/2}$, with $\rho_{M}$ the total local matter density. The
densities are written separately here because $\rho_{M}(1-\epsilon_{M})$ is
constant in these simulations and in general may evolve differently than
$\rho_{\rm DDM}$. The isobaric transition condition can then be transformed
into a temperature threshold $T_{iso}$, where

$T_{iso}\propto\epsilon_{M}^{2/3}\frac{\rho_{M}(1-\epsilon_{M})+\rho_{\rm
DDM}}{\rho_{\rm DDM}^{2/3}},$ (7)

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
$\epsilon$ 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-$\epsilon$ panels of Figure 4.
The $\epsilon=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.

|
---|---
|
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 $\epsilon$, 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)$\mathrm{K}$$. Significantly, no trajectories exhibit the high-
temperature quasi-isothermal behaviour seen there, due to the greater initial
$\text{H}_{\rm D,2}$ abundance. Varying $\epsilon$ most directly affects the
initial total number density, with low temperature trajectories at lower
$\epsilon$ values also demonstrating the alternating adiabatic/isobaric
evolution that results in low-density cooling. The threshold value of
$\epsilon=0.356$, below which the isobaric phase of cooling is significant, is
derived from Eq.(7).

Lastly, in Figures 5a-5b we consider dark matter that is fully dissipative
($\epsilon=1$) and demonstrate how varying the values of $m$, $M$, $\alpha$,
and $\xi$ can drastically change the evolution of the halos. In Figure 5a,
with $m=$250\text{\,}\mathrm{keV}$$, $M=$20\text{\,}\mathrm{GeV}$$,
$\alpha=2/137$, and $\xi=0.02$ four different behaviors emerge. For these
parameters, atomic collisional excitation cooling becomes efficient at
approximately $3000\text{\,}\mathrm{K}$, so for halos with initial
temperatures 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
$\text{H}_{\rm D,2}$ production prevents efficient molecular cooling. In a
(relatively) small temperature range (in magenta), however, enough
$\text{H}_{\rm 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-density-efficient high-density rovibrational
cooling (which scales as $n$), around $n_{\rm
tot}\approx$1\text{\,}{\mathrm{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
$\text{H}_{D,3}$ reactions, the behavior above $n_{\rm
tot}\approx${10}^{8}\text{\,}{\mathrm{cm}}^{-3}$$, i.e. in the high density
and opacity regime, is uncertain. Lastly, in the lowest temperature halos, the
trace $\text{H}_{\rm 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\text{\,}\mathrm{Gyr}$(black with diamond) or more than
$10\text{\,}\mathrm{Gyr}$(black with star), at time of simulation termination.

In Figure 5b, with $m=$1\text{\,}\mathrm{MeV}$$,
$M=$0.1\text{\,}\mathrm{GeV}$$, $\alpha=137^{-1}$, and $\xi=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\text{\,}\mathrm{K}$, the halo collapses and heats adiabatically until
$n_{\rm tot}\approx$600\text{\,}{\mathrm{cm}}^{-3}$$, at which point halos
begin alternating between adiabatic and isobaric evolution and the amount of
generated $\text{H}_{\rm D,2}$ becomes critical. At the lowest temperatures,
the halos enter pseudo-equilibrium, as before. However, for the starting
temperatures in the approximately
$150\text{\,}\mathrm{K}800\text{\,}\mathrm{K}$ range, the halo does not form
sufficient $\text{H}_{\rm 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
rotational transition than in the Standard Model, as the rovibrational cooling
channels (which re-scale inversely to the dark proton mass, see Equations
B1-B3) have significantly increased magnitude. At higher starting
temperatures, halos achieve the high-temperature quasi-isothermal 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 $\text{H}_{\rm D}^{-}$ and $\text{H}_{\rm D,2}^{+}$
paths contribute to $\text{H}_{\rm D,2}$ formation, rapidly forming large
amounts of $\text{H}_{\rm D,2}$, which, combined with the increased atomic
cooling, drops the halo once again to the point of low-temperature quasi-
isothermal evolution.

(a) $m=$250\text{\,}\mathrm{keV}$$, $M=$20\text{\,}\mathrm{GeV}$$,
$\alpha=2/137$ and $\xi=0.02$

(b) $m=$1\text{\,}\mathrm{MeV}$$, $M=$0.1\text{\,}\mathrm{GeV}$$,
$\alpha=137^{-1}$ and $\xi=0.05$

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 $\text{H}_{\rm D,2}$ production (mustard) and low-temperature
quasi-isothermal cooling from rapid $\text{H}_{\rm D,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_{\rm vib}/T_{\rm C.E.peak}\approx 0.4\sqrt{($1836$\,m)/M}\approx
1.7$, the lowest vibrational energy transition occurs at a higher temperature
than the atomic collisional-excitation-cooling peak.

##  4 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.

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.

##  Appendix 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\text{\,}\mathrm{keV}$$,
$M=m_{p}=$0.938\text{\,}\mathrm{GeV}$$, and $\alpha_{D}=\alpha=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 $\Lambda$ has units of
$\mathrm{erg}\text{\,}{\mathrm{cm}}^{-3}\text{\,}{\mathrm{s}}^{-1}$ and we
assume chemical equilibrium in the figures. We use

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

In addition, we define a temperature re-scaled by the atomic energy scale as

$\tilde{T}_{a}=\frac{T}{r_{\alpha}^{2}r_{m}}\,.$ (A2)

In Figure 6, we demonstrate how the individual components contribute to the
overall cooling rate, $\Lambda$, for Standard Model values and the two sets of
dark parameters used in Section 3.1,
$(m,M,\alpha,\xi)=\\{($250\text{\,}\mathrm{keV}$,$20\text{\,}\mathrm{GeV}$,2/137,0.02),\,($1\text{\,}\mathrm{MeV}$,$0.1\text{\,}\mathrm{GeV}$,137^{-1},0.05)\\}$.

(a)

(b)

(c)

Figure 6: Plots demonstrating the contributions to the net dark atomic cooling
process and their parameter dependence. The parameter values are
$(m/$\mathrm{keV}$,M/$\mathrm{GeV}$,\alpha,\xi)=\\{(511,0.938,137^{-1},1),(250,20,2/137,0.02),(10^{3},0.1,137^{-1},0.05)\\}$,
corresponding to panels {6a,6b,6c} and the parameter values used in Figures
{2,5a,5b}. In all cases, the dark gas is assumed to be in chemical
equilibrium, the total particle density is $n_{\rm
tot}=$1\text{\,}{\mathrm{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 $\xi$
dependence of inverse Compton scattering.

###  A.1 Recombination

The dark atomic recombination ($\text{H}_{\rm D}^{+}+\text{e}_{\rm
D}\rightarrow\text{H}_{\rm D}+\gamma_{D}$) cooling rate is simply given as the
thermal average of the collision kinetic energy summed over all energy levels
(Rosenberg & Fan, 2017), or

$\Lambda_{\rm rec}=\sum_{n=1}^{\infty}\langle(\text{K.E.})\,\sigma_{{\rm
rec},n}v\rangle n_{\text{H}_{\rm D}^{+}}n_{\text{e}_{\rm D}}.$ (A3)

To compute the re-scaled rate then, we need the $\sigma_{{\rm rec},n}$ re-
scaling, which, from Ryan et al. (2022), is

$\sigma_{{\rm rec},n,{\rm DM}}=r_{\alpha}^{5}r_{\Delta E}^{-2}\sigma_{{\rm
rec},n}.$ (A4)

The kinetic energy term introduces an additional factor of $T$, so the final
re-scaling is

$\displaystyle\Lambda_{\rm rec,DM}(T)$ $\displaystyle=\sqrt{\frac{r_{\Delta
E}}{r_{m}}}r_{\Delta E}\left(r_{\alpha}^{5}r_{\Delta
E}^{-2}\right)\Lambda_{\rm rec}\left(\frac{T}{r_{\Delta E}}\right)$ (A5)
$\displaystyle=r_{\alpha}^{5}\,r_{m}^{-1/2}r_{\Delta E}^{-1/2}\Lambda_{\rm
rec}\left(\frac{T}{r_{\Delta E}}\right)$ (A6)
$\displaystyle=r_{\alpha}^{4}\,r_{m}^{-1}\Lambda_{\rm rec}(\tilde{T}_{a}),$
(A7)

where the first term in the first line comes from the $\sqrt{T/m}$ term in the
thermal average and the second from the collision kinetic energy, and in the
third line we used $\Delta E\propto m\alpha^{2}$. For the analytic rate, we
avoid the full integral calculation from Rosenberg & Fan (2017) to decrease
computation time and simply use their high and low limits (defining
$y^{2}=(m\,\alpha^{2})/(2k_{B}T)$, where $k_{B}$ is the Boltzmann constant),

$\frac{\Lambda_{\rm rec}}{n_{\text{e}_{\rm D}}n_{\text{H}_{\rm
D}^{+}}}=\begin{cases}4.7\times
10^{-25}\,r_{\alpha}^{3}\,r_{m}^{-3/2}\left(\frac{T}{10^{5}}\right)^{1/2}(0.74+\log{y^{2}}+\frac{1}{3y^{2}})&y\gg
1\\\ 1.1\times
10^{-25}\,r_{\alpha}^{5}\,r_{m}^{-1/2}\left(\frac{10^{6}}{T}\right)^{1/2}(5+y^{2}(2.860+14/3\log{y^{2}}))&y\ll
1\end{cases}.$ (A8)

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, $\Delta E\propto
m\alpha^{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 $\sigma_{\rm
ci,D}=r_{\alpha}^{-2}\,r_{m}^{-2}\sigma_{\rm ci}$ and the final re-scaling is
then

$\displaystyle\Lambda_{\rm ci,DM}(T)$ $\displaystyle=\sqrt{\frac{r_{\Delta
E}}{r_{m}}}r_{\Delta E}\left(r_{\alpha}^{-2}\,r_{m}^{-2}\right)\Lambda_{\rm
ci}\left(\frac{T}{r_{\Delta E}}\right)$ (A9)
$\displaystyle=r_{\alpha}\,r_{m}^{-1}\Lambda_{\rm ci}(\tilde{T}_{a}).$ (A10)

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
(${\rm Ei}(z)=-\int_{-z}^{\infty}e^{-t}/t\,dt$),

$\displaystyle\frac{\Lambda_{\rm ci}}{n_{\text{e}_{\rm D}}n_{\text{H}_{\rm
D}}}$
$\displaystyle=$3.9\text{\times}{10}^{-18}$r_{\alpha}^{2}\,r_{m}^{-1/2}\sqrt{\frac{${10}^{5}\text{\,}\mathrm{K}$}{T}}f(y^{2})$
(A11) $\displaystyle f(y^{2})$
$\displaystyle=\frac{1}{2}\left(e^{-y^{2}}+y^{2}{\rm Ei}(-y^{2})\right).$
(A12)

As a well-known special function, ${\rm 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 $\mathcal{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 $\alpha$.

Figure 7: Comparison of the DARKATOM, ADARKATOM, and binary-encounter-Bethe
collisional ionization rates for $m=$511\text{\,}\mathrm{keV}$$ and
$\alpha=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, $\sigma_{\rm
ce,D}=r_{\alpha}^{-2}\,r_{m}^{-2}\sigma_{\rm 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,

$\displaystyle\Lambda_{\rm ce,DM}(T)$
$\displaystyle=r_{\alpha}\,r_{m}^{-1}\Lambda_{\rm ce}(\tilde{T}_{a}).$ (A13)

Analytically, the $1s\rightarrow 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\alpha^{2})/(2k_{B}T)$),

$\displaystyle\frac{\Lambda_{\rm ce}}{n_{\text{H}_{\rm D}^{+}}n_{\text{H}_{\rm
D}}}$
$\displaystyle=$3.9\text{\times}{10}^{-18}$\,r_{\alpha}^{2}\,r_{m}^{-1/2}\sqrt{\frac{10^{5}}{T}}g(y^{2})$
(A14) $\displaystyle g(y^{2})$
$\displaystyle=\int_{\frac{\sqrt{3}}{2}y}^{\infty}\text{d}u\frac{u\text{e}^{-u^{2}}}{1+\frac{7y^{2}}{4u^{2}}}\log{\frac{4u}{y}}$
(A15)

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

$\displaystyle\Lambda_{\rm ff,D}(T)$ $\displaystyle=\sqrt{\frac{2\pi
k_{B}T}{3m}}\frac{16\hbar^{2}\alpha^{3}c}{3m}n_{\text{e}_{\rm
D}}n_{\text{H}_{\rm D}^{+}}g_{ff,D}$ (A16)
$\displaystyle=$1.4\text{\times}{10}^{-27}$\,r_{\alpha}^{3}\,r_{m}^{-3/2}\left(\frac{T}{1\text{K}}\right)^{1/2}n_{\text{e}_{\rm
D}}n_{\text{H}_{\rm D}^{+}}g_{ff,D}$ (A17)
$\displaystyle=r_{\alpha}^{4}\,r_{m}^{-1}\Lambda_{\rm ff}(\tilde{T}_{a}).$
(A18)

That is, the analytic expression from Rosenberg & Fan (2017) and the re-scaled
rate are identical. The default value of $g_{ff,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),

$\displaystyle\Lambda_{\rm ics,D}$
$\displaystyle=\frac{4k_{B}(T_{\text{e}_{\rm
D}}-T_{\gamma_{D}})}{mc^{2}}c\left[\frac{8\pi}{3}\left(\frac{\alpha\hbar}{mc}\right)^{2}\right]n_{\text{e}_{\rm
D}}a_{r}T_{\gamma_{D}}^{4}$ (A19)
$\displaystyle=$1.0\text{\times}{10}^{-37}$\left(\frac{T-T_{\gamma_{D}}}{$1\text{\,}\mathrm{K}$}\right)r_{\alpha}^{2}\,r_{m}^{-3}\left(\frac{T_{\gamma_{D}}}{$1\text{\,}\mathrm{K}$}\right)^{4}$
(A20) $\displaystyle=r_{\alpha}^{2}\,r_{m}^{-3}\Lambda_{\rm
ics}(T,T_{\gamma_{D}}),$ (A21)

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_{\gamma_{D}}=(1+z)T_{\gamma,D}=(1+z)\,\xi\,T_{\rm CMB}$ for redshift $z$.

##  Appendix B Dark Molecular Processes

For molecular processes DarkKROME uses the rates derived in Ryan et al.
(2022). The added molecular reactions are listed in Table 4 and we describe
our implementation of the thermal processes here. We also give a brief
derivation of the $H_{2}+H\rightarrow 3H$ reaction scaling. Both $H_{2}$
collisional cooling and the chemical reaction cooling processes are included
by selecting -cooling=DARKMOL.

###  B.1 $H_{2}$ Collisional Cooling

$H_{2}$ collisional cooling, also referred to as line cooling, involves the
collisional excitation of $H_{2}$ and subsequent radiative cooling. Our basic
approach to implementing the dark version involves re-scaling the line cooling
rates for $\text{H}_{\rm D,2}-\\{\text{H}_{\rm D},\text{H}_{\rm
D,2},\text{e}_{\rm D},\text{H}_{\rm D}^{+}\\}$ collisions from Glover (2015)
using the process described in Ryan et al. (2022). Effectively, this can be
condensed to re-scaling the high-density rotational and vibrational cooling
rates, $H_{R},H_{V}$ (using the KROME notation Grassi et al. (2014))

$\displaystyle H_{R,D}(T)$
$\displaystyle=r_{\alpha}^{9}\,r_{m}^{8}\,r_{M}^{-6}H_{R}(\tilde{T}_{R})$ (B1)
$\displaystyle H_{V,D}(T)$
$\displaystyle=r_{\alpha}^{9}\,r_{m}^{5}\,r_{M}^{-3}H_{V}(\tilde{T}_{V}),$
(B2)

where $\tilde{T}_{R}=r_{\alpha}^{2}\,r_{m}^{2}\,r_{M}^{-1}\,T$ and
$\tilde{T}_{V}=r_{\alpha}^{2}\,r_{m}^{3/2}\,r_{M}^{-1/2}\,T$, the low-density
rovibrational cooling rate, $\Lambda_{\text{H}_{\rm D,2},\\{\text{H}_{\rm
D},\text{H}_{\rm D,2}\\}}$ as

$\Lambda_{\text{H}_{\rm D,2},\\{\text{H}_{\rm D},\text{H}_{\rm
D,2}\\}}=\begin{cases}r_{\alpha}\,r_{m}\,r_{M}^{-2}\;\begin{cases}\Lambda_{\text{H}_{2},\\{\text{H},\text{H}_{2}\\}}\left(\tilde{T}_{r}\right)\;&T\leq
T_{0,r}\\\
\text{lerp}\left(\Lambda_{\text{H}_{2},\\{\text{H},\text{H}_{2}\\}}\left(\tilde{T}_{r}\right),T_{0,r},\tilde{T}_{r}\right)&T_{0,r}\leq
T\leq T_{0}\end{cases}\\\
r_{\alpha}\,r_{m}^{1/4}\,r_{M}^{-5/4}\;\begin{cases}\text{lerp}\left(\Lambda_{\text{H}_{2},\\{\text{H},\text{H}_{2}\\}}\left(\tilde{T}_{v}\right),T_{0,v},\tilde{T}_{v}\right)&T_{0}<T\leq
T_{0,v}\\\
\Lambda_{\text{H}_{2},\\{\text{H},\text{H}_{2}\\}}\left(\tilde{T}_{v}\right)&T_{0,v}\leq
T\end{cases}\end{cases},$ (B3)

the $\Lambda_{\text{H}_{\rm D,2},\text{e}_{\rm D}}$ cooling rate as

$\Lambda_{\text{H}_{\rm D,2},\text{e}_{\rm
D}}=r_{\alpha}^{1}\,r_{M}^{-1}\Lambda_{\text{H}_{2},\text{e}}(\tilde{T}_{r}),$
(B4)

and the $\Lambda_{\text{H}_{\rm D,2},\text{H}_{\rm D}^{+}}$ cooling rate as

$\Lambda_{\text{H}_{\rm D,2},\text{H}_{\rm
D}^{+}}=r_{\alpha}^{1}\,r_{m}^{1/2}\,r_{M}^{-3/2}\Lambda_{\text{H}_{2},\text{H}^{+}}(\tilde{T}_{r}),$
(B5)

Here we have defined $T_{0}$ as $855.833\text{\,}\mathrm{K}$ for $H_{2}-H$ and
$5402.44\text{\,}\mathrm{K}$ for $H_{2}-H_{2}$,
$T_{0,r}=r_{\alpha}^{2}\,r_{m}^{2}\,r_{M}^{-1}\,T_{0}$,
$T_{0,v}=r_{\alpha}^{2}\,r_{m}^{3/2}\,r_{M}^{-1/2}\,T_{0}$, and
$\text{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
$\text{H}_{\rm 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/$\mathrm{eV}$ | Notes |  | Reaction | Energy/$\mathrm{eV}$ | Notes
---|---|---|---|---|---|---
$\text{H}_{\rm D,2}+\text{H}_{\rm D}\rightleftharpoons 3\text{H}_{\rm D}$ | $-4.48$ | 1 |  | $\text{H}_{\rm D}^{-}+H\rightarrow H_{2}+\text{e}_{\rm D}$ | $3.53$ |
$\text{H}_{\rm D,2}+\text{e}_{\rm D}\rightarrow 2\text{H}_{\rm D}+\text{e}_{\rm D}$ | $-4.48$ | 2 |  | $\text{H}_{\rm D,2}^{+}+\text{H}_{\rm D}\rightarrow\text{H}_{\rm D,2}+\text{H}_{\rm D}^{+}$ | $1.83$ |
$2\text{H}_{\rm D,2}\rightleftharpoons\text{H}_{\rm D,2}+2\text{H}_{\rm D}$ | $-4.48$ |  1, 3 |  |  |  |

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
$\text{e}_{\rm D}=r_{\alpha}^{2}\,r_{m}E_{SM}$.

  * 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

$\Lambda_{j}=E_{j}\,k_{j}n(R_{j1})n(R_{j2}),$ (B6)

where $E_{j}$ the energy consumed/produced and $k_{j}$ is the rate for
reaction $j$, and the net rate is just $\Lambda_{\rm
chem}=\sum_{j}\Lambda_{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=\left(1+n_{\rm cr}/n_{\rm tot}\right)$,
however, following Hollenbach & McKee (1979), whose parametric dependence must
still be determined. The critical density is approximated in the Standard
Model rate as

$n_{\rm cr}=\frac{A_{\rm vib}}{\gamma_{20}^{\rm H}x_{\rm H}+\gamma_{10}^{\rm
H_{2}}x_{2}},$ (B7)

with $A_{\rm vib}$ the Einstein $A$ coefficient, $x_{\rm H},x_{2}$ the
relative abundances of $\rm H$ and $\rm H_{2}$, and $\gamma_{\Delta v}^{\rm
H,H_{2}}$ the collisional de-excitation rate coefficients for the $\Delta v$
vibrational transition. From Ryan et al. (2022), these quantities all re-scale
with the dark parameters with the re-scalings described therein, and so the
critical density re-scales as

$\displaystyle n_{\rm cr,D}$
$\displaystyle=r_{\alpha}^{8}\,r_{m}^{19/4}\,r_{M}^{-7/4}n_{\rm
cr}(\tilde{T}_{v})$ (B8) $\displaystyle\tilde{T}_{v}$
$\displaystyle=r_{\alpha}^{-2}\,r_{m}^{-3/2}\,r_{M}^{1/2}\,T.$

Since the network used here from Ryan et al. (2022) only contains a minimal
set of reactions, only the reactions $\text{H}_{\rm D}^{-}+\text{H}_{\rm
D}\rightarrow\text{H}_{\rm D,2}+\text{e}_{\rm D}$, $\text{H}_{\rm
D,2}^{+}+\text{H}_{\rm D}\rightarrow\text{H}_{\rm D,2}+\text{H}_{\rm D}^{+}$,
$\text{H}_{\rm D,2}+2\text{H}_{\rm D}\rightarrow 2\text{H}_{\rm D,2}$, and
$\text{H}_{\rm D,2}+\text{H}_{\rm D}\rightleftharpoons 3\text{H}_{\rm D}$
contribute to the cooling and heating in the simulations presented here. The
cooling and heating reactions are included if either -cooling=DARKMOL or
-heating=DARKMOL option is specified, as the dark analog of the equivalent
options -cooling=CHEM or -heating=CHEM.

###  B.3 $H_{2}$ Dissociation Scaling

From Ryan et al. (2022), the overall rate re-scaling for the collisional
reaction $\text{H}_{2}+\text{H}\rightarrow 3\text{H}$ is given by

$k_{\rm diss,D}(T)=r_{\alpha}^{-1}\,r_{m}^{-3/2}\,r_{M}^{-1/2}\;k_{\rm
diss}(\tilde{T}_{a}),$ (B9)

The Standard Model rate used in both DarkKROME and KROME is based on a master
rate formulation from Martin et al. (1996), with the general form

$\log(k_{\rm
diss}(T))=\log(k_{h}(T))-\frac{\log(k_{h}(T))-\log(k_{l}(T))}{1+(n_{H}/n_{\rm
cr}(T))^{p}},$ (B10)

where $k_{h,l}(T)$ are the high(low)-density limits, $n_{H}$ is the $\rm H$
number density, and $n_{\rm cr}$ is the critical density at which downwards
energy transitions switch from radiative to collisional. Thus, for the dark
rate, we also need to re-scale the critical density using both the overall and
temperature re-scaling from Equation B8.

##  Appendix C Reactions file

DarkKROME includes a basic reactions file, react_dark containing several
example reactions and parameter definitions. The included reactions are listed
in Table 4, while the parameter definitions include setting $\xi=0.01$ and the
dark free-free gaunt factor $g_{ff,D}=1.5$, used in dark bremsstrahlung. The
format of the actual reactions is identical to that of KROME, using the $Q$
notation as described in Section 2.

| Reaction | Source |  |  | Reaction | Source
---|---|---|---|---|---|---
1 | $\text{H}_{\rm D}^{+}+\text{e}_{\rm D}\rightarrow\text{H}_{\rm D}+\gamma_{D}$ | Rosenberg & Fan (2017) |  | 10 | $\text{H}_{\rm D,2}^{+}+\text{H}_{\rm D}\rightarrow\text{H}_{\rm D,2}+\text{H}_{\rm D}^{+}$ | Ryan et al. (2022)
2 | $\text{H}_{\rm D}+\gamma_{D}\rightarrow\text{H}_{\rm D}^{+}+\text{e}_{\rm D}$ | Rosenberg & Fan (2017) |  | 15 | $\text{H}_{\rm D,2}+\text{H}_{\rm D}^{+}\rightarrow\text{H}_{\rm D,2}^{+}+\text{H}_{\rm D}$ | Ryan et al. (2022)
1 | $\text{H}_{\rm D}+\text{e}_{\rm D}\rightarrow\text{H}_{\rm D}^{+}+2\text{e}_{\rm D}$ | Rosenberg & Fan (2017) |  | 18 | $\text{H}_{\rm D,2}+\gamma_{D}\rightarrow\text{H}_{\rm D,2}^{+}+\text{e}_{\rm D}$ | Ryan et al. (2022)
3 | $\text{H}_{\rm D}+\text{e}_{\rm D}\rightarrow\text{H}_{\rm D}^{-}+\gamma_{D}$ | Ryan et al. (2022) |  | 1 | $\text{H}_{\rm D,2}+\text{H}_{\rm D}\rightarrow 3\text{H}_{\rm D}$ |  Ryan et al. (2022)2
4 | $\text{H}_{\rm D}^{-}+\gamma_{D}\rightarrow\text{H}_{\rm D}+\text{e}_{\rm D}$ | Ryan et al. (2022) |  | 1 | $3\text{H}_{\rm D}\rightarrow\text{H}_{\rm D,2}+\text{H}_{\rm D}$ | Ryan et al. (2022)
5 | $\text{H}_{\rm D}^{-}+\text{H}_{\rm D}\rightarrow\text{H}_{\rm D,2}+\text{e}_{\rm D}$ | Ryan et al. (2022) |  | 1 | $\text{H}_{\rm D,2}+2\text{H}_{\rm D}\rightarrow 2\text{H}_{\rm D,2}$ | Ryan et al. (2022)
7 | $\text{H}_{\rm D}^{-}+\text{H}_{\rm D}^{+}\rightarrow 2\text{H}_{\rm D}$ | Ryan et al. (2022) |  | 1 | $2\text{H}_{\rm D}+\text{H}_{\rm D}^{+}\rightarrow\text{H}_{\rm D,2}+\text{H}_{\rm D}^{+}$ | Ryan et al. (2022)
8 | $\text{H}_{\rm D}+\text{H}_{\rm D}^{+}\rightarrow\text{H}_{\rm D,2}^{+}+\gamma_{D}$ | Ryan et al. (2022) |  | 1 | $2\text{H}_{\rm D}+\text{H}_{\rm D}^{+}\rightarrow\text{H}_{\rm D,2}^{+}+\text{H}_{\rm D}$ | Ryan et al. (2022)
9 | $\text{H}_{\rm D,2}^{+}+\gamma_{D}\rightarrow\text{H}_{\rm D}+\text{H}_{\rm D}^{+}$ | Ryan et al. (2022) |  |  |  |

Table 4: Dark particle reactions included in the DarkKROME react_dark file.
Reactions are numbered here according to Galli & Palla (1998).

  * 1

Not included in Galli & Palla (1998).

  * 2

See also Appendix B.3.

## References

  * Abbott et al. (2005) Abbott, B. P., et al. 2005, Phys. Rev. D, 72, 082002, doi: 10.1103/PhysRevD.72.082002
  * Abbott et al. (2008) —. 2008, Phys. Rev. D, 77, 062002, doi: 10.1103/PhysRevD.77.062002
  * Abbott et al. (2018) —. 2018, Phys. Rev. Lett., 121, 231103, doi: 10.1103/PhysRevLett.121.231103
  * Abbott et al. (2019) —. 2019, Phys. Rev. Lett., 123, 161102, doi: 10.1103/PhysRevLett.123.161102
  * Abel et al. (1997) Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New A, 2, 181, doi: 10.1016/s1384-1076(97)00010-9
  * Ackerman et al. (2009) Ackerman, L., Buckley, M. R., Carroll, S. M., & Kamionkowski, M. 2009, Phys. Rev. D, 79, doi: 10.1103/physrevd.79.023519
  * Aghanim et al. (2020) Aghanim, N., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Agrawal et al. (2017) Agrawal, P., Cyr-Racine, F.-Y., Randall, L., & Scholtz, J. 2017, J. Cosmology Astropart. Phys, 2017, 021, doi: 10.1088/1475-7516/2017/08/021
  * Agrawal et al. (2021) Agrawal, P., et al. 2021, EPJC, 81, doi: 10.1140/epjc/s10052-021-09703-7
  * Arkani-Hamed et al. (2016) Arkani-Hamed, N., Cohen, T., D’Agnolo, R. T., et al. 2016, Phys. Rev. Lett., 117, 251801, doi: 10.1103/PhysRevLett.117.251801
  * Berezhiani et al. (2001) Berezhiani, Z., Comelli, D., & Villante, F. L. 2001, PhLB, 503, 362, doi: 10.1016/s0370-2693(01)00217-9
  * Billard et al. (2021) Billard, J., et al. 2021, Direct Detection of Dark Matter – APPEC Committee Report.  https://arxiv.org/abs/2104.07634
  * Bird et al. (2016) Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, Phys. Rev. Lett., 116, 201301, doi: 10.1103/PhysRevLett.116.201301
  * Boddy et al. (2016) Boddy, K. K., Kaplinghat, M., Kwa, A., & Peter, A. H. 2016, Phys. Rev. D, 94, 123017
  * Bramante & Elahi (2015) Bramante, J., & Elahi, F. 2015, Phys. Rev. D, 91, 115001, doi: 10.1103/PhysRevD.91.115001
  * Bramante & Linden (2014) Bramante, J., & Linden, T. 2014, Phys. Rev. Lett., 113, 191301, doi: 10.1103/PhysRevLett.113.191301
  * Bramante et al. (2018) Bramante, J., Linden, T., & Tsai, Y.-D. 2018, Phys. Rev. D, 97, 055016, doi: 10.1103/PhysRevD.97.055016
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Callaway (1994) Callaway, J. 1994, ADNDT, 57, 9, doi: 10.1006/adnd.1994.1009
  * Capelo et al. (2018) Capelo, P. R., Bovino, S., Lupi, A., Schleicher, D. R. G., & Grassi, T. 2018, MNRAS, 475, 3283, doi: 10.1093/mnras/stx3355
  * Cen (1992) Cen, R. 1992, ApJS, 78, 341, doi: 10.1086/191630
  * Chacko et al. (2018) Chacko, Z., Curtin, D., Geller, M., & Tsai, Y. 2018, JHEP, 2018, 163, doi: 10.1007/jhep09(2018)163
  * Cline et al. (2014) Cline, J. M., Liu, Z., Moore, G. D., & Xue, W. 2014, Phys. Rev. D, 90, 015023
  * Collins et al. (2021) Collins, M. L. M., Read, J. I., Ibata, R. A., et al. 2021, MNRAS, 505, 5686, doi: 10.1093/mnras/stab1624
  * Cyr-Racine et al. (2014) Cyr-Racine, F.-Y., de Putter, R., Raccanelli, A., & Sigurdson, K. 2014, Phys. Rev. D, 89, 063517, doi: 10.1103/physrevd.89.063517
  * Cyr-Racine et al. (2021) Cyr-Racine, F.-Y., Ge, F., & Knox, L. 2021, A Symmetry of Cosmological Observables, and a High Hubble Constant as an Indicator of a Mirror World Dark Sector.  https://arxiv.org/abs/2107.13000
  * Cyr-Racine & Sigurdson (2013) Cyr-Racine, F.-Y., & Sigurdson, K. 2013, Phys. Rev. D, 87, doi: 10.1103/physrevd.87.103515
  * D’Amico et al. (2017) D’Amico, G., Panci, P., Lupi, A., Bovino, S., & Silk, J. 2017, MNRAS, doi: 10.1093/mnras/stx2419
  * de Lavallaz & Fairbairn (2010) de Lavallaz, A., & Fairbairn, M. 2010, Phys. Rev. D, 81, 123521, doi: 10.1103/PhysRevD.81.123521
  * De Luca et al. (2020) De Luca, V., Franciolini, G., Pani, P., & Riotto, A. 2020, J. Cosmology Astropart. Phys, 06, 044, doi: 10.1088/1475-7516/2020/06/044
  * de Martino et al. (2020) de Martino, I., Chakrabarty, S. S., Cesare, V., et al. 2020, Univ, 6, 107, doi: 10.3390/universe6080107
  * Diego (2020) Diego, J. M. 2020, Phys. Rev. D, 101, 123512, doi: 10.1103/PhysRevD.101.123512
  * Fan et al. (2013) Fan, J., Katz, A., Randall, L., & Reece, M. 2013, PDU, 2, 139, doi: 10.1016/j.dark.2013.07.001
  * Feng et al. (2009) Feng, J. L., Kaplinghat, M., Tu, H., & Yu, H.-B. 2009, J. Cosmology Astropart. Phys, 2009, 004, doi: 10.1088/1475-7516/2009/07/004
  * Foot & Vagnozzi (2015) Foot, R., & Vagnozzi, S. 2015, Phys. Rev. D, 91, 023512, doi: 10.1103/physrevd.91.023512
  * Foot & Vagnozzi (2016) —. 2016, J. Cosmology Astropart. Phys, 2016, 013, doi: 10.1088/1475-7516/2016/07/013
  * Galli & Palla (1998) Galli, D., & Palla, F. 1998, A&A, 335, 403.  https://arxiv.org/abs/astro-ph/9803315
  * Ghalsasi & McQuinn (2018) Ghalsasi, A., & McQuinn, M. 2018, Phys. Rev. D, 97, 123018, doi: 10.1103/PhysRevD.97.123018
  * Glover (2012) Glover, S. 2012, in The First Galaxies (Springer Berlin Heidelberg), 103–174, doi: 10.1007/978-3-642-32362-1_3
  * Glover (2015) Glover, S. C. O. 2015, MNRAS, 451, 2082, doi: 10.1093/mnras/stv1059
  * Glover & Abel (2008) Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627, doi: 10.1111/j.1365-2966.2008.13224.x
  * Goldberg & Hall (1986) Goldberg, H., & Hall, L. J. 1986, PhLB, 174, 151, doi: 10.1016/0370-2693(86)90731-8
  * Goldman & Nussinov (1989) Goldman, I., & Nussinov, S. 1989, Phys. Rev. D, 40, 3221, doi: 10.1103/PhysRevD.40.3221
  * Gow et al. (2020) Gow, A. D., Byrnes, C. T., Hall, A., & Peacock, J. A. 2020, J. Cosmology Astropart. Phys, 01, 031, doi: 10.1088/1475-7516/2020/01/031
  * Grassi et al. (2014) Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386, doi: 10.1093/mnras/stu114
  * Gresham & Zurek (2019) Gresham, M. I., & Zurek, K. M. 2019, Phys. Rev. D, 99, 083008, doi: 10.1103/PhysRevD.99.083008
  * Gurian et al. (2022) Gurian, J., Jeong, D., Ryan, M., & Shandera, S. 2022, ApJ, 934, 121, doi: 10.3847/1538-4357/ac75e4
  * Hollenbach & McKee (1979) Hollenbach, D., & McKee, C. F. 1979, ApJ, 41, 555, doi: 10.1086/190631
  * Kaplan et al. (2010) Kaplan, D. E., Krnjaic, G. Z., Rehermann, K. R., & Wells, C. M. 2010, J. Cosmology Astropart. Phys, 2010, 021, doi: 10.1088/1475-7516/2010/05/021
  * Kaplan et al. (2011) —. 2011, J. Cosmology Astropart. Phys, 2011, 011, doi: 10.1088/1475-7516/2011/10/011
  * Kim (2001) Kim, Y.-K. 2001, Phys. Rev. A, 64, 032713, doi: 10.1103/physreva.64.032713
  * Kouvaris & Tinyakov (2011) Kouvaris, C., & Tinyakov, P. 2011, Phys. Rev. D, 83, 083512, doi: 10.1103/PhysRevD.83.083512
  * Kouvaris et al. (2018) Kouvaris, C., Tinyakov, P., & Tytgat, M. H. G. 2018, Phys. Rev. Lett., 121, 221102, doi: 10.1103/PhysRevLett.121.221102
  * Latif et al. (2019) Latif, M. A., Lupi, A., Schleicher, D. R. G., et al. 2019, MNRAS, 485, 3352, doi: 10.1093/mnras/stz608
  * Lin (2019) Lin, T. 2019, POS, 333, 009, doi: 10.22323/1.333.0009
  * Martin et al. (1996) Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461, 265, doi: 10.1086/177053
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge New York: Cambridge University Press)
  * Nitz & Wang (2021a) Nitz, A. H., & Wang, Y.-F. 2021a, Phys. Rev. Lett., 126, 021103, doi: 10.1103/PhysRevLett.126.021103
  * Nitz & Wang (2021b) —. 2021b, ApJ, 915, 54, doi: 10.3847/1538-4357/ac01d9
  * Nitz & Wang (2021c) —. 2021c, Phys. Rev. Lett., 127, 151101, doi: 10.1103/PhysRevLett.127.151101
  * Omukai (2000) Omukai, K. 2000, ApJ, 534, 809, doi: 10.1086/308776
  * Peterkops (1977) Peterkops, R. 1977, Theory of ionization of atoms by electron impact (Boulder: Colorado Associated University Press)
  * Petraki & Volkas (2013) Petraki, K., & Volkas, R. R. 2013, IJMPA, 28, 1330028, doi: 10.1142/s0217751x13300287
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Prieto et al. (2015) Prieto, J., Jimenez, R., Haiman, Z., & González, R. E. 2015, MNRAS, 452, 784, doi: 10.1093/mnras/stv1234
  * Randall & Scholtz (2015) Randall, L., & Scholtz, J. 2015, J. Cosmology Astropart. Phys, 2015, 057, doi: 10.1088/1475-7516/2015/09/057
  * Rosenberg & Fan (2017) Rosenberg, E., & Fan, J. 2017, Phys. Rev. D, 96, doi: 10.1103/physrevd.96.123001
  * Ryan et al. (2022) Ryan, M., Gurian, J., Shandera, S., & Jeong, D. 2022, ApJ, 934, 120, doi: 10.3847/1538-4357/ac75ef
  * Rybicki & Lightman (1985) Rybicki, G. B., & Lightman, A. P. 1985, Radiative Processes in Astrophysics (John Wiley & Sons Inc.)
  * Sasaki et al. (2016) Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. 2016, Phys. Rev. Lett., 117, 061101, doi: 10.1103/PhysRevLett.117.061101
  * Schiff (1968) Schiff, L. 1968, Quantum mechanics (New York: McGraw-Hill)
  * Schumann (2019) Schumann, M. 2019, JPhG, 46, 103003, doi: 10.1088/1361-6471/ab2ea5
  * Shandera et al. (2018) Shandera, S., Jeong, D., & Gebhardt, H. S. G. 2018, Phys. Rev. Lett., 120, doi: 10.1103/physrevlett.120.241102
  * Singh et al. (2021) Singh, D., Ryan, M., Magee, R., et al. 2021, Phys. Rev. D, 104, 044015, doi: 10.1103/physrevd.104.044015
  * Sofue & Rubin (2001) Sofue, Y., & Rubin, V. 2001, ARA&A, 39, 137, doi: 10.1146/annurev.astro.39.1.137
  * Strigari (2013) Strigari, L. E. 2013, Phys. Rep., 531, 1, doi: 10.1016/j.physrep.2013.05.004
  * Suazo et al. (2019) Suazo, M., Prieto, J., Escala, A., & Schleicher, D. 2019, ApJ, 885, 127, doi: 10.3847/1538-4357/ab45eb
  * Yoshida et al. (2006) Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6, doi: 10.1086/507978
  * Zurek (2014) Zurek, K. M. 2014, Phys. Rep., 537, 91, doi: 10.1016/j.physrep.2013.12.001