††thanks<EMAIL_ADDRESS>https://orcid.org/0009-0009-0290-5755††thanks<EMAIL_ADDRESS>http://orcid.org/0000-0001-8837-4127

# Kaluza-Klein Graviton Freeze-In and Big Bang Nucleosynthesis

Mathieu Gross1 Dan Hooper2,3,4 1Université Paris-Saclay, CNRS/IN2P3, IJCLab,
91405 Orsay, France  2University of Chicago, Department of Astronomy and
Astrophysics, Chicago, Illinois 60637, USA  3University of Chicago, Kavli
Institute for Cosmological Physics, Chicago, Illinois 60637, USA  4Fermi
National Accelerator Laboratory, Theoretical Astrophysics Group, Batavia,
Illinois 60510, USA

###### Abstract

In models featuring extra spatial dimensions, particle collisions in the early
universe can produce Kaluza-Klein gravitons. Such particles will later decay,
potentially impacting the process of Big Bang nucleosynthesis. In this paper,
we consider scenarios in which gravity is free to propagate throughout $n$
flat, compactified extra dimensions, while the fields of the Standard Model
are confined to a 3+1 dimensional brane. We calculate the production and decay
rates of the states that make up the Kaluza-Klein graviton tower and determine
the evolution of their abundances in the early universe. We then go on to
evaluate the impact of these decays on the resulting light element abundances.
We identify significant regions of previously unexplored parameter space that
are inconsistent with measurements of the primordial helium and deuterium
abundances. In particular, we find that for the case of one extra dimension
(two extra dimensions), the fundamental scale of gravity must be
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}2\times 10^{13}\,{\rm GeV}$
($M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}10^{10}\,{\rm GeV}$) unless the temperature of the early universe
was never greater than $T\sim 2\,{\rm TeV}$ ($T\sim 1\,{\rm GeV}$). For larger
values of $n$, these constraints are less stringent. For the case of $n=6$,
for example, our analysis excludes all values of $M_{\star}$ less than $\sim
10^{6}\,{\rm GeV}$, unless the temperature of the universe was never greater
than $T\sim 3\,{\rm TeV}$. The results presented here severely limit the
possibility that black holes were efficiently produced through particle
collisions in the early universe’s thermal bath.

††preprint: FERMILAB-PUB-24-0354-T

##  I Introduction

Measurements of the primordial light element abundances provide us with our
earliest probe of cosmic history, allowing us to constrain the expansion rate
and overall composition of our universe as early as $\sim 1\,{\rm s}$ after
the Big Bang. In particular, this information indicates that the early
universe was radiation dominated and at least as hot as a few MeV [1, 2, 3, 4,
5]. Little is known, however, about the thermal history of our universe prior
to the onset of Big Bang nucleosynthesis.

At extremely high temperatures, particles in the thermal plasma can scatter to
produce gravitational excitations, leading to the production of a stochastic
gravitational wave background [6, 7]. Such considerations allow us to
constrain the maximum temperature of the very early universe, $T_{\rm
max}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm] $\sim$}}~{}M_{\rm
Pl}\sim 10^{19}\,{\rm GeV}$ [8, 9]. If those gravitational waves were later
diluted, such as through inflation, this constraint would instead apply to the
temperature of subsequent reheating, $T_{\rm
RH}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm] $\sim$}}~{}M_{\rm
Pl}\sim 10^{19}\,{\rm GeV}$.

If our universe has extra spatial dimensions, gravitational excitations can be
be produced more efficiently and at lower temperatures. As a result, we can
potentially place much more stringent constraints on the maximum temperature
of the universe at early times in such scenarios. In this study, we will focus
on the model proposed by Arkani-Hamed, Dimopoulos, and Dvali (ADD) [10, 11,
12], which features $n$ extra dimensions that are flat and compactified on a
torus of radius, $R$. Unlike gravity, all of the Standard Model fields are
confined to a $3$-dimensional brane, the volume of which constitutes the $3+1$
dimensional spacetime that we experience.

This class of models was originally proposed as a possible solution to the
electroweak hierarchy problem [10]. In particular, the effective 4-dimensional
Planck scale, $M_{\rm Pl}\approx 1.22\times 10^{19}\,{\rm GeV}$, is related to
the fundamental $n+4$ dimensional Planck scale, $M_{\star}$, according to the
following [10, 12]:

$\displaystyle M^{2}_{\rm Pl}=(2\pi)^{n}\,R^{n}\,M_{\star}^{2+n}.$ (1)

Thus, for the appropriate values of $R$ and $n$, the fundamental scale of
gravity could be similar to electroweak scale, $M_{\star}\sim{\rm TeV}$. In
such a scenario, the apparent hierarchy between the Planck scale and the
electroweak scale would be a consequence of the Standard Model’s localization
on the $3+1$ dimensional brane.

Since the ADD model was proposed more than two decades ago, stringent
constraints have been placed on this class of scenarios. In particular, data
from the Large Hadron Collider require that $M_{\star}$ must be greater than
several TeV [13, 14]. Tests of the gravitational force law at sub-millimeter
distances further constrain $R\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\
[-0.07cm] $\sim$}}~{}30\,\mu{\rm m}$ [15, 16, 17], corresponding to
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}5.4\times 10^{8}\,{\rm GeV}$ for $n=1$ and $M_{\star}>3.6\times
10^{3}\,{\rm GeV}$ for $n=2$. The requirement that neutron stars are not
overly heated by Kaluza-Klein graviton decays further requires
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}1.7\times 10^{5}\,{\rm GeV}$ for $n=2$ and
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}7.6\times 10^{4}\,{\rm GeV}$ for $n=3$ [18] (for further
discussion, see the entry on Extra Dimensions in the Particle Data Group’s
Review of Particle Physics [19]).

In light of the constraints from the Large Hadron Collider, we will not
attempt to motivate this study by appealing to the electroweak hierarchy
problem, but rather by the broader possibility of extra spatial dimensions,
such as within the context of string theory [11, 20]. With this in mind, we
will consider values of $M_{\star}$ that range from several TeV up to the
Planck scale.

In this class of models, particle collisions in the early universe can result
in the efficient production of Kaluza-Klein gravitons. For gravitons lighter
than $m_{\rm KK}\sim 10^{6}\,{\rm GeV}$, such states will decay during or
after the era of Big Bang Nucleosynthesis, producing energetic Standard Model
particles that can break apart helium nuclei through the processes of
photodissassociation or hadrodissociation. Such decays can reduce of the
primordial helium abundance and increase the abundance of primordial
deuterium. We explore the impact of these decays and use measurements of the
light element abundances to place constraints on this class of scenarios.

The remainder of this paper is organized as follows. In Sec. II, we evaluate
the evolution of the Kaluza-Klein graviton abundance in the early universe,
including their production via freeze-in and their subsequent decays. In Secs.
III and IV, we discuss the impact of these decays on the primordial light
element abundances and use this information to place constraints on this class
of models. In Sec. V, we consider the impact of these constraints on the
possibility that black holes could be efficiently produced through particle
collisions in the early universe. We summarize our main results in Sec. VI.

##  II Kaluza-Klein Graviton Freeze In

In the ADD scenario, all of the Standard Model fields are restricted to
propagate within the $3+1$ dimensional brane. In constrast, gravitons are free
to propagate throughout the $n+4$ dimensional bulk. To observers on the brane,
these massless spin-$2$ gravitons appear as massive Kaluza-Klein states. More
specifically, for each level of the Kaluza-Klein graviton tower, there exists
one spin-$2$ state, $\tilde{h}_{m}$, $(n-1)$ spin-$1$ states, and $n(n-1)/2$
spin-$0$ states, $\tilde{\phi}_{m}$, all with masses given by
$m_{\tilde{h}_{m}}=m_{\tilde{\phi}_{m}}=m/R$, where $m$ is the level of the
Kaluza-Klein tower. The spin-$1$ states are entirely decoupled and will play
no role in the calculations performed here [21].

The spin-2 Kaluza-Klein gravitons decay to Standard Model fields with the
following partial widths [21]:

$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow\gamma\gamma}$
$\displaystyle=\frac{m^{3}_{\tilde{h}_{m}}}{80\pi M^{2}_{\rm Pl}}$ (2)
$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow ZZ}$
$\displaystyle=\frac{13m^{3}_{\tilde{h}_{m}}}{960\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{Z}^{2}}{m^{2}_{\tilde{h}_{m}}}\bigg{)}^{1/2}\bigg{(}1+\frac{56m_{Z}}{169m_{\tilde{h}_{m}}}+\frac{36m_{Z}^{2}}{169m^{2}_{\tilde{h}_{m}}}\bigg{)}$
$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow WW}$
$\displaystyle=\frac{13m^{3}_{\tilde{h}_{m}}}{480\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{W}^{2}}{m^{2}_{\tilde{h}_{m}}}\bigg{)}^{1/2}\bigg{(}1+\frac{56m_{W}}{169m_{\tilde{h}_{m}}}+\frac{36m_{W}^{2}}{169m^{2}_{\tilde{h}_{m}}}\bigg{)}$
$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow gg}$
$\displaystyle=\frac{m^{3}_{\tilde{h}_{m}}}{10\pi M^{2}_{\rm Pl}}$
$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow HH}$
$\displaystyle=\frac{m^{3}_{\tilde{h}_{m}}}{960\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{H}^{2}}{m^{2}_{\tilde{h}_{m}}}\bigg{)}^{5/2}$
$\displaystyle\Gamma_{\tilde{h}_{m}\rightarrow f\bar{f}}$
$\displaystyle=\frac{g_{f}\,m^{3}_{\tilde{h}_{m}}}{640\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{f}^{2}}{m^{2}_{\tilde{h}_{m}}}\bigg{)}^{3/2}\bigg{(}1+\frac{8m_{f}}{3m_{\tilde{h}_{m}}}\bigg{)},$

where $g_{f}$ is the number of internal degrees-of-freedom of the fermion
species (4 for each charged lepton, 2 for each neutrino, and 12 for each
quark).

For $m_{\tilde{h}_{m}}\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm]
$\sim$}}~{}{\rm TeV}$, the sum of these partial widths is given by111Note that
this quantity is sometimes calculated including decays into right-handed
neutrinos, in which case the numerical prefactor is instead given by
$292/960\pi=73/240\pi$.

$\displaystyle\Gamma_{\tilde{h}_{m}}\approx\frac{283m^{3}_{\tilde{h}_{m}}}{960\pi
M^{2}_{\rm Pl}}.$ (3)

Figure 1: The lifetimes of Kaluza-Klein gravitons as a function of their mass,
for the case of $M_{\star}=10^{6}\,{\rm GeV}$ and $n=6$ (for which the
lightest Kaluza-Klein states have a mass of $1/R\approx 273\,{\rm GeV}$). The
black dots represent the spin-2 states ($\tilde{h}_{m}$), while the red
squares correspond to the spin-0 states ($\tilde{\phi}_{m}$). We expect those
Kaluza-Klein modes with $\tau\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\
[-0.07cm] $\sim$}}~{}1\,{s}$ to decay after the onset of Big Bang
nucleosynthesis, potentially impacting the primordial light element
abundances.

The partial widths of the spin-0 Kaluza-Klein graviton states are given by
[21]:

$\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow\gamma\gamma}$
$\displaystyle=0$ (4) $\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow ZZ}$
$\displaystyle=\frac{m^{3}_{\tilde{\phi}_{m}}}{(n+2)48\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{Z}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}^{1/2}\bigg{(}1-\frac{4m_{Z}}{m_{\tilde{\phi}_{m}}}+\frac{12m_{Z}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}$
$\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow WW}$
$\displaystyle=\frac{m^{3}_{\tilde{\phi}_{m}}}{(n+2)24\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{W}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}^{1/2}\bigg{(}1-\frac{4m_{W}}{m_{\tilde{\phi}_{m}}}+\frac{12m_{W}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}$
$\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow gg}$ $\displaystyle=0$
$\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow HH}$
$\displaystyle=\frac{m^{3}_{\tilde{\phi}_{m}}}{(n+2)48\pi M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{H}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}^{1/2}\bigg{(}1+\frac{2m_{H}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}^{2}$
$\displaystyle\Gamma_{\tilde{\phi}_{m}\rightarrow f\bar{f}}$
$\displaystyle=\frac{g_{f}\,m^{2}_{f}m_{\tilde{\phi}_{m}}}{(n+2)48\pi
M^{2}_{\rm
Pl}}\bigg{(}1-\frac{4m_{f}^{2}}{m^{2}_{\tilde{\phi}_{m}}}\bigg{)}^{1/2}\bigg{(}1-\frac{2m_{f}}{m_{\tilde{\phi}_{m}}}\bigg{)},$

which for $m_{\tilde{\phi}_{m}}\gg{\rm TeV}$ sums to

$\displaystyle\Gamma_{\tilde{\phi}_{m}}\approx\frac{m^{3}_{\tilde{\phi}_{m}}}{(n+2)12\pi
M^{2}_{\rm Pl}}.$ (5)

In Fig. 1, we plot of the lifetimes of the Kaluza-Klein gravitons as a
function of their mass, for the case of $M_{\star}=10^{6}\,{\rm GeV}$ and
$n=6$.

Figure 2: The production rate of Kaluza-Klein gravitons divided by three
powers of temperature, $P/T^{3}$, for the case of $M_{\star}=10^{6}\,{\rm
GeV}$ and $n=6$. We show results for Kaluza-Klein gravitons with levels of
$m=10^{4},10^{3},10^{2},10$, and 1. The solid red lines correspond to the
spin-2 states, $\tilde{h}_{m}$, while the dashed red lines correspond to the
scalar states, $\tilde{\phi}_{m}$. For comparison, we also show the value of
the hubble rate, $H$, which is much larger than the production rate of Kaluza-
Klein gravitons, ensuring that equilibrium is never achieved.

Kaluza-Klein gravitons can be produced in the early universe through the
inverse decays of Standard Model particles, such as
$e^{+}e^{-}\rightarrow\tilde{h}_{m}$, for example. It follows from the
principle of detailed balance that, in equilibrium, the production rate of a
given particle species will be equal to the rate of its destruction. We can
use this information to calculate the production rates of the various Kaluza-
Klein graviton modes from their decay widths [22, 23, 24]:

$\displaystyle P_{\tilde{h}_{m}}$ $\displaystyle=n^{\rm
Eq}_{\tilde{h}_{m}}\,\langle\Gamma_{\tilde{h}_{m}}\rangle$ (6)
$\displaystyle\approx n^{\rm
Eq}_{\tilde{h}_{m}}\,\Gamma_{\tilde{h}_{m}}\,\frac{K_{1}(m_{\tilde{h}_{m}}/T)}{K_{2}(m_{\tilde{h}_{m}}/T)},$
$\displaystyle P_{\tilde{\phi}_{m}}$ $\displaystyle=n^{\rm
Eq}_{\tilde{\phi}_{m}}\,\langle\Gamma_{\tilde{\phi}_{m}}\rangle$
$\displaystyle\approx n^{\rm
Eq}_{\tilde{\phi}_{m}}\,\Gamma_{\tilde{\phi}_{m}}\,\frac{K_{1}(m_{\tilde{\phi}_{m}}/T)}{K_{2}(m_{\tilde{\phi}_{m}}/T)},$

where $\langle\Gamma_{\tilde{h}_{m}}\rangle$ and
$\langle\Gamma_{\tilde{\phi}_{m}}\rangle$ are the thermally averaged decay
widths [25], $T$ is the temperature of the decaying particle population, and
$K_{1}$ and $K_{2}$ are modified Bessel functions of the second kind. The
equilibrium number densities of Kaluza-Klein gravitons are given by

$\displaystyle n^{\rm Eq}_{\tilde{h}_{m}}$
$\displaystyle=\frac{5}{2\pi^{2}}\int_{m_{\tilde{h}_{m}}}^{\infty}\frac{(E^{2}-m^{2}_{\tilde{h}_{m}})^{1/2}}{e^{E/T}-1}EdE$
(7) $\displaystyle n^{\rm Eq}_{\tilde{\phi}_{m}}$
$\displaystyle=\frac{n(n-1)}{4\pi^{2}}\int_{m_{\tilde{\phi}_{m}}}^{\infty}\frac{(E^{2}-m^{2}_{\tilde{\phi}_{m}})^{1/2}}{e^{E/T}-1}EdE,$

where the factor of $n(n-1)/2$ in the second expression accounts for the
multiplicity of spin-0 states at each level of the Kaluza-Klein tower.

For $T\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm] $\sim$}}~{}{\rm
TeV}$, these production rates reduce to

$\displaystyle P_{\tilde{h}_{m}}$
$\displaystyle\approx\frac{5\zeta(3)T^{3}}{\pi^{2}}\,\frac{283m^{3}_{\tilde{h}_{m}}}{960\pi
M^{2}_{\rm
Pl}}\,\frac{K_{1}(m_{\tilde{h}_{m}}/T)}{K_{2}(m_{\tilde{h}_{m}}/T)}$ (8)
$\displaystyle P_{\tilde{\phi}_{m}}$
$\displaystyle\approx\frac{n(n-1)\zeta(3)T^{3}}{2\pi^{2}}\,\frac{m^{3}_{\tilde{\phi}_{m}}}{(n+2)12\pi
M^{2}_{\rm
Pl}}\,\frac{K_{1}(m_{\tilde{\phi}_{m}}/T)}{K_{2}(m_{\tilde{\phi}_{m}}/T)}.$

Note that the production rate of spin-2 Kaluza-Klein gravitons exceeds that of
spin-0 states by a factor of $\sim 283(n+2)/[16n(n-1)]\approx 4.7-35$ (for
$n=6-2$). The lifetimes of the scalar modes, however, are longer than those of
the spin-2 modes by a factor of $\sim 283/70$, leading to complementary
impacts on the light element abundances.

We are now in a position to calculate the evolution of the abundances of the
Kaluza-Klein gravitons by solving the following coupled set of differential
equations (for each level of the Kaluza-Klein tower, $m$):

$\displaystyle\frac{dn_{\tilde{h}_{m}}}{dt}$
$\displaystyle=-3Hn_{\tilde{h}_{m}}+P_{\tilde{h}_{m}}-n_{\tilde{h}_{m}}\langle\Gamma_{\tilde{h}_{m}}\rangle$
(9) $\displaystyle\frac{dn_{\tilde{\phi}_{m}}}{dt}$
$\displaystyle=-3Hn_{\tilde{\phi}_{m}}+P_{\tilde{\phi}_{m}}-n_{\tilde{\phi}_{m}}\langle\Gamma_{\tilde{\phi}_{m}}\rangle,$

where $H=(8\pi\rho/3M^{2}_{\rm Pl})^{1/2}$ is the rate of Hubble expansion and
$\rho=\rho_{\rm
SM}+\sum_{m}\rho_{\tilde{h}_{m}}+\sum_{m}\rho_{\tilde{\phi}_{m}}$ is the total
energy density of the universe, $\rho_{\rm SM}=\pi^{2}g_{\star}T^{4}/30$ is
the energy density in Standard Model particles, and $g_{\star}$ is the number
of relativistic degrees-of-freedom. We evolve the temperature of the Standard
Model bath by applying entropy conservation, $T\propto
a^{-1}g_{\star,S}^{-1/3}$, where $g_{\star,S}$ is the number of relativistic
degrees-of-freedom in entropy.

In Fig. 2, we plot the production rate of Kaluza-Klein gravitons (divided by
three powers of the temperature) for several selected levels of the Kaluza-
Klein tower, and for the case of the $M_{\star}=10^{6}\,{\rm GeV}$ and $n=6$.
Since these production rates are many order of magnitude below the rate of
Hubble expansion, the Kaluza-Klein graviton abundances never reach their
equilibrium values, placing us safely within the regime of thermal freeze-in.
Unlike more typical freeze-in scenarios, however, the large multiplicity of
Kaluza-Klein states can greatly enhance the total energy density of these
particles that is produced. The case shown in Fig. 2, for example, is
effectively that of the simultaneous freeze-in of $\sim 2\times 10^{4}$
different particle species, each of which contributes to the total resulting
abundance.

Figure 3: The accumulated non-relativistic densities of Kaluza-Klein
gravitons, $\rho_{i}=n_{i}m_{i}$, divided by the total energy density in
Standard Model particles for the case of $M_{\star}=10^{6}\,{\rm GeV}$ and
$n=6$, and for two values of the initial temperature, $T_{\rm max}$. We show
results for the Kaluza-Klein modes with levels of $m=10^{3},10^{2},10$, and 1,
as well as for the sum of all modes. The red lines represent the abundances of
the spin-2 states, while the black lines include both spin-$2$ and scalar
states. Note that the $m=10^{3}$ states begin to appreciably decay prior to
the onset of Big Bang nucleosynthesis.

In Fig. 3, we plot the evolution of the non-relativistic densities of Kaluza-
Klein gravitons, $\rho_{i}=n_{i}m_{i}$, compared to the total energy density
in Standard Model particles, for the case of $M_{\star}=10^{6}\,{\rm GeV}$ and
$n=6$, and for two values of the initial temperature, $T_{\rm max}$. Results
are shown several selected values of the Kaluza-Klein levels, $m$, as well as
for the sum of all modes. In the left frame, one can notice that the
$m=10^{3}$ states begin to appreciably decay prior to the onset of Big Bang
nucleosynthesis.

##  III Kaluza-Klein Graviton Decays During Big Bang Nucleosynthesis

Measurements of the primordial deuterium [26, 27] and helium [28, 29, 30]
abundances provide us with the earliest probe of our universe’s thermal
history, confirming that our universe was radiation dominated and generally
well-described by $\Lambda$CDM cosmology throughout the era of Big Bang
nucleosynthesis, which began $t\sim 1\,{s}$ after the Big Bang [31, 32, 33,
34, 35, 36]. Such measurements allow us to place stringent constraints on the
expansion history of our universe, as well as on any energy injection that may
have taken place during or after this era [37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60].

Figure 4: Constraints on a particle species decaying to quarks from
measurements of the primordial element abundances, as presented in Ref. [47].
From top-to-bottom, the curves correspond to the upper limits on the abundance
of particles of mass $10^{6}\,{\rm GeV}$, $10^{5}\,{\rm GeV}$, $10^{4}\,{\rm
GeV}$, $10^{3}\,{\rm GeV}$, $10^{2}\,{\rm GeV}$, and $30\,{\rm GeV}$. These
constraints are presented in terms of the number density of decaying particles
per unit entropy (prior to their decays), $Y_{X}=n_{X}/s$.

The presence of Kaluza-Klein gravitons in the early universe could have
potentially impacted the light element abundances in a number of ways. In
particular, their decay products could have broken up helium nuclei through
the processes of photodissociation and hadrodissociation, reducing the
abundance of primordial helium while enhancing that of deuterium. These and
other such processes have been modeled in detail in a number of publicly
available codes [61, 62, 63]. In this study, we make use of the results of
Kawasaki et al. (2018), who evaluated the impact of decaying particles on the
resulting light element abundances [47] (see also Refs. [44, 46, 64, 48, 49,
50, 51, 52, 53, 54, 55, 56, 58, 59]). In particular, the authors of Ref. [47]
derived constraints on the lifetime and abundance of a decaying particle
species for various values of the particle’s mass and dominant decay modes.
These constraints were presented in terms of the mass of the decaying
particles multiplied by the number of such particles per unit entropy, $MY$,
as evaluated prior to their decays, $t\ll\tau$. Whereas that study considered
the impact of only one decaying particle species at a time, we are concerned
here with the decays of the entire tower of Kaluza-Klein gravitons. To recast
the results of Kawasaki et al. for the case at hand, we treat all Kaluza-Klein
states with lifetimes within a given decade (for example,
$\tau=10^{2}-10^{3}\,{\rm s}$) as a single particle species, with a lifetime
and mass equal to that of the median Kaluza-Klein state within that group. We
then compare this to Fig. 12 of Kawasaki et al. in order to determine whether
a given scenario is consistent with the measured helium and deuterium
abundances. We repeat this procedure for each decade of lifetime, allowing us
to produce conservative constraints on the values of $M_{\star}$ and $n$, as a
function of the initial temperature of the universe, $T_{\rm max}$. We include
in our calculations only those graviton decays that proceed to quarks or
gluons, and apply Kawasaki’s constraints on the $u\bar{u}$ channel (which are
almost entirely indistinguishable from those found in the cases of other
hadronic final states). The constraints are presented by Kawasaki et al. are
shown in Fig. 4.

Figure 5: The upper limits derived in this study on the maximum temperature of
the universe, $T_{\rm max}$, as a function of $M_{\star}$ and $n$. The dashed
portions of the lines represent values of $M_{\star}$ that are ruled out by
other considerations (see Sec. I). We only consider values of the temperatures
that are below the fundamental Planck scale, $T<M_{\star}$.

##  IV Results

The main results of this study are present in Fig. 5, where we plot the
maximum temperature of the universe that is consistent with the measured light
element abundances, as a function of $M_{\star}$ and $n$. Where the curves are
solid in this figure, the light element abundances provide the strongest
constraint on this class of models. In contrast, whereas the curves are
dashed, the constraints derived here are less restrictive than others that
have been presented in the literature (as summarized in Sec. I).

To understand the results that appear in Fig. 5, note that we are applying
constraints that apply to particles with lifetimes ranging from $\tau\sim
0.1\,{\rm s}$ to $10^{12}\,{\rm s}$. From Eq. 3, we find that this range
corresponds to Kaluza-Klein gravitons with masses in the range of
$m_{\tilde{h}_{m}}\sim 10\,{\rm GeV}$ to $\sim 2\times 10^{5}\,{\rm GeV}$.
Thus, if the universe was never at a temperature greater than
$\mathcal{O}({\rm GeV})$, the only Kaluza-Klein gravitons that could be
produced would be too long-lived to be constrained by the present analysis
(although other constraints, such as those from the cosmic microwave
background could still be potentially restrictive). On the other hand, in
order for the Kaluza-Klein gravitons to not decay prior to 0.1 s, (or
equivalently, for those states to be lighter than $m_{\tilde{h}_{m}}\sim
2\times 10^{5}\,{\rm GeV}$), the value of $M_{\star}$ must not be too high.
For the $n=1$ case, for example, the first Kaluza-Klein mode has a mass of
$m_{\tilde{h}_{1}}\sim 2\times 10^{5}\,{\rm GeV}$ for $M_{\star}\sim
10^{14}\,{\rm GeV}$, explaining why our constraints do not extend to larger
values of the fundamental Planck scale. Alternatively, for $n=6$,
$m_{\tilde{h}_{1}}\sim 2\times 10^{5}\,{\rm GeV}$ corresponds to a value of
$M_{\star}\sim 10^{8}\,{\rm GeV}$. In this later case, you may notice that our
constraints only extend up to $M_{\star}\sim 10^{6}\,{\rm GeV}$. The reason
that they do not extend to higher values of $M_{\star}$ is that the number of
decaying particles per unit entropy is not particular large in this case and
the decays occur around $\tau\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\
[-0.07cm] $\sim$}}~{}0.1\,{\rm s}$, where the constraints are rather weak (see
Fig. 4). If we instead consider $M_{\star}\sim 10^{6}\,{\rm GeV}$, the
lightest Kaluza-Klein mode has a mass of $m_{\tilde{h}_{1}}\sim 300\,{\rm
GeV}$ and a lifetime of $\tau\sim 4\times 10^{7}\,{\rm s}$. Constraints in
this lifetime range are much more stringent, allowing us to exclude values of
$T_{\rm max}$ that are comparable to or larger than the mass of the lightest
Kaluza-Klein mode.

##  V Implications for Black Hole Production

It has long been appreciated that small-scale inhomogeneties in the early
universe may have led to the formation of primordial black holes [65, 66].
Alternatively, in models with extra spatial dimensions, black holes could have
been produced through the collisions of particles in the thermal bath, in
particular if the total energy of a collision exceeds the fundamental Planck
scale, $M_{\star}$ [67] (see also, Refs. [68, 69]).

In the context of the flat and compactified extra dimensions that we have
considered in this study, the differential production rate of black holes in a
thermal bath of temperature, $T$, is given by [70]

$\displaystyle\frac{d\Gamma_{\rm BH}}{dM_{\rm
BH}}=\frac{g_{\star}(T)^{2}}{8\pi^{4}}\bigg{(}\frac{8\Gamma(\frac{n+3}{2})}{n+2}\bigg{)}^{\frac{2}{n+1}}M_{\rm
BH}\,T^{2}\,\bigg{(}\frac{M_{\rm
BH}}{M_{\star}}\bigg{)}^{\frac{2n+4}{n+1}}\bigg{[}\frac{M_{\rm
BH}}{T}K_{1}(M_{\rm BH}/T)+2K_{2}(M_{\rm BH}/T)\bigg{]}\,\Theta(M_{\rm
BH}-M_{\star}),$ (10)

where $\Gamma$ is the gamma function, $K_{1}$ and $K_{2}$ are again the
modified Bessel functions of the second kind, and $\Theta$ is the Heavyside
step function. For $T\sim M_{\star}$, black holes can be produced at a very
high rate. In particular, after dropping order one factors, this expression
integrates to $\Gamma_{\rm BH}\sim M_{\star}^{4}$ for temperatures near the
fundamental Planck scale. This indicates that the black hole production rate
per particle could potentially be larger than the Hubble rate by a huge
factor, roughly $\sim(\Gamma_{\rm BH}/n_{\rm SM})/H\sim
4(M_{\star}^{4}/T^{3})/(T^{2}/M_{\rm Pl})\sim 4M_{\rm Pl}/M_{\star}$. In
contrast, the black hole production rate is dramatically suppressed at
temperatures lower than $M_{\star}$. This behavior is confirmed in Fig 6,
where we plot the integrated production rate of black holes per Standard Model
particle per Hubble time, for several values of $M_{\star}$ and for $n=1,2,4$
and 6 (this result depends only very weakly on the value of $n$).

Figure 6: The integrated production rate of black holes per Standard Model
particle per Hubble time, for several values of $M_{\star}$ and for $n=1,2,4$
and 6 (this result depends only very weakly on the value of $n$). For
comparison, we plot as a dashed curve the value at which each Standard Model
particles produces an average of one black hole per Hubble time, $\Gamma_{\rm
BH}/(n_{\rm SM}H)=1$.

If a large number of microscopic black holes had been generated in the early
universe, these objects could have had a number of potentially observable
impacts. In particular, the products of their Hawking evaporation could have
included particles that would act as dark radiation (and contribute to the
effective number of neutrino species, $N_{\rm eff}$), or contribute to the
dark matter density [71] (see also, Refs. [72, 73, 74]). Alternatively, in
some scenarios microscopic black holes could grow rapidly through accretion,
significantly delaying their evaporation [67, 75, 76].

The constraints obtained in this study severely limit the rate at which black
holes could have been produced through particle collisions in the early
universe. This can be seen by comparing the black hole production rates shown
in Fig. 6 to the constraints we have presented in Fig. 5. For the case of
$M_{\star}=10^{6}\,{\rm GeV}$, for example, our constraints derived from the
primoridal light element abundances require $T_{\rm
max}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}3.2\,{\rm TeV}$ for any value of $n\leq 6$. For such temperatures,
the black hole production rate is vanishingly small, $\Gamma_{\rm BH}/(n_{\rm
SM}H)\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}10^{-120}$. Thus, for this value of $M_{\star}$, particle
collisions will not produce any significant abundance of black holes in the
early universe. If we consider larger values of $M_{\star}$, black hole
production could still be potentially important. For $M_{\star}=10^{9}\,{\rm
GeV}$ ($10^{12}\,{\rm GeV}$), for example, large black hole production rates
are still possible, provided that $n\geq 3$ ($n\geq 2$).

##  VI Summary and Discussion

In this study, we have considered the impact of Kaluza-Klein gravitons on the
primordial light element abundances, as established during the era of Big Bang
nucleosynthesis. In particular, we have focused on the ADD-scenario, which
features $n=1-6$ extra spatial dimensions which are flat and compactified
around a torus of radius, $R$. We have taken the fields of the Standard Model
to be confined to a 3+1 dimensional brane, allowing the fundamental scale of
gravity to be much lower than the effective 4-dimensional Planck scale,
$M_{\star}\ll M_{\rm Pl}$ [10, 10, 12].

In this model, massless spin-2 gravitons propagating in the $n+4$ dimensional
bulk appear to observers on the brane as a tower of massive Kaluza-Klein
states, with an evenly-spaced series of masses, $m_{\rm KK}=m/R$, where $m$ is
the level of the Kaluza-Klein tower. These Kaluza-Klein gravitons can be
produced in the early universe and, for $m_{\rm
KK}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}10^{6}\,{\rm GeV}$, will subsequently decay during or after the era
of Big Bang nucleosynthesis. Such decays can produce energetic Standard Model
particles which can break up helium nuclei, significantly altering the
predicted abundances of primoridal helium and deuterium.

We have calculated the production and decay rates of Kaluza-Klein gravitons in
the early universe, and evaluated the impact of these decays on the primordial
light element abundances. These abundances, in turn, allow us to place
constraints on Kalzua-Klein gravitons with masses in the range of $m_{\rm
KK}\sim 10\,{\rm GeV}$ to $\sim 2\times 10^{5}\,{\rm GeV}$, as summarized in
Fig. 5. For one extra dimension, $n=1$, we exclude all values of
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}2\times 10^{13}\,{\rm GeV}$
($\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}10^{10}\,{\rm GeV}$) unless the maximum temperature of the universe
was less than $\sim 2\,{\rm TeV}$ ($\sim 1\,{\rm GeV}$). For larger values of
$n$, our constraints are somewhat less stringent. For the case of $n=6$, for
example, our analysis excludes all values of
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}(1-2)\times 10^{6}\,{\rm GeV}$ unless the maximum temperature of
the universe was less than $\sim 3\,{\rm TeV}$.

The results presented here severely limit the possibility that black holes may
have been produced in significant numbers through particle collisions in the
early universe’s thermal bath [70]. For the case of
$M_{\star}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm]
$\sim$}}~{}10^{6}\,{\rm GeV}$, we can rule out the possibility that any
appreciable abundance of black holes was formed through such collisions, for
any value of $n\leq 6$. For larger values of $M_{\star}$, thermal black hole
production may still have been potentially important. For
$M_{\star}=10^{9}\,{\rm GeV}$, for example, efficient black hole production
could have occurred in the early universe, provided that $n\geq 3$.

Before closing, we note that the models we are considering here are similar
to, but not the same as, the “Dark Dimension” scenario that has been proposed
within the context of the Swampland program of string theory [77, 78]. In
particular, whereas Kaluza-Klein gravitons can decay into lighter KK modes in
the Dark Dimension picture, extra-dimensional momentum conservation prevents
such decays in the ADD model (decays into Standard Model particles are
allowed, as the brane is a dynamical object which can recoil to take on the
momentum of the decaying graviton). We leave the consideration of this class
of models to future work.

Acknowledgments

We would like to thank Gordan Krnjaic, Yann Mambrini, Jamie Law-Smith, and Joe
Lykken for helpful discussions. DH is supported by the Fermi Research
Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of
Energy, Office of High Energy Physics. This project has also received funding
from the European Union’s Horizon Europe research and innovation program under
the Marie Skłodowska-Curie Staff Exchange grant agreement No. 101086085.

## References

  * [1] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram, and S. Hannestad, “MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles,” JCAP 12 (2019) 012, arXiv:1908.10189 [hep-ph].
  * [2] S. Hannestad, “What is the lowest possible reheating temperature?,” Phys. Rev. D 70 (2004) 043506, arXiv:astro-ph/0403291.
  * [3] G. F. Giudice, E. W. Kolb, and A. Riotto, “Largest temperature of the radiation era and its cosmological implications,” Phys. Rev. D 64 (2001) 023508, arXiv:hep-ph/0005123.
  * [4] M. Kawasaki, K. Kohri, and N. Sugiyama, “MeV scale reheating temperature and thermalization of neutrino background,” Phys. Rev. D 62 (2000) 023506, arXiv:astro-ph/0002127.
  * [5] M. Kawasaki, K. Kohri, and N. Sugiyama, “Cosmological constraints on late time entropy production,” Phys. Rev. Lett. 82 (1999) 4168, arXiv:astro-ph/9811437.
  * [6] J. Ghiglieri, G. Jackson, M. Laine, and Y. Zhu, “Gravitational wave background from Standard Model physics: Complete leading order,” JHEP 07 (2020) 092, arXiv:2004.11392 [hep-ph].
  * [7] J. Ghiglieri and M. Laine, “Gravitational wave background from Standard Model physics: Qualitative features,” JCAP 07 (2015) 022, arXiv:1504.02569 [hep-ph].
  * [8] B. X. Hu and A. Loeb, “An Upper Limit on the Initial Temperature of the Radiation-Dominated Universe,” JCAP 01 (2021) 041, arXiv:2004.02895 [astro-ph.CO].
  * [9] A. Ringwald, J. Schütte-Engel, and C. Tamarit, “Gravitational Waves as a Big Bang Thermometer,” JCAP 03 (2021) 054, arXiv:2011.04731 [hep-ph].
  * [10] N. Arkani-Hamed, S. Dimopoulos, and G. R. Dvali, “The Hierarchy problem and new dimensions at a millimeter,” Phys. Lett. B 429 (1998) 263–272, arXiv:hep-ph/9803315.
  * [11] I. Antoniadis, N. Arkani-Hamed, S. Dimopoulos, and G. R. Dvali, “New dimensions at a millimeter to a Fermi and superstrings at a TeV,” Phys. Lett. B 436 (1998) 257–263, arXiv:hep-ph/9804398.
  * [12] N. Arkani-Hamed, S. Dimopoulos, and G. R. Dvali, “Phenomenology, astrophysics and cosmology of theories with submillimeter dimensions and TeV scale quantum gravity,” Phys. Rev. D 59 (1999) 086004, arXiv:hep-ph/9807344.
  * [13] CMS Collaboration, A. M. Sirunyan et al., “Search for new physics in final states with an energetic jet or a hadronically decaying $W$ or $Z$ boson and transverse momentum imbalance at $\sqrt{s}=13\text{ }\text{ }\mathrm{TeV}$,” Phys. Rev. D 97 no. 9, (2018) 092005, arXiv:1712.02345 [hep-ex].
  * [14] ATLAS Collaboration, M. Aaboud et al., “Search for dark matter and other new phenomena in events with an energetic jet and large missing transverse momentum using the ATLAS detector,” JHEP 01 (2018) 126, arXiv:1711.03301 [hep-ex].
  * [15] E. G. Adelberger, J. H. Gundlach, B. R. Heckel, S. Hoedl, and S. Schlamminger, “Torsion balance experiments: A low-energy frontier of particle physics,” Prog. Part. Nucl. Phys. 62 (2009) 102–134.
  * [16] J. Murata and S. Tanaka, “A review of short-range gravity experiments in the LHC era,” Class. Quant. Grav. 32 no. 3, (2015) 033001, arXiv:1408.3588 [hep-ex].
  * [17] W.-H. Tan, S.-Q. Yang, C.-G. Shao, J. Li, A.-B. Du, B.-F. Zhan, Q.-L. Wang, P.-S. Luo, L.-C. Tu, and J. Luo, “New Test of the Gravitational Inverse-Square Law at the Submillimeter Range with Dual Modulation and Compensation,” Phys. Rev. Lett. 116 no. 13, (2016) 131101.
  * [18] S. Hannestad and G. G. Raffelt, “Supernova and neutron star limits on large extra dimensions reexamined,” Phys. Rev. D 67 (2003) 125008, arXiv:hep-ph/0304029. [Erratum: Phys.Rev.D 69, 029901 (2004)].
  * [19] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [20] G. Shiu and S. H. H. Tye, “TeV scale superstring and extra dimensions,” Phys. Rev. D 58 (1998) 106007, arXiv:hep-th/9805157.
  * [21] T. Han, J. D. Lykken, and R.-J. Zhang, “On Kaluza-Klein states from large extra dimensions,” Phys. Rev. D 59 (1999) 105006, arXiv:hep-ph/9811350.
  * [22] L. J. Hall and D. Tucker-Smith, “Cosmological constraints on theories with large extra dimensions,” Phys. Rev. D 60 (1999) 085008, arXiv:hep-ph/9904267.
  * [23] A. de Giorgi and S. Vogl, “Dark matter interacting via a massive spin-2 mediator in warped extra-dimensions,” JHEP 11 (2021) 036, arXiv:2105.06794 [hep-ph].
  * [24] A. de Giorgi and S. Vogl, “Warm dark matter from a gravitational freeze-in in extra dimensions,” JHEP 04 (2023) 032, arXiv:2208.03153 [hep-ph].
  * [25] M. Escudero, D. Hooper, G. Krnjaic, and M. Pierre, “Cosmology with A Very Light Lμ $-$ Lτ Gauge Boson,” JHEP 03 (2019) 071, arXiv:1901.02010 [hep-ph].
  * [26] R. Cooke, M. Pettini, R. A. Jorgenson, M. T. Murphy, and C. C. Steidel, “Precision measures of the primordial abundance of deuterium,” Astrophys. J. 781 no. 1, (2014) 31, arXiv:1308.3240 [astro-ph.CO].
  * [27] S. Riemer-Sørensen and E. S. Jenssen, “Nucleosynthesis Predictions and High-Precision Deuterium Measurements,” Universe 3 no. 2, (2017) 44, arXiv:1705.03653 [astro-ph.CO].
  * [28] Y. I. Izotov, T. X. Thuan, and N. G. Guseva, “A new determination of the primordial He abundance using the He i $\lambda$10830 Å emission line: cosmological implications,” Mon. Not. Roy. Astron. Soc. 445 no. 1, (2014) 778–793, arXiv:1408.6953 [astro-ph.CO].
  * [29] E. Aver, K. A. Olive, and E. D. Skillman, “The effects of He I $\lambda$10830 on helium abundance determinations,” JCAP 07 (2015) 011, arXiv:1503.08146 [astro-ph.CO].
  * [30] A. Peimbert, M. Peimbert, and V. Luridiana, “The primordial helium abundance and the number of neutrino families,” Rev. Mex. Astron. Astrofis. 52 no. 2, (2016) 419–424, arXiv:1608.02062 [astro-ph.CO].
  * [31] D. N. Schramm and M. S. Turner, “Big Bang Nucleosynthesis Enters the Precision Era,” Rev. Mod. Phys. 70 (1998) 303–318, arXiv:astro-ph/9706069.
  * [32] G. Steigman, “Primordial Nucleosynthesis in the Precision Cosmology Era,” Ann. Rev. Nucl. Part. Sci. 57 (2007) 463–491, arXiv:0712.1100 [astro-ph].
  * [33] F. Iocco, G. Mangano, G. Miele, O. Pisanti, and P. D. Serpico, “Primordial Nucleosynthesis: from precision cosmology to fundamental physics,” Phys. Rept. 472 (2009) 1–76, arXiv:0809.0631 [astro-ph].
  * [34] R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, “Big Bang Nucleosynthesis: 2015,” Rev. Mod. Phys. 88 (2016) 015004, arXiv:1505.01076 [astro-ph.CO].
  * [35] C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, “Precision big bang nucleosynthesis with improved Helium-4 predictions,” Phys. Rept. 754 (2018) 1–66, arXiv:1801.08023 [astro-ph.CO].
  * [36] B. D. Fields, K. A. Olive, T.-H. Yeh, and C. Young, “Big-Bang Nucleosynthesis after Planck,” JCAP 03 (2020) 010, arXiv:1912.01132 [astro-ph.CO]. [Erratum: JCAP 11, E02 (2020)].
  * [37] S. Sarkar, “Big bang nucleosynthesis and physics beyond the standard model,” Rept. Prog. Phys. 59 (1996) 1493–1610, arXiv:hep-ph/9602260.
  * [38] M. Pospelov and J. Pradler, “Big Bang Nucleosynthesis as a Probe of New Physics,” Ann. Rev. Nucl. Part. Sci. 60 (2010) 539–568, arXiv:1011.1054 [hep-ph].
  * [39] K. Jedamzik and M. Pospelov, “Big Bang Nucleosynthesis and Particle Dark Matter,” New J. Phys. 11 (2009) 105028, arXiv:0906.2087 [hep-ph].
  * [40] M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, “BBN constraints on MeV-scale dark sectors. Part I. Sterile decays,” JCAP 02 (2018) 044, arXiv:1712.03972 [hep-ph].
  * [41] C. Keith, D. Hooper, N. Blinov, and S. D. McDermott, “Constraints on Primordial Black Holes From Big Bang Nucleosynthesis Revisited,” Phys. Rev. D 102 no. 10, (2020) 103512, arXiv:2006.03608 [astro-ph.CO].
  * [42] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev. D 81 (2010) 104019, arXiv:0912.5297 [astro-ph.CO].
  * [43] G.-y. Huang, T. Ohlsson, and S. Zhou, “Observational Constraints on Secret Neutrino Interactions from Big Bang Nucleosynthesis,” Phys. Rev. D 97 no. 7, (2018) 075009, arXiv:1712.04792 [hep-ph].
  * [44] L. Forestell, D. E. Morrissey, and G. White, “Limits from BBN on Light Electromagnetic Decays,” JHEP 01 (2019) 074, arXiv:1809.01179 [hep-ph].
  * [45] M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, “BBN constraints on MeV-scale dark sectors. Part II. Electromagnetic decays,” JCAP 11 (2018) 032, arXiv:1808.09324 [hep-ph].
  * [46] P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, “BBN constraints on the annihilation of MeV-scale dark matter,” JCAP 04 (2019) 029, arXiv:1901.06944 [hep-ph].
  * [47] M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, “Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles,” Phys. Rev. D 97 no. 2, (2018) 023502, arXiv:1709.01211 [hep-ph].
  * [48] V. Poulin and P. D. Serpico, “Nonuniversal BBN bounds on electromagnetically decaying particles,” Phys. Rev. D 91 no. 10, (2015) 103007, arXiv:1503.04852 [astro-ph.CO].
  * [49] R. H. Cyburt, J. Ellis, B. D. Fields, F. Luo, K. A. Olive, and V. C. Spanos, “Nucleosynthesis Constraints on a Massive Gravitino in Neutralino Dark Matter Scenarios,” JCAP 10 (2009) 021, arXiv:0907.5003 [astro-ph.CO].
  * [50] M. Kusakabe, T. Kajino, T. Yoshida, T. Shima, Y. Nagai, and T. Kii, “New Constraints on Radiative Decay of Long-Lived X Particles in Big Bang Nucleosynthesis with New Rates of Photodisintegration Reactions of 4He,” Phys. Rev. D 79 (2009) 123513, arXiv:0806.4040 [astro-ph].
  * [51] K. Jedamzik, “Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles,” Phys. Rev. D 74 (2006) 103509, arXiv:hep-ph/0604251.
  * [52] K. Kohri, T. Moroi, and A. Yotsuyanagi, “Big-bang nucleosynthesis with unstable gravitino and upper bound on the reheating temperature,” Phys. Rev. D 73 (2006) 123511, arXiv:hep-ph/0507245.
  * [53] K. Jedamzik, “Did something decay, evaporate, or annihilate during Big Bang nucleosynthesis?,” Phys. Rev. D 70 (2004) 063524, arXiv:astro-ph/0402344.
  * [54] M. Kawasaki, K. Kohri, and T. Moroi, “Big-Bang nucleosynthesis and hadronic decay of long-lived massive particles,” Phys. Rev. D 71 (2005) 083502, arXiv:astro-ph/0408426.
  * [55] M. Kawasaki, K. Kohri, and T. Moroi, “Hadronic decay of late - decaying particles and Big-Bang Nucleosynthesis,” Phys. Lett. B 625 (2005) 7–12, arXiv:astro-ph/0402490.
  * [56] R. H. Cyburt, J. R. Ellis, B. D. Fields, and K. A. Olive, “Updated nucleosynthesis constraints on unstable relic particles,” Phys. Rev. D 67 (2003) 103521, arXiv:astro-ph/0211258.
  * [57] K. Kohri, “Primordial nucleosynthesis and hadronic decay of a massive particle with a relatively short lifetime,” Phys. Rev. D 64 (2001) 043515, arXiv:astro-ph/0103411.
  * [58] M. Kawasaki, K. Kohri, and T. Moroi, “Radiative decay of a massive particle and the nonthermal process in primordial nucleosynthesis,” Phys. Rev. D 63 (2001) 103502, arXiv:hep-ph/0012279.
  * [59] E. Holtmann, M. Kawasaki, K. Kohri, and T. Moroi, “Radiative decay of a longlived particle and big bang nucleosynthesis,” Phys. Rev. D 60 (1999) 023506, arXiv:hep-ph/9805405.
  * [60] M. Kawasaki and T. Moroi, “Electromagnetic cascade in the early universe and its application to the big bang nucleosynthesis,” Astrophys. J. 452 (1995) 506, arXiv:astro-ph/9412055.
  * [61] A. Arbey, J. Auffinger, K. P. Hickerson, and E. S. Jenssen, “AlterBBN v2: A public code for calculating Big-Bang nucleosynthesis constraints in alternative cosmologies,” Comput. Phys. Commun. 248 (2020) 106982, arXiv:1806.11095 [astro-ph.CO].
  * [62] C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, “Precision Big Bang Nucleosynthesis with the New Code PRIMAT,” JPS Conf. Proc. 31 (2020) 011034, arXiv:1909.12046 [astro-ph.CO].
  * [63] O. Pisanti, A. Cirillo, S. Esposito, F. Iocco, G. Mangano, G. Miele, and P. D. Serpico, “PArthENoPE: Public Algorithm Evaluating the Nucleosynthesis of Primordial Elements,” Comput. Phys. Commun. 178 (2008) 956–971, arXiv:0705.0290 [astro-ph].
  * [64] M. Kawasaki and T. Moroi, “Gravitino production in the inflationary universe and the effects on big bang nucleosynthesis,” Prog. Theor. Phys. 93 (1995) 879–900, arXiv:hep-ph/9403364.
  * [65] B. J. Carr and S. W. Hawking, “Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168 (1974) 399–415.
  * [66] B. J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201 (1975) 1–19.
  * [67] J. A. Conley and T. Wizansky, “Microscopic Primordial Black Holes and Extra Dimensions,” Phys. Rev. D 75 (2007) 044006, arXiv:hep-ph/0611091.
  * [68] P. C. Argyres, S. Dimopoulos, and J. March-Russell, “Black holes and submillimeter dimensions,” Phys. Lett. B 441 (1998) 96–104, arXiv:hep-th/9808138.
  * [69] R. C. Myers and M. J. Perry, “Black Holes in Higher Dimensional Space-Times,” Annals Phys. 172 (1986) 304.
  * [70] A. Friedlander, K. J. Mack, S. Schon, N. Song, and A. C. Vincent, “Primordial black hole dark matter in the context of extra dimensions,” Phys. Rev. D 105 no. 10, (2022) 103508, arXiv:2201.11761 [hep-ph].
  * [71] A. Friedlander, N. Song, and A. C. Vincent, “Dark matter from higher-dimensional primordial black holes,” Phys. Rev. D 108 no. 4, (2023) 043523, arXiv:2306.01520 [hep-ph].
  * [72] D. Hooper, G. Krnjaic, and S. D. McDermott, “Dark Radiation and Superheavy Dark Matter from Black Hole Domination,” JHEP 08 (2019) 001, arXiv:1905.01301 [hep-ph].
  * [73] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott, and R. Petrossian-Byrne, “Hot Gravitons and Gravitational Waves From Kerr Black Holes in the Early Universe,” arXiv:2004.00618 [astro-ph.CO].
  * [74] A. Friedlander, N. Song, and A. C. Vincent, “Dark matter from hot big bang black holes,” Phys. Rev. D 108 no. 8, (2023) L081301, arXiv:2303.07372 [hep-ph].
  * [75] R. Guedens, D. Clancy, and A. R. Liddle, “Primordial black holes in brane world cosmologies: Accretion after formation,” Phys. Rev. D 66 (2002) 083509, arXiv:astro-ph/0208299.
  * [76] A. S. Majumdar, “Domination of black hole accretion in brane cosmology,” Phys. Rev. Lett. 90 (2003) 031303, arXiv:astro-ph/0208048.
  * [77] E. Gonzalo, M. Montero, G. Obied, and C. Vafa, “Dark dimension gravitons as dark matter,” JHEP 11 (2023) 109, arXiv:2209.09249 [hep-ph].
  * [78] M. Montero, C. Vafa, and I. Valenzuela, “The dark dimension and the Swampland,” JHEP 02 (2023) 022, arXiv:2205.12293 [hep-th].


# Towards a FAIR Documentation of Workflows and Models in Applied Mathematics

Marco Reidelbach Corresponding author Mathematics of Complex Systems, Zuse
Institute Berlin, 14195 Berlin, Germany<EMAIL_ADDRESS>Björn
Schembera  Institute of Applied Analysis and Numerical Simulation, University
of Stuttgart, 70569 Stuttgart, Germany<EMAIL_ADDRESS>stuttgart.de Marcus Weber  Mathematics of Complex Systems, Zuse Institute
Berlin, 14195 Berlin, Germany<EMAIL_ADDRESS>
###### Abstract

Modeling-Simulation-Optimization workflows play a fundamental role in applied
mathematics. The Mathematical Research Data Initiative, MaRDI, responded to
this by developing a FAIR and machine-interpretable template for a
comprehensive documentation of such workflows. MaRDMO, a Plugin for the
Research Data Management Organiser, enables scientists from diverse fields to
document and publish their workflows on the MaRDI Portal seamlessly using the
MaRDI template. Central to these workflows are mathematical models. MaRDI
addresses them with the MathModDB ontology, offering a structured formal model
description. Here, we showcase the interaction between MaRDMO and the
MathModDB Knowledge Graph through an algebraic modeling workflow from the
Digital Humanities. This demonstration underscores the versatility of both
services beyond their original numerical domain.

##  1 Introduction

Mathematical research data holds a pivotal role in advancing scientific
understanding across various disciplines, ranging from core mathematical
sciences to applied fields such as Engineering, Physics, and Digital
Humanities. Mathematical research data is not limited to numerical or symbolic
data contained in datasets, but also to data about the model, the solution
algorithm and many more appearing in the Model-Simulation-Optimization (MSO)
workflow [11, 17]. In response to the intricate nature of mathematical
workflows and the need for standardized documentation, the Mathematical
Research Data Initiative (MaRDI [17]), a project under the German National
Research Data Infrastructure (NFDI [8]), aimed at developing a robust research
data infrastructure for mathematics, developed a workflow documentation
template adhering to the FAIR (Findable, Accessible, Interoperable, and
Reusable) principles [20]. This template provides a comprehensive
documentation of MSO workflows in applied mathematics, including mathematical
models, methods, software, hardware, input, and output data addressing
specific research objectives [3]. The template is related to the developments
in the NFDI4ING consortium [16], where the Metadata4Ing [1] ontology was
created to capture engineering workflows. To facilitate the widespread
adoption of this documentation template, MaRDI introduced the MaRDMO Plugin
[13], a tool integrated in the Research Data Management Organiser (RDMO [4]).
Leveraging the popularity of RDMO as the most widely used software for the
creation of data management plans (DMPs) in Germany [5], MaRDMO streamlines
the process of MSO workflow documentation by fetching additional information
from various sources such as Wikidata [18], swMath [7] and zbMath [10] and
making the complete documentation accessible via the MaRDI
Portal111https://portal.mardi4nfdi.de, a wikibase for Mathematics. This semi-
automatic approach enhances efficiency while ensuring the completeness and
accuracy of the documented workflows.

At the heart of these MSO workflows lies the proper documentation of
mathematical models. Recognizing this fundamental aspect, MaRDI has developed
the Mathematical Models Database (MathModDB) ontology [15, 14], meticulously
designed to capture essential elements associated with mathematical models,
including research fields, problems, formulations, quantities, and tasks.
Numerical models documented through this scheme will be made accessible to the
public through the MathModDB Knowledge Graph (KG), ensuring widespread
availability and utilization of these valuable resources.

In this paper, we aim to establish a connection between MaRDMO and the
MathModDB KG, allowing researchers to access a standardized collection of
mathematical models through MaRDMO. Additionally, MaRDMO can serve as an
additional interface for MathModDB to gather new models from scientists across
disciplines via RDMO. Through a practical example of a Logical Data Analysis
algebraic workflow from the Digital Humanities, we demonstrate how this
integration facilitates comprehensive documentation and analysis, enhancing
reproducibility, transparency, collaboration, and interdisciplinary innovation
in the scientific community. In addition, we show the transferability of
MaRDMO and MathModDB, initially developed for capturing MSO workflows and
their numerical models, to algebraic modeling applications.

##  2 MaRDMO, MathModDB and their Connection

###  2.1 MaRDMO

The MaRDMO plugin builds upon the capabilities of RDMO, which enables the
creation of DMPs through customizable questionnaires. By now, general
questionnaires exist reflecting the requirements of individual funding
organizations, as well as subject-specific
catalogs.222https://github.com/rdmorganiser/rdmo-catalog. With MaRDMO we
introduced a questionnaire333https://github.com/MarcoReidelbach/MaRDMO-
Questionnaire initially tailored for the documentation of MSO workflows and a
customized export plugin444https://github.com/MarcoReidelbach/MaRDMO-Export-
Plugin. Such plugins extend the functionality of RDMO, e.g. facilitating an
export to Zenodo555https://github.com/rdmorganiser/rdmo-plugins. In turn,
MaRDMO enables direct export to the MaRDI Portal. As far as we know, this is
the first connection of RDMO with a KG.

The MaRDMO questionnaire guides researchers through the documentation process,
facilitating the capture of crucial information about their workflows. Divided
into four sections, the questionnaire covers various aspects of the research
process, including general aspects (1), models, variables and parameters (2),
process information (3), and reproducibility considerations (4). Through a
series of structured questions, researchers can input details about their
research objectives, mathematical models, software and hardware used, input
and output data, methods, and more (c.f. complete questionnaire66footnotemark:
6).

The guiding principle of MaRDMO is reusing existing information wherever
possible from established data sources. This approach ensures proper
integration into the existing research data landscape while minimizing
duplication of effort. The latest version of MaRDMO has streamlined the
retrieval of additional information, now seamlessly conducted in the
background. This enhancement, coupled with dynamic on-the-fly queries of
repositories enhances usability and functionality of the plugin. For instance,
when the DOI of a publication related to the workflow is provided, MaRDMO
automatically retrieves all relevant citation and author information. This
information is then presented for validation and completion, ensuring proper
documentation of the workflow’s context.

Upon completion of the questionnaire, MaRDMO facilitates the export of the
documented workflow to the MaRDI Portal. Therefore, a comprehensive summary is
generated following the MaRDI template and published as wikipage. Essential
details (related publication, scientific fields and mathematical areas,
applied mathematical model, methods, software, input and output data) are
integrated into the KG of the MaRDI Portal, making them accessible to other
researchers for search and exploration. Thereby, MaRDMO empowers researchers
to efficiently capture and share their MSO workflows, fostering collaboration,
reproducibility, and transparency within the scientific community.

###  2.2 MathModDB

The MathModDB ontology is designed to specifically address the intricate
nature of mathematical models in the domain of MSO. Within the vast landscape
of scientific reasoning, mathematical models play a fundamental role, serving
as indispensable tools for abstraction, formalization, analysis, and
comprehension across diverse disciplines. However, the complexity and
diversity inherent in mathematical models necessitate a unified semantic
framework to comprehensively capture their essence.

Recognizing this imperative, MathModDB emerges as a sophisticated ontology
crafted to meet the multifaceted demands of mathematical modeling. By
structuring mathematical knowledge into a coherent framework, MathModDB
enhances the semantic representation of mathematical models as they appear in
the MSO workflow, while fostering interoperability and accessibility across
domains.

The ontology is structured around eight classes essential for describing
mathematical models comprehensively, a first result of iterative (and still
ongoing) development [14, 15] shaped by project internal discussions and
valuable feedback from the mathematics community. These classes have evolved
over time to reflect the diverse needs and perspectives within the
mathematical modeling domain, ensuring a comprehensive representation of
mathematical models. Moreover, this iterative process underscores the
ontology’s adaptability, with the potential for further evolution in response
to emerging trends and evolving requirements within the mathematics community.

The classes of MathModDB include the Mathematical Model itself, the Research
Field within which the model operates, the Research Problem it addresses, the
Mathematical Formulations formalizing the model, the Quantities and Quantity
Kinds involved in the Mathematical Formulations, the Computational Tasks
associated with the model, and the Publications inventing, studying,
surveying, or using the model. The Computational Task class was introduced to
MathModDB only recently to accommodate the diverse range of tasks or questions
posed to a model, resulting in different formulations, inputs and outputs.

###  2.3 MaRDMO and MathModDB

The connection of MaRDMO and MathModDB marks a significant advancement in the
documentation of workflows and mathematical models, so far only stemming from
the domain of MSO. Prior to the integration, MaRDMO’s documentation of
mathematical models, not yet present in established repositories, was limited,
as researchers could only provide basic information such as model name,
description, main subject, defining formulas, and identifiers. This often fell
short of providing a comprehensive understanding of the model, especially when
considered in isolation from the contextual workflow.

To address this limitation, the MathModDB ontology has been integrated into
MaRDMO, facilitating the generation of comprehensive and insightful
documentation for the models used in MSO workflows. This integration
introduces additional question sets aligned with the classes of MathModDB.
Following MaRDMO’s guiding principle of utilizing existing resources, these
question sets leverage the MathModDB KG. Now, mathematical models can be
documented by connecting documented aspects of other models using the
MathModDB vocabulary. In cases where no suitable entity exists, new ones can
be created, with each domain featuring mandatory and optional fields. To
integrate new mathematical models and their components properly into the
existing data landscape, MaRDMO is also able to establish model
interconnections through the MathModDB vocabulary, e.g. using generalizations,
specifications or combinations. While these interconnections are crucial for
maximizing the utility of the MathModDB KG, they might pose challenges for
individual researchers, leaving it an optional asset. If researchers opt not
to interconnect their model, the responsibility falls on domain experts to
curate the MathModDB KG. This dual approach ensures that the KG remains a
valuable asset while also accommodating the varying needs and perspectives of
researchers.

##  3 A semantic Representation of an Algebraic Modeling Workflow

We now show a first proof of concept of how a workflow from algebraic modeling
can be semantically represented by MaRDMO and MathModDB.

###  3.1 Introduction

The Logical Data Analysis algebraic modeling workflow, rooted in the Digital
Humanities, particularly Egyptology, intertwines the disciplines of
Mathematics and Archaeology to unearth hidden patterns within a database.
Originating from the Cachette de Karnak, an archaeological repository
unearthed in 1903 by G. Legrain [6], the workflow aims to discern underlying
rules governing the destruction patterns observed in ancient Egyptian objects.
Legrain’s findings revealed common destruction patterns such as missing heads,
amputated body parts, and fragmented pieces, sparking the inquiry into whether
specific rules governed these occurrences.

###  3.2 Object Data and Encoding

An online database666https://www.ifao.egnet.net/bases/cachette/ cataloging
statues discovered in the Cachette provides the foundation for the analysis.
Experts from Egyptology identified 16 potentially significant properties
inherent in a subset of 333 artifacts, which were then encoded into binary
numbers, denoting the presence or absence of each property.

###  3.3 Interdisciplinary Collaboration

The encoded dataset was then passed from Egyptology to Mathematics, marking
the interdisciplinary exchange of data. Mathematicians then applied an object
comparison model [19] (c.f. Fig 1), utilizing a boolean ring, to unravel the
underlying rules governing the destruction patterns. Using the “Rules and
Pattern” algorithm777https://github.com/pynoor/The-RAP-Algorithm, written in
Julia [2] and leveraging the OSCAR package [12], the generator of the ideal
was computed as a comprehensible Gröbner Basis comprising 172 logical rules.

Figure 1: Schematic representation of the object comparison model within the
field of Egyptology, designed to identify destruction rules. Developed by
Weber et al. (2022) [19], the model employs a boolean ring as its formulation,
with object properties serving as quantities, represented as booleans. The
model is associated with a task, featuring its own formulation, input, and
output quantities.

###  3.4 Interpretation and Validation

The computational output sheds light on potential underlying rules. The
validation of these expressions for scientific and statistical relevance falls
upon Egyptology experts, marking the return of data to its disciplinary
origin. This interdisciplinary exchange ensures a comprehensive evaluation of
the identified rules.

###  3.5 Reproducibility and Future Applications

The algebraic modeling workflow presented here follows a structured process,
encompassing data selection and encoding, computational analysis, and
interpretation of results. The methodology is reproducible, underscoring its
reliability for future studies. Beyond Archaeology, the Object Comparison
Model holds promise for solving diverse computational tasks and addressing
distinct problems across various domains, e.g. literary analysis [19].
Computational tasks, other than the “Extraction of Logical Rules” include
object sortings and feature extractions, broadening the applicability of the
model beyond its initial context.

The documentation of the Logical Data Analysis workflow, following the MaRDI
template, can be accessed through the MaRDI Portal888Documentation Video,
Wikipage and Knowledge Graph Entry:
https://portal.mardi4nfdi.de/wiki/MaRDMO
https://portal.mardi4nfdi.de/wiki/Logical\\_Data\\_Analysis\\_for\\_Egyptian\\_Objects
https://portal.mardi4nfdi.de/wiki/Item:Q6032641. The underlying “Object
Comparison Model” is visualized in Fig 1 since the MathModDB KG is not yet
publicly accessible. All aspects of the algebraic modeling workflow described
before could be documented through the MaRDMO Plugin and the integrated
MathModDB Ontology.

##  4 Conclusion and Outlook

The integration of MathModDB into MaRDMO significantly enhances the quality of
the workflow documentation, providing a more comprehensive and insightful
portrayal of applied mathematical models. Through RDMO, we expect, once again,
to reach scientists from diverse disciplines and facilitate their access to
standardized and enriched model documentation.

We have shown how a research project from a seemingly completely different
field, namely the humanities, can benefit from our approach. This involves the
differentiation of destruction patterns in ancient objects, which can be
modeled using mathematical/algebraic methods and semantically enriched and
represented with the help of the tools presented (MaRDMO and MathModDB).
Moreover, we have demonstrated that our solutions work beyond the domain of
the classic numerical MSO workflow as they were applied to an algebraic
modeling workflow. All the relevant information was captured by tools
initially developed for MSO.

Looking ahead, our next objective is to connect MaRDMO with the Mathematical
Algorithms Database (MathAlgoDB [9]) KG999https://algodata.mardi4nfdi.de/,
also developed within MaRDI, to bolster the method component of the workflow
documentation. This integration could leverage existing connections between
MathAlgoDB and MathModDB, offering users a curated selection of algorithms to
address specific tasks associated with particular models, thus further
enhancing usability. Moreover, given the interdisciplinary nature of these
workflows, we aim to collaborate with the other consortia of the NFDI to
explore potential connections with non-mathematical services, expanding the
scope and utility of MaRDMO in the future. In addition, further research is
needed to show that the connection between MaRDMO and MathModDB also works for
more advanced workflows in the field of algebra.

#### Acknowledgements

Marco Reidelbach, Björn Schembera and Marcus Weber are supported by MaRDI,
funded by the Deutsche Forschungsgemeinschaft (DFG), project number 460135501,
NFDI 29/1 “MaRDI – Mathematische Forschungsdateninitiative”.

## References

  * [1] Arndt, S., Farnbacher, B., Fuhrmans, M., Hachinger, S., Hickmann, J., Hoppe, N., et al.: Metadata4Ing: An Ontology for Describing the Generation of Research Data within a Scientific Activity. (Feb 2022). https://doi.org/10.5281/zenodo.7706017
  * [2] Bezanson, J., Edelman, A., Karpinski, S., Shah, V.B.: Julia: A fresh approach to numerical computing. SIAM review 59(1), 65–98 (2017), https://doi.org/10.1137/141000671
  * [3] Boege, T., Fritze, R., Görgen, C., Hanselman, J., Iglezakis, D., Kastner, L., et al.: Research-data management planning in the German mathematical community. European Mathematical Society Magazine (2023). https://doi.org/10.4171/MAG/152
  * [4] Engelhardt, C., Enke, H., Klar, J., Ludwig, J., Neuroth, H.: Research data management organiser. In: Proceedings of the 14th International Conference on Digital Preservation. pp. 25–29 (2017)
  * [5] Enke, H., Hausen, D., Henzen, C., Jagusch, G., Krause, C., Schönau, S., et al.: Data management planning: Concept for setting up a working group in the nfdi section common infrastructures. Zenodo (2023). https://doi.org/10.5281/zenodo.7540682
  * [6] G. Legrain: Les récentes découvertes de Karnak. Bulletin de l’institut d’Égypte 5, 109–120 (1904)
  * [7] Greuel, G.M., Sperber, W.: swMATH – An Information Service for Mathematical Software. In: Mathematical Software – ICMS 2014. pp. 691–701. Springer Berlin Heidelberg (2014)
  * [8] Hartl, N., Wössner, E., Sure-Vetter, Y.: Nationale Forschungsdateninfrastruktur (NFDI). Informatik Spektrum 44, 370–373 (2021)
  * [9] Himpe, C., Kleikamp, H., Fritze, R., Rave, S.: MaRDI Task Area 2 - Scientific Computing @ WWU Münster. AlgoData - Algorithm Knowledge Graph - Ontology (Version 0.1) (2022), https://mardi4nfdi.de/algodata/0.1
  * [10] Hulek, K., Teschke, O.: Die transformation von zbmath zu einer offenen plattform für die mathematik. Mitteilungen der Deutschen Mathematiker-Vereinigung 28(2), 108–111 (2020). https://doi.org/doi:10.1515/dmvm-2020-0031
  * [11] Koprucki, T., Kohlhase, M., Tabelow, K., Müller, D., Rabe, F.: Model Pathway Diagrams for the Representation of Mathematical Models. Optical and Quantum Electronics 50, 1–9 (2018). https://doi.org/10.1007/s11082-018-1321-7
  * [12] Oscar – open source computer algebra research system, version 0.12.1-dev (2023), https://www.oscar-system.org
  * [13] Reidelbach, M., Ferrer, E., Weber, M.: MaRDMO Plugin - Document and Retrieve Workflows Using the MaRDI Portal. In: Proceedings of the 1st Conference on Research Data Infrastructure (CoRDI) - Connecting Communities (2023). https://doi.org/10.52825/cordi.v1i.254
  * [14] Schembera, B., Wübbeling, F., Kleikamp, H., Biedinger, C., Fiedler, J., Reidelbach, M., et al.: Ontologies for models and algorithms in applied mathematics and related disciplines. arXiv preprint arXiv:2310.20443 (2023). https://doi.org/10.48550/arXiv.2310.20443
  * [15] Schembera, B., Wübbeling, F., Koprucki, T., Biedinger, C., Reidelbach, M., Schmidt, B., et al.: Building Ontologies and Knowledge Graphs for Mathematics and its Applications. In: Proceedings of the 1st Conference on Research Data Infrastructure (CoRDI) - Connecting Communities (2023). https://doi.org/10.52825/cordi.v1i.255
  * [16] Schmitt, R.H., Anthofer, V., Auer, S., Başkaya, S., Bischof, C., Bronger, T., et al.: NFDI4Ing-the National Research Data Infrastructure for Engineering Sciences. Zenodo (2020). https://doi.org/10.5281/zenodo.4015201
  * [17] The MaRDI consortium: MaRDI: Mathematical Research Data Initiative Proposal (2022). https://doi.org/10.5281/zenodo.6552436
  * [18] Vrandečić, D.: Wikidata: a new platform for collaborative data collection. In: Proceedings of the 21st International Conference on World Wide Web. p. 1063–1064. Association for Computing Machinery (2012). https://doi.org/10.1145/2187980.2188242
  * [19] Weber, M. and Fackeldey, K.: The Mathematics of Comparing Objects. Zenodo (2022). https://doi.org/doi:10.48550/arXiv.2201.07032
  * [20] Wilkinson, M.D., Dumontier, M., Aalbersberg, I.J., Appleton, G., Axton, M., Baak, A., et al.: The FAIR Guiding Principles for scientific data management and stewardship. Scientific data 3(1), 1–9 (2016). https://doi.org/10.1038/s41597-019-0009-6.


# Multi-task Learning for Source Attribution and Field Reconstruction for
Methane Monitoring

Arka Daw § Dept. of Computer Science
Virginia Tech
Blacksburg, United States
<EMAIL_ADDRESS>Kyongmin Yeo  IBM T.J. Watson Research Center
Yorktown Heights, United States
<EMAIL_ADDRESS>Anuj Karpatne  Dept. of Computer Science
Virginia Tech
Blacksburg, United States
<EMAIL_ADDRESS>Levente Klein  IBM T.J. Watson Research Center
Yorktown Heights, United States
<EMAIL_ADDRESS>
###### Abstract

Inferring the source information of greenhouse gases, such as methane, from
spatially sparse sensor observations is an essential element in mitigating
climate change. While it is well understood that the complex behavior of the
atmospheric dispersion of such pollutants is governed by the Advection-
Diffusion equation, it is difficult to directly apply the governing equations
to identify the source location and magnitude (inverse problem) because of the
spatially sparse and noisy observations, i.e., the pollution concentration is
known only at the sensor locations and sensors sensitivity is limited. Here,
we develop a multi-task learning framework that can provide high-fidelity
reconstruction of the concentration field and identify emission
characteristics of the pollution sources such as their location, emission
strength, etc. from sparse sensor observations. We demonstrate that our
proposed framework is able to achieve accurate reconstruction of the methane
concentrations from sparse sensor measurements as well as precisely pin-point
the location and emission strength of these pollution sources.

###### Index Terms:

Physics-informed Machine Learning, Uncertainty Quantification, Climate and
Sustainability.

§§footnotetext: Work done during internship at IBM T.J.Watson Research Center.

##  I Introduction

Methane is one of the potent greenhouse gasses [6, 5, 16] that is emitted into
the atmosphere through leakages in natural gas systems, raising livestocks, or
via natural sources such as wetlands. These methane emissions caused by human
activities have further been identified as a major contributor to climate
change [3, 15, 14]. The global warming potential of methane is 25 times larger
than CO2, over a 100 year, and most emission are preventable with monitoring
and sensing solutions [2].Thus, inferring the location and emission strength
of pollution sources such as methane leaks are both essential in monitoring
the air-quality as well as mitigating the climate change. However, one of the
major challenges in localization of these emission sources is the
unavailability of high-resolution methane concentration maps, and these
individual pollution sources are to be estimated from the limited
concentration measurements from a sparse sensor network.

Digital twins are information system that create digital replication of the
physical state and temporal evolution constrained by observations and the
governing laws of physics[12]. Planetary digital twins are proposed to be
developed to better understand the long term impact of climate change [1].
Digital twins can model the complexity of the forward propagation of multi
source dispersion [8] and enable a digital playground to study the best
measurement approach to reconstruct information and identify sources. The
digital twin can be informed by sensors and satellite measurements as
realistic constrain of methane dispersion[18]. Here we demonstrate how a
dynamic, two dimensional propagation of pollution plumes [20, 19], under
realistic scenario can be used to reconstruct the propagation field and
identify sources using a data driven approach.

Here, we consider atmospheric inverse models that aim to either reconstruct
the concentration field [10] or identify pollution source information [7] of
airborne pollutants. While atmospheric inverse models have been studied
extensively, due to the ill-posed nature of the inverse problem, it still
remains as a challenging topic. Recently, potential of deep learning
approaches for the inverse problem has been demonstrated [4] [11]. We are
interested in developing a deep-learning-based atmospheric inverse model,
utilizing a digital twin to generate a range of possible emission scenarios.

We develop a multi-task learning framework that can provide a high-fidelity
reconstruction of the spatio-temporal concentration field and identify the
emission characteristics of the pollution sources such as their location, and
emission strengths from a time series of sensor measurements. We propose a
novel 3D diffusive masked convolution to gradually propagate the information
from the sensor locations to the unobserved field using a diffusion process.
We demonstrate that our multi-task model is able to achieve accurate
reconstruction of the pollution concentrations from sparse sensor measurements
as well as precisely pin-point the location and emission strength of these
pollution sources.

Figure 1: Neural Network Architecture of our proposed multi-task learning
framework.

##  II Background and Problem Setup

Here, we assume that there is a digital twin for the domain of interest,
$\Omega$. The digital twin consists of a numerical weather simulation model,
that can generate a meteorological condition, such as the wind velocity field,
$\bm{u}$, a computational model for atmospheric dispersion, and prior
information, such as the locations of the potential emission sources, possible
emission scenarios, and the locations of the sensors. We assume a two-
dimensional domain, which is discretized by equispaced $N_{x}\times N_{y}$
grid points in the east-west ($x$) and the north-south ($y$) directions.

In the digital twin, the concentration for a particular emission scenario can
be computed by solving the following advection-diffusion equation,

$\displaystyle\frac{\partial\phi}{\partial
t}+\mathbf{u\cdot}\nabla\phi-K\nabla^{2}\phi=\sum_{i=1}^{N_{s}}q_{i}(\mathbf{x},t);\>\>\mathbf{x}\in\Omega,t\in\mathbb{R}^{+}$
(1)

where $\phi$ is the concentration, $\mathbf{u}$ is the wind velocity field,
$K$ is the turbulent diffusivity, $q_{i}(\bm{x},t)$ is the emission strength
of the $i-th$ source, and $N_{s}$ denote the total number of sources. The
emission sources are point sources;
$q_{i}(x,t)=c_{i}\delta(\bm{x}-\bm{x}_{i}^{s})$, where $c_{i}$ is the emission
strength, $\bm{x}^{s}_{i}$ is the location of the source, and $\delta(.)$ is
the Dirac delta function.

Let the time-series of measurements from the $N_{o}$ sensors be
$\bm{\Phi}=\\{\bm{\Phi}_{1},...,\bm{\Phi}_{N_{o}}\\}$, where $\Phi_{i}$
denotes the time series of the measurements of the $i$-th sensor, i.e.,
$\bm{\Phi}_{i}=[\Phi_{i}^{1},\Phi_{i}^{2},...,\Phi_{i}^{N_{t}}]$, with $N_{t}$
denoting the number of time steps. The location of the $i$-th sensors are
denoted by $\mathbf{x}_{i}^{\Phi}$.

_Problem Statement_ : Given the time series of the sensor measurements
$\Phi\in\mathbb{R}^{N_{o}\times N_{t}}$, their locations
$\mathbf{x}^{\Phi}=\\{\mathbf{x}^{\Phi}_{1},\mathbf{x}^{\Phi}_{2},...,\mathbf{x}^{\Phi}_{N_{o}}\\}$,
and the wind velocity field $\mathbf{u}$, we aim to estimate the concentration
$\phi$ on the $N_{t}\times N_{x}\times N_{y}$ spatio-temporal grid.
Additionally, we are also interested in estimating the emission
characteristics of the $N_{s}$ potential emission sources, such as their
constant emission strengths $\mathbf{c}=[c_{1},c_{2},...,c_{N_{s}}]$, and
their locations
$\mathbf{x^{s}}=\\{\mathbf{x}^{s}_{1},\mathbf{x}^{s}_{2},...,\mathbf{x}^{s}_{N_{s}}\\}$.
We assume that the digital twin platform can be used to generate the data on
the $N_{t}\times N_{x}\times N_{y}$ grid for all possible emission scenarios.

##  III Proposed Method

###  III-A Multi-task learning framework

Multi-task learning [21, 17] is a learning paradigm, where the knowledge from
one task can be utilized to improve the performance of the model on other
similar tasks. Here, we aim to leverage the strong correlation between the two
tasks of identifying the source information and reconstructing the
concentration field. Note that, through the advection-diffusion equation, the
concentration field is strongly coupled with the source characteristics given
the wind velocity field $\mathbf{u}$. We propose to learn a shared encoder
$g(\boldsymbol{\theta}):[\Phi,\mathbf{x}^{\Phi},\mathbf{u}]\rightarrow\mathbf{z}$
which aims to learn a latent representation $\mathbf{z}$. Then, this latent
representation $\mathbf{z}$, is fed into task-specific decoders
$f_{1}(\boldsymbol{\omega}):\mathbf{z}\rightarrow\phi$ and
$f_{2}(\boldsymbol{\psi}):\mathbf{z}\rightarrow[\mathbf{c},\mathbf{x}^{s}]$ to
find the solution of the inverse problem and obtain the reconstructed field.

The multi-task learning framework can then be optimized using the objective
function
$\mathcal{L}(\boldsymbol{\theta},\boldsymbol{\omega},\boldsymbol{\psi})=\mathcal{L}^{\text{recon}}(\boldsymbol{\theta},\boldsymbol{\omega})+\mathcal{L}^{\text{inverse}}(\boldsymbol{\theta},\boldsymbol{\psi})$,
where $\mathcal{L}^{\text{recon}}(\boldsymbol{\theta},\boldsymbol{\omega})$ is
the loss function for reconstructing the spatio-temporal concentration field,
and $\mathcal{L}^{\text{inverse}}(\boldsymbol{\theta},\boldsymbol{\psi})$ is
the loss function for the inverse problem. The schematic of our proposed
multi-task learning framework is shown in Figure 1.

###  III-B Field Reconstruction using Diffusive Masked Convolution

We employ a series of masked convolutions for the reconstruction of the
spatio-temporal concentration field from the sparse observations. We first
define a “masked” convolution operation by restricting the convolution
operation on the masked region centered around the sensor locations where we
actually have the observations. The masked convolution operation is defined
as,

$\displaystyle x^{\prime}=\begin{cases}W^{T}(X\odot
M)\frac{\text{sum}(\mathds{1})}{\text{sum}(M)}+b,&\text{if sum}(M)>0\\\
0,&\text{otherwise}\end{cases}$ (2)

where $W$ and $b$, respectively, denote the weights and bias of a
convolutional filter, $X$ is the current feature (pixel) values for the
convolutional sliding window, $M$ is a corresponding real-valued mask between
0 and 1 (i.e., $M_{ij}\in[0,1])$, and $\odot$ represents element-wise
multiplication operation. Note that the masked convolution is defined similar
to the Partial Convolutions [9] introduced for an inpainting problem.

Then, the mask, $M$, is updated using a diffusion process after each masked
convolution layer. The diffusion process is modeled by a Spatial Gaussian
Convolutional Kernel of size $k\times k$ denoted by
$\tilde{W_{k}}\in\mathbb{R}_{+}^{k\times k}$ as follows:

$\displaystyle\tilde{W_{k}}(i,j)=\exp(-\frac{1}{2\sigma_{k}^{2}}\Big{[}(i-\frac{k-1}{2})^{2}+(j-\frac{k-1}{2})^{2}\Big{]});$
(3) $\displaystyle i,j\in\\{0,1,..,k-1\\},$

where $i,j$ are the $i$-th row and $j$-th column of the $k\times k$ kernel,
$\sigma_{k}$ is the only learnable parameter and represents the standard
deviation of the spatial kernel. In other words, if $\sigma_{k}$ is large then
the spatial-convolution kernel $\tilde{W_{k}}$ would diffuse the mask values
over a larger spatial region. Thus, the learnable-parameter $\sigma_{k}$
controls the rate of diffusion between the different layers.

The spatial kernel $\tilde{W}_{k}$ can be repeated across time once (since the
size of the kernel along time is one to ensure time-invariance) and the input
$C_{in}$ and output $C_{out}$ channels to form the time-invariant Diffusion
Kernel $W^{\text{diff}}_{k}$, and the mask is updated using the diffusion as
follows:

$\displaystyle M^{\prime}=(W^{\text{diff}}_{k})^{T}M$ (4)

After the mask update we additionally clip the mask values greater than one.
We first initialize the mask at the input layer $M_{0}$ by the sparse sensor-
network locations, as shown in Equation 5.

$\displaystyle
M_{0}(\mathbf{x})=\begin{cases}1,&\text{if}\>\>\mathbf{x}\in\mathbf{x^{\Phi}}\\\
0,&\text{otherwise}\end{cases}$ (5)

Then, we stack multiple Diffusive Masked Convolution layers so that the mask
$M_{0}$ gradually grows after each convolution layer and ultimately we can
reconstruct the full field from the sparse sensor measurements. Note that this
formulation using the Diffusive Masked Convolution can encode arbitrarily
placed sensors as the input mask. Thus, by training such a model with
different sparse sensor measurements, it is possible to predict the
concentration fields on unseen sensor-network configurations.

###  III-C Gaussian Negative Log-likelihood Formulation for Uncertainty
Quantification

In this problem setup, we are also interested in obtaining the uncertainty
estimates of the predicted concentration fields. Here, we employ a Gaussian
model, where $\phi\sim\mathcal{N}(\phi_{\mu},\sigma^{2})$ to quantify the
uncertainty in the predicted concentration field. The negative log-likelihood
loss function [13] is given as,

$\displaystyle\mathcal{L}^{\text{recon}}=\frac{1}{2}\mathbb{E}_{x\sim
p_{data}}[\log\sigma^{2}+\frac{(\phi_{true}-\phi_{\mu})^{2}}{\sigma^{2}}]$ (6)

where $\phi_{\mu}$ and $\sigma^{2}$ and denote the predicted mean and variance
of the concentration field respectively. Later, in Section IV we demonstrate
that the Gaussian negative log-likelihood loss formulation prevents over-
fitting on the ground truth concentration field and offers a smooth estimation
of the predicted mean concentration $\phi_{\mu}$.

###  III-D Emission Characteristics Estimation

Next, we consider the task of solving the inverse problem of estimating the
emission characteristics. One of the challenges is that the number of
potential sources can vary, thus, the decoder $f_{2}$ should be able to handle
this varying output size. We propose to divide the 2D spatial domain
($N_{x}\times N_{y}$) into a $S\times S$ grid. For each cell in the $S\times
S$ grid we predict the following source characteristic vector
$[p_{i},x_{i},y_{i},c_{i}]$, where $i\in\\{1,2,...,S^{2}\\}$, $p_{i}$ is the
probability of containing a source in the grid cell $i$, $x_{i}$ and $y_{i}$
represent the relative location of the source with respect to the top left
corner of the cell, and $c_{i}$ represents the emission strength of the
source.

The objective function for estimating the emission characteristics can be
evaluated as follows:

$\displaystyle\mathcal{L}^{\text{inverse}}=$
$\displaystyle\lambda_{src}\sum_{i=1}^{S^{2}}\mathds{1}_{ij}^{src}[(x_{i}^{s}-\hat{x}_{i})^{2}+(y_{i}^{s}-\hat{y}_{i})^{2}]$
$\displaystyle+\lambda_{src}\sum_{i=1}^{S^{2}}\mathds{1}_{ij}^{src}[(\sigma(\hat{p}_{i})-1)^{2}]$
$\displaystyle+\lambda_{nosrc}\sum_{i=1}^{S^{2}}\mathds{1}_{ij}^{nosrc}[(\sigma(\hat{p}_{i}))^{2}]$
$\displaystyle+\lambda_{src}\sum_{i=1}^{S^{2}}\mathds{1}_{ij}^{src}[(c_{i}-\text{Softplus}(\hat{c}_{i}))^{2}]$
(7)

where, $\mathds{1}_{i}^{src}$ denotes if the source occurs in the grid cell
$i$. Note that we use a SoftPlus activation function on the predicted emission
strengths $\text{Softplus}(\hat{c}_{i})$ to enforce that the emission
strengths are always positive, i.e., $c_{i}\in\mathbb{R}^{+}$.

##  IV Results

Figure 2: Demonstrating the reconstruction of the global field from sparse
sensor measurements on a representative test example. The white crosses denote
the position of the sensors. Top row: shows the ground truth concentration
fields at various time stamps, Bottom row: the predicted mean of the
concentration fields

###  IV-A Experiment and Evaluation Setup

_Experimental Setup_ : We performed our experiments on 4000 simulations of the
forward problem of the advection-diffusion equation, with varying source
locations and their emission strengths under stochastically generated wind
conditions. In each simulation, the number of sources is randomly selected
between one to four, and then randomly placed within the domain. These set of
simulations were there divided into training and test sets to train our
proposed multi-task learning framework. The proposed multi-task learning
framework was trained on $80\%$ of the simulations and tested on the rest.
During every epoch, for each simulation in the training set we first randomly
select the number of sensors, such that $N_{s}$ varies between 30 and 100, and
then placed them randomly inside the domain. This sampling strategy allows us
to generalize on the field reconstruction task from any arbitrary sensor-
network during the inference stage.

|  | Reconstruction
---
RMSE
Precision | Recall |  | Source Mag.
---
MSE
| Rel. Loc.
---
MSE
w/o G-NLL | 7.96 | 70.6 | 94.2 | 0.178 | 1.38e-7
w. G-NLL | 11.29 | 73.3 | 95.6 | 0.209 | 1.32e-7
TABLE I: Comparing the performance of the model trained with and without the
Gaussian Negative Log-likelihood (G-NLL) loss on different evaluation metrics.

_Evaluation Metrics_ : We define the following metrics to evaluate the
performance of our model.

  * •

_Reconstruction RMSE_ : It is defined as the Root Mean-Squared Error (RMSE)
between the predicted concentration field and the ground truth concentration
field on the test set. A lower value of RMSE indicates higher fidelity in
reconstruction of the concentration field.

  * •

_Precision and Recall_ : To evaluate the reliability of our model in detecting
the emission sources, we utilize the precision and recall metrics. For an
ideal model, both the precision and recall should be equal to one. Note that
the precision and recall are computed on the detections made w.r.t. the
$S\times S$ grid.

  * •

_Relative Location MSE_ : It is defined as the Mean Squared Error (MSE)
between the relative locations of the ground truth emission sources and the
predicted source locations inside each grid cell. Thus, the Rel. Loc. MSE
represents the efficacy of the model in localizing the sources inside the each
grid. Therefore, by combining the precision, recall and the Rel. Loc. MSE we
are able to fully evaluate the model’s capability in localizing the emission
sources.

  * •

_Source Magnitude MSE_ : It is defined as the MSE between the predicted source
magnitudes and the ground truth magnitude of the emission sources. It is used
as a metric to quantify the emission strength identification ability of the
model.

Figure 3: Localization of the pollution sources on a representative test
example. Figure 4: Emission strength estimation on a representative test
example.

###  IV-B Comparing Model Performance

Table I shows the comparison of the performance of the model trained with and
without the Gaussian Negative Log-Likelihood (G-NLL) Loss on the different
evaluation metrics.

Figure 5: Sensitivity to the number of sensors with varying training set sizes
(500-3500 simulations) Figure 6: Sensitivity to the number of pollution with
varying training set sizes (500-3500 simulations)

Figure 2 shows an example of the reconstructed concentration field on a test
dataset. The proposed diffusive masked convolution models have a
Reconstruction RMSE of 7.96 when trained without the Gaussian NLL loss, but
the Reconstruction RMSE increases to 11.29 when trained with the Gaussian NLL
loss. The larger RMSE of G-NLL is mainly due to the misfit at the emission
source location. In G-NLL, the estimated concentration field becomes smooth
and the delta-function-like behavior at the source location is treated as a
noise. On the other hand, without G-NLL, the model tries to follow the delta-
function-like behavior, which results in an artificial oscillation in the
solution.

We also demonstrate a fairly high recall of around $95\%$, and a precision of
$70\%$ for both models. This suggests that the model is able to recover $95\%$
of the emission sources that were present in the simulations. However, the low
precision suggests that the model also generates about $30\%$ false positives
in its predictions. Further, we observe that once a grid cell detects a
emission source, both models can exactly pin-point the location of the source
inside the grid with a Rel. Loc. MSE of around $1.3e-7$.

We also show an accurate estimation of the emission strength for the sources
as two models display a Source Mag. MSE of 0.18 and 0.21 respectively. Figure
3 and 4 shows the ability of our proposed framework to accurately identify the
pollution sources as well as approximate their emission strengths on an
example case.

###  IV-C Analyzing Sensitivity to the Number of Sensors

Next, we analyze the sensitivity of our model performance on the number of
sensors. We vary the number of sensors from 20 to 110 and compare the
performance of our model on various training sizes ranging from 500
simulations to 3500 simulations. From Figure 5, it can be seen that all of the
evaluation metrics gradually improve as the size of the training set
increases. The reconstruction RMSE decreases as the number of sensors
increase. This trend is expected as having more sensors observations make the
reconstruction problem less ill-posed, thereby making it easier to obtain
higher fidelity field reconstructions. We also observe that the with lower
number of sensors the precision decreases significantly while the recall
remains the same. This suggests that with less sensors the model is unable to
exactly identify the location of the sources and ends up predicting a large
number of false-positives that decrease the precision. Additionally, having
more sensors improves the Source Mag. MSE, suggesting that the model is able
to better estimate the source magnitudes (i.e., solve the inverse problem
better) with more sensor measurements.

Figure 7: Analyzing the empirical coverage of the 95% Prediction Interval
(shown as the shaded region) for the model trained on Gaussian Negative Log-
likelihood (NLL) loss for a particular spatial slice of $y=30$ and $y=37$ (a
line passing through 3 pollution sources, that can be identified using the
peaks in the ground truth) on an example case.

###  IV-D Analyzing Sensitivity to the Number of Emission Sources

Similarly, we analyze the sensitivity of the trained model on the number of
emission sources by varying them from one to five sources. From Figure 6, we
can see that the Reconstruction RMSE linearly increases with the number of
sources. This is because the magnitude of the concentration field is
proportional to the number of emission sources, and thus the same is reflected
on the reconstruction error. We also observe that the precision improves as
the number of sources increases and the recall remains almost unaffected.
Further, the source magnitude MSE increases with the number of sources, as it
becomes difficult to disentangle the contributions from multiple emission
sources.

###  IV-E Effects of the NLL formulation

We further qualitatively analyze the predictive intervals of the model trained
with the Gaussian Negative Log-likelihood loss. For a fixed value of the
y-coordinate ($y=30$ and $y=37$) we plot the variations in the concentration
fields w.r.t. the x-axis at various times in Figure 7. We observe that the
predictive intervals estimated by our model always engulf the ground truth
concentration. Additionally, we show that the Gaussian-Negative Loglikelihood
(NLL) formulation for estimating uncertainties allows the estimate of the mean
concentration field to be a smooth function, as the large fluctuations near
the pollution sources are taken care by the predicted variance, thereby
preventing an oscillatory solution. This behavior can also be seen in Figure
8, where we see that the absolute error without the NLL shows rectangular
oscillatory patterns, which is not observed for the model trained with NLL.

Figure 8: Comparing the Absolute Errors of the model trained with and without
the Gaussian Negative Log-likelihood (NLL) loss at different time intervals on
an example test case.

##  V Conclusion

We present a multi-task learning framework for identifying potential emission
sources and obtaining reconstruction of spatio-temporal concentration fields
from sparse sensor measurements. We also propose a novel diffusive-masked
convolution operations that employs the diffusion process in performing
masked-convolutions. The diffusive-masked convolution is realized by a
spatial-Gaussian convolution kernel, followed by diffusing the mask to nearby
regions. Thus, by stacking multiple such layers we are able to iteratively
diffuse the information from a sparse sensor measurements (represented using
the initial mask) in a principled manner until the learned representations
spread over the entire spatio-temporal domain. We also demonstrate precise
reconstruction of the concentration field along with accurate localization and
emission strength estimation of the pollution sources on the test simulations
using our proposed framework.

## References

  * [1] Peter Bauer, Bjorn Stevens, and Wilco Hazeleger.  A digital twin of earth for the green transition.  Nature Climate Change, 11(2):80–83, 2021.
  * [2] Olivier Boucher, Pierre Friedlingstein, Bill Collins, and Keith P Shine.  The indirect global warming potential and global temperature change potential due to methane oxidation.  Environmental Research Letters, 4(4):044007, 2009.
  * [3] Mingkui Cao, Keith Gregson, and Stewart Marshall.  Global methane emission from wetlands and its sensitivity to climate change.  Atmospheric environment, 32(19):3293–3299, 1998.
  * [4] Kai Fukami, Romit Maulik, Nesar Ramachandra, Koji Fukagata, and Kunihiko Taira.  Global field reconstruction from sparse sensors with voronoi tessellation-assisted deep learning.  Nature Machine Intelligence, 3(11):945–951, 2021.
  * [5] Robert W Howarth.  A bridge to nowhere: methane emissions and the greenhouse gas footprint of natural gas.  Energy Science & Engineering, 2(2):47–60, 2014.
  * [6] Robert W Howarth, Renee Santoro, and Anthony Ingraffea.  Methane and the greenhouse-gas footprint of natural gas from shale formations.  Climatic change, 106(4):679–690, 2011.
  * [7] Youngdeok Hwang, Hang J Kim, Won Chang, Kyongmin Yeo, and Yongku Kim.  Bayesian pollution source identification via an inverse physics model.  Computational Statistics & Data Analysis, 134:76–92, 2019.
  * [8] Anna Karion, Thomas Lauvaux, Israel Lopez Coto, Colm Sweeney, Kimberly Mueller, Sharon Gourdji, Wayne Angevine, Zachary Barkley, Aijun Deng, Arlyn Andrews, et al.  Intercomparison of atmospheric trace gas dispersion models: Barnett shale case study.  Atmospheric chemistry and physics, 19(4):2561–2576, 2019.
  * [9] Guilin Liu, Fitsum A Reda, Kevin J Shih, Ting-Chun Wang, Andrew Tao, and Bryan Catanzaro.  Image inpainting for irregular holes using partial convolutions.  In Proceedings of the European conference on computer vision (ECCV), pages 85–100, 2018.
  * [10] Xiao Liu, Kyongmin Yeo, and Siyuan Lu.  Statistical modeling for spatio-temporal data from stochastic convection-diffusion processes.  Journal of the American Statistical Association, 117(539):1482–1499, 2022.
  * [11] Mirco Milletarı, Sara Malvar, Yagna D Oruganti, Leonardo O Nunes, Yazeed Alaudah, and Anirudh Badam.  Source attribution and emissions quantification for methane leak detection: A non-linear bayesian regression approach.  2020\.
  * [12] Steven A. Niederer, Michael S. Sacks, Mark Girolami, and Karen Willcox.  A digital twin of earth for the green transition.  Nature Computational Science, 1(5), 2021.
  * [13] David A Nix and Andreas S Weigend.  Estimating the mean and variance of the target probability distribution.  In Proceedings of 1994 ieee international conference on neural networks (ICNN’94), volume 1, pages 55–60. IEEE, 1994.
  * [14] Fiona M O’Connor, O Boucher, N Gedney, CD Jones, GA Folberth, R Coppell, P Friedlingstein, WJ Collins, J Chappellaz, J Ridley, et al.  Possible role of wetlands, permafrost, and methane hydrates in the methane cycle under future climate change: A review.  Reviews of Geophysics, 48(4), 2010.
  * [15] Dave Reay and Pete Smith.  Methane and climate change.  Routledge, 2010.
  * [16] ED Schulze, S Luyssaert, P Ciais, A Freibauer, IA Janssens, Jean-François Soussana, P Smith, John Grace, I Levin, B Thiruchittampalam, et al.  Importance of methane and nitrous oxide for europe’s terrestrial greenhouse-gas balance.  Nature geoscience, 2(12):842–850, 2009.
  * [17] Trevor Standley, Amir Zamir, Dawn Chen, Leonidas Guibas, Jitendra Malik, and Silvio Savarese.  Which tasks should be learned together in multi-task learning?  In International Conference on Machine Learning, pages 9120–9132. PMLR, 2020.
  * [18] Theodore G Van Kessel, Levente J Klein, Muralidhar Ramachandran, Eric J Zhang, and Hendrik Hamann.  Satellite guided mobile wireless methane detection for oil and gas operations.  In 2020 IEEE 6th World Forum on Internet of Things (WF-IoT), pages 1–4. IEEE, 2020.
  * [19] Chulin Wang, Eloisa Bentivegna, Wang Zhou, Levente Klein, and Bruce Elmegreen.  Physics-informed neural network super resolution for advection-diffusion models.  preprint arXiv:2011.02519, 2020.
  * [20] Chulin Wang, Kyongmin Yeo, Xiao Jin, Andres Codas, Levente J. Klein, and Bruce Elmegreen.  S3rp: Self-supervised super-resolution and prediction for advection-diffusion process.  arXiv preprint arXiv:2111.04639, 2021.
  * [21] Yu Zhang and Qiang Yang.  A survey on multi-task learning.  IEEE Transactions on Knowledge and Data Engineering, 2021.