#  Primordial black holes, early galaxies, and antimatter in the Milky Way
Plenary ralk presented at $6^{th}$ International Conference on Particle
Physics and Astrophysics (ICCPA-2022)

A. D. Dolgova,b

aDepartment of Physics, Novosibirsk State University,

Pirogova st. 2, Novosibirsk, 630090 Russia

bBogolyubov Laboratory of Theoretical Physics, Joint Institute for Nuclear
Research,

Joliot-Curie st. 6, Dubna, Moscow region, 141980 Russia

Astronomical observations strongly incompatible with the canonical
cosmological model are reviewed. In particular too early formation of
galaxies, as discovered by HST and JWST, are discussed in detail. Other data
revealing highly dense population of the very young universe with plethora of
other different types of objects are presented. It is demonstrated that
similar or maybe even more pronounced problems can be seen also in the present
day universe. It is argued that all of the above mentioned problems can be
nicely fixed by assumption that the universe is filled with primordial black
holes in wide mass interval from a fraction of the solar mass up to
supermassive BH. The mechanism of PBH formation presented in 1993 is
described. The predicted by this mechanism log-normal mass spectrum of such
PBH is shown to agree very well with the data. Finally possible rich
population of our Galaxy by antimatter is discussed and new ways of its
identification are presented

##  1 Introduction

Discoveries of the last decades, especially by Hubble Space Telescope (HST)
and James Webb Space Telescope (JWST), created strong confusion among the
astronomical community. The early universe at redshifts $z\sim 10$ and the age
of a few hundred million years was found to be densely populated by galaxies,
quasars, gamma-bursters, supernovae, and contained a very high level of dust.
According to the common understanding the available universe age was by far
less than the necessary time for the creation of such rich universe
population. The frustrating feeling came to life that the standard
cosmological model, well supported by the theory and ”experiment”, was on the
verge of breaking.

Not less striking and equally puzzling evidence, that does not fit the
conventional frameworks, came from the observations of the contemporary
universe. The review of the problems, we encounter both in the early as well
as in the present day universe, at the stage of art which existed 4 years ago,
can be found in ref. [2]. But since that time much higher amount of
controversies was accumulated thanks to active and precise astronomical
observations.

All those surprising phenomena were in fact anticipated thirty years ago
according to our paper [3], where a new mechanism of formation of highly
massive primordial black hole (PBH) was proposed. This mechanism was further
elaborated in ref. [4]. A very simple log-normal mass spectrum of PBH was
predicted that, as later verified, very well agrees with observational data.
The abundant PBH population in very wide mass interval can eliminate the
tension between theory and observations, in particular, because supermassive
PBH could seed galaxy formation as it was envisaged in refs. [3, 4].

In addition to high mass PBH formation, the mechanism of refs [3, 4] could
lead to noticeable antimatter population of galaxies. In particular
antimatter, including antistars, may exist in the Milky Way. This exciting
predition seems to be confirmed by the recent studies.
Outline of this talk.
1\. Discovery by JWST of the dense population of the early universe with
galaxies, which according to the conventional cosmology cannot be there.
2\. More data contradicting traditional cosmology and astrophysics.
3\. Anticipated resolution of the problems by PBH pre-suggested in papers [3,
4].
4\. Prediction and test of Log-normal mass spectrum of PBHs. The only mass
spectrum confirmed by observatons.
5\. Black dark matter.
6\. Gravitational waves and PBH.
7\. Antimatter in the Galaxy (antistars, anit-nuclej, positrons).
8\. Basics of the mechanism of PBH and antimatter creation .

##  2 Strong blow to conventional cosmology by JWST

Observatkions of the several recent months made by JWST created almost panic
among traditional cosmologists and astrophysicists. It was discovered that the
pretty young universe with the age 200-300 million years contains plenty of
bright galaxies, see e.g. [5]-[13], which simply cannot be there according to
the accepted faith.

As is stated in the JWST publications: an unexpectedly large density (stellar
mass density $\varrho_{*}\gtrsim 10^{6}M_{\odot}$ Mpc-3) of massive galaxies
(stellar masses $M_{*}\geq 10^{10.5}M_{\odot}$) are discovered at extremely
high redshifts $z\gtrsim 10$. The following galaxies with record redshifts and
the corresponding age of the universe are observed as reported by: CEERS
(Cosmic Evolution Early Release Science): $z=14.3\pm 0.4$, $t_{U}=264$ Myr;
$z=16.7$, $t_{U}=235$ Myr Such unbelievably early observed galaxies forced the
JWST team to some retreat reported as: ”Bit of panic: Astronomers forced to
rethink early JWST findings. Revised instrument calibrations are bedevilling
work on the distant Universe”, according to
https://www.nature.com/articles/d41586-022-…

Several examples of the JWST data with very high redshifts are presented in
Fig. 1 and compared with theoretical expectations by the standard $\Lambda$CDM
model.

Figure 1: Comparison of some JWST data with theoretical expectation

In Fig. 2 the results of JWST are compared with those by HST for two events
for which both telescopes registered the same objects. An impressive agreement
is demonstrated.

Figure 2:

According to the canonical theory of large scale structure (LSS) formation the
density contrast $\Delta\equiv\delta\varrho/\varrho$ started to rise at the
onset of the matter dominated stage at $z=10^{4}$. After that $\Delta$ evolved
as the cosmological scale factor. Since initially $\Delta_{in}\lesssim
10^{-4}$, by the present time it may reach unity and after that fast LSS
formation takes place (violent relaxation - strong rising of the gravitational
field of the inhomogeneity) leading to the observed highly inhomogeneous
universe at the galactic and galaxy clusters scales.

In a simple way the process of structure foirmation can be understood as
follows. The velocity of the Hubble runaway at distance $r$ is $v_{H}=Hr$ and
the virial velocity in the gravitational field of the inhomogeneity is

$\displaystyle v^{2}_{grav}=\frac{4\pi r^{3}}{rm_{Pl}^{2}}\delta\varrho$ (1)

Using $H^{2}=(8\pi\varrho)/(3m_{pl}^{2})$ we find $v_{grav}\geq v_{H}$ if
$\Delta>1$. The probability of such huge density fluctuation for the flat
spectrum of perturbations is quite low.

There are two effects operating in the same directions. Firstly the available
time is constrained by the universe age, which is essentially equal to
$t_{U}=1/H$ and is quite short. In addition, a large value of $H$ means
expansion is very fast. That suppresses the efficiency of the structure
formation.

##  3 Problems preceding JWSP

Similar serious problems are known already for several years. The Hubble space
telescope discovered that the early universe, at $z=6-7$ was too densely
populated with quasars, alias SMBH, supernovae, gamma-bursters and happened to
be very dusty. No understanding is found in conventional cosmology how all
these creature were given birth to in such a short time. Moreover great lots
of phenomena in the $\sim 15$ billion year old, present day universe are in
strong tension with the conventional cosmological expectations.

HST sees the universe up to $z=6-7$, but accidentally a galaxy at $z\approx
12$ has been discovered for which both Hubble and Webb are in good agreement,
as we have already mentioned in the previous section. Still, despite the
earlier discoveries by HST, only after publications of JWST data the
astronomical establishment became seriously worried.

To summarise: observational data of the last decades present more and more
evidence indicating existence of the objects contradicting conventional
astrophysics and cosmology in the present day and in quite young universe.
Rephrasing Marcellus from ”The Tragedy of Hamlet, Prince of Denmark” we can
say ”Something is rotten in the state of Denmark the Universe”. However, all
the problems can be neatly solved if the universe is populated by primordial
black holes.

##  4 BH types by formation mechanisms

There three known types of BH depending upon the mechanism of their creation:

  1. 1.

Astrophysical black holes.
This BHs are created by the collapse of a star which exhausted its nuclear
fuel. The expected masses should start immediately above the neutron star
mass, i.e. about ${3M_{\odot}}$, but noticeably below $100M_{\odot}$. Instead
we observe that the BH mass spectrum in the galaxy has maximum at ${M\approx
8M_{\odot}}$ with the width $\sim(1-2)M_{\odot}$. The result is somewhat
unexpected but an explanations in the conventional astrophysical frameworks is
not excluded.
Recently LIGO/Virgo discovered BHs with masses close to $100M_{\odot}$. Their
astrophysical origin was considered to be complitely impossible. Now some,
though quite exotic, formation mechanisms have been suggested.

  2. 2.

Accretion created BHs.
Such BHs are formed by the accretion of matter on the mass excess in galactic
centres. It is known that in any large galaxy at the centre there exists a
supermassive black holes (SMBH) with masses varying from a few millions
$M_{\odot}$ (e,g, Milky Way) up to almost hundred billions $M_{\odot}$.
However, the conventional accretion mechanisms are not efficient enough to
create such monsters during the universe life-time, ${t_{U}\approx 14.6}$ Gyr.
At least 10-fold longer time is necessary, some references can be found in
[2], to say nothing about SMBH in 10 times younger universe.

  3. 3.

Primordial black holes (PBH).
PBH are supposed to be formed in the very early universe during pre-stellar
epoch, i.e. prior to star formation .The idea of primordial black holes and a
possible mechanism of their creation was pioneered by Zeldovich and Novikov
[14]. According to their idea, the density contrast in the early universe
inside the bubble radius, essentially equal to the cosmological horizon, might
accidentally happen to be large, ${\delta\varrho/\varrho\approx 1}$, then that
piece of volume would be inside its gravitational radius i.e. it became a PBH,
that decoupled from the cosmological expansion.
The mechanism was elaborated later by Hawking [15], and by Carr and Hawking
[16].

##  5 BH types by masses

Rather arbitrary black holes are separated into three groups depending on
their mass:
1\. Supermassive black holes (SMBH): $M=(10^{6}-10^{10})M_{\odot}$.
2\. Intermediate mass black holes (IMBH): $M=(10^{2}-10^{5})M_{\odot}$.
3\. Solar mass black holes: masses from a fraction of $M_{\odot}$ up to
$100M_{\odot}$.
4\. There can be also very light black holes, not yet observed, with masses in
the region $\sim 10^{20}$ g; they may be the ”particles” of the cosmological
dark matter.
The origin of most of these BHs is unclear, except maybe of the BHs with
masses of a few solar masses, which may be astrophysical.
Extremely unexpected was very high abundance of IMBH which are appearing
during last several years in huge numbers.
The assumption that (almost) all these black holes in the universe are
primordial except possibly the very light ones, strongly reduces or even
eliminates the tension between their observed abundances and possible
mechanisms of their formation.

##  6 Problems of the contemporary universe. Summary.

1\. SMBH in all large galaxies. The universe age is too short for their
formation through the commonly accepted accretion mechanism.

2\. Several SMBH are found in very small galaxies and even in (almost) empty
space, where not only the time duration but also an amount of material is
insufficient. An interesting recent observation was made by the Hobby-Eberly
Telescope at Texas’s McDonald Observatory suggesting the presence of a black
hole with a mass of about 17 billion $M_{\odot}$ equivalent to 14% of the
total stellar mass of the galaxy. Usually the mass of the central BH is about
0.1 % of the galaxy mass. This SMBH was observed by the analysis of the
motions of the stars near the center of the galaxy.

There appeared recently fresh evidence [17] indicating to supermassive BH with
the mass $3\times 10^{6}M_{\odot}$ in dwarf galaxy Leo 1. Much more new data
are presented practically today [18]. Six dwarf galaxies are identified that
have X-ray AGN. They are presumbly powered by SMBHs of $M>10^{7}M_{\odot}$.

It is not excluded, that such SMBHs, that are not hosted by a large galaxy,
might be pushed out of large galaxies in the process of galaxy collisions.
Such catastrophic event may even create plenty of wandering single
supermassive black holes. However, taking into account a large number of such
exotics, much more natural seems that all SMBH in small galaxies are
primordial. Simply they were unlucky not to acquire their own large galaxy,
since there was not enough matter around to build large galaxies.

3\. Problems with the BH mass spectrum in the Galaxy. According to the data
the masses are concentrated in the narrow interval ${(7.8\pm 1.2)M_{\odot}}$.

4\. Origin and properties of the sources of the observed gravitational waves,
encounter considerable difficulties, if one tries to explain them assuming
astrophysical formation of back hole binaries emitting the observed
gravitational radiation.

5\. IMBH, with ${M\sim(10^{3}-10^{5})M_{\odot}}$ are unexpectedly discovered
in dwarfs and globular clusters. Their origin is unclear, if they are not
primordial.

6\. Invisible Massive Astrophysical Compact Halo Objects (MACHOs), non-
luminous objects with masses ${\sim 0.5M_{\odot}}$ observed through
microlensing. It is unknown what are they and how they were created.

7\. Existence of very unusual stars in the Galaxy, among which there are too
fast moving stars and stars with unusual chemistry. Moreover, too old stars,
are found. Many of them look older than the Galaxy and maybe one is even older
that the universe (sic!?).

An assumption, that the black holes mentioned in the list above, are
primordial eliminates all the problems. The mechanism of PBH formation
suggested in papers [3, 4] predicts also existence of the unusual stars
mentioned in point 7.

##  7 Observations of black holes

The ancient point of view: BH are objects with so strong gravitational field
that nothing can escape it. According to Mitchell (1784): there may be bodies
for which the second cosmic velocity is larger than the speed of light. They
do not shine and do not reflect light, i.e. are absolutely dark, invisible.

However, the truth is quite the opposite, black holes are very well seen. BH
can emit all kind of radiation through the Hawking evaporation (though nobody
has yet seen it). The most powerful sources of radiation in the universe are
SMBH - quasars, point-like objects shining as a thousands of galaxies.
The methods of the BH observation include:
1\. Central mass estimate through analysis of stellar motion around the
supposed BH as e.g. discovery of BH in the center of the Millki Way.
2\. Gravitational lensing (MACHO and some other BHs).
3\. Electromagnetic radiation from the accreting matter; it is the mechanism
of quasar central engine, but mush smaller BH are also observed that way.
However, all these methods allow only to establish that there is a large mass
inside small volume. We need theory to conclude that there should be a black
hole inside. But the following method is free from this restriction:
4\. Registration of gravitational waves from coalescing double systems of
black holes. The data directly show that there are exactly coalescence of BHs.
This is the first test of General Relativity for strong fields and the first
observational proof of existence of the Schwarzschild solution.

##  8 Galaxy seeding by SMBH

We see thar a large amount of observational data are at odds with the
conventional model but nicely agrees with the model of creation of primordial
black holes and primordial stars suggested in refs. [3, 4]. According to the
standard approach the SMBH in galactic centres are formed after galaxies were
created by accretion mechanism. In papers. [3, 4] the validity of the opposite
scenario was conjectured, namely, SMBHs were formed first and subsequently
seeded galaxy formation. The hypothesis advocated in these works allows to
explain presence of SMBH in all large and several small galaxies accessible to
observation and resolves the problem of very early existence of galaxies
observed by HST and JWST.

This statement was recently rediscovered in ref. [19] according to which
”Recent observations with JWST have identified several bright galaxy
candidates at $z\gtrsim 10$, some of which appear unusually massive (up to
$\sim 10^{11}\ \rm M_{\odot}$)… Such early formation of massive galaxies is
difficult to reconcile with standard $\Lambda\rm CDM$ predictions, …we show
that the observed massive galaxy candidates can be explained… if structure
formation is accelerated by massive ($\gtrsim 10^{9}\ \rm M_{\odot}$)
primordial black holes that enhance primordial density fluctuations.”

Very recently, in December, 2022, there appeared another paper on the
possibility of SMBH impact on JWST-galaxy formation [20].

##  9 PBH and inflation

The mechanism suggested in ref. [3, 4] introduced some new features which were
later explored in a series of subsequent works. The proposed there scenarios
are heavily based on the Affleck-Dine [21] model of baryogenesis, that permits
to create very interesting features of the PBH population or some other
macroscopic compact objects, see below, sec. 15

In paper [3] inflationary mechanism was first implied for PBH formation. It
allowed to create PBH with huge masses, much larger than those in the
previously studied models. A year later inflationary creation of PBH was
explored in ref. [22] soon after in ref.[23], and two years later in [24].
Nowadays there exploded an avalanche of papers on inflationary formation of
PBH.

Presently inflationary mechanism of PBH production is commonly used. However,
except for predicting large masses of PBH, the models do not have much
predictive power because the mass spectra of the created PBHs are quite
complicated and strongly parameter dependent. No simple analytic expressions
have been presented. The only exception is the mechanism of refs [3, 4], which
predicts extremely simple log-normal mass spectrum of PBH:

$\displaystyle\frac{dN}{dM}=\mu^{2}\exp{[-\gamma\ln^{2}(M/M_{0})]}.$ (2)

The central value mass can be predicted theoretically [25]: ${M_{0}\sim
10M_{\odot}}$. It is equal to the horizon mass at QCD phase transition from
the free quark-gluon phase to the confinement phase. To be more precise the
horizon mass is approximately equal to $10M_{\odot}$ for the cosmic plasma
with vanishingly small chemical potential $\mu$. In our case $\mu$ is supposed
to be large of the order of the plasma temperature. In this case the phase
tranasion was probably delayed and the horizon mass could be 2-3 times larger.

An impressive feature of the the log-normal mass spectrum with the predicted
value of $M_{0}$ is that it is the only known spectrum tested by ”experiment”
in very good agreement with the observed densities of black holes in all mass
intervals from the solar mass BH, up to black holes with intermediates masses,
and further up to supermassive black holes. In partcular, the mechanism
developed in [3, 4] allows to explain the presence of SMBH in all large and
several small galaxies accessible to observation. Especially impressive is the
confirmation of the model by the chirp mass binaries measured by LIGO/Virgo
which is discussed in the following section.

##  10 Gravitational waves from BH binaries

There is general agreement between several groups that the gravitational waves
discovered by LIGO/Virgo interferometers originated from PBH binaries. We
discuss this issue here following our paper [26]. There are three problems
which indicate that the sources of GWs are most naturally primordial black
holes:
1\. Origin of heavy BHs (with masses ${\sim 30M_{\odot}}$). To form so heavy
BHs, the progenitors should have ${M>100M_{\odot}}$ and a low metal abundance
to avoid too much mass loss during the evolution. Such heavy stars might be
present in young star-forming galaxies but they are not observed in the
necessary amount. Recently there emerged much more striking problem because of
the observation of BH with ${M\sim 100M_{\odot}}$. Formation of such black
holes in the process of stellar collapse was considered to be strictly
forbidden. Some exotic mechanisms might be possibly allowed, such as e.g. BH
formation in the process of collapse of supermassive star heated by dark
matter annihilation inside [27].
On the other hand, primordial black holes with the observed by LIGO masses may
be easily created with sufficient density.
2\. Formation of BH binaries from the original stellar binaries. Stellar
binaries are formed from common interstellar gas clouds and are quite frequent
in galaxies. If BH is created through stellar collapse, small non-sphericity
results in a huge velocity of the BH and the binary is destroyed. BH formation
from PopIII stars and subsequent formation of BH binaries with tens of
${M_{\odot}}$ is estimated to be small. The problem of the binary formation is
simply solved if the observed sources of GWs are the binaries of primordial
black holes. They were at rest in the comoving volume, when inside horizon
they are gravitationally attracted and may loose energy due to dynamical
friction in the early universe. The probability for them to become
gravitationally bound is significant.
The conventional astrophysical scenario is not excluded but less natural.
3\. Low spins of the coalescing BHs . The low values of the BH spins sae
observed in GW150914 and in almost all (except for three) other events. It
strongly constrains astrophysical BH formation from close binary systems.
Astrophysical BHs are expected to have considerable angular momentum,
nevertheless the dynamical formation of double massive low-spin BHs in dense
stellar clusters is not excluded, though difficult. On the other hand, PBH
practically do not rotate because vorticity perturbations in the early
universe are vanishingly small. Still, individual PBH forming a binary
initially rotating on elliptic orbit could gain collinear spins about 0.1 -
0.3, rising with the PBH masses and eccentricity [28, 29]. This result is in
agreement with the GW170729 LIGO event produced by the binary with masses
${50M_{\odot}}$ and ${30M_{\odot}}$ and and GW190521.

To summarize: each of the mentioned problems may be solved in the conventional
frameworks but it looks much simpler to assume that the LIGO/Virgo sources are
primordial.

##  11 Chirp mass distribution

Two rotating gravitationally bound massive bodies are known to emit
gravitational waves. In quasi-stationary inspiral regime, the radius of the
orbit and the rotation frequency are approximately constant and the GW
frequency is twice the rotation frequency. The luminosity of the GW radiation
is:

$\displaystyle
L=\frac{32}{5}\,m_{Pl}^{2}\left(\frac{M_{c}\,\omega_{orb}}{m_{Pl}^{2}}\right)^{10/3}\,,$
(3)

where $M_{1}$, $M_{2}$ are the masses of two bodies in the binary system and
${M_{c}}$ is the so called chirp mass:

$\displaystyle M_{c}=\frac{(M_{1}\,M_{2})^{3/5}}{(M_{1}+M_{2})^{1/5}}\,,$ (4)

and

$\displaystyle\omega^{2}_{orb}=\frac{M_{1}+M_{2}}{m_{Pl}^{2}R^{3}}\,.$ (5)

The available data on the chirp mass distribution of the black holes in the
coalescing binaries in O1-O3 LIGO/Virgo runs are analyzed and compared with
theoretical expectations based on the hypothesis that these black holes are
primordial with log-normal mass spectrum [30]. The inferred best-fit mass
spectrum parameters are: $M_{0}=17M_{\odot}$ and ${\gamma=0.9}$, fall within
the theoretically expected range and show excellent agreement with
observations. On the opposite, binary black hole formation based on massive
binary star evolution require additional adjustments to reproduce the observed
chirp mass distribution. The results are presented in Figs. 3 and 4.

Figure 3: Model distribution $F_{PBH}(<M)$ with parameters $M_{0}$ and
${\gamma}$ for two best Kolmogorov-Smirnov tests. EDF= empirical distribution
function. Figure 4: Cumulative distributions $F(<M)$ for several astrophysical
models of binary BH coalescences.

So we can conclude that PBHs with log-normal mass spectrum perfectly fit the
data, while astrophysical BHs seem to be disfavoured.

New recent dats on GW observations were analysed by K.A. Postnov in his talk
at XXXIV International Workshop on High Energy Physics ”From Quarks to
Galaxies: Elucidating Dark Sides” are depicted in fig. 5.

In this talk is also presented an approximate fitting of the observed chirp-
mass distribution in the O1-O3 LVK GW compact binary coalescences (from GWTC-3
catalog) by two independent PBH populations with initial log-normal mass
distributions $M_{0}^{(1)}=5M_{\odot}$ and $M_{0}^{(2)}=30M_{\odot}$, see fig.
6.

Figure 5:  Figure 6: Approximation of the observed chirp-mass distribution in
the O1-O3 LVK GW compact binary coalescences (from GWTC-3 catalog) by two
independent PBH populations, from K.Postnov talk

In fig. 7 one can see the same approximation but by the simplest model of
astrophysical BH formation from the collapse of the CO-core of a massive star
and standard common-envelope parameter, with taking into account evolution of
star-formation rate in the universe with redshift (Postnov & Kuranov, 2019),
plus a population of PBHs with log-normal initial mass distribution with
$M_{0}^{(2)}=33M_{\odot}$ can be seen in Fig. 7.

Figure 7: The simplest model of astrophysical BH formation from the collapse
of the CO-core of a massive star, from K.Postniv talk

The picture looks more complicated than the earlier one described at the
beginning of this section. One possible interpretation is that there are two
populations of PBH with log-normal mass spectrum each but with different
values of $M_{0}$. The model can be modified to allow that and it was even
discussed earlier in the author’s papers, but aesthetically it looks not so
nice. Another option is that there a mixture of primordial and astrophysical
black holes observed by LIGO/Virgo. and even a possibility that some binaries
consist of a pair of primordial and astrophysical black holes. Hopefully
future data will help to resolve all the controverses.

##  12 Black Dark Matter

The first suggestion that PBH might be dark matter (DM) ”particles” was made
by S. Hawking in 1971 [15] and later by G. Chapline in 1975 [31], who noticed
that low mass PBHs might be abundant in the present-day universe with the
density comparable to the density of dark matter. The scale independent
spectrum of cosmological perturbations was assumed and thus the flat mass
spectrum of PBH in log mass interval was derived:

$\displaystyle{dN=N_{0}(dM/M)}$ (6)

with maximum mass ${M_{max}\lesssim 10^{22}}$ g, which hits the allowed for DM
mass range, see figure 8 below.

Chronologically next paper on DM consisting of PBHs [3] predicted extended
(log-normal) mass spectrum and permitted PBHs to have masses as high as
$10^{6}M_{\odot}$ and possibly even higher.

Upper limits on the cosmological mass density of PBHs are summarised in
reviews [32, 33] and presented here in fig. 8. The result are derived for
monochromatic mass spectrum of PBHs and generally model-dependent.

Figure 8: Constraints on $f(M)$ for a monochromatic mass function, from
evaporations (red), lensing (blue), gravitational waves (GW) (gray), dynamical
effects (green), accretion (light blue), CMB distortions (orange) and large-
scale structure (purple). Evaporation limits from the extragalactic gamma-ray
background (EGB), the Voyager positron flux (V) and annihilation-line
radiation from the Galactic centre (GC). Lensing limits from microlensing of
supernovae (SN) and of stars in M31 by Subaru (HSC), the Magellanic Clouds by
EROS and MACHO (EM) and the Galactic bulge by OGLE (O). Dynamical limits from
wide binaries (WB), star clusters in Eridanus II (E), halo dynamical friction
(DF), galaxy tidal distortions (G), heating of stars in the Galactic disk (DH)
and the CMB dipole (CMB). Large scale structure constraints(LSS). Accretion
limits from X-ray binaries (XB) and Planck measurements of CMB distortions
(PA). The incredulity limits (IL) correspond to one PBH per relevant
environment (galaxy, cluster, Universe). There are four mass windows (A, B, C,
D) in which PBHs could have an appreciable density.

The presented bounds are criticised in ref. [34] where it is argued that the
obtained in the quoted paper results i”…. reopens the possibility for dark
matter in the form of LIGO-mass PBHs.”

In ref. [35] it is stated that ”The most questionable step in this chain of
arguments is the use of overly simplified accretion models. We compare how the
same accretion models apply to X-ray observations from supermassive black
holes SMBHs, M87 and Sgr $A^{*}$. The comparison of these two SMBHs with
intermediate mass MACHOs suggests that the latter could, after all, provide a
significant constituent of all the dark matter.”

So the presented bounds on the cosmological density of PBH should be taken
with a portion of caution.

##  13 Dwarfs and globular clusters

A large number of intermediate mass black holes, that wee discovered during
last decades, hardly fits the narrow frameworks of the standard cosmological
model. However, if they are primordial with the determined above parameters of
the log-normal mass spectrum, their number is just what is necessary to
explain the data and in particular to understand the mechanism of the origin
of dwarf galaxies and globular clusters which is not well understood in the
conventional cosmology. Recently possible discovery of SMBH in dwarf galaxy
Leo 1 was announced [17]. Such SMBH surely could not be created by the
galactic matter accretion to the galaxy centre.

Much more new data are presented practically today [18]. Six dwarf galaxies
are identified that have X-ray AGN, powered by SMBHs of $M>10^{7}M_{\odot}$.
Most probably these dwarfs were seeded by primordial supermassive black holes
in accordance with ref [36].

As argued in ref. [36] IMBHs with masses of a few thousand solar mass, or
higher, can seed formation of globular clusters (GCs) and dwarf galaxies. In
the last several years such IMBH inside GSs are observed confirming this
suggestion. For example the BH with the mass $M\sim 10^{5}M_{\odot}$ is
discovered almost yesterday in the dwarf galaxy SDSS J1521+1404.

The astrophysical origin of IMBH encounter serious problems for all mass
values, but nevertheless huge number of them are discovered with all possible
masses. On the contrary, the described above model of PBH formation
excellently resolves all the inconsistencies.

##  14 Intermediate summary and antimatter in the Galaxy

The mechanism of PBH formation suggested in refs [3, 4] neatly cure all the
problem related to the observed population of the universe at high redshifts
as well as of the present day universe; The predicted log-normal spectrum of
PBH is tested and confirmed by the observations (the only one existing in the
literature). The predicted existence of IMBH in globular clusters is
confirmed.
So the model works great. Thus the seemingly crazy by-product of refs [3, 4],
namely prediction of antimatter in the Galaxy can come true as well. Probably
it is indeed the case. A surprisingly huge flux of cosmic positrons, of He-
antinuclei, and possibly even a population of antistars seem to be observed.

###  14.1 Anti-evidence: cosmic positrons

The observation of intense 0.511 line see refs [37, 38, 39], and earlier
references therein, presents a strong proof of abundant positron population in
the Galaxy. In the central region of the Galaxy electron–positron annihilation
proceeds at a surprisingly high rate, creating the flux:

$\displaystyle\Phi_{511\;{\rm keV}}={1.07\pm 0.03\cdot 10^{-3}}\;{\rm
photons\;cm^{-2}\,s^{-1}}.$ (7)

The width of the line is about 3 keV. It proves that the annihilation takes
place at rest. The emission mostly goes from the Galactic bulge and at much
lower level from the disk, The source of 0.511 MeV line in the Galactic bulge.
even got the name ”Great Annihihilator”.

Until recently the commonly accepted explanation was that ${e^{+}}$ are
created in the strong magnetic fields of pulsars but the recent results of AMS
probably exclude this mechanism, since the spectra of $\bar{p}$ and $e^{+}$ at
high energies are identical [40, 41]. It means that their origin is the same
and since antiprotons are not created in magnetic fields of pulsars, the
conclusion might follow that positrons are not produced by pulsars as well.
However, this might be true only for positrons with very high energies, about
or above tera-electronvolts.

###  14.2 Anti-evidence: cosmic antinuclei

The striking result reported by AMS team is the registration of anti-Helium-3
and anti-Helium-4 nuclei. Namely in 2018 AMS-02 announced possible observation
of six ${\overline{He}^{3}}$ and two ${\overline{He}^{4}}$ [42] [40]. In 2022
already 7 $\overline{D}$ ($E\lesssim 15$ GeV) and 9 $\overline{He}$, ($E\sim
50$ GeV) were observed. Surprisingly high fraction $\overline{He}/He\sim
10^{-9}$ was registered [41]. It is not excluded that the flux of anti-helium
is even much higher because low energy ${\overline{He}}$ may escape
registration in AMS.

Secondary production of different antinuclei in cosmic ray was estimated in
ref. [43]. According to this work anti-deuterium could be tmost efficiently
produced in the collisions ${\bar{p}\,p}$ or ${\bar{p}\,He}$ which can create
the flux ${\sim 10^{-7}/m^{2}/s^{-1}}$/steradian/GeV/neutron), i.e. 5 orders
of magnitude below the observed flux of antiprotons. Antihelium could be
created in the similar reactions and the fluxes of ${\overline{He}^{3}}$ and
${\overline{He}^{4}}$, that could be created in cosmic rays would respectively
be 4 and 8 orders of magnitude smaller than the flux of anti-D.

After the AMS announcements of observations of anti-$He^{4}$ there appeared
theoretical attempts to create anti-$He^{4}$ through dark matter annihilation.
This possibility does not look natural. Moreover, DM annihilation would
presumably create strong cosmic ray flux of other particles, which is not
observed.

In accordance with our model the observed antinuclej can be signatures of
primordial antimatter. However if they are synthesised in the standard big
bang anti-nucleosynthesis (anti-BBN) one would naturally expect the same
abundances of light elements as those created by the canonical BBN. According
to the latter the abundances of deuterium and helium-3 are much smaller then
that of helium-4, approximately by 4 orders of magnitude, while the relative
fraction of these antinuclej are approximately equal. There might be some
astrophysical explanation of that or this anomaly is related to the fact that
in our model antismatter is created in bubbles with unusually high baryon-to-
photon ratio $\beta$. In the canonical BBN $\beta\sim 10^{-9}$, while in our
case it may be as large as unity. However if $\beta\sim 1$ there is no
primordial D. On the other hand in our scenario the formation of primordial
elements takes place inside non-expanding compact stellar-like objects with
fixed temperature. If the temperature is sufficiently high, this so called BBN
may stop before abundant He formation and ends with almost equal abundances of
D and He. One can see that looking at abundances of light elements at a
function of temperature. If it is so, antistars may have equal amount of
$\overline{D}$ and $\overline{He}$.

###  14.3 Anti-evidence: antistars in the Galaxy

Last year a sensational announcement has been made [44] about possible
discovery of 14 antistars in our Galaxy, Milky Way. Quoting the authios: ”We
identify in the catalog 14 antistar candidates not associated with any objects
belonging to established gamma-ray source classes and with a spectrum
compatible with baryon-antibaryon annihilation”. The map of the observed anti-
sources is presented in fig. 9.

Figure 9:  Positions and energy flux in the 100 MeV - 100 GeV range of
antistar candidates selected in 4FGL-DR2. Galactic coordinates. The background
image shows the Fermi 5-year all-sky photon counts above 1 GeV

An additional possible method for antistar detection in the Galaxy or in its
halo has been proposed in ref. [45]. In astrophysically plausible cases of the
interaction of neutral atmospheres or winds from antistars with ionised
interstellar gas, the hadronic annihilation will be preceded by the formation
of excited ${p\bar{p}}$ and $He{\bar{p}}$ atoms. These atoms rapidly cascade
down to low levels prior to annihilation giving rise to a series of narrow
lines which can be associated with the hadronic annihilation gamma-ray
emission. The most significant are L (3p-2p) 1.73 keV line (yield more than
90%) from ${p\bar{p}}$ atoms, and M (4-3) 4.86 keV (yield $\sim 60$%) and L
(3-2) 11.13 keV (yield about 25%) lines from $He^{4}\bar{p}$ atoms. These
lines can be probed in dedicated observations by forthcoming sensitive X-ray
spectroscopic missions XRISM and Athena and in wide-field X-ray surveys like
SRG/eROSITA all-sky survey.

Bounds on the possible density if antistars in the Galaxy were studied in
several our papers [46, 47, 48]. It was shown that the restrictions are rather
mild and an abundant density of compact anti-stars in the universe even in the
Galaxy does not violate existing observations. The reason is that the
annihilation proceeds on a thin surface layer with a very short depth by the
order equal to the proton mean free path in the dense stellar medium. On the
other hand if the there were disperse antimatter clouds, the annihilation
would be by far more efficient and if so, the anticlouds did not survive to
our time, though they might have been existed in the early universe

A vey impressive would be star-antistar collision, which may even be a quasi-
periodic process of a star-antistar direct contact, explosion forcing them
apart and return to each other by gravitatinal attraction, etc…

##  15 PBH and anti-creation mechanism

The model of PBH and antimatter creation of refs. [3, 4] is essentially based
on the supersymmetry (SUSY) motivated baryogenesis, proposed by Affleck and
Dine (AD) [21], though the full extend SUSY is not necessary. SUSY predicts
existence of scalar field $\chi$ (or several such fields) with non-zero baryon
number, $B\neq 0$. Another important feature of the scenario is existence of
the flat directions in the self-potential of $\chi$, i.e. the directions along
which the potential does not rise. Simple examples of such quadratic and
quartic potentials with flat directions are the following, for quartic self-
interaction

$\displaystyle U_{\lambda}(\chi)=\lambda|\chi|^{4}\left(1-\cos 4\theta\right)$
(8)

and of the quadratic mass-like term:

$\displaystyle U_{m}(\chi)=m^{2}|\chi|^{2}[[1-\cos(2\theta+2\alpha)]$ (9)

where ${\chi=|\chi|\exp(i\theta)}$ and ${m=|m|e^{\alpha}}$. If ${\alpha\neq
0}$, C and CP are broken. In GUT SUSY baryonic number is naturally non-
conserved, due non-invariance of ${U(\chi)}$ with respect to phase rotation,
$\chi\rightarrow\chi\exp(i\theta)\chi$.

In the course of inflation $\chi$ quantum-fluctuates along flat directions
with increasing amplitude due to quantum instability of massless fields at De
Sitter stage [49, 50]. Thus taking into account that the wave length of the
fluciuations exponentially rises, $\chi$ could acquire a large classical
value.

When inflation is over and the symmetry maintaining flat directions brakes
down, $\chi$ starts to evolve to the equilibrium point, ${\chi=0}$, according
to the equation of the Newtonian mechanics with the liquid friction term:

$\displaystyle\ddot{\chi}+3H\dot{\chi}+U^{\prime}(\chi)=0.$ (10)

Due to quantum fluctuations orthogonal to the flat directions $\chi$ obtains
momentum in this orthogonal direction. This is how baryonic number of $\chi$
is generated:

$\displaystyle B_{\chi}=\dot{\theta}|\chi|^{2},$ (11)

$B$ is analogous to mechanical angular momentum in the two dimensional complex
plane $[{\Re}e\chi,{\Im}m\,\chi]$.

Decays of $\chi$ into quarks and antiquarks is supposed to conserve baryonic
number and transforms the baryonic number of $\chi$ into the baryonic number
of quarks. In this process a huge cosmological baryon asymmetry can be
generated, much larger than the observed one, $\beta\approx 10^{-9}$. If
$m\neq 0$, the angular momentum, B, could be generated due to possible
different direction of the quartic and quadratic valleys at low ${\chi}$. In
this case orthogonal quantum fluctuation are unnecessary. If CP-odd phase
${\alpha}$ is small but non-vanishing, both baryonic and antibaryonic domains
might be formed with possible dominance of one of them. Matter and antimatter
domains might be created with possible global ${B\neq 0}$.

In the model of ref. [3, 4] the AD, scenario of baryogenesis was essentially
modified due to addition of interaction of the Affleck-Dine field $\chi$ with
the inflaton $\Phi$. The interaction potential is taken in the form:

$\displaystyle
U=g|\chi|^{2}(\Phi-\Phi_{1})^{2}+\lambda|\chi|^{4}\,ln\left(\frac{|\chi|^{2}}{\sigma^{2}}\right)+(\lambda_{1}\chi^{4}+h.c.)+(m^{2}\chi^{2}+h.c.).$
(12)

The first term in this expression is the new type of interaction potential
between $\Phi$ and $\chi$. $\Phi_{1}$ is a constant, It is the amplitude of
the inflaton field, taken by $\Phi$ in the process of inflation. The remaining
duration of inflation after $\Phi$ passed $\Phi_{1}$ should secure the number
of e-foldings about 30-40 to allow for formation of sufficiently massive PBH.
Though the interaction between $\chi$ and $\Phi$ looks rather artificial, it
is not so. This is the general renormalizable coupling of two scalar fields.

The second term in eq. (12) is the Coleman-Weinberg potential [51] which
arises as a resukt of one-loop quantum corrections to $\lambda|\chi|^{4}$
interaction.

The remaining two terms are the toy model ones describing flat directions.

The constants $\lambda_{1}$ and $m$ are generally speaking complex. This may
lead to C and CP violation. However the charge symmetry would be broken only
if the relative phase of $\lambda_{1}$ and $m$ is nonzero. Coupling of $\chi$
to fermions (quarks) can break C and CP as well.

When $\Phi>\Phi_{1}$, potential (12) has deep minimum near $\chi=0$ and $\chi$
classically stays there. In the course of inflation $\Phi$ drops down and at
some moment reaches $\Phi_{1}$. So the barrier disappears and the window to
the flat direction opens. During this period, when $\Phi$ stays close to
$\Phi_{1}$, the field ${\chi}$ started to diffuse away from the old minimum,
according to quantum diffusion equation derived by Starobinsky which we have
generalised to a complex field ${\chi}$. At some stage for sufficiently large
$chi$ the diffusion turns into classical motion.

If the window to the flat direction, when ${\Phi\approx\Phi_{1}}$ is open only
during relatively short period, cosmologically small but possibly
astronomically large bubbles with high $\beta$ could be created, occupying a
small fraction of the universe volume. Indeed, when $\Phi$ passes the value
$\Phi_{1}$ sufficiently far, the old minimum at $\chi=0$ reappears and $\chi$
goes back to zero. While $\chi$ is large it propagates along the flat
direction of the quartic potential. When finally $\chi$ becomes small, it
starts to feel quadratic potential and in the process of motion from quartic
to quadratic potentials $\chi$ acquires a large angular momentum, that is a
large baryonic number.

If the probability of $\chi$ to reach a large value is not too big,
cosmologically small but possibly astrophysically large bubbles with high
baryon density would be formed while the rest of the universe has normal
${{\beta\approx 6\cdot 10^{-10}}}$, created by small ${\chi}$.

Initially large isocurvature perturbations are generated in chemical content
of massless quarks, while, density perturbations stay practically zero.
Density perturbations are generated rather late after the QCD phase
transition, when massless quarks turns into heavy nucleons with masses about 1
GeV, much larger than the temperature of the phase transition, $T_{qcd}\sim
100$ MeV, The emerging universe looks like a piece of Swiss cheese, where
holes are high baryonic density objects occupying a minor fraction of the
universe volume. These Hi-B Bubbles (HBB) mostly turn into primordial black
holes The evolution of the effective potential of $\chi$ is presented in Fig.
10.

Figure 10:  Schematic of the effective potential for $\chi$ and its evolution

This mechanism of massive PBH formation quite different from all others. The
fundament of PBH creation is build at inflation by making large isocurvature
fluctuations at relatively small scales, with practically vanishing density
perturbations, which appeared only much later. The mass spectrum of PBH
reflects the distribution of the bubble by size during inflation. Its log-
normal form is a general feature of diffusion processes.

##  16 Summary of the results and conclusion

$\bullet$ A large lot of outstanding problems of the canonical cosmology can
be nicely resolved if the universe is populated by primordial black holes with
masses in the interval from the solar mass up to supermassive BHs with masses
of the order of billion solar masses.
$\bullet$ The inverted mechanism of galaxy formation is proposed when firstly
a SMBH was formed and later it seeded the galaxy formation by gravitational
attraction of the surrounding matter. Thus an existence of SMBH in almost
empty space can be understood.
$\bullet$ The early galaxies observed by HST and JWST, when the universe was
only several hundred million years old, could be created if seeded by SMBHs.
$\bullet$ Existence of a noticeable number of galaxies with masses that are is
too small to allow for creation the observed SMBHs inside, can be explained if
these SMBH are primordial.
$\bullet$ Creation of SMBH in large contemporary galaxies by conventional
accretion mechanism demands time larger than the universe age. The problem
disappears if the central SPBH is primordial.
$\bullet$ Observations of several dwarf galaxies with SMBH in their centres,
confirmed our prediction made several years ago.
$\bullet$ The theoretically predicted log-normal mass spectrum of PBH is
verified by the chirp mass distribution of the gravitational waves observed by
LIGO/Virgo. The agreement between observation and theory is impressively good.
$\bullet$ PBHs formed according to our scenario explain the peculiar features
of the sources of GWs observed by LIGO/Virgo, e.g. an existence of BH with
$M=100M_{\odot}$
$\bullet$ The density of the intermediate mass black holes (IMBH),
$M=(10^{2}-10^{5})M_{\odot}$ well agrees with their primordial origin.
Assumption of the astrophysical formation of IMBH enoounters problems.
$\bullet$ Predicted by the model extremely old stars seem to exist even, one
”older than universe star” is found; the older age is mimicked by the unusual
initial chemistry. The model also predicts too fast moving stars, which are
also observed.
$\bullet$ Natural consequence of the suggested model of PBH creation leads to
noticeable population of our Galaxy by antimatter. This striking consequence
seems to be confirmed by recent observations.

## Acknowledgement

The work was supported by the Russian Federation Ministry of Science and
Higher Education (project FSUS-2022-0015).

## References

  * [1]
  * [2] A.D. Dolgov, Massive and supermassive black holes in the contemporary and early Universe and problems in cosmology and astrophysics, Usp. Fiz. Nauk 188 (2018) 2, 121-142, Phys. Usp. 61 (2018) 2, 115-132.
  * [3] A. Dolgov, J. Silk, Baryon isocurvature fluctuations at small scale and baryonic dark matter, PRD 47 (1993) 4244-4256.
  * [4] A.D. Dolgov, M. Kawasaki, N. Kevlishvili Inhomogeneous baryogenesis, cosmic antimatter, and dark matter, Nucl. Phys.B 807 (2009) 229-250, e-Print: 0806.2986 [hep-ph]
  * [5] M. Castellano et al, ”Early results from GLASS-JWST. III: Galaxy candidates at ${z\sim(9-15)}$”, arXiv:2207.09436;
  * [6] S.L. Finkelstein, et al ”Long Time Ago in a Galaxy Far, Far Away: A Candidate ${z\sim 14}$ Galaxy…” arXiv:2207.12474.
  * [7] N. J. Adams, et al, ”Discovery and properties of ultra-high redshift galaxies (${9<z<12}$)…” arXiv:2207.11217.
  * [8] H. Yan et al, ”First Batch of Candidate Galaxies at Redshifts 11 to 20… ” arXiv:2207.11558.
  * [9] N. Leethochawalit, et al, Early results from GLASS-JWST. X: Rest-frame UV-optical properties of galaxies at ${7<z<9}$”, arXiv:2207.11135;
  * [10] C. T. Donnan, et al, ”JWST Geers CEERS-93316, existed just 235 million years after the Big Bang,” arXive 2207.12356.
  * [11] Y. Harikane, et al ”Comprehensive Study on Galaxies at ${z=9-17}$ Found in the Early JWST Data: UV Luminosity Functions and Cosmic Star-Formation History at the Pre-Reionization Epoch” arXiv:2208.01612.
  * [12] A. Ferrara, A. Pallottini, P. Dayal, On the stunning abundance of super-early, massive galaxies revealed by JWST, arXiv:2208.00720
  * [13] M. Haslbauer, P. Kroupa, A.H. Zonoozi, H. Haghi, Has JWST already falsified dark-matter-driven galaxy formation?, arXiv:2210.14915 [astro-ph.GA].
  * [14] Ya. B. Zeldovich, I.D. Novikov, The Hypothesis of Cores Retarded During Expansion and the Hot Cosmological Model, Astronomicheskij Zhurnal, 43 (1966) 758: Soviet Astronomy, 10(4), 602-603 (1967).
  * [15] S. Hawking, ”Gravitationally collapsed objects of very low mass”.
  * [16] B. J. Carr and S. W. Hawking, ”Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * [17] F. Pacucci, A. Loeb, ”Accretion from Winds of Red Giant Branch Stars May Reveal the Supermassive Black Hole in Leo I”, ApJL 940 L33, 2022.
  * [18] M. Mezcua, M. Siudek, H. Suh, et al, ”Overmassive black holes in dwarf galaxies out to $z\sim 0.9$ in the VIPERS survey”, arXiv:2212.14057 [astro-ph.GA]
  * [19] B. Liu, V. Bromm, ”Accelerating early galaxy formation with primordial black holes”, arXiv:2208.13178,
  * [20] M. Volonteri, M Habouzit, M. Colpi, ”What if young z¿9 JWST galaxies hosted massive black holes?” e-Print: 2212.04710.
  * [21] I. Affleck, M. Dine, A New Mechanism for Baryogenesis, Nucl. Phys. B 249 (1985) 361-380.
  * [22] B.J. Carr, J.H. Hilbert, J.E. Lidsey, ”Black hole relics and inflation: Limits on blue perturbation spectra”, Phys.Rev.D 50 (1994) 4853, astro-ph/9405027
  * [23] P. Ivanov, P. Naselsky, I. Novikov, Inflation and primordial black holes as dark matter, PRD 50 (1994) 7173.
  * [24] J. Garcia-Bellido(, A.D. Linde, Density perturbations and black hole formation in hybrid inflation Phys.Rev.D 54 (1996) 6040-6058, e-Print: astro-ph/9605094 [astro-ph]
  * [25] A. Dolgov, K. Postnov, Why the mean mass of primordial black hole distribution is close to $10M_{\odot}$”. JCAP 07 (2020) 063.
  * [26] S. Blinnkov, A. Dolgov, N. Porayko, K. Postnov, ”Solving puzzles of GW150914 by primordial black holes,” CAP 1611 (2016), 036.
  * [27] J. Ziegler, K. Freese, ”Filling the Black Hole Mass Gap: Avoiding Pair Instability in Massive Stars through Addition of Non-Nuclear Energy”, Phys. Rev. D 104, 043015 (2021), arXiv:2010.00254:
  * [28] K. Postnov, N. Mitichkin, Spins of primordial binary black holes before coalescence, JCAP 1906 (2019) no.06, 044, arXiv 1904.00570v5 [astro-ph.HE]
  * [29] K. Postnov, A. Kuranov, N. Mitichkin, ”Spins of black holes in coalescing compact binaries”, Physics-Uspekhi vol. 62, No. 11, (2019), arXiv:1907.04218) .
  * [30] A.D. Dolgov, A.G. Kuranov, N.A. Mitichkin, S. Porey, K.A. Postnov, O.S. Sazhina, I.V. Simkine On mass distribution of coalescing black holes, JCAP 12 (2020) 017, e-Print: 2005.00892.
  * [31] G.F. Chapline, ”Cosmological effects of primordial black holes”. Nature, 253, 251 (1975).
  * [32] B. Carr, K. Kohri, Y. Sendouda, J. Yokoyama Constraints on Primordial Black Holes, arXiv:2002.12778 [astro-ph.CO]
  * [33] B. Carr, F. Kuhnel, Primordial black holes as dark matter candidates, SciPost Phys.Lect.Notes 48 (2022), e-Print: 2110.02821 [astro-ph.CO]
  * [34] C. Boehm, et al, Eliminating the LIGO bounds on primordial black hole dark matter, arXiv:2008.10743
  * [35] C. Corianò, P.H. Frampton, Does CMB Distortion Disfavour Intermediate Mass Dark Matter? arXiv:2012.13821 [astro-ph.GA]
  * [36] A. Dolgov, K. Postnov, ”Globular Cluster Seeding by Primordial Black Hole Population”, JCAP 04 (2017) 036, e-Print: 1702.07621 [astro-ph.CO].
  * [37] G. Weidenspointner et al., ”The sky distribution of positronium annihilation continuum emission measured with SPI/INTEGRAL”, Astron. Astrophys. 450, 1013 (2006); astro-ph/0601673.
  * [38] J. Knodlseder et al., ”The sky distribution of positronium annihilation continuum emission measured with SPI/INTEGRAL”, Astron. Astrophys. 441, 513 (2005); [arXiv:astro-ph/0506026].

  * [39] P. Jean et al., ”Spectral analysis of the Galactic e+e- annihilation emission”, Astron. Astrophys. 445, 579 (2006). [arXiv:astro-ph/0509298].
  * [40] S. Ting, L’Aquila Joint Astroparticle Colloquium, 10th November, 2021.
  * [41] S. Ting, COSPAR 2022, 16-24 July.Choutko
  * [42] A. Choutko, AMS-02 Collaboration, AMS Days at La Palma, La Palma, Canary Islands, Spain, (2018);
S. Ting, Latest Results from the AMS Experiment on the International Space
Station. Colloquium at CERN, May, 2018.

  * [43] R. Duperray, B. Baret, D. Maurin, et al, ”Flux of light antimatter nuclei near Earth, induced by cosmic rays in the Galaxy and in the atmosphere”, Phys.Rev.D 71 (2005) 083013 • e-Print: astro-ph/0503544 [astro-ph].
  * [44] S. Dupourqué, L. Tibaldo, P. von Ballmoos, ”Constraints on the antistar fraction in the Solar System neighbourhood from the 10-year Fermi Large Area Telescope gamma-ray source catalog”, Phys Rev D.103.083016 103 (2021) 083016.
  * [45] A.E. Bondar, S.I. Blinnikov, A.M. Bykov, A.D. Dolgov, K.A. Postnov, ”X-ray signature of antistars in the Galaxy”, JCAP 03 (2022) 03, 009, e-Print: 2109.12699 [astro-ph.HE].
  * [46] C. Bambi, A.D. Dolgov, Antimatter in the Milky Way, Nucl.Phys. B 784 (2007) 132-150 e-Print: astro-ph/0702350,
  * [47] A.D. Dolgov, S.I. Blinnikov, Stars and Black Holes from the very Early Universe, Phys.Rev.D 89 (2014) 2, 021301 e-Print: 1309.3395.
  * [48] S.I. Blinnikov, A.D. Dolgov, K.A. Postnov, Antimatter and antistars in the universe and in the Galaxy, Phys.Rev.D 92 (2015) 2, 023516 e-print: 1409.5736
  * [49] A.D. Linde, Scalar field fluctuations in the expanding universe and the new inflationary universe scenario, Phys. Lett. B116 (1982) 335-339
  * [50] A.A. Starobinski, Dynamics of phase transitions in the new inflationary universe scenario and generation of perturbations, Phys.Lett, B117 (1982) 175-178;
  * [51] S.R. Coleman, E.J. Weinberg, ’Radiative Corrections as the Origin of Spontaneous Symmetry Breaking”, Phys.Rev.D 7 (1973) 1888-1910.


# Dual-Frequency Quantum Phase Estimation Mitigates the Spectral Leakage of
Quantum Algorithms

Yifeng Xiong, , Soon Xin Ng, , Gui-Lu Long, ,
and Lajos Hanzo The authors are with School of Electronics and Computer
Science, University of Southampton, SO17 1BJ, Southampton (UK).L. Hanzo would
like to acknowledge the financial support of the Engineering and Physical
Sciences Research Council projects EP/P034284/1 and EP/P003990/1 (COALESCE) as
well as of the European Research Council’s Advanced Fellow Grant QuantCom
(Grant No. 789028). This work is also supported in part by China Scholarship
Council (CSC).

###### Abstract

Quantum phase estimation is an important component in diverse quantum
algorithms. However, it suffers from spectral leakage, when the reciprocal of
the record length is not an integer multiple of the unknown phase, which
incurs an accuracy degradation. For the existing single-sample estimation
scheme, window-based methods have been proposed for spectral leakage
mitigation. As a further advance, we propose a dual-frequency estimator, which
asymptotically approaches the Cramér-Rao bound, when multiple samples are
available. Numerical results show that the proposed estimator outperforms the
existing window-based methods, when the number of samples is sufficiently
high.

###### Index Terms:

Quantum phase estimation, spectral leakage, dual-frequency estimator,
algorithmic error mitigation.

##  I Introduction

Quantum phase estimation [1, 2, 3] is a key enabler of quantum computational
speedup over classical computers. It is a widely used component of quantum
algorithms providing substantial acceleration over their best classical
counterpart, including Shor’s factoring algorithm [4], the quantum clock
synchronization algorithm [5, 6], the Harrow-Hassidim-Lloyd (HHL) algorithm
[7, 8] conceived for solving linear systems, and the quantum counting
algorithm [9, 10]. Recently, iterative quantum phase estimators have also been
proposed for potential application in near-term noisy intermediate-scale
quantum (NISQ) computers [11, 12, 13].

The circuit diagram of quantum phase estimation is portrayed in Fig. 1, where
the initial state $\ket{\psi}$ of the data register is an eigenstate of the
unitary operator ${\mathcal{U}}$ satisfying
${\mathcal{U}}\ket{\psi}=e^{j\varphi}\ket{\psi}$, and $\mathrm{QFT}_{M}^{-1}$
denotes the inverse quantum Fourier transform [14, Sec. 5.1] applied to $M$
qubits, which is a quantum-domain version of the classical discrete Fourier
transform [15, Sec. 8]. This algorithm aims for estimating the phase
$\varphi$. Broadly speaking, the $2^{M}-1$ controlled unitaries produce a
$2^{M}$-point sinusoidal signal in the control register, whose frequency
corresponds exactly to the desired phase $\varphi$. This may be viewed as the
quantum-domain version of the power method of eigendecomposition [16]. When
efficient implementations of the power of the unitary operator ${\mathcal{U}}$
are available, exponential speedup over classical algorithms is possible, as
seen in Shor’s algorithm [4].

Quantum phase estimation yields the exact result, when the recording length
$2^{M}$ is an integer multiple of the period of the sinusoidal signal [14,
Sec. 5.2.1]. When this is not the case, the spectral leakage problem [17]
arises, which is an important topic in classical signal processing methods
related to the discrete Fourier transform. The quantum-domain version of this
problem is somewhat more grave, since the measurement outcome would assume
erroneous values at non-negligible probabilities, hence the information about
the correct phase is lost.

One of the conventional solutions to the spectral leakage issue is to multiply
the time-domain signal by a smooth “window function” [17]. This idea has been
applied to the quantum phase estimation problem in [18], where the authors
show that the cosine window [15, Sec. 7] is optimal in terms of the mean-
squared error (MSE) of single-sample-based estimation. An efficient quantum
circuit-based implementation of the cosine window has later been proposed in
[19], also showing that it has satisfactory performance in terms of the error
probability.

Figure 1: The circuit diagram of quantum phase estimation.

While existing treatises focus on single-sample-based estimation, in this
treatise, we consider the spectral leakage mitigation attained by multiple
samples (i.e. measurement outcomes). This is motivated by the fact that the
coherence time of NISQ computers is limited [20], hence classical computing
power may be harnessed, by appropriately fusing multiple samples. Explicitly,
our contributions are as follows:

  * •

We derive the Cramér-Rao bound (CRB) of the quantum phase estimation problem,
which sheds light upon the asymptotic multiple-sample performance of practical
estimators.

  * •

We propose a dual-frequency estimator that asymptotically attains the CRB upon
increasing the number of samples. We term this as the asymptotic regime.

  * •

Using numerical simulations, we demonstrate that the proposed dual-frequency
estimator outperform the cosine window-based solution, when the number of
samples is sufficiently large.

The rest of this treatise is organized as follows. We first formulate the
spectral leakage problem and discuss the existing countermeasures in Section
II. Next, we derive the CRB and propose the dual-frequency estimator in
Section III. The numerical results are then presented in Section IV, and we
conclude in Section V.

##  II Problem Formulation and Related Works

Let us commence by deriving of the output state of the quantum circuit used
for phase estimation and shown in Fig. 1. After the Hadamard gates, the joint
quantum state of both the control register and the target register becomes

$\ket{+}^{\otimes
M}\ket{\psi}=\frac{1}{\sqrt{2^{M}}}\sum_{n=0}^{2^{M}-1}\ket{n}\ket{\psi},$ (1)

which will be referred to as the initial state of the quantum phase estimation
algorithm in the rest of this chapter. The subsequent
controlled-${\mathcal{U}}$ gates transform the initial state into the
following state

$\displaystyle\ket{\phi}$
$\displaystyle=\frac{1}{\sqrt{2^{M}}}\sum_{n=0}^{2^{M}-1}\ket{n}{\mathcal{U}}^{n}\ket{\psi}$
(2)
$\displaystyle=\frac{1}{\sqrt{2^{M}}}\bigotimes_{k=0}^{M-1}\left(\ket{0}+e^{j2^{k}\varphi}\ket{1}\right)\ket{\psi}$
$\displaystyle=\frac{1}{\sqrt{2^{M}}}\sum_{n=0}^{2^{M}-1}e^{jn\varphi}\ket{n}\ket{\psi}.$

Note that the inverse quantum Fourier transform $\mathrm{QFT}_{M}^{-1}$ may be
expressed as follows111In classical signal processing theory, the transform
having the minus sign “$-$” in the exponents is typically referred to as the
discrete Fourier transform, and the transform having the plus sign is referred
to as the inverse discrete Fourier transform. However, the quantum computing
community is using a different convention, which refers the transform (3) as
the inverse quantum Fourier transform [14, Section 5.1].

$\mathrm{QFT}_{M}^{-1}=\frac{1}{\sqrt{N}}\sum_{m=0}^{N-1}e^{-j\frac{2\pi
m}{N}}\ket{m}\bra{n},$ (3)

where $N=2^{M}$ is the record length. Hence the quantum state of the control
register at the output of the $\mathrm{QFT}_{M}^{-1}$ may be expressed as

$\ket{\phi}_{\rm
out}=\frac{1}{N}\sum_{m=0}^{N-1}\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
m}{N}\right)}\ket{m}.$ (4)

The probability of observing the outcome $\ket{y}$ is thus given by

$\displaystyle f(y;\varphi)$ $\displaystyle=|\bra{y}\phi\rangle_{\rm
out}|^{2}$ (5)
$\displaystyle=\frac{1}{N^{2}}\left|\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
y}{N}\right)}\right|^{2}.$

When $\varphi=2\pi k/N,~{}k\in\mathbb{Z}$, we see that the probability
$f(y;\varphi)$ takes the maximum value of $1$ at $y=k$, hence the phase
estimation yields the exact solution. However, when $\varphi$ is not an
integer multiple of $2\pi/N$, $f(y;\varphi)$ is non-zero for almost all $0\leq
y\leq N-1$, causing large estimation errors. This phenomenon is referred to as
“spectral leakage” in the literature of classical signal processing [17, 21].

The spectral leakage issue can be mitigated using the classic windowing method
[22], multiplying the input signal by a “window” function defining a weighting
vector $\bm{\alpha}$. This results in the following observation probability
distribution

$f(y;\varphi)=\frac{1}{N}\left|\sum_{n=0}^{N-1}\alpha_{n}e^{jn\left(\varphi-\frac{2\pi
y}{N}\right)}\right|^{2},$ (6)

where $\alpha_{n}$ is the $n$-th entry of $\bm{\alpha}$, and
$\|\bm{\alpha}\|=1$. Upon comparing (5) to (6), we see that the original
quantum phase estimation corresponds to a rectangular window of
$\bm{\alpha}=\frac{1}{\sqrt{N}}\bm{1}$. This weighting procedure may be
implemented upon replacing the Hadamard gates in the dashed box of Fig. 1 by
customized state preparation circuits [19].

It has been shown that in terms of the lowest MSE, the optimal window is the
cosine window [18, 23], given by

$\alpha_{n}^{(\cos)}=\sqrt{\frac{2}{N}}\sin\left[\frac{\pi(n-1)}{N}\right].$
(7)

The authors of [19] have designed an efficient state preparation circuit for
the cosine window. The effect of the window may be intuitively interpreted by
considering the corresponding observation probability distribution, as
portrayed in Fig. 2. Observe that the “sidelobes” of the cosine window are
much lower than those of the rectangular window, hence the probability that
extremely large errors occur is substantially reduced. However, it is also
seen that the cosine window has a wider mainlobe. If the sidelobes can be
suppressed without widening the mainlobe, we may achieve better performance
than that of the cosine window. This motivates our dual-frequency estimator,
which will be detailed in Section III.

When multiple samples are available, it is possible to further improve the
accuracy by simply taking the average of these samples, which we refer to as
the “sample-mean estimator”. The sample-mean estimator has a time complexity
on the order of $O(N_{\rm s})$, where $N_{\rm s}$ is the number of samples. In
fact, when the sidelobes are effectively suppressed, the sample-mean estimator
is sufficient for obtaining a near-optimal performance (in the sense of CRB),
as will be detailed in Section III.

Figure 2: Illustration of the observation probability distribution
$f(y;\varphi)$ for rectangular and cosine windows, respectively.

##  III Proposed Estimator

In this section, we first derive the CRB of the quantum phase estimation
problem, which will be used as a performance benchmark in the numerical
simulations of Section IV. Moreover, it also motivates us to propose the dual-
frequency estimator detailed in Section III-B, which asymptotically approaches
the CRB upon increasing the number of samples. The generic CRB of the quantum
phase estimation problem has been derived in the literature [24, 25]. In this
section, we take into account the effect of the window function on the phase
estimation performance.

###  III-A CRB Analysis

The CRB is a lower bound on the mean-squared error of estimators [26], which
is tight, when the noise is weak or the number of observations is large. In
light of this, it may be used as the metric for determining the optimal window
function for a large number of samples.

In general, given the likelihood function $f(y;\theta)$ of an observation $y$,
the CRB of the parameter $\theta$ is given by

$\mathbb{E}\\{(\theta-\hat{\theta})^{2}\\}\geq\frac{1}{\mathrm{FI}(\theta)},$
(8)

where

$\mathrm{FI}(\theta)=\mathbb{E}\left\\{\left[\frac{\partial\ln
f(y;\theta)}{\partial\theta}\right]^{2}\right\\}.$ (9)

is the Fisher information of $\theta$ [26]. When there are multiple
independent and identical distributed observations, the total Fisher
information is the sum of the Fisher information of all individual
observations.

In the context of quantum phase estimation, the likelihood function of a
single observation is given by (6), which may be rewritten in a more compact
form as

$\displaystyle f(y;\varphi)$ $\displaystyle=\frac{1}{N}|\bm{e}_{y}^{\rm
H}\bm{\alpha}|^{2}$ (10) $\displaystyle=\frac{1}{N}\bm{e}_{y}^{\rm
H}\bm{\alpha}\bm{\alpha}^{\rm H}\bm{e}_{y},$

where

$\bm{e}_{y}:=\left[1,~{}e^{j\left(\varphi-\frac{2\pi
y}{N}\right)},\dotsc,~{}e^{j(N-1)\left(\varphi-\frac{2\pi
y}{N}\right)}\right]^{\rm T}$

denotes the vector containing all the phase terms in (5). For the simplicity
of further derivation, we rewrite the Fisher information in (9) in the
following alternative form

$\displaystyle\mathrm{FI}(\varphi)$
$\displaystyle=\mathbb{E}\left\\{\left[\frac{\partial\ln
f(y;\theta)}{\partial\theta}\right]^{2}\right\\}$ (11)
$\displaystyle=\sum_{y}f(y;\theta)\left[\frac{\partial\ln
f(y;\theta)}{\partial\theta}\right]^{2}$
$\displaystyle=\sum_{y}\frac{1}{f(y;\varphi)}\cdot\left(\frac{\partial}{\partial\varphi}f(y;\varphi)\right)^{2},$

where the partial derivative may be expressed as

$\displaystyle\frac{\partial}{\partial\varphi}f(y;\varphi)$
$\displaystyle=\frac{1}{N}\left[\left(\frac{\partial\bm{e}_{y}}{\partial\varphi}\right)^{\rm
T}\bm{\alpha}^{*}\bm{\alpha}^{\rm
T}\bm{e}_{y}^{*}+\left(\frac{\partial\bm{e}_{y}^{*}}{\partial\varphi}\right)^{\rm
T}\bm{\alpha}\bm{\alpha}^{\rm H}\bm{e}_{y}\right]$ (12)
$\displaystyle=\frac{j}{N}(\bm{e}^{T}\bm{N}\bm{\alpha}^{*}\bm{\alpha}^{T}\bm{e}^{*}-\bm{e}^{H}\bm{N}\bm{\alpha}\bm{\alpha}^{H}\bm{e})$
$\displaystyle=\frac{j}{N}\bm{e}_{y}^{\rm H}(\bm{\alpha}\bm{\alpha}^{\rm
H}\bm{N}-\bm{N}\bm{\alpha}\bm{\alpha}^{\rm H})\bm{e}_{y},$

where we have

$\bm{N}=\mathrm{diag}\\{0,1,\dotsc,N-1\\}.$

Upon introducing $\bm{A}=\bm{\alpha}\bm{\alpha}^{\rm H}$, the Fisher
information may be expressed as

$\displaystyle\mathrm{FI}(\varphi)$
$\displaystyle=\frac{1}{N}\sum_{y=0}^{N-1}\frac{|\bm{e}_{y}^{\rm
H}(\bm{\alpha}\bm{\alpha}^{\rm H}\bm{N}-\bm{N}\bm{\alpha}\bm{\alpha}^{\rm
H})\bm{e}_{y}|^{2}}{\bm{e}_{y}^{\rm H}\bm{\alpha}\bm{\alpha}^{\rm
H}\bm{e}_{y}}$ (13)
$\displaystyle=\frac{1}{N}\sum_{y=0}^{N-1}\frac{|\bm{e}^{H}[\bm{A},\bm{N}]\bm{e}|^{2}}{\bm{e}^{H}\bm{A}\bm{e}}$
$\displaystyle=\frac{4}{N}\sum_{y=0}^{N-1}\frac{\mathrm{Im}^{2}\\{\bm{e}_{y}^{\rm
H}\bm{A}\bm{N}\bm{e}_{y}\\}}{\bm{e}_{y}^{\rm H}\bm{A}\bm{e}_{y}}.$

When there are $N_{\rm s}$ samples, due to their mutual independence, the
total Fisher information is simply formulated as $N_{\rm
s}\cdot\mathrm{FI}(\varphi)$. This implies that the optimal window in the
asymptotic regime of infinite samples may actually be determined by
considering the Fisher information of the window applied to a single sample.

Figure 3: The square-root CRB vs. number of samples of various window
functions, compared with the RMSE of corresponding sample-mean estimators.

We may use (13) to obtain an intuition about the optimal window in the
asymptotic regime. In Fig. 3, we compare the square-root CRB of various
windows, where we set the record length to $N=128$, corresponding to $7$
control qubits. Observe that the rectangular window (equivalent to no
windowing at all over a finite interval) corresponds to the lowest CRB. This
suggests that it should yield the best asymptotic performance among all window
functions. However, we see that the naive sample-mean-based estimator operates
far from the CRB. By contrast, the sample-mean estimators of cosine and
Bartlett windows exhibit near-CRB performances. We will further address this
issue in Section III-B.

###  III-B The Dual-Frequency Estimator

Let us commence our discussion from constructing an approximate maximum
likelihood estimator based on the rectangular window with the aid of $N_{\rm
s}$ samples. We first denote the samples as
$\bm{y}=[y_{1},~{}y_{2},\dotsc,~{}y_{N_{\rm s}}]^{\rm T}$. From (5), the exact
maximum likelihood estimator is given by

$\hat{\varphi}_{\rm ML}=\mathop{\rm
argmax}_{\varphi\in[0,2\pi)}~{}\sum_{i=1}^{N_{s}}\ln\left|\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
y_{i}}{N}\right)}\right|^{2}.$ (14)

Since the samples only take integer values between $0$ and $N-1$, we may use
an alternative parametrization of $\bm{z}\in\mathbb{R}^{N}$, whose $n$-th
entry $z_{k}$ represents the number of times that $n-1$ occurs in the samples
$\bm{y}$. Thus we may rewrite (14) as

$\hat{\varphi}_{\rm ML}=\mathop{\rm
argmax}_{\varphi\in[0,2\pi)}~{}\sum_{k=1}^{N}z_{k}\ln\left|\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
k}{N}\right)}\right|^{2}.$ (15)

A naive strategy of solving (15) is an exhaustive search of $\varphi$ in
$[0,2\pi)$, which is computationally expansive. To simplify the problem, we
first obtain a rough estimate of $\varphi$ as

$\hat{\varphi}_{\rm
rough}=\frac{2\pi}{N}\left[(\mathop{\mathrm{argmax}}_{k}z_{k})-1\right].$ (16)

Note that the worst-case complexity of this step is on the order of $O(N_{\rm
s})$ when $N_{\rm s}\ll N$. Next we conduct a refined search within
$[\hat{\varphi}_{\rm rough}-2\pi/N,\hat{\varphi}_{\rm rough}+2\pi/N]$.
Although more sophisticated optimization techniques may have better
performance, here we consider the simple approach of a uniform grid search
over $O(\sqrt{N_{\rm s}})$ grid points, inspired by the fact that the Fisher
information is on the order of $O(N_{\rm s})$. To avoid the summation over $n$
in (15), we consider the following approximation for large $N$:

$\frac{1}{N}\left|\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
k}{N}\right)}\right|\approx\left|\mathrm{sinc}\left(\frac{N\varphi}{2\pi}-k\right)\right|.$
(17)

Using the above approximation, the worst-case complexity of the grid search is
reduced to $O(N_{\rm s}^{1.5})$. We denote the final estimate as
$\hat{\varphi}_{\rm AML}$, which may be obtained as

$\hat{\varphi}_{\rm
AML}=\mathop{\mathrm{argmax}}_{\begin{subarray}{c}\varphi\in{\mathscr{F}}\\\
\varphi=\hat{\varphi}_{\rm rough}+\frac{4k\pi}{NN_{\rm
g}},~{}k\in\mathbb{Z}\end{subarray}}\sum_{k=1}^{N}z_{k}\ln\left|\mathrm{sinc}\left(\frac{N\varphi}{2\pi}-k\right)\right|,$

where $N_{\rm g}$ is the number of grid points. The correction term is denoted
as $e_{\rm AML}=\hat{\varphi}_{\rm AML}-\hat{\varphi}_{\rm rough}$.

This estimator, however, does not in general yield satisfactory performance.
We may develop some further intuition concerning this issue by revisiting Fig.
2. Observe that most of the information about $\varphi$ is contained in the
pair of sample points within the mainlobe. When $\varphi$ is near $(2\pi
k+1)/N,~{}k\in\mathbb{Z}$, the largest and the second largest entries in
$\bm{z}$ may be used as reliable estimates of the two sample points. However,
when $\varphi$ is close to $2\pi k/N,~{}k\in\mathbb{Z}$, it is likely that the
second largest entry in $\bm{z}$ corresponds to the first sidelobe, causing an
unexpected estimation error. As portrayed in Fig. 4, the difficulty of
distinguishing the highest sidelobe from the sample point in the mainlobe
having smaller likelihood incurs an ambiguity problem for the maximum
likelihood estimator.

Figure 4: Illustration of the ambiguity problem associated with the
approximate maximum likelihood estimator when $\varphi$ is close to $2\pi
k/N,~{}k\in\mathbb{Z}$.

To understand this ambiguity problem in a more precise manner, let us consider
a concrete example. Specifically, we set the record length to $N=100$, the
number of samples to $N_{\rm s}=30$, and then plot the corresponding
estimation error based on the approximate maximum likelihood estimator
characterized in Fig. 5.222In practice, the record length can only be an
integer power of 2. Here we choose $N=100$ only for the purpose of a clearer
illustration. Observe that the errors are sometimes large when $\varphi$ is
close to an integer multiple of $1/2\pi N$. In particular, the magnitude of
the error grows linearly with the distance to the closest integer multiple of
$1/2\pi N$. As we have discussed previously, this phenomenon originates from
the fact that two possible interpolations of the sample points are hardly
distinguishable, hence there is a non-zero probability that the estimator
yields the erroneous result corresponding to the “mirror point” of the correct
one across the line of $\varphi=k/2\pi N$, where $k$ is the closest integer
multiple of $1/2\pi N$.

Figure 5: Scatter plot of the errors of the approximate maximum likelihood
estimator for $\varphi$ uniformly distributed over a fixed interval.

To tackle this problem, we split the samples into two sets, each having
$N_{\rm s}/2$ samples. For the first set, we compute the aforementioned
approximate maximum likelihood estimator, and obtain $\hat{\varphi}_{{\rm
rough},1}$ and $e_{{\rm AML},1}$. For the second set, we apply a frequency
offset of $1/2N$ (one half of a frequency resolution unit) to the control
register, resulting in the following maximum likelihood problem

$\displaystyle\hat{\varphi}_{{\rm ML},2}$ $\displaystyle=\mathop{\rm
argmax}_{\varphi\in[0,2\pi)}~{}\sum_{k=1}^{N}z_{k}\ln\left|\sum_{n=0}^{N-1}e^{jn\left(\varphi-\frac{2\pi
k}{N}+\frac{\pi}{N}\right)}\right|^{2}$ (18) $\displaystyle\hskip
8.53581pt+\frac{\pi}{N}.$

Similarly, we may obtain $\hat{\varphi}_{{\rm rough},2}$ and $e_{{\rm
AML},2}$. Based on these results, we construct four candidate intermediate
estimates, contained in a vector as follows:

$\displaystyle\bm{u}$ $\displaystyle\\!=\\!\Big{[}\hat{\varphi}_{{\rm
rough},1}\\!+\\!e_{{\rm AML},1},~{}\hat{\varphi}_{{\rm rough},1}\\!-\\!e_{{\rm
AML},1},$ (19) $\displaystyle\hskip 19.91692pt\\!\hat{\varphi}_{{\rm
rough},2}\\!+\\!e_{{\rm AML},2},~{}\hat{\varphi}_{{\rm rough},2}\\!-\\!e_{{\rm
AML},2}+\frac{2\pi}{N}\Big{]}^{\rm T}.$

Finally, we find the pair of entries in $\bm{u}$ that are the closest to each
other, meaning that

$(i,j)=\mathop{\mathrm{argmax}}_{(i,j)}~{}(u_{i}-u_{j})^{2},$ (20)

and the final estimate is given by

$\hat{\varphi}_{\rm DF}=\frac{1}{2}(u_{i}+u_{j}).$ (21)

The overall complexity of the dual-frequency estimator is at most $O(N_{\rm
s}^{1.5})$, which is comparable to the $O(N_{\rm s})$ complexity of the simple
sample-mean estimator. In practice, the complexity can be lower than $O(N_{\rm
s}^{1.5})$. To elaborate, note that the actual complexity of the dual-
frequency estimator is proportional to

$C_{\rm DF}=N_{\rm bin}\sqrt{N_{\rm s}},$ (22)

where $N_{\rm bin}$ denotes the number of non-empty “bins” $z_{k}\neq 0$. In
fact, most of the bins would be empty, since the probability of observing
those outcomes would be extremely low. Moreover, the bins containing a small
number of samples would provide less information about the desired phase
$\varphi$, since they are much noisier than the bins containing more samples.
For the numerical simulations presented in this chapter, we choose $8$ bins
having the largest number of samples. In light of this, the actual complexity
of the dual-frequency estimator becomes $O(N_{\rm s}+\sqrt{N})$, where the
$O(N_{\rm s})$ comes from the fact that re-organizing the samples into bins
would require $O(N_{\rm s})$ operations.

Figure 6: Scatter plot of the errors of the dual-frequency estimator for
$\varphi$ uniformly distributed over a fixed interval.

Next we show that the dual-frequency estimator indeed resolves the ambiguity
problem using a concrete example. Similar to Fig. 5, we set the record length
to $N=100$ and the number of samples to $N_{\rm s}=30$. The estimation error
of the dual-frequency estimator is portrayed in Fig. 6. As seen from the
figure, the dual-frequency estimator only produces small errors, and the
linearly increasing trend of errors seen in Fig. 5 does not appear. This
implies that the ambiguity point becomes distinguishable from the correct
solution.

The rationale of the dual-frequency estimator may be intuitively understood as
follows. To circumvent the difficulty of distinguishing the correct solution
from the ambiguity point, we include both points into the set of candidate
estimates, and also apply the same technique to the set of samples associated
with the $(1/2N)$-frequency offset. Since the ambiguity points of both sets
are unlikely to be the same due to the frequency offset, we may then identify
the correct solution by finding the matching pair of candidate solutions using
(20).

Finally, we show that the $(1/2N)$-frequency offset exploited in the estimator
may be implemented using single-qubit phase rotation gates, hence the
computational overhead of state preparation is negligible. Note that this
offset may be obtained by initializing the input state of the control register
as

$\ket{\phi}_{\rm in}=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}e^{\frac{j\pi
n}{N}}\ket{n}.$ (23)

This state admits the following simple tensor-product form

$\ket{\phi}_{\rm
in}=\bigotimes_{m=1}^{M}\frac{1}{\sqrt{2}}\left(\ket{0}+e^{\frac{j\pi
2^{m-1}}{N}}\ket{1}\right),$ (24)

which may be implemented by applying the phase rotation gate
${\mathcal{R}}_{\rm z}\left(\pi N^{-1}2^{m-1}\right)$ to each $m$-th qubit in
the control register. This method bears some resemblance with the
controlled-${\mathcal{U}}$ gates in the quantum phase estimation circuit,
which also constructs a sinusoidal signal (in other words, a frequency shift)
on the control register. This implies that the qubit complexity of the dual-
frequency estimator is $\lceil\log_{2}N\rceil$, which is the same as that of
the windowing methods.

##  IV Numerical Results

In this section, we characterize the performance of the proposed estimator
using numerical simulations. We first consider the relationship between the
root-mean-squared error (RMSE) of the estimators and the number of samples. In
this example, we set the number of control qubits to $M=7$, corresponding to
$N=128$. The number of samples $N_{\rm s}$ varies from $2$ to $100$. The
unknown phase $\varphi$ is randomly chosen from $(0,2\pi)$ in each of the
$10^{5}$ Monte Carlo trails. The RMSE of the dual-frequency estimator as well
as that of the sample-mean estimator based on the cosine window is portrayed
in Fig. 7, where the corresponding CRBs are also incorporated as benchmarks.

Figure 7: The phase estimation RMSE vs. the number of samples of both the
proposed dual-frequency estimator and the cosine window method, compared to
the corresponding CRBs.

Observe that all curves in Fig. 7 exhibit the same linear trend in the
asymptotic regime on the log-log scale. This implies that both estimators have
the same asymptotic error scaling as the CRB, which scales on the order of
$O(1/\sqrt{N_{\rm s}})$ (in terms of RMSE). When $N_{\rm s}$ is insufficiently
large, the dual-frequency estimator exhibits a performance similar to that of
the sample-mean estimator. But as $N_{\rm s}$ increases, the performance of
the dual-frequency estimator “switches” to near-CRB operation. Observe
furthermore that it outperforms the cosine window-based method for $N_{\rm
s}\geq 16$. The relatively low accuracy of the dual-frequency estimator in the
small-$N_{\rm s}$ regime originates from the fact that the number of samples
is to small for us to construct any reliable candidate estimate (as given in
(19)). Consequently, the sample points that are far away from the true value
of $\varphi$ (i.e., the “outliers”) cannot be reliably identified, hence would
cause large estimation errors.

Figure 8: The phase estimation RMSE vs. the record length $N$ of the proposed
dual-frequency estimator, the approximate maximum likelihood estimator, and
the cosine window method, compared to the corresponding CRBs.

Next we consider the dependence of the RMSE on the record length $N$. We set
the number of samples to $N_{\rm s}=30$. The number of control qubits varies
from $6$ to $10$, corresponding to $N=64,128,256,512,1024$. The corresponding
results are portrayed in Fig. 8. We also incorporate the approximate maximum
likelihood estimator in Fig. 8 for a better illustration.

Observe that the sample-mean estimator based on the rectangular window
exhibits an $O(1/\sqrt{N})$ scaling, while the others exhibit $O(1/N)$
scaling. This is a phenomenon that has also been observed in [18], which
suggests that the rectangular window does not provide a substantial quantum
speedup in the sense of RMSE scaling, since the $O(1/N)$ scaling (i.e., the
“Heisenberg limit [27]) is an important characteristic of quantum algorithms
conceived for phase estimation.

We also observe that the dual-frequency estimator does not exhibit the
aforementioned bimodal phenomenon with respect to $N$. This suggests that for
large $N$, we may still use a constant number of samples (for example, $N_{\rm
s}\geq 20$ as indicated by Fig. 7) to achieve a near-CRB performance.

Figure 9: The RMSE of the dual-frequency estimator vs. the record length, for
$N_{\rm s}$ values in different regions.

In Fig. 9, we portray the dependence of the RMSE on the record length $N$ for
different values of $N_{\rm s}$. Specifically, we pick three values of $N_{\rm
s}$ from the ambiguity region ($N_{\rm s}\leq 10)$, the transition region
($10<N_{\rm s}\leq 20$), and the asymptotic region ($N_{\rm s}>20$). We
observe that in the asymptotic region ($N_{\rm s}=30$), the dual-frequency
estimator exhibits a near-CRB performance, and has an order of $O(1/N)$ RMSE
scaling. By contrast, in the ambiguity region, the estimator exhibits an order
of $O(1/\sqrt{N})$ RMSE scaling. Finally, in the transition region, the order
of RMSE scaling is between $O(1/\sqrt{N})$ and $O(1/N)$.

##  V Conclusions

A dual-frequency phase estimator was proposed for mitigating the spectral
leakage-induced error in the quantum phase estimation algorithm based on
multiple samples. We have presented its CRB analysis in the asymptotic regime.
This CRB also inspires the design of our dual-frequency estimator. Compared to
the naive sample-mean estimator, the proposed estimator attains the Heisenberg
limit in the sense of exhibiting RMSE scaling on the order of $O(1/N)$.
Furthermore, the estimator is capable of outperforming the cosine window,
which is shown to be optimal for single-sample estimation, when the number of
samples is sufficiently large (but constant with respect to $N$). Our method
may inspire future research on the algorithmic error mitigation for a broader
range of quantum algorithms.

## References

  * [1] A. Y. Kitaev, “Quantum measurements and the Abelian stabilizer problem,” _arXiv preprint_ , 1995. [Online]. Available: https://arxiv.org/abs/quant-ph/9511026
  * [2] D. S. Abrams and S. Lloyd, “Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors,” _Phys. Rev. Lett._ , vol. 83, no. 24, p. 5162, Dec. 1999.
  * [3] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, “Quantum algorithms revisited,” _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , vol. 454, no. 1969, pp. 339–354, Jan. 1998.
  * [4] P. W. Shor, “Algorithms for quantum computation: Discrete logarithms and factoring,” in _Proc. 35th Annual Symp. Foundations of Computer Science_ , Santa Fe, New Mexico, USA, Nov. 1994, pp. 124–134.
  * [5] I. L. Chuang, “Quantum algorithm for distributed clock synchronization,” _Phys. Rev. Lett._ , vol. 85, pp. 2006–2009, Aug. 2000.
  * [6] J. Zhang, G. L. Long, Z. Deng, W. Liu, and Z. Lu, “Nuclear magnetic resonance implementation of a quantum clock synchronization algorithm,” _Phys. Rev. A_ , vol. 70, p. 062322, Dec. 2004.
  * [7] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” _Phys. Rev. Lett._ , vol. 103, no. 15, Oct. 2009.
  * [8] S.-J. Wei, Z.-R. Zhou, D. Ruan, and G.-L. Long, “Realization of the algorithm for system of linear equations in duality quantum computing,” in _Proc. 2017 IEEE 85th Veh. Technol. Conf. (VTC Spring)_ , Sydney, Australia, Jun. 2017, pp. 1–4.
  * [9] G. Brassard, P. HØyer, and A. Tapp, “Quantum counting,” in _Proc. 25th International Colloquium on Automata, Languages, and Programming_ , Aalborg, Denmark, Jul. 1998, pp. 820–831.
  * [10] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation,” _Contemporary Mathematics_ , vol. 305, pp. 53–74, 2002\.
  * [11] M. Dobšíček, G. Johansson, V. Shumeiko, and G. Wendin, “Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark,” _Phys. Rev. A_ , vol. 76, no. 3, p. 030306, Sep. 2007.
  * [12] N. Wiebe and C. Granade, “Efficient Bayesian phase estimation,” _Phys. Rev. Lett._ , vol. 117, no. 1, p. 010503, Jun. 2016.
  * [13] A. D. Corcoles, M. Takita, K. Inoue, S. Lekuch, Z. K. Minev, J. M. Chow, and J. M. Gambetta, “Exploiting dynamic quantum circuits in a quantum algorithm with superconducting qubits,” _arXiv preprint_ , 2021. [Online]. Available: https://arxiv.org/abs/2102.01682
  * [14] M. A. Nielsen and I. L. Chuang, _Quantum Computation and Quantum Information_ , 2nd ed. New York, NY, USA: Cambridge University Press, 2011.
  * [15] A. V. Oppenheim and R. W. Schafer, _Discrete-Time Signal Processing_ , 3rd ed. USA: Prentice Hall Press, 2009\.
  * [16] G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 3rd ed. Baltimore, MD, USA: Johns Hopkins University Press, 1996.
  * [17] F. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” _Proc. IEEE_ , vol. 66, no. 1, pp. 51–83, Jan. 1978\.
  * [18] W. van Dam, G. M. D’Ariano, A. Ekert, C. Macchiavello, and M. Mosca, “Optimal quantum circuits for general phase estimation,” _Phys. Rev. Lett._ , vol. 98, p. 090501, Mar. 2007.
  * [19] G. Rendon, T. Izubuchi, and Y. Kikuchi, “Effects of cosine tapering window on quantum phase estimation,” _arXiv preprint_ , 2021. [Online]. Available: https://arxiv.org/abs/2110.09590
  * [20] J. Preskill, “Quantum Computing in the NISQ era and beyond,” _Quantum_ , vol. 2, pp. 1–21, Aug. 2018.
  * [21] S. Ye, J. Sun, and E. Aboutanios, “On the estimation of the parameters of a real sinusoid in noise,” _IEEE Signal Process. Lett._ , vol. 24, no. 5, pp. 638–642, May 2017.
  * [22] S. Djukanović, “An accurate method for frequency estimation of a real sinusoid,” _IEEE Signal Process. Lett._ , vol. 23, no. 7, pp. 915–918, Jul. 2016.
  * [23] H. Shu, J. Wu, L. Senhadji, and L. Luo, “New fast algorithm for modulated complex lapped transform with sine windowing function,” _IEEE Signal Process. Lett._ , vol. 16, no. 2, pp. 93–96, Feb. 2009.
  * [24] U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, “Optimal quantum phase estimation,” _Phys. Rev. Lett._ , vol. 102, p. 040403, Jan. 2009.
  * [25] N. Wiebe and C. Granade, “Efficient Bayesian phase estimation,” _Phys. Rev. Lett._ , vol. 117, p. 010503, Jun. 2016.
  * [26] S. M. Kay, _Fundamentals of statistical signal processing: Estimation theory_ , 1st ed. Prentice-Hall, Inc., 1993\.
  * [27] W. Górecki, R. Demkowicz-Dobrzański, H. M. Wiseman, and D. W. Berry, “$\pi$-corrected Heisenberg limit,” _Phys. Rev. Lett._ , vol. 124, p. 030501, Jan. 2020.

  *[NISQ]: noisy intermediate-scale quantum
  *[MSE]: mean-squared error
  *[CRB]: Cramér-Rao bound
  *[RMSE]: root-mean-squared error