# Reconstruction of Primordial Power Spectrum of curvature perturbation from
the merger rate of Primordial Black Hole Binaries

Rampei Kimura1, Teruaki Suyama2, Masahide Yamaguchi2 and Ying-li Zhang3,4
1Waseda Institute for Advanced Study, Waseda University, 1-6-1 Nishi-Waseda,
Shinjuku, Tokyo 169-8050, Japan
2Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8551, Japan
3School of Physics Science and Engineering, Tongji University, Shanghai
200092, China
4Institute for Advanced Study of Tongji University, Shanghai 200092, China

###### Abstract

The properties of primordial curvature perturbations on small scales are still
unknown while those on large scales have been well probed by the observations
of the cosmic microwave background anisotropies and the large scale structure.
In this paper, we propose the reconstruction method of primordial curvature
perturbations on small scales through the merger rate of binary primordial
black holes, which could form from large primordial curvature perturbation on
small scales.

###### pacs:

04.50.Kd,95.30.Sf,98.80.-k

##  I Introduction

The idea of Primordial Black Holes (PBHs) can be traced back to 1960s, which
was firstly investigated by Zel’dovich and Novikov Zeldovich:1967 . Later, in
1971, Hawking pointed out that in the primordial universe, black holes (BHs)
can be formed by the gravitational collapse of a highly overdense region of
inhomogeneities Hawking:1971ei . Since then, the studies of PBHs have
attracted a lot of attention mainly in two aspects. Firstly, since PBHs only
interact via gravity, without introducing new physics beyond the standard
model, they are natural candidates for dark matter (DM) which comprises 25$\%$
of the critical density (the earliest idea can be found in 1970
Chapline:1975ojl ). Most of the constraints on the fraction of PBHs in dark
matter, come from the gravitational lensing and dynamical constraints.
However, recent reanalysis opened a parameter window where PBHs with mass
$10^{-16}-10^{-10}M_{\odot}$ can constitute all of the DM Niikura:2017zjd ;
Katz:2018zrn ; Smyth:2019whb ; Montero-Camacho:2019jte . This refreshes the
interest of the investigation of PBHs as the origin of DM Pi:2017gih ;
Cai:2018tuh . Secondly, the detected binary black hole (BBH) mergers during
the past O1, O2 and first half of O3 runs of the Laser Interferometer
Gravitational-wave Observatory (LIGO) Abbott:2016blz ; Abbott:2016nmj ;
TheLIGOScientific:2016pea ; Abbott:2017vtc ; Abbott:2017oio ; Abbott:2017gyy ;
LIGOScientific:2018mvr ; Venumadhav:2019lyq ; Abbott:2020niy  have shown two
unexpected characteristics of binary black hole (BBH): the high masses (larger
than 20 $M_{\odot}$) and low effective spin. The origin of these heavy BHs and
the formation of such BBHs which merge within the age of the Universe are
still in debate Belczynski:2010tb ; Dominik:2012kk ; Dominik:2013tma ;
Dominik:2014yma ; Belczynski:2016obo  (for a comprehensive summarization, see
e.g. TheLIGOScientific:2016htt ; Miller:2016krr ). One of the possible
explanations is that the mergers could be primordial in origin Bird:2016dcv ;
Clesse:2016vqa ; Sasaki:2016jop ; Sasaki:2018dmp ; Garcia-Bellido:2020pwq .

Up to now, most of works concentrate on predicting the PBH construction and
its associated BBH mergers for (given) primordial curvature perturbations.
However, in near future, we expect to detect much more BBH mergers. In this
new era, the following (opposite) question will become more and more
important: How can we reconstruct the features of primordial curvature
perturbations such as power spectra and bispectra from the observations of BBH
mergers? In fact, if (some of) these BBH mergers could be attributed to PBHs,
we will be able to reconstruct the features of primordial curvature
perturbations which would predict PBHs and then such PBHs are responsible for
the detected BBH mergers. This reconstruction program is quite important
because, over the past two decades, on cosmological scales (larger than Mpc
scale), the observations of the cosmic microwave background (CMB) anisotropies
and the large scale structure indicate an almost scale-invariant primordial
curvature power spectrum with accurately measured amplitude and tilt
Spergel:2003cb ; Aghanim:2018eyx ; Akrami:2018odb . However, on small scales
(smaller than Mpc scale), we still lack the information of primordial
curvature power spectrum though some upper bounds are obtained Chluba:2019kpb
; Jeong:2014gna ; Nakama:2014vla ; Inomata:2016uip ; Allahverdi:2020bys . This
reconstruction program could directly determine the information of primordial
curvature perturbations on small scales. As far as we know, no such attempt
has yet been done albeit its importance.

In this paper, we will organize such a reconstruction program. Our finding is
that the reconstruction program must go through three distinct steps to
successfully reconstruct the primordial power spectrum. In step (i), one needs
to relate the feature of primordial curvature perturbations (such as power
spectrum) to the (smoothed) density contrast (characterized by the mass
variance for example). In step (ii), the PBH mass function needs to be related
to the (smoothed) density contrast. In step (iii), the PBH mass function needs
to be related to the (would-be) observed rates of BBH mergers associated with
PBHs. In each step, there should be one-to-one correspondence for this
reconstruction program to work. As we will explain in details later, each step
suffers from its own theoretical uncertainties. They stem from i)very rare and
highly non-linear nature of the PBH formation which defeats a first-principle
calculation and forces us to employ some phenomenological modeling where some
ambiguities come into play, and ii)our ignorance of the model of the early
universe. As a result, it is inevitable to adopt some assumptions to achieve
the reconstruction of the primordial power spectrum. In the next section, we
will explain those assumptions one by one. For the technical assumption, we
will also discuss how we should change the reconstruction procedure when the
assumption is relaxed.

This paper is organized as follows. In Section II, we state the underlying
four assumptions for the calculation and describe the motivations for these
assumptions. In Section III, we divide the reconstruction procedure into three
steps. In Section III.1, we relate the power spectrum of the Gaussian
primordial curvature perturbations to the variance of the density perturbation
smoothed by the top-hat window function. In Section III.2, we discuss the
correspondence between the variance and the PBH mass function with the
critical collapse phenomena being taken into account. In Section III.3, by
adopting the most updated formula of the merger rate, we finally show how to
relate the PBH mass function to the (would-be) observed rates of BBH mergers
associated with PBHs, assuming some of the observed binary BH merger events
are attributed to PBH mergers, and show the PBH mass function can be
reconstructed by putting these three steps together once the merger rate with
sufficient mass resolution is determined. As an example, in Section IV, we
consider the simplest case where the effects of critical collapse and
suppression factor are neglected and then find the one-to-one correspondence
between the mass function and the merger rate in an analytic way. In Section
V, we draw the conclusion of this paper and give some discussions for further
extension. In Appendix A, for completeness, we make a general mathematical
discussion on the one-to-one correspondence between the variance and the
primordial power spectrum with an arbitrary window function.

##  II Underlying assumptions

In this section, we state the following four assumptions that will be used in
the next section:

(0)

In the future, some of the observed binary BH merger events are attributed to
mergers of PBHs formed, and the merger rate distribution of the PBHs in the
mass plane is observationally determined.

(1)

PBHs were formed out of rare high-$\sigma$ peaks of the primordial curvature
perturbations in radiation dominated era.

(2)

The window function takes the top-hat form in $k$-space;

(3)

The primordial curvature perturbation follows Gaussian distribution and the
effect of non-linearity between the curvature and the density perturbations is
ignored;

About the assumption (0)
The origin of the BHs detected by LIGO/Virgo experiments remains an open
question. Broadly, there are two astrophysical channels of the BBH formation
which could possibly explain the observed merger events (for instance, see
Mapelli:2018uds ). One is the isolated binary formation channel in which the
massive stars in a binary evolve into the BH binary after experiencing the
common envelope phase through which the binary shrinks its orbit and ends up
with a tightly bound system. The other is the dynamical formation channel in
which the two BHs form a binary by dynamical encounters which occurs in dense
environments such as globular clusters. Another intriguing possibility is the
PBH scenario, which could explain not only the observed binary black hole
(BBH) merger events, but also the fact that many observed BBHs have large mass
and low spin Abbott:2020niy . Given the large theoretical uncertainties in
predicting the merger rate distribution of the individual astrophysical
channels, it is a challenging task to determine how much the PBHs contribute
to the observed merger rate distribution from the data (see Kocsis:2017yty ;
Liu:2018ess ; Chen:2019irf ; Garriga:2019vqu ; Clesse:2020ghq ; Hall:2020daa ;
Wong:2020yig  as recent works along this direction). In order to achieve this,
we would ultimately need to reduce the theoretical uncertainties of the
astrophysical channels and subtract the astrophysical contributions from the
observationally determined merger rate distribution to identify the PBH
signal. Whether the merger rate distribution obtained in such a way is
consistent with the mergers of PBHs can be checked by using the consistency
relation, which is briefly discussed in the last section in this paper (see
Kocsis:2017yty ; Liu:2018ess ; Garriga:2019vqu  as relevant previous studies).
Alternatively, the issue of discriminating from the astrophysical BHs can be
circumvented if merger events at cosmological redshifts $z\gtrsim 20$
Nakamura:2016hna  or of sub-solar mass BHs are detected since astrophysical
BHs are thought to be very rare at such high redshifts and heavier than the
solar mass (see Authors:2019qbw  which conducted search for sub-solar PBH
mergers.). In these cases, it would be straightforward to observationally
determine the PBH merger rate distribution. In this paper, we simply adopt the
assumption (0) and suppose that one could know the merger rate of BBH
originating from the PBHs because the main purpose of this paper is to give
the method of the reconstruction of primordial curvature perturbations from
such merger rates.

About the assumption (1)
Though various mechanism to produce PBHs such as cosmic topological defects
Hawking:1982ga ; Crawford:1982yz ; Hawking:1987bn ; Moss:1994iq , interacting
dark matter clumps Shandera:2018xkn  have been proposed so far, one of the
most interesting possibilities is large curvature perturbations generated
during inflationary expansion Ivanov:1994pa ; GarciaBellido:1996qt ;
Bullock:1996at ; Ivanov:1997ia ; Cai:2019jah ; Palma:2020ejf ; Cheong:2019vzl
; Zhou:2020kkf ; Chen:2020uhe . In fact, through this mechanism of large
curvature perturbations, a lot of works have been done to “predict” PBHs,
which might account for dark matter, MACHOs, (unidentified) cosmic rays, and
so on. In these works, given a specific model of inflation, primordial
curvature perturbations were predicted, and then, the formation of PBHs
originating from their designed power spectra has been discussed. In
principle, it is possible that early matter dominated era existed prior to the
radiation dominated and PBHs were formed in such matter dominated era.
Formation process of the PBHs in such era is very different from that in the
radiation dominated era Khlopov:1980mg ; Harada:2016mhb . Given that the PBH
formation in the radiation dominated era is the most standard possibility
considered in the literature, we assume that PBHs originate from the
primordial curvature perturbations in the radiation dominated era.

About the assumption (2)
The assumption (2) is imposed in order to obtain the analytic relation between
the primordial power spectrum and the variance (see (7)). Of course there is
no physical principle to choose a certain form of the window function. It is
evident that different form of window function will change the final result,
but it has been shown in Ando:2018qdb ; Young:2019osy  that the effect of the
choice of window function will cause the uncertainty in the amplitude of the
power spectrum only up to $\mathcal{O}(10\%)$. Hence, we concentrate on the
top-hat form as a typical example in the main text. A numerical strategy to
reconstruct the primordial power spectrum for other types of the window
function is also briefly discussed in the appendix.

About the assumption (3)
The primordial curvature perturbation is non-Gaussian in general. In the
standard single-field slow-roll inflation models, the non-Gaussianity is
suppressed by the slow-roll parameters. For other classes of inflation models,
the level of non-Gaussianity as well as the types of non-Gaussianity are
model-dependent. If the primordial curvature perturbations are non-Gaussian,
the PBH mass function and the PBH abundance depend not only on the power
spectrum but also on higher-order statistics. In other words, the
reconstruction of the power spectrum from the observations of the PBH mergers
does not work unless we make specific assumptions on the higher-order
statistics. The higher-order statistics will become important if the
probability to realize the overdensity to form a PBH, which corresponds the
tail of the probability density in realistic situation, differs significantly
from the one given by the Gaussisan distribution. For instance, for the
curvature perturbation with the local-type non-Gaussianity parametrized by the
$f_{\rm NL}$ parameter, the non-Gaussianity will change the PBH abundance from
the Gaussian case if $|f_{\rm NL}|\gtrsim\sigma^{2}/\delta_{\rm th}^{3}$,
where $\sigma^{2}$ is the variance of the perturbation and $\delta_{\rm th}$
is the threshold for the PBH formation. Since there is no standard non-
Gaussian shape which is well-motivated than the other non-Gaussian shapes, we
assume that the curvature perturbation is Gaussian. Moreover, even if we
assume the Gaussianity of the curvature perturbation, the non-linearity
intrinsic to GR makes the density perturbation, which is relevant to the PBH
formation, non-Gaussian. We also neglect the effect of this non-Gaussianity
for simplicity, which is the assumption (3). 111The effect of non-Gaussianity
can be found in literatures, for example, Ref Bullock:1996at ; Ivanov:1997ia ;
PinaAvelino:2005rm ; Lyth:2012yp ; Byrnes:2012yx ; Shandera:2012ke ;
Young:2013oia ; Franciolini:2018vbk ; Cai:2018dig ; Yoo:2019pma ; Atal:2019erb
; Ezquiaga:2019ftu ; Young:2014oea ; Young:2015cyn ; DeLuca:2019qsy .

The assumption (0) is mandatory (almost by definition) for the reconstruction
program to be achieved successfully. In this sense, this assumption is crucial
in this paper. The assumption (2) is not essential to the reconstruction
although this is practically convenient from the technical point of view. It
is not clear if the violation of the assumption (1) and(or) the assumption (3)
can spoil the reconstruction program. Answer may depend on the early universe
model for the PBH formation and the type of non-Gaussianity of the primordial
curvature perturbations that we employ. In that case, one needs to check the
presence of the one-to-one correspondence in each step given below. Addressing
this issue is beyond the scope of this paper.

##  III The principle of a reconstruction program of primordial curvature
perturbations

Under the assumptions (0)–(3), we describe the principle of the reconstruction
program. The logic chain of the reconstruction program is given as follows:

$\displaystyle\mathcal{P}_{\mathcal{R}}(k)\stackrel{{\scriptstyle\rm(i)}}{{\Longleftrightarrow}}\sigma^{2}(R)\stackrel{{\scriptstyle\rm(ii)}}{{\Longleftrightarrow}}f(m)\stackrel{{\scriptstyle\rm(iii)}}{{\Longleftrightarrow}}\mathcal{R}(m_{1},m_{2},t)\,,$
(1)

where $\mathcal{P}_{\mathcal{R}}(k)$ is the power spectrum of the primordial
curvature perturbations, $\sigma^{2}(R)$ is the variance of the density
perturbation smoothed over a comoving length scale $R$, $f(m)$ is the PBH mass
function, and $\mathcal{R}(m_{1},m_{2},t)$ is the merger rate density of the
PBH binaries with the masses $m_{1}$ and $m_{2}$ at the cosmic time $t$. Their
detailed definitions will be given later. In what follows, we explain these
three steps one by one.

###  III.1 Step (i)

In radiation dominated universe, an initially overdense region collapses to a
BH right after the horizon reentry if the average amplitude of the density
contrast of the overdense region evaluated at the time of the horizon reentry
is larger than the threshold. In the literature, the density contrast $\Delta$
on the comoving slice is sometimes used as a useful quantity to measure the
PBH formation Young:2014ana . The density contrast smoothed over the comoving
scale $R$ is defined by

$\displaystyle\Delta_{R}(t,{\bf x})=\int W({\bf y}-{\bf x};R)\Delta(t,{\bf
y})d{\bf y}.$ (2)

Here $W({\bf y};R)$ is the window function. Notice that $\Delta$ evolves as
$\propto a^{2}$ on super-Hubble scales. What is relevant to the formation of
PBH with mass $M$ is the variance of $\Delta_{R}$,

$\displaystyle\sigma^{2}(R)=\langle\Delta_{R}^{2}(t_{R},{\bf
x})\rangle=\int_{0}^{\infty}W^{2}(kR)\mathcal{P}_{\Delta}(t_{R},k)~{}\mathrm{d}(\ln
k),$ (3)

where $R$ is determined as a function of $M$ by requiring that $M$ is the
horizon mass evaluated at the time $t=t_{R}$ (i.e. $a(t_{R})R=1/H(t_{R})$)
when the scale $R$ reenters the Hubble horizon, and $W(kR)$ is the Fourier
transformation of the window function. In order to connect $\sigma^{2}(R)$ to
${\cal P}_{\cal R}(k)$, we use the relation between the power spectrum for the
density contrast $\mathcal{P}_{\Delta}$ and the one for the curvature
perturbation during the radiation-dominated era given by
$\mathcal{P}_{\mathcal{R}}$ as (see e.g., Green:2004wb )

$\displaystyle\mathcal{P}_{\Delta}(t,k)=\frac{16}{81}\left(\frac{k}{aH}\right)^{4}\mathcal{P}_{\mathcal{R}}(k).$
(4)

The right-hand side of this relation is the leading order in the Taylor
expansion in the power of $k/(aH)$ and higher order terms are ignored. We use
this relation in evaluating the right-hand side of Eq. (3), which will be a
good approximation since the unphysical contributions from the sub-Hubble
modes are suppressed thanks to the window function. Then, we have

$\displaystyle\sigma^{2}(R)=\frac{16}{81}\int_{0}^{\infty}W^{2}(kR){(kR)}^{4}{\cal
P}_{\cal R}(k)~{}\mathrm{d}(\ln k).$ (5)

Based on the assumption (2), we adopted the top-hat form window function in
$k$-space

$\displaystyle W(kR)=\left\\{\begin{matrix}\displaystyle\quad
1\quad\,;&0<k<1/R\,,&~{}\\\ \\\ \displaystyle\quad
0\quad;&\mathrm{otherwise}\,,\end{matrix}\right.$ (6)

Taking the derivative of (5) with respect to $1/R$, we then obtain

$\displaystyle\frac{81}{16}R^{4}\frac{{\rm d}}{{\rm
d}R}\left(\frac{\sigma^{2}(R)}{R^{4}}\right)\bigg{|}_{1/R=k}=-k\mathcal{P}_{\cal
R}(k).$ (7)

This is the desired relation which enables us to directly reconstruct the
power spectrum of the curvature perturbation $\mathcal{P}_{\mathcal{R}}(k)$
from (the derivative of) the variance $\sigma^{2}(R)$. This reconstruction
method works in an ideal case where $\sigma^{2}(R)$ is determined by
observations as a continuous function of $R$. In realistic observations,
however, $\sigma^{2}(R)$ will be measured only for some discrete values of $R$
unless we adopt a specific fitting function for which free parameters entering
the fitting function are fixed by observations. Furthermore, it is not clear
how to extend the reconstruction method given above to other types of the
window functions. Although the top-hat form in $k$-space is technically
conveninent in the sense that it is straightforward to solve the integral
equation (5) for ${\cal P}_{\cal R}$, this window function may not be
considered as the most natural one from the physical point of view given that
the top-hat window function, when Fourier-transformed into the real space, is
spatially extended beyond the Hubble horizon and thus includes contribution of
the perturbations outside the Hubble horizon which should not affect the PBH
formation. In the appendix, we discuss about how to reconstruct ${\cal
P}_{\cal R}$ when $\sigma^{2}(R)$ is given only discretely for a fixed window
function which is not necessarily a top-hat form in $k$-space.

On the other hand, the variance $\sigma^{2}(R)$ can be uniquely determined
from $\mathcal{P}_{\mathcal{R}}(k)$ by (5). Therefore, the step (i) in (1) is
reversible.

###  III.2 Step (ii)

Our next task is to reconstruct the variance $\sigma^{2}(R)$ from the PBH mass
function $f(M)$. In this paper, the PBH mass function $f(M)$ is defined such
that the quantity $f(M){\rm{d}}M$ represents the probability that a randomly
chosen PBH is in the mass range $\left(M,M+{\rm{d}}M\right)$,

$\displaystyle\int_{0}^{\infty}f(M){\rm{d}}M=1.$ (8)

Computing the PBH mass function from the primordial perturbations has been a
long-standing topic. Since the formation of PBHs is a very rare and highly
non-linear phenomena, there is no first principle derivation of the PBH mass
function 222One may think the PBH mass function can be obtained by conducting
the cosmological simulations which evolve the random primordial fluctuations
initially defined on super-Hubble scales until all the interesting scales have
reentered the Hubble horizon and by counting the number of the produced PBHs
and measuring the mass distribution. However, since the PBH formation is very
rare, the size of the simulation box must be extraordinarily large to contain
enough number of PBHs at the end of simulation. This is impossible at the
current computational resource.. In the literature, semi-analytic formulations
of the PBH mass function based on the numerical relativity computations of the
PBH formation in the spherically symmetric system and its phenomenological
extrapolation to the more general random configuration without spherical
symmetry have been addressed and refined repeatedly. In this paper, we adopt
the mass function given in Byrnes:2018clq ; Wang:2019kaf  which connects the
variance $\sigma^{2}(R)$ with the PBH mass function $f(M)$ with the effect of
the critical phenomena being included such that
$M=KM_{R}\left(\Delta-\Delta_{\rm th}\right)^{\gamma}$. It is given by

$\displaystyle f(M)$ $\displaystyle=\frac{1}{f_{\rm
PBH}M}\int_{-\infty}^{\infty}\frac{K}{\sqrt{2\pi}\gamma\sigma(R)}\left(\frac{M_{\rm
eq}}{M_{R}}\right)^{\frac{1}{2}}{\left(\frac{M}{KM_{R}}\right)}^{1+\frac{1}{\gamma}}\exp\left[-\frac{1}{2\sigma^{2}(R)}{\left(\Delta_{\rm
th}+{\left(\frac{M}{KM_{R}}\right)}^{\frac{1}{\gamma}}\right)}^{2}\right]\frac{{\rm
d}M_{R}}{M_{R}}\,,$ (9)

where $f_{\rm PBH}$ is the mass fraction of PBHs in dark matter, $M_{\rm eq}$
is the horizon mass at the matter-radiation equality epoch, and the
integration variable $M_{R}$ is the horizon mass at $t=t_{R}$,

$M_{R}=\frac{4\pi}{3}\rho(t_{R})H^{-3}(t_{R})=\frac{1}{2GH(t_{R})},$ (10)

by which $R$ is related to $M_{R}$. This integration arises since the
overdense region with fixed $R$ can collapse to the PBH of any mass in the
range $0<M\lesssim M_{R}$ when the critical phenomena is taken into account.
$\gamma\approx 0.36$ is the critical exponent of the critical collapse,
$\Delta_{\rm th}$ is the threshold of the PBH formation, which depends both on
the perturbation profile and the adopted window function Young:2019osy , and
$K=3.3$ is a numerical constant.

Eq. (9) is a non-linear integral equation for $\sigma^{2}(R)$ and it is
impossible to find analytic solution. Furthermore, unfortunately, we were not
able to find the numerical method to systematically reconstruct
$\sigma^{2}(R)$ for any given $f(M)$. This is a purely mathematical problem
and providing a full solution to this problem is beyond the scope of this
paper. Instead, here we sketch one ad hoc approach which may work sufficiently
in some cases. To this end, we first notice from Eq. (9) that, from the
physical point of view, the mass function $f(M)$ is expected not to decay
faster than $M^{1/\gamma}$ as $M$ is decreased. In order to see this, let us
assume an extreme case where $\sigma^{2}(R)=0$ for $R<R_{*}$. In this case,
the range of integration in Eq. (9) is effectively restricted to
$M_{R}>M_{R_{*}}$. Thus, for $M\ll M_{R_{*}}$, Eq. (9) becomes

$f(M)\approx\frac{M^{\frac{1}{\gamma}}}{f_{\rm PBH}}\int_{\ln
M_{R_{*}}}^{\infty}\frac{K}{\sqrt{2\pi}\gamma\sigma(R)}\left(\frac{M_{{\rm
eq}}}{M_{R}}\right)^{\frac{1}{2}}{\left(KM_{R}\right)}^{-1-\frac{1}{\gamma}}\exp\left(-\frac{\Delta_{\rm
th}^{2}}{2\sigma^{2}(R)}\right){\rm d}\ln M_{R}.$ (11)

The integrand is independent of $M$ and thus $f(M)\propto
M^{\frac{1}{\gamma}}$ (For the case of the monochromatic power spectrum, this
feature was derived in Niemeyer:1997mt  for the Gaussian perturbation and was
shown to hold for the general non-Gaussian perturbation as well in
Yokoyama:1998xd ). This tail $f(M)\propto M^{\frac{1}{\gamma}}$ reflects the
fact that small-mass PBHs ($M\ll M_{R_{*}}$) are inevitably produced by the
critical phenomena. In this sense, this gives the minimum amount of PBHs even
in the absence of the primordial perturbations at the scale corresponding to
$M(\ll M_{R_{*}})$. In reality, $\sigma^{2}(R)$ will not be exactly zero for
$R<R_{*}$, and such additional perturbations will produce extra PBHs to the
ones mentioned above, which would only make $f(M)$ decrease slower than
$M^{\frac{1}{\gamma}}$. Thus, $f(M)$ is not expected to decay faster than
$M^{\frac{1}{\gamma}}$ as we decrease $M$. This consideration tells us that if
the observed mass function turns out to decay faster than this limit, it
strongly disfavors the scenario of the PBH origin, which itself provides one
universal test of the PBH scenario. Since our main purpose in this subsection
is the reconstruction of $\sigma^{2}(R)$ from $f(M)$ under the assumption that
the observed BHs have been robustly identified as PBHs and their mass function
$f(M)$ has been well measured, we hereafter suppose that the observed $f(M)$
is consistent with this limit.

If the observed $f(M)$ shows the scaling behavior $\propto
M^{\frac{1}{\gamma}}$ in a certain mass range, then the above argument
suggests that $\sigma^{2}(R)$ in the corresponding range of $R$ is negligibly
small compared to that for larger $R$. If, on the other hand, $f(M)$ is
enhanced than the lower limit $M^{\frac{1}{\gamma}}$ for a wide mass range, it
implies that such PBHs were produced by nearly scale-independent
$\sigma^{2}(R)$ in the same range since $f(M)$ is exponentially sensitive to
$\sigma^{2}(R)$ (In any realistic situation, $\sigma^{2}(R)$ must be much
smaller than $\Delta_{\rm th}^{2}$ in order to avoid the overproduction of
PBHs). Fig. 1 shows the integrand of Eq. (9) as a function of
$\frac{KM_{R}}{M}$. We find that the integrand is sharply peaked at

$\frac{KM_{R}}{M}\approx\left(\frac{2\Delta_{\rm
th}}{(2+5\gamma)\sigma^{2}(R)}\right)^{\gamma}.$ (12)

In the figure, the rapid fall-off on the left side is explained by the fact
that the production of PBHs on super-Hubble scales does not occur, and the
decay on the right side is due to the critical phenomena. Then, applying the
method of the saddle point to Eq. (9), we can approximately evaluate the
integral as

$f(M)\approx f_{\rm
PBH}^{-1}{\left(\frac{K}{M}\right)}^{\frac{3}{2}}\left(\frac{M_{\rm
eq}}{\Delta_{\rm
th}}\right)^{\frac{1}{2}}{\left[\left(1+\frac{5}{2}\gamma\right)\frac{\sigma^{2}(R)}{\Delta_{\rm
th}}\right]}^{\frac{1}{2}+\frac{3}{2}\gamma}\exp\left(-\frac{\Delta_{\rm
th}^{2}}{2\sigma^{2}(R)}\right).$ (13)

If $\sigma^{2}(R)$ does not vary significantly, this gives an approximate
relation between $f(M)$ and $\sigma^{2}(R)$ (The scale $R$ for a fixed $M$ is
determined by Eq. (12)). At this level of approximation, it is clear that
there is one-to-one correspondence between $f(M)$ and $\sigma^{2}(R)$, and
$\sigma^{2}(R)$ can be uniquely reconstructed from $f(M)$.

On the other hand, the non-critical collapse case corresponds to setting the
constants $\gamma=0$ and $K\approx 0.3$. In such case, starting from the first
line of (9), the mass function can be expressed as

$\displaystyle f(M)\bigg{|}_{\rm non-critical}$
$\displaystyle=\frac{\sqrt{M_{\rm{eq}}}}{2f_{\rm
PBH}}\left(\frac{K}{M}\right)^{\frac{3}{2}}~{}{\rm
erfc}\left(\frac{\Delta_{\rm th}}{\sqrt{2}\sigma(R)}\right)$
$\displaystyle\approx\sqrt{\frac{M_{\rm eq}}{2\pi}}\frac{\sigma(R)}{f_{\rm
PBH}\Delta_{\rm
th}}{\left(\frac{K}{M}\right)}^{\frac{3}{2}}\exp\left(-\frac{\Delta_{\rm
th}^{2}}{2\sigma^{2}(R)}\right)\,,$ (14)

where ${\rm erfc}(x)$ is the complementary error function defined by

$\displaystyle\mathrm{erfc}\left(x\right)\equiv\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-z^{2}}{\rm{d}}z\,,$
(15)

and we used the asymptotic behavior of error function ${\rm
erfc}(x)\approx(\sqrt{\pi}x)^{-1}\exp(-x^{2})$ for $x\gg 1$ in the second
step. It should be noted that (13) with $\gamma=0$ differs from (III.2) by a
$\sqrt{2\pi}\approx 2.5$ pre-factor. This is because (13) is obtained by the
saddle point method which has assumed an approximate Gaussian peak. It may
differ from the exact integrand by a factor, which is within $\mathcal{O}(1)$
when we make the comparison between (13) and (III.2).

Therefore, we conclude that $\sigma^{2}(R)$ can be uniquely reconstructed from
$f(M)$ at least within the approximations made above.

Figure 1: Integrand of Eq. (9) as a function of $\frac{KM_{R}}{M}$ when
$\sigma^{2}(R)$ is independent of $R$. We use $\gamma=0.36,\Delta_{\rm
th}=0.4,\sigma^{2}=0.01$. Normalization of the vertical axis is chosen
arbitrarily and only relative height is meaningful.

###  III.3 Step (iii)

Our final step is to reconstruct the mass function $f(M)$ from the merger rate
density ${\cal R}(m_{1},m_{2},t)$ in the mass plane. It is known that PBHs can
form binaries either in radiation dominated era or at low redshift (e.g.
Sasaki:2018dmp ). In the former case, two neighboring PBHs which are
accidentally closer than the mean distance between PBHs become gravitationally
bounded and the tidal force from the surrounding PBHs gives the initial
angular momentum of the PBH binary. In the latter case, two PBHs randomly
moving inside the halo have a near miss at which they becomes gravitationally
bounded due to the energy loss by GW radiation. If the fraction of PBHs in
dark matter $f_{\rm PBH}$ is much smaller than unity, which we assume in this
paper, the former channel dominates the merger rate over the latter one. The
merger rate density in the former channel for the case of the extended PBH
mass function has been studied in Raidal:2017mfl ; Kocsis:2017yty ;
Raidal:2018bbj . In this paper, we adopt the most updated one Raidal:2018bbj
given by

$\displaystyle\mathcal{R}(m_{1},m_{2},t)=\frac{1.6\times 10^{6}}{{\rm
Gpc}^{3}{\rm yr}}f_{\rm
PBH}^{\frac{53}{37}}{\left(\frac{t}{t_{0}}\right)}^{-\frac{34}{37}}{\left(\frac{m_{t}}{M_{\odot}}\right)}^{-\frac{32}{37}}{\left(\frac{m_{t}^{2}}{m_{1}m_{2}}\right)}^{\frac{34}{37}}S[f|f_{\rm
PBH},m_{t}]m_{1}m_{2}f(m_{1})f(m_{2}).$ (16)

Here $m_{t}=m_{1}+m_{2}$ and $S$, which is the functional of $f(m)$, is the
so-called suppression factor which accounts for the distruption of the PBH
binaries in later Universe. Now, our problem is to obtain $f(m)$ out of
$\mathcal{R}(m_{1},m_{2},t)$ which is supposed to have been determined by
observations. If the dependence of the suppression factor on $m_{t}$ can be
ignored, one can easily estimate the ratio of $f(m_{1})$ and $f(m_{2})$ as

$\frac{f(m_{1})}{f(m_{2})}={\left(\frac{m_{2}}{m_{1}}\right)}^{\frac{21}{37}}\sqrt{\frac{\mathcal{R}(m_{1},m_{1},t)}{\mathcal{R}(m_{2},m_{2},t)}}.$
(17)

The normalization condition (8) gives the overall factor of $f(m)$. Finally,
$f_{\rm PBH}$ is obtained by substituting the obtained $f(m)$ into Eq. (16).
If the dependence of the suppression factor on $m_{t}$ is non-negligible,
obtaining $f(m)$ requires more procedures. From Eq. (16), we have

$\displaystyle\frac{f(m_{1})f(m_{2})}{f^{2}\left(\frac{m_{1}+m_{2}}{2}\right)}={\left(\frac{{\left(m_{1}+m_{2}\right)}^{2}}{4m_{1}m_{2}}\right)}^{\frac{3}{37}}\frac{\mathcal{R}(m_{1},m_{2},t)}{\mathcal{R}(\frac{m_{1}+m_{2}}{2},\frac{m_{1}+m_{2}}{2},t)}.$
(18)

In the following, we mainly concentrate on the case where the dependence of
the suppression factor on $m_{t}$ is non-negligible.

We can construct $f(m)$ satisfying (18) in the following way. We first choose
two masses $m_{1*}$ and $m_{2*}$. ($m_{1*}<m_{2*}$) arbitrarily and assign
values given by hand to $f(m_{1*})$ and $f(m_{2*})$. Then, the above equation
gives $f(\frac{m_{1*}+m_{2*}}{2})$. From the values of $f(m_{1*})$ and
$f(\frac{m_{1*}+m_{2*}}{2})$, we can obtain $f(\frac{3m_{1*}+m_{2*}}{4})$. By
repeating this procedure for $n$ times, we obtain the following recurrence
relation

$\displaystyle\frac{f(m_{1*})f(M_{n-1})}{f^{2}\left(M_{n}\right)}=\left(\frac{M_{n}^{2}}{m_{1*}M_{n-1}}\right)^{3/37}\frac{\mathcal{R}(m_{1*},M_{n-1},t)}{\mathcal{R}\left(M_{n},M_{n},t\right)}\,,$
(19) $\displaystyle M_{n}\equiv m_{1*}+\frac{\Delta m}{2^{n}}\,,\qquad\Delta
m\equiv m_{2*}-m_{1*}\,,\qquad n=0,1,2,...$ (20)

where $M_{0}=m_{2*}$. From this relation, $f(M_{n})$ can be further expressed
as

$\displaystyle f(M_{n})$
$\displaystyle=f(m_{1*})\left[\left(\frac{m_{2*}}{m_{1*}}\right)^{\frac{3}{37}}\frac{f(m_{2*})}{f(m_{1*})}\right]^{\frac{1}{2^{n}}}\left(\frac{m_{1*}}{M_{n}}\right)^{\frac{3}{37}}\prod_{i=1}^{n}\left[\frac{\mathcal{R}\left(M_{n-i+1},M_{n-i+1},t\right)}{\mathcal{R}(m_{1*},M_{n-i},t)}\right]^{\frac{1}{2^{i}}}\,.$
(21)

Figure 2: Illustration of the construction of the mass function $f(M_{n})$.
According to (20), we divide the mass between $m_{1*}$ and $m_{2*}$ into $n$
pieces such as $M_{1}$, $M_{2}$, …. The corresponding mass functions,
connected by the dashed blue lines, are then calculated by (21). For mass
value $m\in(M_{i+1},M_{i})$, we can set $f(m_{1})=f(M_{i+1})$ and
$f(m_{2})=f(M_{i})$ in (18), respectively, and repeat the above procedure to
obtain the mass function for a more narrow mass range.

As illustrated in Fig. 2, for mass value $m\in(M_{i+1},M_{i})$, we can set
$f(m_{1})=f(M_{i+1})$ and $f(m_{2})=f(M_{i})$ in (18), respectively, and
repeat the above procedure to obtain the mass function in a more narrow mass
range. Therefore, provided with enough data of merger rate for different mass
range, principally speaking, we can obtain the mass function in a sufficiently
narrow mass range. In this case, the normalization condition (8) becomes the
following form:

$\displaystyle\sum_{n=0}^{\infty}f(M_{n})\frac{\Delta m}{2^{n+1}}=1\,.$ (22)

The value of $f(m)$ for $m$ outside $(m_{1*},m_{2*})$ can be obtained in a
similar manner by plugging $m_{1}\to
2m_{1*}-m_{2*}<m_{1*},\frac{m_{1}+m_{2}}{2}\to m_{1*},m_{2}\to m_{2*}$ (or
$m_{1}\to m_{1*},\frac{m_{1}+m_{2}}{2}\to m_{2*},m_{2}\to
2m_{2*}-m_{1*}>m_{2*}$) in Eq. (18). The solution $f(m)$ obtained in this way
contains two free parameters (i.e. $f(m_{1*}),f(m_{2*})$) which corresponds to
the invariance of Eq. (18) under $f(m)\to Ae^{Bm}f(m)$ ($A,B$:constant). One
of the two parameters (we take it to be $A$) is fixed by the normalization
condition (22). Thus, the general solution of Eq. (18) can be expressed as
$f(m)=c(B)f_{0}(m)e^{Bm}$, where $f_{0}(m)$ is a particular solution and
$c(B)$ is the normalization constant. Substituting this solution to Eq. (16)
will enable us to fix $f_{\rm PBH}$ and $B$ simultaneously since they enter
differently in Eq. (16).

##  IV A simple example: neglecting the effect of critical collapse and
suppression factor

In this section, we consider the simplest case where the effects of the
critical collapse and suppression factor are neglected. In such case, the
primordial power spectrum can be analytically expressed in terms of the merger
rate and other parameters, which we will show below. We start with the
assumption (0) that the function of merger rate density $\mathcal{R}(m,m,t)$
has been obtained from observation data. We then use step (iii) to determine
the mass function. Let us choose a fixed mass point $m_{*}$ and design a value
for $f(m_{*})$. Since the effect of the suppression factor is neglected, using
(17), the mass function can be written as

$\displaystyle f(m)=g(t)m^{-\frac{21}{37}}\sqrt{\mathcal{R}(m,m,t)}\,,\qquad
g(t)\equiv
f(m_{*})m_{*}^{\frac{21}{37}}\left[\mathcal{R}(m_{*},m_{*},t)\right]^{-\frac{1}{2}}\,.$
(23)

It should be noted that $f(m)$ is independent from $t$ because of the
normalization condition (8). Secondly, from step (ii), since the effect of
critical collapse is neglected here, the variance can be expressed in terms of
the mass function by inverting (III.2),

$\displaystyle\sigma(m(R))=\frac{\Delta_{\rm th}}{\sqrt{2}}\left[{\rm
erfc}^{-1}\left(\frac{2f_{\rm PBH}\Omega_{\rm CDM}}{\left(M_{\rm
eq}K^{3}\right)^{\frac{1}{2}}}m^{\frac{3}{2}}f(m)\right)\right]^{-1}\,,$ (24)

where ${\rm erfc}^{-1}(x)$ is the inverse function of the complementary error
function defined by ${\rm erfc}\left[{\rm erfc}^{-1}(x)\right]=x$. We notice
that in this simple case, the mass of PBH is approximated as the fraction of
horizon mass such that $m=KM_{R}$. Hence, from (10), the mass of PBH is
related to the horizon scale $k=1/R$ as Green:2004wb

$\displaystyle m(k)=\left(\frac{k_{\rm eq}}{k}\right)^{2}M_{\rm
eq}K\left(\frac{g_{*,\rm eq}}{g_{*}}\right)^{\frac{1}{3}}\,,$ (25)

where $g_{*}$ is the number of relativistic degrees of freedom which is
expected to be of order $10^{2}$ in the early universe. At matter radiation
equality, we have $g_{*,\rm eq}\approx 3$ and $k_{\rm eq}=0.01{\rm Mpc}^{-1}$.
Inserting (23) and (25) into (24), we obtain

$\displaystyle\sigma(m(k))=\frac{\Delta_{\rm th}}{\sqrt{2}}\left[{\rm
erfc}^{-1}(A(m))\right]^{-1}\bigg{|}_{m=m(k)}\,,\qquad A(m)\equiv\frac{2f_{\rm
PBH}\Omega_{\rm CDM}}{\left(M_{\rm
eq}K^{3}\right)^{\frac{1}{2}}}\left(g(t)m^{\frac{69}{74}}\sqrt{\mathcal{R}(m,m,t)}\right)\,.$
(26)

Finally, we take step (i), i.e., insert (26) into (7), then the primordial
power spectrum can be obtained as

$\displaystyle\mathcal{P}_{\cal R}(k)$
$\displaystyle=\frac{81}{16}\left(4\sigma^{2}+k\frac{{\rm d}\sigma^{2}}{{\rm
d}k}\right)$
$\displaystyle=\frac{81}{8}\sigma^{2}\bigg{[}2-\sqrt{\pi}A(m)\left(\frac{69}{74}+\frac{m}{2\mathcal{R}(m,m,t)}\frac{{\partial}\mathcal{R}(m,m,t)}{{\partial}m}\right)\times$
$\displaystyle\qquad\times\left[{\rm
erfc}^{-1}\left(A(m)\right)\right]^{{\color[rgb]{0,0,1}-1}}\exp\left\\{\left[{\rm
erfc}^{-1}\left(A(m)\right)\right]^{2}\right\\}\bigg{]}_{m=m(k)}\,,$ (27)

where we have used the relation ${\rm d}\left[{\rm erfc}^{-1}(x)\right]/{\rm
d}x=-(\sqrt{\pi}/2){\rm exp}\left(\left[{\rm
erfc}^{-1}\left(x\right)\right]^{2}\right)$ in the last step. We notice that
in such a simple case, once the parameters such as $f_{\rm PBH}$, $K$ and
$\Delta_{\rm th}$ are known, the primordial powerspectrum is solely determined
by the merger rate of the PBH binaries with the same mass.

##  V Conclusion and discussion

In this paper, we have proposed the reconstruction method of primordial
curvature perturbations on small scales through the merger rate of binary
primordial black holes, which could form from large primordial curvature
perturbation on small scales during radiation dominated era. This is motivated
by the current observations of gravitational wave events by LIGO/VIRGO
collaboration, which indicate high mass and low spin of the BBHs. These
observations refresh the interest that PBHs may be one of the candidates of
the corresponding BBHs.

Adopting the scenario that PBHs originate from the primordial curvature
perturbations on small scales, we discussed how to reconstruct the power
spectrum of primordial curvature perturbations from observations of (possible)
PBH merger events, based on the following assumptions : (0) In the future,
some of the observed binary BH merger events are attributed to mergers of PBHs
formed, and the merger rate distribution of the PBHs in the mass plane is
observationally determined, (1) PBHs were formed out of rare high-$\sigma$
peaks of the primordial curvature perturbations in radiation dominated era,
(2) The window function takes the top-hat form in $k$-space, (3) The
primordial curvature perturbation follows Gaussian distribution and the effect
of non-linearity between the curvature and the density perturbations is
ignored. The reconstruction program is based on three steps. In the step (i),
one could relate the primordial power spectrum $\mathcal{P}_{\cal R}(k)$ and
the variance of the density perturbation $\sigma$ via (7). In the step (ii),
we found the relation (13) between the variance $\sigma$ and the PBH mass
function $f(M)$ where the effect of the critical phenomena is included, given
in Byrnes:2018clq ; Wang:2019kaf . In the step (iii), we presented the way to
obtain the PBH mass function $f(M)$ from the merger rate ${\cal
R}(m_{1},m_{2},t)$ of binary PBHs, adopting the theoretical expression of the
merger rate given in Raidal:2018bbj . By taking all steps, the power spectrum
for the curvature perturbations can be uniquely determined by the PBH merger
rate of binary black holes, provided with enough number of observed PBH merger
events with various mass.

It should be noted that the violation of the assumptions (1)–(3) might not
necessarily invalidate our methodology of the reconstruction. For example, in
the case of generic window function, the reconstruction of the primordial
power spectrum from the variance is still possible as presented in the
Appendix A. Although the reconstruction procedures in this case becomes more
complicated than in the case of the top-hat window function, the methodology
for step (i) still holds. In this sense, the principle of our methodology is
quite general though, once one relaxes each assumption, one needs to check
whether there is one-to-one correspondence in each step of the reconstruction
program.

On the other hand, there has been a long debate on which formalism should be
used to calculate the PBH abundance and the mass function. Besides the
formalism we have adopted in the paper, the peak theory formalism
Bardeen:1985tr  is another one of the mostly used formalism which accounts the
high-$\sigma$ peaks of the primordial perturbations for the PBH formation
Green:2004wb ; Suyama:2019npc ; Germani:2019zez ; Yoo:2020dkz . It remains to
be addressed whether the reconstruction can be achieved in this case.

Finally, it is interesting to notice that the merger rate of BBHs with a same
mass is needed in the reconstruction procedure, and one with different mass
can be used for the consistency check. Therefore, this reconstruction of the
power spectrum could be also a test of the PBH scenarios. In fact, we can
derive a consistency relation by using (16) as follows:

$\displaystyle\mathcal{R}(m_{1},m_{2},t)$
$\displaystyle=\left(\frac{m_{1}+m_{2}}{2\sqrt{m_{1}m_{2}}}\right)^{\frac{36}{37}}S[f|f_{\rm
PBH},m_{t}]\sqrt{\frac{\mathcal{R}(m_{1},m_{1},t)\mathcal{R}(m_{2},m_{2},t)}{S[f|f_{\rm
PBH},2m_{1}]S[f|f_{\rm PBH},2m_{2}]}}\,.$ (28)

When the mass dependence of the suppression factor can be neglected, (28)
reduces to the consistency relation derived in Ref. Kocsis:2017yty . If the
obtained mass function $f(m)$ determined only through $\mathcal{R}(m,m,t)$
cannot reproduce the merger rate for different masses adequately, that is, it
does not satisfy the above consistency relation (28), it suggests us to
include additional effects which are not included in Raidal:2018bbj . Finally,
even if we could include all of the relevant effects and the reconstruction
program would not work, it might rule out the PBH scenarios. Thus, this
reconstruction program is useful even to prove or to disprove the assumption
(0).

## Acknowledgements

T. S. is supported by the MEXT Grant-in-Aid for Scientific Research on
Innovative Areas No. 17H06359, and No. 19K03864. M. Y. is supported in part by
JSPS Grant-in-Aid for Scientific Research Numbers 18K18764, and JSPS Bilateral
Open Partnership Joint Research Projects. Y. Z. is supported by the
Fundamental Research Funds for the Central Universities; the National Natural
Science Foundation of China under Grants No. 12035011, No. 11975167, No.
11761161001, No. 11535004, and No. 11961141003; the National Key Research and
Development Program of China under Grant No. 2018YFA0404403. This work was
supported by Mitsubishi Foundation.

##  Appendix A Reconstruction of ${\cal P}_{\cal R}(k)$ from $\sigma^{2}(R)$

The purpose of this appendix is to study how to reconstruct ${\cal P}_{\cal
R}$ from $\sigma^{2}(R)$ which is given by observations only for discrete
values of $R$. Our starting point is Eq. (5) which gives $\sigma^{2}(q)$ in
terms of ${\cal P}_{\cal R}$ (Here, $q$ is wavenumber and we adopt the abuse
notation $\sigma^{2}(q)\equiv\sigma^{2}(R=1/q)$ in stead of $\sigma^{2}(R)$
for notational simplicity.). We suppose observations have determined the
values of $\sigma^{2}$ at discrete wavenumbers ${\bf
Q}=(q_{1},q_{2},\cdots,q_{N})$. Since there are only $N$ measured values, we
can determine the values of ${\cal P}_{\cal R}$ at maximally $N$ different
wavenumbers. At this point, there arises an ambiguity as to what combination
of the $N$ different wavenumbers we choose for determining ${\cal P}_{\cal R}$
since, in general, the reconstructed ${\cal P}_{\cal R}$ will depend on the
choice of the wavenumbers. Here, we consider the natural choice that the
wavenumbers for ${\cal P}_{\cal R}$ are also ${\bf Q}$ with the expectation
that the ambiguities mentioned above become negligible for suficiently large
$N$. Under this prescription, the discretized version of Eq. (5) can be
written as

$\sigma^{2}(q_{i})=\frac{16}{81}\sum_{j=1}^{N}\frac{\Delta
q_{j}}{q_{j}}W^{2}(q_{j}/q_{i}){(q_{j}/q_{i})}^{4}{\cal P}_{\cal R}(q_{j}),$
(29)

where $\Delta q_{j}$ is a suitable measure that approximates the integration
whose explicit expression depends on ${\bf Q}$. These are the linear algebraic
equations for $N$ variables ${\bf x}=({\cal P}_{\cal R}(q_{1}),\cdots,{\cal
P}_{\cal R}(q_{N}))$. Formally, they can be expressed by $N\times N$ matrix
$M$ and the vector ${\bf S}$ as

$M{\bf x}={\bf S},$ (30)

with the following identifications

$M_{ij}=\frac{16}{81}\frac{\Delta
q_{j}}{q_{j}}W^{2}(q_{j}/q_{i}){(q_{j}/q_{i})}^{4},~{}~{}~{}~{}~{}{\bf
S}=(\sigma^{2}(q_{1}),\cdots,\sigma^{2}(q_{N})).$ (31)

Thus, we can reconstruct ${\cal P}_{\cal R}$ within the discretization
approximation if ${\rm dim}~{}{\rm Ker}M=0$. In such as case, the
reconstructed ${\cal P}_{\cal R}$ is given by

${\bf x}=M^{-1}{\bf S}.$ (32)

Whether $M^{-1}$ exists or not will in general depend on the shape of the
window function as well as the discretization of the wavenumbers, and in
reality we have to check the existence of $M^{-1}$ in each case.

## References

  * (1) Y. B. Zel’dovich and I. D. Novikov, Soviet Astronomy 10, 602 (1967).
  * (2) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * (3) G. F. Chapline, Nature 253, no. 5489, 251 (1975).
  * (4) H. Niikura, M. Takada, N. Yasuda, R. H. Lupton, T. Sumi, S. More, T. Kurita, S. Sugiyama, A. More and M. Oguri, et al. Nature Astron. 3, no.6, 524-534 (2019) [arXiv:1701.02151 [astro-ph.CO]].
  * (5) A. Katz, J. Kopp, S. Sibiryakov and W. Xue, JCAP 1812, 005 (2018) [arXiv:1807.11495 [astro-ph.CO]].
  * (6) N. Smyth, S. Profumo, S. English, T. Jeltema, K. McKinnon and P. Guhathakurta, Phys. Rev. D 101, no.6, 063005 (2020) [arXiv:1910.01285 [astro-ph.CO]].
  * (7) P. Montero-Camacho, X. Fang, G. Vasquez, M. Silva and C. M. Hirata, JCAP 1908, 031 (2019) [arXiv:1906.05950 [astro-ph.CO]].
  * (8) S. Pi, Y. l. Zhang, Q. G. Huang and M. Sasaki, JCAP 05, 042 (2018) [arXiv:1712.09896 [astro-ph.CO]].
  * (9) Y. F. Cai, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. Lett. 121, no.8, 081306 (2018) [arXiv:1805.03639 [astro-ph.CO]].
  * (10) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016) [arXiv:1602.03837 [gr-qc]].
  * (11) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 24, 241103 (2016) [arXiv:1606.04855 [gr-qc]].
  * (12) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. X 6, no. 4, 041015 (2016) Erratum: [Phys. Rev. X 8, no. 3, 039903 (2018)] [arXiv:1606.04856 [gr-qc]].
  * (13) B. P. Abbott et al. [LIGO Scientific and VIRGO Collaborations], Phys. Rev. Lett. 118, no. 22, 221101 (2017) Erratum: [Phys. Rev. Lett. 121, no. 12, 129901 (2018)] [arXiv:1706.01812 [gr-qc]].
  * (14) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 119, no. 14, 141101 (2017) [arXiv:1709.09660 [gr-qc]].
  * (15) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Astrophys. J. 851, no. 2, L35 (2017) [arXiv:1711.05578 [astro-ph.HE]].
  * (16) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. X 9, no. 3, 031040 (2019) [arXiv:1811.12907 [astro-ph.HE]].
  * (17) T. Venumadhav, B. Zackay, J. Roulet, L. Dai and M. Zaldarriaga, Phys. Rev. D 101, no.8, 083030 (2020) [arXiv:1904.07214 [astro-ph.HE]].
  * (18) R. Abbott et al. [LIGO Scientific and Virgo], arXiv:2010.14527 [gr-qc].
  * (19) K. Belczynski, M. Dominik, T. Bulik, R. O’Shaughnessy, C. Fryer and D. E. Holz, Astrophys. J. 715, L138 (2010) [arXiv:1004.0386 [astro-ph.HE]].
  * (20) M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012) [arXiv:1202.4901 [astro-ph.HE]].
  * (21) M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel and R. O’Shaughnessy, Astrophys. J. 779, 72 (2013) [arXiv:1308.1546 [astro-ph.HE]].
  * (22) M. Dominik et al., Astrophys. J. 806, no. 2, 263 (2015) [arXiv:1405.7016 [astro-ph.HE]].
  * (23) K. Belczynski, D. E. Holz, T. Bulik and R. O’Shaughnessy, Nature 534, 512 (2016) [arXiv:1602.04531 [astro-ph.HE]].
  * (24) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Astrophys. J. 818, no. 2, L22 (2016) [arXiv:1602.03846 [astro-ph.HE]].
  * (25) M. Coleman Miller, Gen. Rel. Grav. 48, no. 7, 95 (2016) [arXiv:1606.06526 [astro-ph.HE]].
  * (26) S. Bird, I. Cholis, J. B. Munoz, Y. Ali-Haimoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli and A. G. Riess, Phys. Rev. Lett. 116, no. 20, 201301 (2016) [arXiv:1603.00464 [astro-ph.CO]].
  * (27) S. Clesse and J. Garcia-Bellido, Phys. Dark Univ. 15, 142 (2017) [arXiv:1603.05234 [astro-ph.CO]].
  * (28) M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Phys. Rev. Lett. 117, no. 6, 061101 (2016) Erratum: [Phys. Rev. Lett. 121, no. 5, 059901 (2018)] [arXiv:1603.08338 [astro-ph.CO]].
  * (29) M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no. 6, 063001 (2018) [arXiv:1801.05235 [astro-ph.CO]].
  * (30) J. Garcia-Bellido, J. F. N. Siles and E. Ruiz Morales, arXiv:2010.13811 [astro-ph.CO].
  * (31) D. N. Spergel et al. [WMAP Collaboration], Astrophys. J. Suppl. 148, 175 (2003) [astro-ph/0302209].
  * (32) Y. Akrami et al. [Planck], Astron. Astrophys. 641, A10 (2020) [arXiv:1807.06211 [astro-ph.CO]].
  * (33) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (34) J. Chluba, A. Kogut, S. P. Patil, M. H. Abitbol, N. Aghanim, Y. Ali-Haïmoud, M. A. Amin, J. Aumont, N. Bartolo and K. Basu, et al. Bull. Am. Astron. Soc. 51, no.3, 184 (2019) [arXiv:1903.04218 [astro-ph.CO]].
  * (35) D. Jeong, J. Pradler, J. Chluba and M. Kamionkowski, Phys. Rev. Lett. 113, 061301 (2014) [arXiv:1403.3697 [astro-ph.CO]].
  * (36) T. Nakama, T. Suyama and J. Yokoyama, Phys. Rev. Lett. 113, 061302 (2014) [arXiv:1403.5407 [astro-ph.CO]].
  * (37) K. Inomata, M. Kawasaki and Y. Tada, Phys. Rev. D 94, no.4, 043527 (2016) [arXiv:1605.04646 [astro-ph.CO]].
  * (38) R. Allahverdi, M. A. Amin, A. Berlin, N. Bernal, C. T. Byrnes, M. Sten Delos, A. L. Erickcek, M. Escudero, D. G. Figueroa and K. Freese, et al. Open J. Astrophys. 4, 2021 [arXiv:2006.16182 [astro-ph.CO]].
  * (39) M. Mapelli, Proc. Int. Sch. Phys. Fermi 200 (2020), 87-121 [arXiv:1809.09130 [astro-ph.HE]].
  * (40) B. Kocsis, T. Suyama, T. Tanaka and S. Yokoyama, Astrophys. J. 854, no. 1, 41 (2018) [arXiv:1709.09007 [astro-ph.CO]].
  * (41) A. Hall, A. D. Gow and C. T. Byrnes, Phys. Rev. D 102, 123524 (2020) [arXiv:2008.13704 [astro-ph.CO]].
  * (42) K. W. K. Wong, G. Franciolini, V. De Luca, V. Baibhav, E. Berti, P. Pani and A. Riotto, [arXiv:2011.01865 [gr-qc]].
  * (43) S. Clesse and J. Garcia-Bellido, arXiv:2007.06481 [astro-ph.CO].
  * (44) J. Garriga and N. Triantafyllou, JCAP 09 (2019), 043 [arXiv:1907.01455 [astro-ph.CO]].
  * (45) Z. C. Chen and Q. G. Huang, JCAP 08 (2020), 039 [arXiv:1904.02396 [astro-ph.CO]].
  * (46) L. Liu, Z. K. Guo and R. G. Cai, Phys. Rev. D 99 (2019) no.6, 063523 [arXiv:1812.05376 [astro-ph.CO]].
  * (47) T. Nakamura, M. Ando, T. Kinugawa, H. Nakano, K. Eda, S. Sato, M. Musha, T. Akutsu, T. Tanaka and N. Seto, et al. PTEP 2016 (2016) no.9, 093E01 [arXiv:1607.00897 [astro-ph.HE]].
  * (48) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 123 (2019) no.16, 161102 [arXiv:1904.08976 [astro-ph.CO]].
  * (49) S. W. Hawking, I. G. Moss and J. M. Stewart, Phys. Rev. D 26, 2681 (1982).
  * (50) M. Crawford and D. N. Schramm, Nature 298, 538 (1982).
  * (51) S. W. Hawking, Phys. Lett. B 231, 237 (1989).
  * (52) I. G. Moss, Phys. Rev. D 50, 676 (1994).
  * (53) S. Shandera, D. Jeong and H. S. G. Gebhardt, Phys. Rev. Lett. 120, no. 24, 241102 (2018) [arXiv:1802.08206 [astro-ph.CO]].
  * (54) P. Ivanov, P. Naselsky and I. Novikov, Phys. Rev. D 50, 7173 (1994).
  * (55) J. Garcia-Bellido, A. D. Linde and D. Wands, Phys. Rev. D 54, 6040 (1996) [astro-ph/9605094].
  * (56) J. S. Bullock and J. R. Primack, Phys. Rev. D 55, 7423 (1997) [astro-ph/9611106].
  * (57) P. Ivanov, Phys. Rev. D 57, 7145 (1998) [astro-ph/9708224].
  * (58) Y. F. Cai, C. Chen, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. D 100, no.4, 043518 (2019) [arXiv:1902.08187 [astro-ph.CO]].
  * (59) D. Y. Cheong, S. M. Lee and S. C. Park, JCAP 01, 032 (2021) [arXiv:1912.12032 [hep-ph]].
  * (60) C. Chen, X. H. Ma and Y. F. Cai, Phys. Rev. D 102, no.6, 063526 (2020) [arXiv:2003.03821 [astro-ph.CO]].
  * (61) G. A. Palma, S. Sypsas and C. Zenteno, Phys. Rev. Lett. 125, no.12, 121301 (2020) [arXiv:2004.06106 [astro-ph.CO]].
  * (62) Z. Zhou, J. Jiang, Y. F. Cai, M. Sasaki and S. Pi, Phys. Rev. D 102, no.10, 103527 (2020) [arXiv:2010.03537 [astro-ph.CO]].
  * (63) M. Y. Khlopov and A. G. Polnarev, Phys. Lett. B 97 (1980), 383-387
  * (64) T. Harada, C. M. Yoo, K. Kohri, K. i. Nakao and S. Jhingan, Astrophys. J. 833 (2016) no.1, 61 [arXiv:1609.01588 [astro-ph.CO]].
  * (65) K. Ando, K. Inomata and M. Kawasaki, Phys. Rev. D 97, no. 10, 103528 (2018) [arXiv:1802.06393 [astro-ph.CO]].
  * (66) S. Young, Int. J. Mod. Phys. D 29, no. 02, 2030002 (2019) [arXiv:1905.01230 [astro-ph.CO]].
  * (67) P. Pina Avelino, Phys. Rev. D 72, 124004 (2005) [arXiv:astro-ph/0510052 [astro-ph]].
  * (68) D. H. Lyth, JCAP 05, 022 (2012) [arXiv:1201.4312 [astro-ph.CO]].
  * (69) C. T. Byrnes, E. J. Copeland and A. M. Green, Phys. Rev. D 86, 043512 (2012) [arXiv:1206.4188 [astro-ph.CO]].
  * (70) S. Shandera, A. L. Erickcek, P. Scott and J. Y. Galarza, Phys. Rev. D 88, no. 10, 103506 (2013) [arXiv:1211.7361 [astro-ph.CO]].
  * (71) S. Young and C. T. Byrnes, JCAP 1308, 052 (2013) [arXiv:1307.4995 [astro-ph.CO]].
  * (72) S. Young and C. T. Byrnes, Phys. Rev. D 91, no.8, 083521 (2015) [arXiv:1411.4620 [astro-ph.CO]].
  * (73) S. Young, D. Regan and C. T. Byrnes, JCAP 02, 029 (2016) [arXiv:1512.07224 [astro-ph.CO]].
  * (74) G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, JCAP 03 (2018), 016 [arXiv:1801.09415 [astro-ph.CO]].
  * (75) R. g. Cai, S. Pi and M. Sasaki, Phys. Rev. Lett. 122, no.20, 201101 (2019) [arXiv:1810.11000 [astro-ph.CO]].
  * (76) V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Unal, JCAP 07 (2019), 048 [arXiv:1904.00970 [astro-ph.CO]].
  * (77) V. Atal, J. Cid, A. Escriva and J. Garriga, JCAP 05, 022 (2020) [arXiv:1908.11357 [astro-ph.CO]].
  * (78) C. M. Yoo, J. O. Gong and S. Yokoyama, JCAP 09, 033 (2019) [arXiv:1906.06790 [astro-ph.CO]].
  * (79) J. M. Ezquiaga, J. García-Bellido and V. Vennin, JCAP 03, 029 (2020) [arXiv:1912.05399 [astro-ph.CO]].
  * (80) S. Young, C. T. Byrnes and M. Sasaki, JCAP 07 (2014), 045 [arXiv:1405.7023 [gr-qc]].
  * (81) A. M. Green, A. R. Liddle, K. A. Malik and M. Sasaki, Phys. Rev. D 70, 041502 (2004) [arXiv:astro-ph/0403181 [astro-ph]].
  * (82) C. T. Byrnes, M. Hindmarsh, S. Young and M. R. S. Hawkins, JCAP 08, 041 (2018) [arXiv:1801.06138 [astro-ph.CO]].
  * (83) S. Wang, T. Terada and K. Kohri, Phys. Rev. D 99 (2019) no.10, 103531 [erratum: Phys. Rev. D 101 (2020) no.6, 069901] [arXiv:1903.05924 [astro-ph.CO]].
  * (84) J. C. Niemeyer and K. Jedamzik, Phys. Rev. Lett. 80 (1998), 5481-5484 [arXiv:astro-ph/9709072 [astro-ph]].
  * (85) J. Yokoyama, Phys. Rev. D 58 (1998), 107502 [arXiv:gr-qc/9804041 [gr-qc]].
  * (86) M. Raidal, C. Spethmann, V. Vaskonen and H. Veermäe, JCAP 02 (2019), 018 [arXiv:1812.01930 [astro-ph.CO]].
  * (87) M. Raidal, V. Vaskonen and H. Veermäe, JCAP 09 (2017), 037 [arXiv:1707.01480 [astro-ph.CO]].
  * (88) J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304 (1986), 15-61
  * (89) T. Suyama and S. Yokoyama, PTEP 2020, no.2, 023E03 (2020) [arXiv:1912.04687 [astro-ph.CO]].
  * (90) C. Germani and R. K. Sheth, Phys. Rev. D 101 (2020) no.6, 063520 [arXiv:1912.07072 [astro-ph.CO]].
  * (91) C. M. Yoo, T. Harada, S. Hirano and K. Kohri, arXiv:2008.02425 [astro-ph.CO].


* [47] X. L. Dong, A. Halevy, C. Yu, Data integration with uncertainty, The VLDB Journal 18 (2) (2009) 469–500.
  * [48] S. Jan, M. Li, G. Al-Sultany, H. Al-Raweshidy, Ontology alignment using rough sets, in: 2011 Eighth Intl. Conference on Fuzzy Systems and Knowledge Discovery (FSKD), Vol. 4, IEEE, 2011, pp. 2683–2686.
  * [49] E. P. Blasch, É. Dorion, P. Valin, E. Bossé, Ontology alignment using relative entropy for semantic uncertainty analysis, in: Proceedings of the IEEE 2010 National Aerospace & Electronics Conference, IEEE, 2010, pp. 140–148.
  * [50] Y. Zhang, A. V. Panangadan, V. K. Prasanna, UFOM: Unified fuzzy ontology matching, in: Proceedings of the 2014 IEEE 15th International Conference on Information Reuse and Integration (IEEE IRI 2014), IEEE, 2014, pp. 787–794.
  * [51] U. Bharambe, S. S. Durbha, R. L. King, Geospatial ontologies matching: An information theoretic approach, in: 2012 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), IEEE, 2012, pp. 2918–2921.
  * [52] D. Faria, C. Pesquita, E. Santos, M. Palmonari, I. F. Cruz, F. M. Couto, The AgreementMakerLight ontology matching system, in: Proc. of the OTM 2013 Conferences - Confederated International Conferences, Springer, 2013, pp. 527–541.  doi:10.1007/978-3-642-41030-7\\_38.
  * [53] G. Stoilos, D. Geleta, J. Shamdasani, M. Khodadadi, A novel approach and practical algorithms for ontology integration, in: ISWC, Springer, 2018, pp. 458–476.
  * [54] G. J. Klir, B. Yuan, Fuzzy sets and fuzzy logic: Theory and applications, Upper Saddle River (1995) 563.
  * [55] V. Cross, Uncertainty in the automation of ontology matching, in: Fourth Intl. Symposium on Uncertainty Modeling and Analysis, 2003. ISUMA 2003., IEEE, 2003, pp. 135–140.
  * [56] J. David, J. Euzenat, F. Scharffe, C. T. dos Santos, The alignment API 4.0, Semantic Web 2 (1) (2011) 3–10.  doi:10.3233/SW-2011-0028.
  * [57] D. Gale, L. S. Shapley, College admissions and the stability of marriage, The American Mathematical Monthly 69 (1) (1962) 9–15.  doi:10.1080/00029890.1962.11989827.
  * [58] J. Munkres, Algorithms for the assignment and transportation problems, Journal of the society for industrial and applied mathematics 5 (1) (1957) 32–38.  doi:10.1137/0105003.
  * [59] I. F. Cruz, F. P. Antonelli, C. Stroe, Efficient selection of mappings and automatic quality-driven combination of matching methods, in: Proceedings of the 4th International Workshop on Ontology Matching (OM) collocated with the 8th International Semantic Web Conference (ISWC), Vol. 551, Citeseer, CEUR-WS.org, 2009, pp. 49–60.
  * [60] A. Tordai, On combining alignment techniques, Ph.D. thesis, Vrije Universiteit, Amsterdam, The Netherlands (2012).
  * [61] A. Solimando, E. Jimenez-Ruiz, G. Guerrini, Minimizing conservativity violations in ontology alignments: Algorithms and evaluation, Knowledge and Information Systems 51 (3) (2017) 775–819.
  * [62] P. Arnold, E. Rahm, Semantic enrichment of ontology mappings: A linguistic-based approach, in: Proc. of the East European Conference on Advances in Databases and Information Systems, Springer, 2013, pp. 42–55.  doi:10.1007/978-3-642-40683-6\\_4.
  * [63] H. H. Do, S. Melnik, E. Rahm, Comparison of schema matching evaluations, in: Web, Web-Services, and Database Systems, NODe 2002 Web and Database-Related Workshops, Vol. 2593, Springer, 2002, pp. 221–237.  doi:10.1007/3-540-36560-5\\_17.
  * [64] M. Ehrig, J. Euzenat, Relaxed precision and recall for ontology matching, in: Proceedings of the K-CAP 2005 Workshop on Integrating Ontologies, Vol. 156, CEUR-WS.org, 2005.
URL http://ceur-ws.org/Vol-156/paper5.pdf

  * [65] J. Euzenat, Semantic precision and recall for ontology alignment evaluation, in: Proceedings of the 20th International Joint Conference on Artificial Intelligence (IJCAI), 2007, pp. 348–353.
URL http://ijcai.org/Proceedings/07/Papers/054.pdf

  * [66] S. Melnik, H. Garcia-Molina, E. Rahm, Similarity flooding: A versatile graph matching algorithm and its application to schema matching, in: Proceedings of the 18th International Conference on Data Engineering, IEEE Computer Society, 2002, pp. 117–128.  doi:10.1109/ICDE.2002.994702.
  * [67] G. C. D. Lab, Our World in Data, Last accessed 06 Oct 2020.
URL https://ourworldindata.org

  * [68] F. Lasso, Countries of the World: Country names linked to region, population, area, size, GDP, mortality, and more, Kaggle, Last accessed 21 Oct 2020 (2018).
URL https://www.kaggle.com/fernandol/countries-of-the-world

  * [69] SRK, Country Statistics – UNData: Dataset of economic, social, infra and environmental indicators, Kaggle, Last accessed 21 Oct 2020 (2018).
URL https://www.kaggle.com/sudalairajkumar/undata-country-profiles

  * [70] FAO, Statistics – Food security indicators, Food and Agriculture Organization of the United Nations, Last accessed 21 Oct 2020 (2014).
URL
http://web.archive.org/web/20141108183923/http://www.fao.org:80/economic/ess/ess-
fs/fs-data/en/#.VF5jZEX7TIU

  * [71] WorldBank, Prevalence of undernourishment (% of population), The World Bank Group, Last accessed 21 Oct 2020 (2020).
URL https://data.worldbank.org/indicator/SN.ITK.DEFC.ZS

  * [72] FAO, Statistics – Food security indicators, Food and Agriculture Organization of the United Nations, Last accessed 21 Oct 2020 (2020).
URL http://www.fao.org/economic/ess/ess-fs/ess-fadata/en/#.Xt04jUVKjIU

  * [73] R. J. Barro, Barro-Ursua Macroeconomic Data, Last accessed 21 Oct 2020 (2010).
URL https://scholar.harvard.edu/barro/publications/barro-ursua-macroeconomic-
data

  * [74] WorldBank, GDP (current US$), The World Bank Group, Last accessed 21 Oct 2020 (2020).
URL https://data.worldbank.org/indicator/NY.GDP.MKTP.CD

  * [75] WorldBank, World Development Indicators (WDI), The World Bank Group, Last accessed 21 Oct 2020 (2020).
URL https://datacatalog.worldbank.org/dataset/world-development-indicators

  * [76] WUP, World urbanization prospects: The 2018 revision, online edition, United Nations, Department of Economic and Social Affairs, Population Division, Last accessed 21 Oct 2020 (2018).
URL https://population.un.org/wup/Download/

  * [77] R. Zijdeman, F. Ribeira da Silva, Life expectancy at birth (Total), IISH Data Collection, V1, Last accessed 21 Oct 2020 (2015).
URL https://hdl.handle.net/10622/LKYT53

  * [78] WorldBank, Life expectancy at birth, total (years), The World Bank Group, Last accessed 21 Oct 2020 (2020).
URL https://data.worldbank.org/indicator/SP.DYN.LE00.IN

  * [79] ClimateWatch, Historical emissions, Data Explorer, Last accessed 21 Oct 2020 (2020).
URL https://www.climatewatchdata.org/data-explorer/historical-emissions

  * [80] FreedomHouse, 2020 list of electoral democracies, Freedom in the World, FIW, Last accessed 21 Oct 2020 (2020).
URL https://freedomhouse.org/report/freedom-world

  * [81] Unicef, Sexual violence in childhood, UNICEF DATA, Last accessed 21 Oct 2020 (June 2020).
URL https://data.unicef.org/topic/child-protection/violence/sexual-violence/

  * [82] J. Euzenat, An API for ontology alignment, in: Proceedings of the third International Semantic Web Conference, ISWC 2004, Vol. 3298, Springer, 2004, pp. 698–712.  doi:10.1007/978-3-540-30475-3\\_48.
  * [83] S. F. Pileggi, J. Hunter, An ontological approach to dynamic fine-grained urban indicators, Procedia Computer Science 108 (2017) 2059 – 2068, international Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland.  doi:https://doi.org/10.1016/j.procs.2017.05.003.
URL http://www.sciencedirect.com/science/article/pii/S187705091730491X

  * [84] S. F. Pileggi, Is the world becoming a better or a worse place? a data-driven analysis, Sustainability 12 (1) (2020).
  * [85] X. Liu, J. Zhang, D. Ma, Y. Bao, Z. Tong, X. Liu, Dynamic risk assessment of drought disaster for maize based on integrating multi-sources data in the region of the northwest of liaoning province, china, Natural hazards 65 (3) (2013) 1393–1409.