#  Modeling early-universe energy injection with Dense Neural Networks

Yitian Sun<EMAIL_ADDRESS>Center for Theoretical Physics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  The NSF AI Institute for
Artificial Intelligence and Fundamental Interactions  Tracy R. Slatyer
<EMAIL_ADDRESS>Center for Theoretical Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  The NSF AI Institute for Artificial
Intelligence and Fundamental Interactions

###### Abstract

We show that Dense Neural Networks can be used to accurately model the cooling
of high-energy particles in the early universe, in the context of the public
code package DarkHistory. DarkHistory self-consistently computes the
temperature and ionization history of the early universe in the presence of
exotic energy injections, such as might arise from the annihilation or decay
of dark matter. The original version of DarkHistory uses large pre-computed
transfer function tables to evolve photon and electron spectra in redshift
steps, which require a significant amount of memory and storage space. We
present a light version of DarkHistory that makes use of simple Dense Neural
Networks to store and interpolate the transfer functions, which performs well
on small computers without heavy memory or storage usage. This method
anticipates future expansion with additional parametric dependence in the
transfer functions without requiring exponentially larger data tables.

††preprint: MIT-CTP/5449

##  I Introduction

Dark matter (DM) constitutes 84% of the matter content in the universe Aghanim
_et al._ (2020) and plays an important role in the evolution of the early
universe. It has so far eluded detection in all channels other than
gravitational interactions. DM annihilation or decay could inject energy in
the form of Standard Model particles, modifying the temperature and ionization
of the intergalactic medium (IGM) and the anisotropies of the cosmic microwave
background (CMB); studies of these observables have placed strong constraints
on such energy injections (e.g. Adams  _et al._ (1998); Chen and Kamionkowski
(2004); Padmanabhan and Finkbeiner (2005); Zhang  _et al._ (2006); Slatyer
_et al._ (2009); Kanzaki  _et al._ (2010); Slatyer (2013); Galli  _et al._
(2013); Slatyer (2016a, b); Liu  _et al._ (2016); Slatyer and Wu (2017);
Poulin  _et al._ (2017); Liu  _et al._ (2021)).

DarkHistory Liu  _et al._ (2020) is a Python package developed to calculate
the evolution of the IGM temperature and ionization in the early universe in
the presence of such exotic energy injections. For an injected spectrum of
Standard Model (SM) particles, it calculates the particle cascade by computing
(1) the production of photons and electrons/positrons by the decay of the
originally injected SM particles; (2) the subsequent secondary particle
cascade and energy deposition arising from this exotic injection of
photons/electrons/positrons, due to interaction with the IGM and the photon
bath; (3) modifications to the IGM’s temperature and ionization from the
secondary particles and their energy deposition, using a simple Three-Level
Atom (TLA) model.

These calculations are carried out in redshift steps starting prior to
recombination (at redshift $1+z=3000$ by default) and ending well after
reionization near the present day (redshift $1+z=4$). In particular, the
particle cascade in step (2) is evaluated using precomputed transfer
functions, which are matrices that take an input spectrum and output the
spectrum of secondary particles (for a given redshift step). DarkHistory
includes the backreaction effects of changes to the ionization level of
matter, which means the transfer functions themselves are functions of the gas
ionization levels, as well as redshift. In the previous version of
DarkHistory, this dependence is realized by interpolating tables of transfer
function matrices on a grid of values for the hydrogen and helium ionization
fractions, as well as a grid of redshift values.

Figure 1: Model schematics of a Dense Neural Network transfer function. A DNN
takes in (the logarithm of) the input and output photon/electron energy and
necessary physical parameters including the redshift $z$, ionized hydrogen
fraction $x_{\text{HII}}$, and singly ionized helium fraction
$x_{\text{HeII}}$, and outputs (the logarithm of) the transfer function value
$P$. The transfer function matrix acting on a given discretized spectrum is
then obtained by evaluating the DNN on the given energy abscissa.

At around $1.5$ Gb per table with 12 tables around this size, the transfer
functions take up significant storage space as well as memory during runtime,
since they are all loaded in a standard run. They will also be difficult to
scale up to include additional parametric dependence, as the expected size
scales exponentially with the number of added parameters. Intuitively, storing
the transfer functions as tables is an over-representation of the information
content within, since the transfer functions, despite not being smooth
globally, can be divided into multiple regions that are relatively smooth,
each of which could plausibly be fitted with analytical functions. However, in
practice, finding such a solution is quite non-trivial, and even if a solution
was found by _ad hoc_ methods, it would be quite specific to individual
transfer functions and potentially difficult to maintain under changes to the
code. In this work, we present a general solution to these issues by replacing
the transfer function tables with trained Dense Neural Networks.

Recently, machine learning and especially Neural Networks (NNs) has seen many
applications in high energy physics and astrophysics LeCun  _et al._ (2015);
Bourilkov (2019); Radovic  _et al._ (2018); Albertsson  _et al._ (2018). NNs
are often used to capture highly nonlinear and complex relations between
inputs and outputs, and in particular, Dense Neural Networks (DNNs), also
known as fully connected Neural Networks, are a class of general function
approximators given sufficient neuron numbers Goodfellow  _et al._ (2016).

Figure 2: DarkHistory flow chart and transfer functions. Starting with
injected photons and electrons, this flow chart illustrates the work flow of
DarkHistory in the evolution of each redshift step. The boxed quantities are
particle spectra and the arrows indicate transfer functions, which act on the
particle spectra. The orange transfer functions are generated or reconstructed
using DNNs, while the blue ones only use tabulated data. (Modified from Fig. 1
in Ref. Liu  _et al._ (2020).)

In the updated version of DarkHistory we present in this work, we use
lightweight DNNs to store and automatically interpolate DarkHistory’s transfer
functions. A transfer function is essentially a multi-dimensional table, with
2 dimensions corresponding to the input and output particle energy, and the
rest corresponding to physical parameters which in this case will be a subset
of redshift $z$, ionized hydrogen fraction $x_{\text{HII}}\equiv
n_{\text{HII}}/n_{\text{H}}$, and singly ionized helium fraction
$x_{\text{HeII}}\equiv n_{\text{HeII}}/n_{\text{H}}$. (As a matter of
convenience, we define the singly ionized helium fraction with the _hydrogen_
number density as the denominator, so that $x_{\text{HeII}}$ and
$x_{\text{HII}}$ can be easily summed.) The DNN-based transfer functions are
shown in schematic form in Fig. 1. We let a DNN take in all of the relevant
parameters on equal footing and predict (the natural logarithm of) the
transfer function value $P$.

Using DNNs as transfer functions has several benefits:

  * •

The DNNs we use are smaller in size compared to the stored transfer function
tables, by a factor of $\sim$400.

  * •

The computed matter temperature history and ionization history match those
calculated using the transfer function tables to within a few percent relative
difference (with sub-percent relative differences in regions when the species
in question are more than 10% ionized), while the spectral distortion due to
upscattered CMB photons (see Sec. II for precise definition) matches to below
the 10 percent level. These errors are small compared to current experimental
uncertainties.

  * •

We expect the DNNs to have improved scaling in size when additional parameters
are added compared to the original tables. (Including a smooth dependence on
an additional parameter might result in a $\mathcal{O}(1)$ increase in the DNN
neuron number to reach similar accuracy due to increased information content.
On the other hand, adding an additional dimension to a data table would
increase its size by a multiplicative factor of the number of bins in the new
parameter.)

  * •

The DNNs automatically interpolate to any the physical parameter values and
input/output particle energies within the trained range. This allows the use
of flexible binning in DarkHistory, and will also allow us to perform
interpolation on sparse training data Duchi  _et al._ (2011, 2013). The latter
may become necessary in future extensions of DarkHistory, when probing
dependence on an increasing number of physical parameters and generating dense
grids of training data becomes computationally expensive.

  * •

The DNNs predict transfer functions quickly, taking up a similar amount of
time to the rest of the evolution routine for injected particles. Thus
compared to retrieving tabular data from memory on a personal computer, the
use of DNNs results in only a $\mathcal{O}(1)$ increase in total runtime.

  * •

Open source building and training tools for NNs and especially simple
architectures like DNNs are readily available. (In this work we use Tensorflow
2.0 Abadi  _et al._ (2015) with Keras Chollet  _et al._ (2015).)

In Sec. II, we introduce DarkHistory and the roles of transfer functions. In
Sec. III, we detail the training and implementation of the DNN transfer
functions. In Sec. IV we present test runs and discuss the performance of DNN
transfer functions compared to baseline DarkHistory. Finally, in Sec. V we
summarize our results, briefly discuss other possible approaches, and outline
some ideas for future expansion with DNN transfer functions in DarkHistory.

##  II Transfer functions in DarkHistory

To better illustrate the role of transfer functions, we briefly introduce the
procedure followed in DarkHistory and sketched in Fig. 2, which is modified
from a flow chart in Ref. Liu  _et al._ (2020). In Fig. 2, boxed quantities
represent particle spectra, and arrows represent transfer functions, which are
functions acting on spectra. DarkHistory stores the free streaming photon
spectrum, the IGM temperature, and the IGM’s ionized hydrogen fraction and
singly ionized helium fraction at each redshift step (these quantities are
assumed to be homogeneous). For each redshift step, DarkHistory:

  1. 1.

Converts injected SM particles at that redshift to injected photons
$\mathbf{N}^{\gamma}_{\text{inj}}$ and electron/positrons
$\mathbf{N}^{e}_{\text{inj}}$ (hereafter referred collectively as electrons).
In DarkHistory, high energy ($>3$ keV) positrons are treated as electrons
since their dominant energy-loss process, Inverse Compton Scattering (ICS) on
the CMB, does not depend on the particle charge. Lower-energy positrons are
assumed to annihilate and the resulting photon spectrum is tracked; their
kinetic energy is approximated as following the same pattern of energy
deposition as that of the electrons (see Ref. Liu  _et al._ (2020) for a more
in-depth discussion).

  2. 2.

Computes any injected electron spectrum’s energy deposition into ionization,
excitation, or heating by applying the transfer function
$\overline{\mathbf{R}}_{c}$. Computes secondary photon and electron spectra
produced from injected electrons due to ICS, positronium formation and decay,
and atomic processes by applying the ICS transfer function
$\overline{\mathbf{T}}_{\text{ICS}}$ and secondary electron transfer function
$\overline{\mathbf{T}}_{e}$, and evaluating the spectrum of gamma rays
produced from positron annihilation
$\overline{\mathbf{N}}^{\gamma}_{\text{pos}}$. These secondary photon spectra
are added to the spectrum of photons $\mathbf{N}^{\gamma}$ propagated from the
previous step, plus any injected photon spectrum.

  3. 3.

Computes the secondary particles and energy deposition produced by a
propagating photon spectrum $\mathbf{N}^{\gamma}$, due to a variety of
processes, including photon-photon scattering, Compton scattering, pair
production, photoionization, and redshifting. The production of secondary
electrons/positrons in the same redshift step, and their subsequent production
of photons via ICS on the CMB or positron annihilation, are also included. The
propagating photon spectrum for the next redshift step is obtained by applying
the high energy photon transfer function $\overline{\mathbf{P}}_{\gamma}$ to
$\mathbf{N}^{\gamma}$. The low energy photon spectrum, which stores photons
below 3 keV that either photoionize within the redshift step or lie below 13.6
eV, is obtained by applying the low energy photon transfer function
$\overline{\mathbf{D}}_{\gamma}$. The low energy electron spectrum, which
stores electrons with kinetic energy below 3 keV where atomic cooling
dominates ICS and is treated separately in the electron cooling module, is
obtained by applying the low energy electron transfer function
$\overline{\mathbf{D}}_{e}$. Finally, the photon’s energy deposition into
ionization, excitation, and heating is obtained by applying the high energy
deposition transfer function $\overline{\mathbf{D}}_{c}^{\text{high}}$.

  4. 4.

Computes the change to IGM temperature and ionization by first calculating the
energy deposition fraction $f_{c}$’s, from the low-energy electron/photon
spectra and direct energy deposition by higher-energy particles, and then
performing the TLA integration (see Ref. Liu  _et al._ (2020) for details).

In general, transfer functions from an input spectrum to another spectrum
(such as $\overline{\mathbf{T}}_{\text{ICS}}$,
$\overline{\mathbf{P}}_{\gamma}$) take up much more space than those
outputting an energy value (such as $\overline{\mathbf{D}}_{c}^{\text{high}}$)
or those that are diagonal (such as $\overline{\mathbf{D}}_{\gamma}$). For
this version of DarkHistory, we focus on replacing the largest spectral
transfer functions with DNNs, but similar procedures can be applied to the
lower-dimension transfer functions in the future. In the following we
introduce the two major types of transfer functions we will replace in more
detail.

###  II.1 ICS transfer functions

As discussed above, the ICS transfer functions describe the spectrum of
scattered photons from the complete cooling of injected electrons. Unlike the
transfer functions applied to the photon spectrum, the ICS transfer functions
$\overline{\mathbf{T}}_{\text{ICS}}$, $\overline{\mathbf{T}}_{e}$ and
$\overline{\mathbf{R}}_{c}$ are not directly interpolated from tables, but
reconstructed from reference ICS transfer function tables. Sec. III.D and
Appendix A of Ref. Liu  _et al._ (2020) describe in detail how this is
achieved, and we only provide a simple summary here: since an electron quickly
deposits all of its energy in one of our redshift steps, in order to obtain
the total secondary spectrum or energy deposition by an electron, we need to
consider multiple interaction events (via ICS and atomic processes). This can
be done recursively: one can reconstruct the full ICS secondary spectrum and
energy deposition for an electron of energy $E$ knowing the same information
for all electrons with $E^{\prime}<E$. With discretized energy abscissa, the
full ICS-induced photon spectrum and energy output can be solved recursively
from the lowest energy bin. As a result, at each redshift, one can solve for
the full electron ICS transfer functions using transfer functions describing a
single ICS scattering event (as well as functions describing the interaction
rates due to atomic processes).

The transfer functions for a single ICS scattering event on the CMB have
simple (approximate) scaling relations with respect to the temperature of the
CMB $T$, as described in detail in Appendix A of Ref. Liu  _et al._ (2020). As
such one can derive ICS transfer functions at different redshifts from a
single transfer function at a fixed redshift, assuming the CMB is the dominant
radiation background (in DarkHistory, $1+z=400$ is used). These reference
transfer functions are interpolated from tables ics_thomson, ics_rel, and
ics_engloss, corresponding respectively to transfer functions for the
secondary photon spectra of nonrelativistic electrons and relativistic
electrons, and the relativistic electron energy loss in a single ICS
scattering event on the CMB. It is these tables that we will fit with DNNs.

###  II.2 Photon transfer function

The high energy photon transfer function $\overline{\mathbf{P}}_{\gamma}$, low
energy photon transfer function $\overline{\mathbf{D}}_{\gamma}$, and low
energy electron transfer function $\overline{\mathbf{D}}_{e}$ are interpolated
from corresponding tables highengphot, lowengphot, and lowengelec. They in
general depend on the CMB temperature through redshift and matter ionization
levels (ionized hydrogen fraction and singly ionized helium fraction). For
each combination of these physical parameters, lowengphot can be represented
as 1-D arrays (with the one dimension being input/output energy) and is a
factor of 500 smaller than the other transfer functions, so it is at present
not replaced with a DNN. We also found that the numerical calculation of
Compton scattering used in the previous version of DarkHistory was inaccurate
in some parts of parameter space (in particular populating kinematically
forbidden regions), and so we have updated the relevant tabular transfer
functions to ensure sufficient accuracy.111This issue primarily affected
secondary electrons in the 10 eV to 3 keV range produced by Compton scattering
of photons with energies between 100 eV and 10s of keV; the effects on
observable quantities were very small for all cases we checked. The updated
tables on Zenodo incorporate this correction. We now describe some special
features of the photon transfer functions:

Figure 3: Comparison of DNN and tabulated transfer functions. The first two
panels show the high energy photon transfer functions generated at redshift
$1+z=300$, ionized hydrogen fraction $x_{\text{HII}}=0.6$ (and singly ionized
helium fraction $x_{\text{HeII}}=0.0$). The transfer function value in the
region encircled by the black dashed line is negative and the $\log_{10}$ of
the absolute value is shown. The third panel shows the error in
$\log_{10}{\mathbf{\overline{P}}_{\gamma}}$ (defined as
$\Delta\log_{10}{\mathbf{\overline{P}}_{\gamma}}\equiv\log_{10}|{\mathbf{\overline{P}}_{\gamma}}_{\text{DNN}}|-\log_{10}|{\mathbf{\overline{P}}_{\gamma}}_{\text{table}}|$),
which is concentrated on the physical features of the transfer functions. The
fourth panel shows the relative error in the energy transfer function (see
more details in Sec. IV.1). The errors are concentrated on the diagonal that
corresponds to photon propagation and redshifting.

##### Redshift regimes and matter ionization dependence.

All photon transfer functions depend on the redshift, as the photon cooling
processes involve interactions with the redshift-dependent photon background
and/or interstellar medium. For late redshifts ($z<40$) encompassing the epoch
of reionization, the transfer functions are allowed to vary with both the
ionized hydrogen fraction $x_{\text{HII}}$ and the singly ionized helium
fraction $x_{\text{HeII}}$, which can be altered by exotic energy injections.
For redshifts between helium recombination and reionization ($40<z<1600$), the
ionized helium fraction can be safely approximated as zero Liu  _et al._
(2020); the transfer functions are pre-computed assuming no helium ionization,
but can depend on the hydrogen ionized fraction. Before helium recombination
($z>1600$), exotic energy injections consistent with current experimental
bounds have little impact on the thermal equilibrium determining hydrogen and
helium ionization levels Liu  _et al._ (2020), so DarkHistory uses RECFAST
Wong  _et al._ (2008) ionization fractions as a baseline to pre-compute the
transfer functions, which only depends on redshift.

For the DNN implementation, the flexibility of the network allows one DNN to
be trained on the entire redshift range $4<z<3000$ with a learned dependence
on the ionization levels, for each transfer function. However, using different
networks for different redshift regimes gives slightly better accuracy, and
the latter approach is chosen in this version of DarkHistory.

##### Redshift step coarsening and energy conservation.

In the previous version of DarkHistory, photon transfer functions are computed
with redshift step $\Delta\log(1+z)=0.001$. One can choose to increase the
(log) step size to multiples of 0.001 to speed up computation. To combine
multiple redshift steps, DarkHistory pre-composes multiple propagation
transfer functions $\overline{\mathbf{P}}_{\gamma}$ and applies them
appropriately to the deposition transfer functions
$\overline{\mathbf{D}}_{\gamma}$, $\overline{\mathbf{D}}_{e}$, and
$\overline{\mathbf{D}}_{c}^{\text{high}}$. Since the DNN implementation of
transfer functions introduces a small amount of error that can accumulate over
many redshift steps, in order to increase numerical stability, we train the
DNNs to reproduce the pre-composed transfer functions, and require the use of
a fixed log redshift step of $\Delta\log(1+z)=0.012$. 222It was shown in Liu
_et al._ (2020) that an increased redshift step size at this level will lead
to errors smaller than $10^{-3}$, and so will not be the primary contributor
to numerical error in our context. Note also that using a $n$-fold coarsened
log redshift step is not equivalent to using every $n^{\text{th}}$ entry in
the transfer function table. The redshift abscissa used for the table are not
directly related to (and in general are coarser than) the redshift steps used
during a run; the table is interpolated to obtain the transfer functions at
each redshift step. To further decrease numerical error, we store the total
fraction of the injected energy entering each type of secondary particle
spectrum, for an injected photon of any given energy, and use these data
(which are of $\sim 1\%$ the size of the transfer functions) to ensure energy
conservation while accounting for energy loss to redshifting. We discuss the
precise procedure of transfer function pre-composing and imposing energy
conservation in Appendix V.1.

Figure 4: The DNN high energy photon transfer function network’s fitting and
interpolation of the table transfer function values. The panels show the DNN
transfer function’s (redshift regime $4<z<40$) dependence on each of the five
input variables, in a limited region around a central value indicated on the
plot. The shaded region is where both the table value and the DNN prediction
are negative (and the log of the absolute value is shown). A fairly consistent
relative error below $10\%$ is observed.

##  III Transfer Functions from Neural Networks

As described earlier and as indicated in Fig. 1, we update the largest
transfer function tables (ics_thomson, ics_rel, ics_engloss, highengphot, and
lowengelec) to DNN networks that take in input and output particle energies,
redshift, and possibly (depending on redshift) the ionized hydrogen fraction
and singly ionized helium fraction, to produce the transfer function value
$P$. In general, $P$ can vary greatly across many orders of magnitude. (For
example, after recombination the probability for a 10 keV photon to free
stream, losing energy only through redshifting, is substantial, while the
chance of it directly producing secondary photons of 1 keV is very close to 0,
since this outcome is not kinematically allowed in a single Compton scattering
event, nor is the scattered electron produced by Compton scattering or
photoionization able to up-scatter other photons to 1 keV.) As such, we train
the networks to output the natural logarithm of the transfer function value
$\log P$. Similarly, since the input/output energy abscissa and our redshift
steps are also binned in log space by default, the networks also take the log
value of these as inputs. The ionized hydrogen fraction and singly ionized
helium fraction are linearly scaled to match the spread of other parameters
before being fed into the DNN.

For high energy photons, the redshift-coarsened transfer functions (see Sec.
II.2) are fitted. Note also that the transfer function value $P$ can be
negative since by convention the CMB spectrum is subtracted from the transfer
function and the negative values (and the positive values at higher energies)
represent CMB photons being upscattered. In this case $\log|P|$ is predicted
by the DNNs, and the negative value region is recovered in a post-processing
step that identifies local minima of $\log|P|$ (which takes up an negligible
time compared to the rest of the DarkHistory routine). Additionally, the
transfer functions values near the diagonal in the input/output energy
dimensions (corresponding to the free-streaming photons) are adjusted to
enforce energy conservation up to redshifting. For details please refer to
Appendix V.1.

After adjusting hyperparameters, we find that for all transfer functions in
question, it is sufficient to use DNNs with 7 hidden layers with 400 neurons
each, making the number of parameters per DNN $\sim 400^{2}\times(7-1)\sim
9.7\times 10^{5}$, and $2.9\times 10^{6}$ for each transfer function built
from 3 such networks.

Training is done with TensorFlow 2.0 Abadi  _et al._ (2015) and Keras Chollet
_et al._ (2015); Adagrad Lydia and Francis (2019) is used as the optimizer,
with the mean squared error of $\log P$ as the loss function. Each DNN is
trained on 2 V100 GPUs for $\mathcal{O}(10)$ hours or equivalent. For each
epoch, training data is generated by interpolating the multi-dimensional
transfer function table on uniformly random sampled inputs. Since training
data is not reused across epoch, there is no concern of overfitting to a fixed
subset of full data set. Trainings are terminated with the evaluation after
each iteration stops improving significantly. To check for systematic offsets
of the table transfer functions and the NNs, we trained multiple NNs with
random initial values and randomly sampled training data. We found no obvious
systematic offsets between the ensemble of NNs and the tables; e.g. at any
given point the different NNs both underpredicted and overpredicted the table
data.

The codes related to the DNN transfer functions are stored in the nntf module
under DarkHistory. A new example file “Example 12: Using Neural Network
transfer functions.ipynb” is provided to demonstrate using the DNN transfer
functions and comparison with the baseline DarkHistory (which will be
available if the appropriate data tables are present).

##  IV Performance

In this section, we describe the accuracy and speed of generation of the DNN
transfer functions, as well as the accuracy of full runs over a range of
simple DM injection scenarios using the DNN transfer functions.

###  IV.1 Transfer function value prediction

In Fig. 3, a high energy photon transfer function generated by a DNN is
compared against one interpolated from tables. As one can see, the errors in
the raw output i.e. logarithm value of the transfer functions are concentrated
near the distinct physical features in the transfer function, such as at the
output photon energy of $\sim 0.5$ MeV corresponding to positronium decay. The
absolute values in logarithm errors $\Delta\log_{10}|P|$ can be interpreted as
relative errors in $|P|$ (up to a $\ln 10$ factor). In this particular slice
through the high energy photon transfer function, the average
$\Delta\log_{10}|P|$ when $|P|>10^{-20}$ is 0.017, corresponding to a relative
error of $\sim 4$%. (The range of $\log_{10}|P|$ is about $\sim-45$ to 6).
Overall $\Delta\log_{10}|P|$ is comparable to this value, for all DNN transfer
functions. A summary of the errors can be found in Tab. 1.

Figure 5: Example evolution of matter temperature and ionization levels. The
evolution of matter temperature and ionization levels, in an example scenario
where 0.1 GeV DM is decaying into $e^{+}e^{-}$ pairs, as predicted using the
original tabulation (black lines) and using the DNN transfer functions (red
dashed lines). The blue line shows the difference between the two, and the
gray dashed line shows the standard evolution without any external injections.

To see the impact of these errors in a DarkHistory evolution run, it is also
useful to look at errors in energy (per bin) transition rates, besides the
particle number transition rates. The energy transfer function de-emphasize
errors at low energies where many particles can be produced with only a small
fraction of the original particle’s energy, and thus such errors have a small
effect on heating and ionization. Since the standard photon transfer functions
are maps between particle number spectra $N_{i}$, the transfer function values
have the physical meaning of number transition rates. The particle energy
spectrum is related to the number spectrum by

$E_{i}=\mathcal{E}_{i}N_{i}=\mathcal{E}_{i}^{2}\left(\frac{\text{d}N}{\text{d}\mathcal{E}}\right)_{i},$
(1)

where $N_{i}$ is the number of particles in the $i$-th bin, and
$\mathcal{E}_{i}$ its central energy. (Note that the energy bins are log-
spaced.) For a particle number transfer function $P$ with

$N^{\text{out}}_{i}=\sum_{j}N^{\text{in}}_{j}P_{ji},$ (2)

the corresponding energy transfer function $P^{E}$ is defined by

$E^{\text{out}}_{i}=\sum_{j}E^{\text{in}}_{j}P_{ji}^{E}\implies
P_{ji}^{E}=P_{ji}\mathcal{E}_{i}/\mathcal{E}_{j}.$ (3)

In the last panel of Fig. 3, we show the relative error in the high energy
photon transfer function. As expected, the errors are concentrated on the
highest output energy for a given input energy. The relative errors are
generally sub-percent.

Fig. 4 shows that the DNN networks interpolate sensibly between the fixed
abscissa values in the transfer function tables, for the high energy photon
transfer function in the lowest redshift regime ($4<z<40$), as an example. The
errors in the transfer function interpolation are consistent with average
values shown in Tab. 1.

The time it takes to generate a photon transfer function is about or less than
1 second on an 8-CPU personal computer. ICS transfer functions that are only
generated once per run takes only slightly longer at $\lesssim 3$ s. The
accuracy and prediction time for all other transfer functions can be seen in
Tab. 1.

Figure 6: Relative errors in temperature, ionization levels, and low energy
photon distortion across a range of DM mass in a scenario with DM decaying to
$e^{+}e^{-}$. The left panel shows the maximum relative error at any point in
the evolution in $4<z<3000$ of matter temperature, ionized hydrogen fraction,
and singly ionized helium fraction across a range of DM exotic electron
injection energies*. (*Due to its very small absolute value, the relative
error for singly ionized helium fraction $x_{\text{HeII}}$ when
$n_{\text{HeII}}/n_{\text{He}}<10^{-3}$ is not included in this plot. Since
$n_{\text{HeII}}/n_{\text{He}}$ changes rapidly between $\mathcal{O}(1)$ and
$<10^{-3}$, we are essentially only counting its relative error when its value
is order unity.) The relative errors are generally below 5%. The right panel
shows the maximum error over the maximum value in the low energy photon
spectral distortion at $z=0$. The spectral distortion errors at the two
dashed-line values of $m_{\text{DM}}=0.1$ GeV and $10$ GeV are shown in Fig.
7.  Figure 7: Example partial low-energy photon spectral distortion in the
present day. The two panels show two examples of low energy photon spectral
distortion outputs (see definition in Sec. II) from two different runs: 0.1
GeV and 10 GeV DM decaying to $e^{-}e^{+}$ pairs. They represent the extremes
of large and small relative errors for runs over the full range of DM masses
we consider, as shown in Fig. 6. The photon number density (per bin) is
normalized against the baryon number density. The black lines represent
outputs generated using tabulated transfer functions while the red dashed
lines represent that using the DNN transfer functions. The blue line shows the
different between the two. Note that the relatively large errors shown in the
right panel are in part due to the error in the location of the zero.

###  IV.2 Performance over a range of scenarios

Fig. 5 shows the evolution of integrated variables in one particular setting:
0.1 GeV DM particles decaying into electron-positron pairs (this is the same
as the example used in Fig. 4 of Ref. Liu  _et al._ (2020)). In this example,
the error introduced by the DNN in the matter temperature history and hydrogen
ionization history is consistently sub-percent, while the error in singly
ionized helium fraction is sub-percent when
$n_{\text{HeII}}/n_{\text{He}}>10^{-3}$. For a scenario where DM decays into
electron-positron pairs, over a range of DM rest masses, the maximum relative
error for temperature and ionized hydrogen fraction over the entire redshift
range is consistently below 2%, as shown in Fig. 6. Taking into account
injection scenarios with DM decaying to photons, and also undergoing $s$-wave
annihilation into $e^{+}e^{-}$ or photons, the relative error is always below
8%. DarkHistory also computes the partial photon spectral distortion from high
energy photon and electron processes, mostly from ICS of electrons or
positrons on the CMB. This distortion is stored in the low energy photon
spectrum. Note that spectral distortions arising from atomic transitions, due
to photons and electrons below 3 keV, are not included in this spectrum (which
is why we label it as “partial” or “incomplete”). The DNN transfer functions
introduce a small amount of error in this spectral distortion, as shown in
Fig. 7 and Fig. 6. While the shape of the spectral distortion is generally
correct, the error in the location of the distortion zero can cause errors
with a magnitude up to 10% of the distortion magnitude due to the errors in
the distortion zero location. In a future update of DarkHistory, we are
anticipating an update to include the correct treatment of the complete photon
spectral distortion including contributions from atomic transitions. A small
photon spectral error would allow the DNN functions to be used simultaneously
with these updates.

In Appendix V.2, we include the errors in some other exotic injection
scenarios, including DM decaying to photon pairs, and DM annihilating to
photon or $e^{+}e^{-}$ pairs. Although these examples cover only a few simple
injection scenarios, they can serve to test the whole range of DarkHistory’s
dependence on transfer functions, since all exotic energy injections are
converted to either $e^{+}e^{-}$ or photon injections. We expect exotic energy
injections with a more complicated injection spectrum (such as e.g.
annihilation to quarks) to have similar amounts of error associated with using
the DNN transfer function.

Figure 8: Relative errors in physical quantities across a range of DM mass, in
a scenario with DM decaying to a pair of photons. The panels have similar
construction to Fig. 6. Note that the matter temperature and ionization
fractions have max relative errors below 2%, and the errors in the low-energy
photon spectral distortion are also below 10%.

##  V Conclusion

In this work, we have made use of simple Dense Neural Networks to approximate
complex and multi-dimensional transfer functions in DarkHistory to reduce
storage and memory usage, as well as to enable the possibility of adding more
parameter dependence to the transfer function. The DNN transfer functions
achieve good accuracy in computing the evolution history of matter temperature
and ionization, as well as the partial CMB spectral distortion evaluated by
the current version of DarkHistory; typical errors are at the few percent
level, and are comparable to or smaller than estimates of systematic
uncertainties in previous studies of constraints on energy injection Galli
_et al._ (2013); Slatyer (2016a); Liu  _et al._ (2020).

The DNN-based functionality is available in the DarkHistory Github repository
at https://github.com/hongwanliu/DarkHistory, and the necessary data files
(one can choose to download the large tables, or DNN and auxiliary files, or
both) are hosted on Zenodo (see the Github repository for details).

While the use of DNNs offers one solution to this challenge, there may well be
other viable solutions. The information in the DNNs still seems likely to be
an over-representation of the piecewise-smooth transfer functions. We have
briefly explored some alternative methods, including fitting to conventional
functions directly, and with the assistance of symbolic regression techniques
Schmidt and Lipson (2009); however, DNN networks stand out as the best
solution (so far) in terms of fitting accuracy and ease of implementation.
There is work ongoing to expand the capabilities of DarkHistory, and we look
forward to exploring NN-based and alternative techniques in this context.

## Appendix

###  V.1 Redshift step coarsening and energy conservation

In this Appendix, we describe how the photon transfer function change with
coarsened redshift step, (expanding on Sec. III.E.3 of Ref. Liu  _et al._
(2020)), and how energy conservation while correctly accounting for photon
energy loss due to redshift is implemented.

####  V.1.1 Transfer functions without coarsening

In DarkHistory’s main.evolve function, evolution is discretized into log-
normal redshift steps with fixed spacing $d$ (dlnz in code) where the next
redshift $z^{\prime}$ is expressed in terms of $z$ such that:

$\log(1+z^{\prime})=\log(1+z)-d.$ (4)

Following DarkHistory’s flow described in Fig. 2, to obtain the propagating
photon spectrum ${\mathbf{N}_{\text{prop}}^{\gamma}}$, low energy photon
spectrum ${\mathbf{N}_{\text{low}}^{\gamma}}$, low energy electron spectrum
${\mathbf{N}_{\text{low}}^{e}}$ and energy deposition array
${E^{\text{high}}_{c}}$, we retrieve the corresponding transfer functions
${\mathbf{\overline{P}}_{\gamma}}$, ${\mathbf{\overline{D}}_{\gamma}}$,
${\mathbf{\overline{D}}_{e}}$, and ${\mathbf{\overline{D}}^{\text{high}}_{c}}$
at a middle redshift $z_{\text{mid}}$ given by

$\log(1+z_{\text{mid}})=\log(1+z)-d/2.$ (5)

The corresponding transfer functions ${\mathbf{\overline{P}}_{\gamma}}$,
${\mathbf{\overline{D}}_{\gamma}}$, and ${\mathbf{\overline{D}}_{e}}$ all take
in the propagating photon spectrum ${\mathbf{N}_{\text{prop}}^{\gamma}}$ at
redshift $z$ and produce secondary spectra at $z^{\prime}$.
(${\mathbf{\overline{D}}^{\text{high}}_{c}}$ produces energy deposition in
this redshift step.) They are implemented as:

$\displaystyle{\mathbf{N}_{\text{prop}}^{\gamma}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\cdot{\mathbf{\overline{P}}_{\gamma}}(z_{\text{mid}})=\sum_{j}{\mathbf{N}_{\text{prop}}^{\gamma}}_{j}{\mathbf{\overline{P}}_{\gamma}}_{ji}$
(6) $\displaystyle{\mathbf{N}_{\text{low}}^{\gamma}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\cdot{\mathbf{\overline{D}}_{\gamma}}(z_{\text{mid}})$
$\displaystyle{\mathbf{N}_{\text{low}}^{e}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\cdot{\mathbf{\overline{D}}_{e}}(z_{\text{mid}})$
$\displaystyle{E^{\text{high}}_{c}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\cdot{\mathbf{\overline{D}}^{\text{high}}_{c}}(z_{\text{mid}})=\sum_{j}{\mathbf{N}_{\text{prop}}^{\gamma}}_{j}{\mathbf{\overline{D}}^{\text{high}}_{c}}_{j},$

where $i$ and $j$ are indices of discretized energy abscissa.

Figure 9: Relative errors in physical quantities across a range of DM mass in
a scenario where DM undergoes $s$-wave annihilation to a electron/positron
pair. The panels have similar construction to Fig. 6.  Figure 10: Relative
errors in physical quantities across a range of DM mass in a scenario where DM
undergoes $s$-wave annihilation to a photon pair. The panels have similar
construction to Fig. 6.

Energy conservation can be enforced straightforwardly: for an injection in any
photon energy bin with central energy ${\mathcal{E}_{i}}$, the total output
energy on the right hand side of the above equations should add up to
${\mathcal{E}_{i}}$ minus the loss to redshifting of the propagating photons.
Photons below a certain energy $E_{\text{relevant}}$ do not contribute to
redshift energy loss because they either dump all of their energies
efficiently within one redshift step or free stream and no longer interact, in
which case they are stored as an array history of low energy photons
${\mathbf{N}_{\text{low}}^{\gamma}}(z^{\prime})$) and not immediately
redshifted.

For the propagating photons that are redshifted by
${\mathbf{\overline{P}}_{\gamma}}$, the energy lost is approximately:

$\displaystyle|E_{\text{redshift},i}|=\left(1-\frac{1+z^{\prime}}{1+z}\right){\mathcal{E}_{i}}$
$\displaystyle=\frac{z-z^{\prime}}{1+z}\,{\mathcal{E}_{i}}\approx
d\,{\mathcal{E}_{i}}$ (7) $\displaystyle\text{for
}{\mathcal{E}_{i}}>E_{\text{relevant}}.$

Let the energy abscissa (log-central energies of each bin) for photons and
electron be ${\mathcal{E}^{\gamma}}_{i}$ and ${\mathcal{E}^{e}}_{i}$
respectively. We can express the above equation as

$|E_{\text{redshift},i}|=d_{i}{\mathcal{E}_{i}}$ (8)

where $d_{i}=d$ when the photon with energy ${\mathcal{E}_{i}}$ should be
redshifted and $d_{i}=0$ otherwise. (This renders $d_{i}$ a function of
redshift and hydrogen and helium ionization levels in general.)

Then the energy conservation constraint can be written as

${\mathcal{E}_{i}}=\sum_{j}\left({\mathbf{\overline{P}}_{\gamma}}_{ij}{\mathcal{E}^{\gamma}}_{j}+{\mathbf{\overline{D}}_{\gamma}}_{ij}{\mathcal{E}^{\gamma}}_{j}+{\mathbf{\overline{D}}_{e}}_{ij}{\mathcal{E}^{e}}_{j}+{\mathbf{\overline{D}}^{\text{high}}_{c}}_{i}+d_{i}{\mathcal{E}_{i}}\right)$
(9)

This relation is imposed by shifting the near diagonal (propagating) part of
the high energy photon transfer function ${\mathbf{\overline{P}}_{\gamma}}$.
Any energy non-conservation due to numerical errors or errors from
approximating the transfer function as DNNs can be absorbed into this shift.

DNN transfer function | $|\Delta\log_{10}|P||$ | Prediction time | Table size | NN size
---|---|---|---|---
high energy photon (compounded) regime 0 | 0.029 | 1.01 s |  |
high energy photon (compounded) regime 1 | 0.043 | 1.04 s | 4.7Gb | 11.4Mb
high energy photon (compounded) regime 2 | 0.016 | 1.02 s |  |
high energy photon (propagator) regime 0 | 0.029 | 1.02 s |  |
high energy photon (propagator) regime 1 | 0.043 | 1.05 s | 4.7Gb | 11.4Mb
high energy photon (propagator) regime 2 | 0.062 | 1.07 s |  |
low energy electron regime 0 | 0.047 | 0.382 s |  |
low energy electron regime 1 | 0.049 | 0.376 s | 4.7Gb | 11.4Mb
low energy electron regime 2 | 0.040 | 0.370 s |  |
ICS Thomson | 0.00199 | 2.82 s | 0.93Gb | 3.8Mb
ICS relativistic | 0.00410 | 2.12 s | 0.93Gb | 3.8Mb
ICS electron energy loss | 0.00250 | 2.30 s | 0.93Gb | 3.8Mb
Table 1: Error, prediction time, and size of DNN transfer functions used in
DarkHistory. The error of the DNN transfer functions in comparison to the
tabular transfer functions is defined as
$|\Delta\log_{10}|P||\equiv|\log_{10}|P_{\text{DNN}}|-\log_{10}|P_{\text{table}}||$
when $|P|>10^{-20}$. Both the errors and prediction times are evaluated from
random draws in the relevant domains for each DNN until the values stabilize.
The prediction times are evaluated on a 8-CPU personal computer. All DNN
transfer functions use 6 hidden layers of 400 neurons each, but the evaluation
time differs due to the length of the input array corresponding to the size of
the evaluated transfer function matrix ($\sim 10^{5}$ for high energy photon,
$\sim 4\times 10^{4}$ for low energy electron, $\sim 2.5\times 10^{5}$ for
ICS).

####  V.1.2 Coarsening

DarkHistory can enlarge the redshift step $d$ to an multiple of a preset step
$d_{0}$. Let the coarsening multiple be $c$, then the next redshift step
$z^{\prime}$ to $z$ is such that

$\log(1+z^{\prime})=\log(1+z)-c\cdot d_{0}.$ (10)

The photon and electron transfer functions are built with log redshift step
$d_{0}$, so we have to reconstruct the transfer functions for $d$. We first
obtain single-step transfer functions (and $d_{i}$) with current ionization
levels and a redshift value $z_{\text{mid}}$ in the middle of the large
redshift step $d$:

$\log(1+z_{\text{mid}})=\log(1+z)-c\cdot d_{0}/2.$ (11)

Then, we approximate the large redshift step as consisting of $c$ single
redshift steps, each with the same transfer function applied. The compounded
transfer functions can be expressed as

$\displaystyle{\mathbf{N}_{\text{prop}}^{\gamma}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)~{}{\mathbf{\overline{P}}_{\gamma}}^{c}$
(12) $\displaystyle{\mathbf{N}_{\text{low}}^{\gamma}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\left(1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}\right){\mathbf{\overline{D}}_{\gamma}}$
$\displaystyle{\mathbf{N}_{\text{low}}^{e}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\left(1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}\right){\mathbf{\overline{D}}_{e}}$
$\displaystyle{E^{\text{high}}_{c}}(z^{\prime})$
$\displaystyle={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\left(1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}\right){\mathbf{\overline{D}}^{\text{high}}_{c}},$

with all transfer function evaluated at $z_{\text{mid}}$. In the DNN
implementation, the compounded transfer functions
${\mathbf{\overline{P}}_{\gamma}}^{c}$,
$\left(1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}\right)$,
and ${\mathbf{\overline{D}}_{e}}$ are learned as DNN networks, and this step
can be carried out without using the value of
${\mathbf{\overline{P}}_{\gamma}}$ itself. (The low energy photon transfer
function ${\mathbf{\overline{D}}_{\gamma}}$ can be quickly reconstructed using
the CMB energy loss information.)

Imposing energy conservation is similar: At each $d_{0}$ step, the redshift
energy loss is ${\mathbf{N}_{\text{prop}}^{\gamma}}(z+nd_{0})\cdot
d_{0i}{\mathcal{E}_{i}}$, where
${\mathbf{N}_{\text{prop}}^{\gamma}}(z+nd_{0})={\mathbf{N}_{\text{prop}}^{\gamma}}(z)\cdot{\mathbf{\overline{P}}_{\gamma}}^{n}$.
So the total redshift energy loss for a ${\mathcal{E}_{i}}$ photon is

$|E_{\text{redshift},i}|=\left[\left(1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}\right)d_{0}\mathcal{E}\right]_{i}.$
(13)

Let
${\mathbf{\overline{S}}_{\gamma,c}}=1+{\mathbf{\overline{P}}_{\gamma}}+{\mathbf{\overline{P}}_{\gamma}}^{2}+\cdots+{\mathbf{\overline{P}}_{\gamma}}^{c-1}$.
Then the energy conservation condition is

$\displaystyle{\mathcal{E}_{i}}={\mathbf{\overline{P}}_{\gamma}}^{c}_{ij}{\mathcal{E}^{\gamma}}_{j}+{\mathbf{\overline{S}}_{\gamma,c}}\Big{(}{\mathbf{\overline{D}}_{\gamma}}_{ij}{\mathcal{E}^{\gamma}}_{j}+{\mathbf{\overline{D}}_{e}}_{ij}{\mathcal{E}^{e}}_{j}$
$\displaystyle+{\mathbf{\overline{D}}^{\text{high}}_{c}}_{i}+d_{i}{\mathcal{E}_{i}}\Big{)}$
(14) $\displaystyle\forall i\text{, sum on }j.$

Again, the propagating photon spectrum can be adjusted to account for energy
non-conservation from numerical and DNN prediction errors.

###  V.2 Performances in other DM injection scenarios

Performance metrics of DarkHistory using DNN transfer functions are further
demonstrated in Fig. 8 (DM decaying to photon pairs), Fig. 9, (DM annihilating
to electron/positron pairs with $s$-wave cross section), and Fig. 10 (DM
annihilating to photon pairs with $s$-wave cross section).

###  V.3 Performances and prediction time of individual DNN transfer
functions.

The prediction accuracy and prediction time of individual DNNs are summarized
in Tab. 1.

## Acknowledgements

It is a pleasure to thank Hongwan Liu, Wenzer Qin, Gregory Ridgway, Jesse
Thaler, Ken van Tilburg, Kathrin Nippel, Nils Schoeneberg, Anastasia Fialkov,
Stefan Heimersheim, Julian Muñoz, Laura Lopez Honorez, Vivian Poulin, and
Yacine Ali-Haimoud for useful conversations pertaining to this work. YS and
TRS were supported by the U.S. Department of Energy, Office of Science, Office
of High Energy Physics of U.S. Department of Energy under grant Contract
Number DE-SC0012567 through the Center for Theoretical Physics at MIT, and the
National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF
AI Institute for Artificial Intelligence and Fundamental Interactions,
http://iaifi.org/). The computations in this paper were run on MIT LNS’s
Erebus machine and the FASRC Cannon cluster supported by the FAS Division of
Science Research Computing Group at Harvard University. This research made use
of the IPython Perez and Granger (2007), Jupyter Kluyver  _et al._ (2016),
matplotlib Hunter (2007), NumPy Harris  _et al._ (2020), Tensorflow Abadi  _et
al._ (2015), Keras Chollet  _et al._ (2015), and SciPy Virtanen  _et al._
(2020) software packages.

## References

  * Aghanim  _et al._ (2020) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, _et al._ , Astronomy & Astrophysics 641, A6 (2020).
  * Adams  _et al._ (1998) J. A. Adams, S. Sarkar, and D. Sciama, Monthly Notices of the Royal Astronomical Society 301, 210 (1998).
  * Chen and Kamionkowski (2004) X. Chen and M. Kamionkowski, Physical Review D 70, 043502 (2004).
  * Padmanabhan and Finkbeiner (2005) N. Padmanabhan and D. P. Finkbeiner, Physical Review D 72, 023508 (2005).
  * Zhang  _et al._ (2006) L. Zhang, X. Chen, Y.-A. Lei, and Z.-g. Si, Physical Review D 74, 103519 (2006).
  * Slatyer  _et al._ (2009) T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, Physical Review D 80, 043526 (2009).
  * Kanzaki  _et al._ (2010) T. Kanzaki, M. Kawasaki, and K. Nakayama, Progress of Theoretical Physics 123, 853 (2010).
  * Slatyer (2013) T. R. Slatyer, Physical Review D 87, 123513 (2013).
  * Galli  _et al._ (2013) S. Galli, T. R. Slatyer, M. Valdes, and F. Iocco, Physical Review D 88, 063502 (2013).
  * Slatyer (2016a) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016a), arXiv:1506.03811 [hep-ph] .
  * Slatyer (2016b) T. R. Slatyer, Phys. Rev. D 93, 023521 (2016b), arXiv:1506.03812 [astro-ph.CO] .
  * Liu  _et al._ (2016) H. Liu, T. R. Slatyer, and J. Zavala, Physical Review D 94, 063507 (2016).
  * Slatyer and Wu (2017) T. R. Slatyer and C.-L. Wu, Physical Review D 95, 023010 (2017).
  * Poulin  _et al._ (2017) V. Poulin, J. Lesgourgues, and P. D. Serpico, Journal of Cosmology and Astroparticle Physics 2017, 043 (2017).
  * Liu  _et al._ (2021) H. Liu, W. Qin, G. W. Ridgway, and T. R. Slatyer, Physical Review D 104, 043514 (2021).
  * Liu  _et al._ (2020) H. Liu, G. W. Ridgway, and T. R. Slatyer, Phys. Rev. D 101, 023530 (2020), arXiv:1904.09296 [astro-ph.CO] .
  * LeCun  _et al._ (2015) Y. LeCun, Y. Bengio, and G. Hinton, nature 521, 436 (2015).
  * Bourilkov (2019) D. Bourilkov, International Journal of Modern Physics A 34, 1930019 (2019).
  * Radovic  _et al._ (2018) A. Radovic, M. Williams, D. Rousseau, M. Kagan, D. Bonacorsi, A. Himmel, A. Aurisano, K. Terao, and T. Wongjirad, Nature 560, 41 (2018).
  * Albertsson  _et al._ (2018) K. Albertsson, P. Altoe, D. Anderson, M. Andrews, J. P. A. Espinosa, A. Aurisano, L. Basara, A. Bevan, W. Bhimji, D. Bonacorsi, _et al._ , in  _Journal of Physics: Conference Series_ , Vol. 1085 (IOP Publishing, 2018) p. 022008.
  * Goodfellow  _et al._ (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Duchi  _et al._ (2011) J. Duchi, E. Hazan, and Y. Singer, Journal of machine learning research 12 (2011).
  * Duchi  _et al._ (2013) J. C. Duchi, M. I. Jordan, and H. B. McMahan, (2013).
  * Abadi  _et al._ (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,”  (2015).
  * Wong  _et al._ (2008) W. Y. Wong, A. Moss, and D. Scott, Monthly Notices of the Royal Astronomical Society 386, 1023 (2008).
  * Lydia and Francis (2019) A. Lydia and S. Francis, Int. J. Inf. Comput. Sci 6, 566 (2019).
  * Schmidt and Lipson (2009) M. Schmidt and H. Lipson, science 324, 81 (2009).
  * Perez and Granger (2007) F. Perez and B. E. Granger, Computing in Science and Engineering 9, 21 (2007).
  * Kluyver  _et al._ (2016) T. Kluyver _et al._ , in  _ELPUB_ (2016).
  * Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
  * Harris  _et al._ (2020) C. R. Harris _et al._ , Nature 585, 357 (2020).
  * Virtanen  _et al._ (2020) P. Virtanen _et al._ , Nature Methods (2020), https://doi.org/10.1038/s41592-019-0686-2.


otherwise numerical instability and round-off errors can occur. Using an
arbitrarily large M value also expands the feasible region of the LP
relaxation unnecessarily and results in longer runtimes, see [10]. We decided
to set `int bigM=k*100` as smaller or larger values led to incorrect values
for some optimization criteria.

The following model $\mathbf{M}_{\mathbf{B}}$ for the Cplex MIP solver uses an
upper triangular matrix of booleans that indicates if a cavity is wired before
another:

    Ψdvar boolean lt[Cavities][Cavities];
    Ψ

If $lt[i,j]=1$, then pfc[i] $<$ pfc[j], otherwise pfc[i] $>$ pfc[j]. Only the
upper triangular matrix is defined, i.e., values are only defined for $i<j$.
To find out if pfc[i] $<$ pfc[j] when $i>j$, the value of $lt[j,i]$ is
determined and inverted. The diagonal of the $lt$ matrix is undefined, because
a cavity cannot be placed before or after itself in the permutation. A
solution can be directly extracted from the $lt$ matrix as the subsequent
example explains.

$lt=\begin{bmatrix}&1&1&1\\\ &&0&1\\\ &&&1\\\ &&&\end{bmatrix}$

$pfc=1,3,2,4$ because

pfc[1]$<$pfc[2], pfc[1]$<$pfc[3], pfc[1]$<$pfc[4], pfc[2]$>$pfc[3],
pfc[2]$<$pfc[4], pfc[3]$<$pfc[4]

Atomic constraints are expressed as entries in the $lt$ matrix:

    Ψforall(a in AtomicConstraints)
    Ψif (a.cbefore < a.cafter) lt[a.cbefore,a.cafter] == 1;
    Ψelse lt[a.cafter,a.cbefore] == 0;
    Ψ

The allDifferent constraint is expressed as inequalities over $lt$ matrix
values stating that either pfc[i] or pfc[j] has to be larger than the other.

    Ψforall(ordered i,j in Cavities)
    Ψpfc[i] - pfc[j] + 1 <= bigM * (1 - lt[i,j]);
    Ψforall(ordered i,j in Cavities)
    Ψpfc[j] - pfc[i] + 1 <= bigM * lt[i,j];
    Ψ

For example, if $lt[i,j]$ is equal to $0$ and therefore pfc[i] $>$ pfc[j], the
first constraint is always satisfied if the big-M constant is chosen big
enough. The right term in the second constraint evaluates to $0$ and the left
term has to be smaller or equal to $0$, which is only possible if the two
values are different. Disjunctive constraints are formulated as inequality
constraints. A satisfied disjunctive constraint is larger or equal to 1,
because we sum up values for each atomic constraint in the disjunction from
the $lt$ matrix. As only the upper triangular matrix of the $lt$ matrix is
defined, case distinctions based on cavity indices have to be made and some
values have to be inverted. With a fully populated matrix $lt$, each
disjunctive constraint is expressed by the following conditional inequations.
Note that for each disjunctive constraint, only one of the if-clauses is
satisfied.

    Ψforall(d in DisjunctiveConstraints){
    Ψif (d.c1before < d.c1after && d.c2before < d.c2after)
    Ψ{lt[d.c1before,d.c1after] + lt[d.c2before,d.c2after]
    Ψ>= 1;}
    Ψif (d.c1before < d.c1after && d.c2before > d.c2after)
    Ψ{lt[d.c1before,d.c1after] + 1-lt[d.c2after,d.c2before]
    Ψ>= 1;}
    Ψif (d.c1before > d.c1after && d.c2before < d.c2after)
    Ψ{1-lt[d.c1after,d.c1before] + lt[d.c2before,d.c2after]
    Ψ>= 1;}
    Ψif (d.c1before > d.c1after && d.c2before > d.c2after)
    Ψ{1-lt[d.c1after,d.c1before] + 1-lt[d.c2after,d.c2before]
    Ψ>= 1;}}
    Ψ

Direct successor constraints also make a case distinction based on cavity
indices. If $i$ is in Directsuccessors, then $i+b$ must follow immediately in
the permutation, which means that the left-hand side of the condition must be
equal to 0 if and only if pfc[i] $<$ pfc[i+b].

    Ψforall(i in Cablestarts:
    Ψi in DirectSuccessors || (i+b) in DirectSuccessors){
    Ψif (i in DirectSuccessors)
    Ψpfc[i+b] - pfc[i] - 1 <= bigM * (1 - lt[i,i+b]);
    Ψif ((i+b) in DirectSuccessors)
    Ψpfc[i] - pfc[i+b] - 1 <= bigM * lt[i,i+b];}Ψ
    Ψ

Minimum and maximum positions of cavities in the permutation are stored in
arrays minimum and maximum. Their values are set by constraints.

    Ψdvar float+ minimum[Cablestarts] in Positions;
    Ψdvar float+ maximum[Cablestarts] in Positions;
    Ψ
    Ψforall(i in Cablestarts) {
    Ψminimum[i] - pfc[i]   <= 0;                        (1)
    Ψminimum[i] - pfc[i+b] <= 0;                        (2)
    Ψpfc[i] - minimum[i]   <= bigM * (1-lt[i,i+b]);     (3)
    Ψpfc[i+b] - minimum[i] <= bigM * lt[i,i+b];         (4)
    Ψmaximum[i] == pfc[i]+pfc[i+b] - minimum[i];}       (5)
    Ψ

Constraints (1) and (2) ensure that the value of minimum[i] is larger or equal
to min(pfc[i],pfc[i+b]). Constraints (3) and (4) ensure that the value of
minimum[i] is smaller or equal to min(pfc[i],pfc[i+b]), depending on which of
the cavities in a job pair occurs first in the permutation. If constraints (1)
to (4) are satisfied, the value of minimum[i] is set correctly. Calculating
the value of maximum[i] is then trivial using constraint (5).

If a job pair is interrupted, its entry in the array cableIsStored is set to
1, which happens if the difference between the minimum and maximum position of
cavities from a job pair in the permutation is larger or equal than 2.

    Ψdvar boolean cableIsStored[Cablestarts];
    Ψ
    Ψforall(p in Cablestarts) {
    Ψ2 - maximum[p] + minimum[p] <= bigM * (1-cableIsStored[p]);
    Ψmaximum[p] - minimum[p] - 1 <= bigM * cableIsStored[p];}
    Ψ

The variable cableIsStoredAtPosition[i,t] is equal to 1 if one cavity from the
job pair $i$ occurs in the permutation before the position $t$, but the other
not. Otherwise cableIsStoredAtPosition[i,t] is equal to 0. The variable
cableIsPluggedBefore[i,t] is equal to 1 if both cavities of a job pair $i$
occur in the permutation before position $t$. The variable
cableIsPluggedAfter[i,t] is equal to 1 if both cavities of a job pair $i$
occur in the permutation after position $t$.

    Ψdvar boolean cableIsStoredAtPosition[Cablestarts,Positions];
    Ψdvar boolean cableIsPluggedBefore[Cablestarts,Positions];
    Ψdvar boolean cableIsPluggedAfter[Cablestarts,Positions];
    Ψ

Because a cable can either be completely inserted at one position in the
permutation or only be inserted with one cable end or not be inserted at all,
the following constraints hold:

    Ψforall(p in Cablestarts, t in Positions) {
    ΨcableIsStoredAtPosition[p,t] + cableIsPluggedBefore[p,t]
    Ψ+ cableIsPluggedAfter[p,t] == 1;
    Ψminimum[p]-t+1 <= bigM*(1-cableIsStoredAtPosition[p,t]);
    Ψt-maximum[p]+1 <= bigM*(1-cableIsStoredAtPosition[p,t]);
    Ψt-minimum[p]   <= bigM*(1-cableIsPluggedBefore[p,t]);
    Ψmaximum[p]-t   <= bigM*(1-cableIsPluggedAfter[p,t]);}
    Ψ

The maximum number of cables $M$ that are simultaneously contained in storage
is set implicitly by a constraint. $M$ is only set correctly, because it is
part of the weighted sum of the four optimization criteria. If not being part
of the weighted sum, the value of $M$ would be set arbitrarily large, such
that it just satisfies the constraint. Similarly, the value $L$ is also set
implicitly by a constraint, whereas the calculation of $N$ and $S$ is
straightforward:

    Ψforall(t in Positions) sum(p in Cablestarts)
    ΨcableIsStoredAtPosition[p,t] <= M;
    Ψ
    Ψforall(p in Cablestarts) maximum[p] - minimum[p] - 1 <= L;
    Ψ
    ΨN == sum(s in SoftAtomicConstraints)
    Ψ(pfc[s.cafter] - pfc[s.cbefore] <= 0);
    Ψ
    ΨS == sum(p in Cablestarts) cableIsStored[p];
    Ψ

The objective function is identical to the one used in the
$\mathbf{M}_{\mathbf{C}}$ model.

    Ψ  ΨΨ  Ψ
    Ψminimize S * pow(k, 3)Ψ+ M * pow(k, 2) + L * pow(k, 1) +  N;
    Ψ

For the Gurobi solver, we manually implemented the models
$\mathbf{M}_{\mathbf{I}}$, $\mathbf{M}_{\mathbf{B}}$, and
$\mathbf{M}_{\mathbf{B}}$ in the C$\sharp$ API of this solver.