††institutetext:  1 Department of Physics and CTC, National Tsing Hua
University, Hsinchu 300, Taiwan

# 511 keV galactic line from first-order phase transitions and primordial
black holes

Po-Yan Tseng1, and Yu-Min Yeh1

###### Abstract

Hawking evaporation of primordial black hole (PBH), with mass $3\times
10^{-17}\lesssim M_{\rm PBH}/M_{\odot}\lesssim 7\times 10^{-17}$ and
fractional abundance $0.01\lesssim f_{\rm PBH}\lesssim 0.5$, well reproduces
511 keV gamma-ray excess from galaxy center. In this work, we investigated the
production mechanism of PBHs base on the first-order phase transition induced
by quartic effective thermal potential of a scalar field in dark sector. We
found the phase transition with vacuum energy, $\mathcal{O}(1)\lesssim
B^{1/4}/{\rm MeV}\lesssim\mathcal{O}(100)$, produces the desired PBH mass and
abundance fraction. Correlated signatures of gravitational wave and
extragalactic gamma-ray from phase transition and black hole evaporation,
respectively, are within $\mu$Ares and AMEGO/e-ASTROGAM projected
sensitivities.

##  1 Introduction

Since the lack of observational evidence of WIMP (Weakly interacting massive
particle), primordial black hole (PBH) gradually becomes one of the promising
dark matter (DM) candidate Hawking:1971ei ; Chapline:1975ojl ; Khlopov:2008qy
; Carr:2016drx ; Carr:2020gox ; Carr:2020xqk ; Green:2020jor  and provides
convincing explanation Clesse:2016vqa ; Bird:2016dcv ; Sasaki:2016jop  for
LIGO/Virgo LIGOScientific:2016aoc ; LIGOScientific:2016sjg ;
LIGOScientific:2017bnn  gravitational wave (GW) signals. Several production
mechanisms have been proposed. After inflation, PBHs can be generated through
the collapse of the overdensity regions, which are developed from the
primordial fluctuation re-entering the horizon Carr:1974nx ; Sasaki:2018dmp .
Other papers discussed PBHs directly form following the bubble collisions of
first-order phase transition (FOPT) Hawking:1982ga ; Moss:1994iq ;
Konoplich:1999qq ; Kodama:1982sf ; Gross:2021qgx ; Baker:2021nyl .

In this work, we consider an alternative PBHs production mechanism from FOPT
in the early Universe. During FOPT, an intermediate state, dubbed as Fermi
balls (FBs), were formed via the expanding bubbles. Then, the attractive
Yukawa interaction triggers a FB collapsing into a PBH Kawana:2021tde . The
success of FB formation requires the fermionic particles $\chi$’s being
trapped and aggregated inside a false vacuum surrounding by the expanding
bubbles of true vacuum. Hong:2020est ; Marfatia:2021twj ; Witten:1984rs ;
Bai:2018dxf . Since the Standard Model (SM) electroweak and QCD (quantum
chromodynamics) phase transitions are smooth crossovers, we therefore invoke
the FOPT being generated in the dark sector. Specifically, a commonly
discussed quartic effective thermal potential of a scalar field $\phi$ is
adopted. In addition, the scalar field couples to the $\chi$ through Yukawa
interaction, which serves two purposes. First, the non-zero VEV (vacuum
expectation value) of the scalar increases the mass $\chi$ in the true vacuum,
due to the energy-momentum conservation at bubble wall, which keeps the
$\chi$’s in the false vacuum if the mass difference is larger than the
critical temperature of the FOPT. The other purpose, the Yukawa interaction
induces attractive force between $\chi$’s and triggers the collapse of a FB
into a PBH. This instability of FB happened when the range of interaction
becomes comparable with the mean separation distance of $\chi$’s in a FB
Kawana:2021tde .

On the astrophysical observations perspective, an excess amount of photons
with energy 511 keV has been confirmed from the central region of Milky Way by
SPI/INTEGRAL Bouchet:2010dj . If someone identified this as a result of
electron-positron annihilation, it requires the injection rate of $2\times
10^{43}$ non-relativistic positrons per second (Weidenspointner:2004my, ;
Churazov:2004as, ; Weidenspointner:2007rs, ; Jean:2005af, ; Prantzos:2005pz,
). Various potential astrophysical sources have been proposed for the excess,
but encountered difficulties explaining the characteristic of the signal
Prantzos:2010wi . In addition, the MeV scale DM annihilation is hardly
consistent with cosmic microwave background observation Wilkinson:2016gsy .
Recently, it has been suggested that the positrons responsible for the 511 keV
line excess might be produced through the Hawking evaporation Hawking:1975vcx
; Gibbons:1977mu  of PBHs Keith:2021guq ; Frampton:2005fk ; Bambi:2008kx ;
Cai:2020fnq ; Laha:2019ssq ; DeRocco:2019fjq . Assuming the PBHs distribution
follows the NFW halo profile and concentrates in the inner Galaxy. The PBH
mass of $M_{\rm PBH}\sim\mathcal{O}(5\times 10^{16})~{}{\rm g}$ and the
$\mathcal{O}(10^{-3})$ fraction of total DM abundance, has been shown, not
only produces the right amount of positron flux but is also consistent with
existing limits of COMPTEL/INTEGRAL gamma-ray observations Keith:2021guq .

In this work, we aim to find the corresponding parameters space under the FOPT
scenario, which generates the PBH mass and relic abundance to explain the 511
keV excess. As a result, the correlated gamma-ray spectra form PBH evaporation
and GW signals from FOPT are predicted. This paper is organized as following:
In section 2, base on NFW and isothermal distributions, we scrutinize the
preferred PBH mass and abundance fraction for 511 keV excess. The realization
of PBH formation and correlated signals production are discussed in section 3
and 4, respectively. Finally, we summarize the results in section 5.

##  2 PBH evaporation

Hawking emission describes a PBH thermally produce primary particles with
masses lighter than the PBH temperature $T_{\rm PBH}=M^{2}_{\rm Pl}/M_{\rm
PBH}$, which numerically formulated as

$\displaystyle T_{\rm PBH}\simeq 5.3~{}{\rm
MeV}\times\left(\frac{10^{-18}M_{\odot}}{M_{\rm PBH}}\right)\,.$ (1)

The secondary particles productions come from the decays or fragmentation of
primary particles (e.g $e^{\pm},\mu^{\pm},\pi^{\pm},\pi^{0}$). A PBH
dominantly emits photons and neutrinos, but for $M_{\rm PBH}/M_{\odot}\lesssim
10^{-16}$, PBH can significantly produces $e^{\pm}$. The emission rate of
primary particle $i$ is given by Hawking:1974rv ; Hawking:1975vcx

$\displaystyle\frac{dN_{i}}{dEdt}=\frac{n^{\rm d.o.f}_{i}\Gamma_{i}(E,M_{\rm
PBH})}{2\pi(e^{E/T_{\rm PBH}}\pm 1)}\,,$ (2)

where $n^{\rm d.o.f}_{i}$ indicates the degrees of freedom of $i$ particle,
and the fermions (bosons) are distinguished by the $+(-)$ in the denominator.
The graybody factor $\Gamma_{i}(E,M)$ varies from particle to particle and is
derived from considering a wave packet scattering in the PBH spacetime
geometry from the PBH horizon to an observer at infinity.

Technically, we use the software package BlackHawk v2.1 Arbey:2019mbc ;
Arbey:2021mbl  to compute the $\gamma$ and $e^{\pm}$ production rates
including both primary and secondary components of PBH evaporation. The
positron emission rate per single PBH can be numerically obtained after
integrating out the evaporation spectrum

$\displaystyle L_{e^{+}}=\int dE\,\frac{dN_{e^{+}}}{dEdt}\,.$ (3)

###  2.1 Galaxy center 511 keV line

The galaxy center 511 keV line observed by SPI/INTEGRAL Bouchet:2010dj
exceeds the astrophysical contributions. In this paper, we focus on the
explanation of PBHs in mass range from $10^{15}\,$g to $2\times 10^{17}\,$g,
which significantly produce particles by Hawking evaporation. Positrons are
abundantly emitted from a PBH and more than 95% of them become non-
relativistic via ionization. Among the non-relativistic positrons, fraction of
$(1-f)$ directly annihilate with an electron to produce two 511 keV photons,
while the others fraction of $f\approx 0.967$ form a positronium bound state
with an electron Keith:2021guq . Among the later case, 25% of positroniums
annihilate to form pair of 511 keV photons, and the rest 75% yield three
photons with energy less than 511 keV. As a Result, the number of 511 keV
photons produced per positron is

$\displaystyle 2(1-f)+\frac{2f}{4}\approx 0.55\,.$ (4)

The flux of 511 keV photons from PBHs near the galactic center, covers solid
angle $\Delta\Omega$, is formulated as

$\displaystyle\Phi_{\mathrm{PBH}}(\Delta\Omega)=\frac{0.55L_{e^{+}}(M_{\textrm{PBH}})f_{\textrm{PBH}}}{4\pi
M_{\textrm{PBH}}}\int_{\Delta\Omega}{\int_{\textrm{l.o.s}}{\rho(\ell,\Omega)d\ell}d\Omega}\,,$
(5)

where $L_{e^{+}}(M_{\textrm{PBH}})$ is the positron production rate from a
PBH, $f_{\textrm{PBH}}$ be the fraction of PBHs abundance to the DM relic
abundance. The integrals are performed over solid angle and along the line-of-
sight. We assume the PBHs distribution follows the Navarro-Frenk-White (NFW)
DM halo profile:

$\displaystyle\rho(r)=\frac{\rho_{0}}{(r/R_{s})^{\gamma}[1+(r/R_{s})]^{3-\gamma}}\,,$
(6)

where $r$ is the distance from the galactic center, and the normalization
parameter $\rho_{0}$ adjusts the DM density near our solar system, $\rho=0.4$
$\textrm{GeV}/\textrm{cm}^{3}$ at $r=8.25$ kpc. Here, we fix $\gamma=1.6$ and
scale radius $R_{s}=20$ kpc in this analysis. Combining with the galactic disk
contribution, the angular distribution of 511 keV line predicted from PBHs
emission is shown in Fig. 1. It distributes as a function of galactic latitude
$b$ after averaging over the longitude profile $-8^{\circ}<\ell<+8^{\circ}$
and is compared with the INTEGRAL data of photon energy from 508.25 keV to
513.75 keV Bouchet:2010dj .

Figure 1: The orange curve is combined 511 keV gamma-ray flux from PBH
evaporation with $M_{\rm PBH}=10^{16}$ g, $f_{\textrm{PBH}}=10^{-4}$ and
astrophysical source Robin:2003uus (blue curve). Comparing to the INTEGRAL
data from 508.25 keV to 513.75 keV Bouchet:2010dj

The best-fit normalization for the 511 keV INTEGRAL data is obtained by the
$\chi^{2}$ test, which is defined by

$\displaystyle\chi^{2}(f_{\mathrm{PBH}})=\sum_{i}{\left[\frac{\Phi_{s}(b_{i},f_{\mathrm{PBH}})+\Phi_{\mathrm{disk}}(b_{i})-\Phi_{\textrm{data}}(b_{i})}{\sigma_{\Phi}(b_{i})}\right]^{2}}\,,$
(7)

and $i$ runs over the INTEGRAL data points in Fig.1. $\Phi_{\mathrm{disk}}$
stems from astrophysical positron emission in the disk. $\Phi_{s}$ represents
the PBHs contribution after taking into account the angular resolution of
INTEGRAL, which is obtained by Gaussian smearing as

$\displaystyle\Phi_{s}(b_{i},f_{\mathrm{PBH}})=\int^{b_{i}+3\sigma_{b_{i}}}_{b_{i}-3\sigma_{b_{i}}}{F(b,b_{i},\sigma_{b_{i}})\Phi_{\mathrm{PBH}}(b,f_{\mathrm{PBH}})db}\,,$
(8)

where $b_{i}$ and $\sigma_{i}$ respectively correspond to the central value
and the horizontal error of each INTEGRAL data point, and then the smearing
function is

$\displaystyle
F(b,b_{i},\sigma_{b_{i}})=\frac{1}{\sqrt{2\pi\sigma^{2}_{b_{i}}}}\exp{\left[-\frac{(b_{i}-b)^{2}}{\sigma^{2}_{b_{i}}}\right]}\,.$
(9)

The yellow region of Fig. 2 shows the best-fit parameter region of $(M_{\rm
PBH},f_{\rm PBH})$, which requires $\chi^{2}-\chi_{\mathrm{min}}^{2}\leq 4$.
For $M_{\rm PBH}\simeq 2\times 10^{17}~{}{\rm g}$, PBHs can serve as 100% DM
abundance and explain the 511 keV line excess simultaneously. For $M_{\rm
PBH}\lesssim 3\times 10^{16}~{}{\rm g}$, PBHs overproduce the extragalactic
gamma-ray, and thus is excluded by current upper bound of COMPTEL.

Figure 2: The yellow and cyan bands, which assume the NFW and Isothermal
profiles respectively, can produce the 511 keV excess observed by INTEGRAL via
the PBH evaporation. The constraints of COMPTEL+EGRET and sensitivities of
eASTROGAM and AMEGO are shown by the solid curves, requiring the extragalactic
photon from PBHs contribution does not exceed the observation data. The
benchmark points (BPs) label the PBHs produced from the FOPT scenario, and
their parameters are listed in Table 1.

###  2.2 Extragalactic gamma-ray

The PBH evaporation also produces extragalactic gamma-ray background, thus the
current gamma-ray observational data can be converted into the upper limit of
the PBHs abundance. To find the experiment limits in the parameter space, we
calculate the PBH extragalactic photon contributions and compare with the
current experimental data of COMPTEL/EGRET/FermiLAT and future sensitivities
of AMEGO/e-ASTROGAM as well. The extragalactic photon contributions due to PBH
is given by Marfatia:2021hcp

$\displaystyle\frac{d\Phi}{dE}=\int^{\mathrm{min}(t_{\mathrm{eva}},t_{0})}_{t_{\mathrm{CMB}}}{c[1+z(t)]\frac{f_{\mathrm{PBH}}\rho_{\mathrm{DM}}}{M_{\mathrm{PBH}}}\frac{d^{2}N_{\gamma}}{d\tilde{E}dt}\bigg{|}_{\tilde{E}=[1+z(t)]E}dt}\,,$
(10)

where $\rho_{\mathrm{DM}}=1.27\mkern 5.0mu\mathrm{GeV}\mathrm{m}^{-3}$ is the
average dark matter density in the extragalactic medium. From this expression,
we can see the extragalactic photon flux is independent from the uncertainties
of DM density profile near the galactic center, thereby it provides more
robust prediction than the galactic one. The lower limit of temporal integral
starts at the time of last scattering, $t_{\mathrm{CMB}}=3.8\times
10^{5}\mkern 5.0mu\mathrm{yr}$, after that photons decoupled from thermal
plasma and can propagate freely. Meanwhile, the upper limit of integral is set
by the shorter period of the lifetime of PBH ($t_{\mathrm{eva}}$) or the age
of Universe ($t_{0}=13.77\times 10^{9}\mkern 5.0mu\mathrm{yr}$). Applying
matter-dominated approximation to the evolution of Universe, the time-redshift
relation is given by $1+z(t)=\left(t_{0}/t\right)^{\frac{2}{3}}$. For
different masses of PBH, we show the upper limits of $f_{\mathrm{PBH}}$ in
Fig. 2, which are compatible with the present COMPTEL/EGRET/FermiLAT upper
bounds (blue curve) and future AMGEGO/e-ASTROGAM gamma-ray sensitivities
(red/black curves).

##  3 First-order phase transition

The PBH can be produced through a two-step processes during FOPT which is
induced by the effective thermal potential of a scalar $\phi$ in the hidden
sector. The FB stems from the aggregation of dark fermions $\chi$’s and plays
the role as intermediate state before PBH formation. Because the $\chi$’s
inside a FB couple through $\phi$, which provides an attractive Yukawa
potential. The length of Yukawa interaction increases as FB temperature
decreases, and thus Yukawa energy eventually dominates the total FB energy
which causing the FB collapses into a PBH Kawana:2021tde .

Specifically, the formation of PBHs can be realized by the Lagrangian
Marfatia:2021hcp

$\displaystyle\mathcal{L}\supset\bar{\chi}(i\not{\partial}-M_{i})\chi-
g_{\chi}\phi\bar{\chi}\chi-V_{\mathrm{eff}}(\phi,T)\,,$ (11)

where $M_{i}$ represents the bare mass of $\chi$ in the false vacuum,
$V_{\mathrm{eff}}$ is the finite-temperature quartic effective potential
commonly appears from theoretical models Dine:1992wr ; Adams:1993zs

$\displaystyle
V_{\mathrm{eff}}(\phi,T)=D(T^{2}-T^{2}_{0})\phi^{2}-(AT+C)\phi^{3}+\frac{\lambda}{4}\phi^{4}\,,$
(12)

which induces the cosmological FOPT. The FOPT happened when the temperature
becomes lower than critical temperature $T_{c}$, which is determined by the
condition $V_{\rm eff}(0,T_{c})=V_{\rm eff}(v_{\phi}(T_{c}),T_{c})$, at this
moment, the false vacuum ($\langle\phi\rangle=0$) was tunnelling to the true
vacuum ($\langle\phi\rangle=v_{\phi}$). We define the parameter $B$ to be the
zero-temperature potential energy density difference between the false and
true vacuum Marfatia:2021twj . Finally, the quartic potential can be described
by these input parameters Marfatia:2021twj

$\lambda,A,B,C,D\,.$

We use the analytical expression for the Euclidean action $S_{3}(T)/T$ of
quartic potentials Adams:1993zs  to compute the bubble nucleation rate per
unit volume $\Gamma(T)$ and then obtain the fraction of space in the false
vacuum $F(t)$ Marfatia:2021twj . The phase transition temperature $T_{\star}$
identified as percolation temperature can be obtained when a fraction $1/e$ of
the space remains in the false vacuum. Such that the corresponding time of
phase transition $t_{\star}$ is given by $F(t_{\star})=1/e\simeq 0.37$.

In general, GW signals accompany the FOPT from bubble collisions, sound waves,
and Magnetohydrodynamic turbulence. Huber:2008hg ; Nakai:2020oit ;
Espinosa:2010hh ; Caprini:2015zlo ; Kehayias:2009tn . In our case, since the
bubble wall reaches a relativistic terminal velocity thus corresponds to a
non-runaway bubble, the sound waves dominates the GW production
Caprini:2015zlo . We follow the semi-analytical approaches to evaluate the GW
spectra Marfatia:2020bcs ; Caprini:2015zlo . GW characteristic frequency
associates with the inverse time duration of the FOPT $\beta/H_{\star}\simeq
T_{\star}\,\left[d(S_{3}/T)/dT\right]$ Nakai:2020oit , which is normalized to
Hubble time scale. In addition to this parameter, the GW spectrum also depends
on $\alpha$, the strength of FOPT, which is defined as the ratio between false
vacuum energy to the total radiation energy.

To find $n_{\rm FB}$ the number density of FB, we need to know the critical
volume $V_{\star}$ which is defined from a volume in false vacuum satisfying
the condition $\Gamma(T_{\star})V_{\star}R_{\star}\simeq v_{w}$ Hong:2020est ,
and thus this volume would not further separate into smaller one during its
shrinking. Here $R_{\star}$ and $v_{w}$ are the radius of critical volume and
the bubble wall velocity, respectively. Under the assumption that each
critical volume corresponds to one FB, the volume fraction of false vacuum
implies the relation $n_{\rm FB}|_{T_{\star}}=F(t_{\star})/V_{\star}$ at
$t_{\star}$.

Inside a FB, to avoid complete annihilation $\bar{\chi}\chi\to\phi\phi$, there
must be a nonzero number density asymmetry $\eta_{\rm
DM}\equiv(n_{\chi}-n_{\bar{\chi}})/s$ (normalized to the entropy density).
Then we define $Q_{\rm FB}$ the total number of $\chi$’s comprising a FB:
$Q_{\rm FB}\equiv\eta_{\rm DM}(s/n_{\rm FB})|_{T_{\star}}$ Hong:2020est .

In general, the dark and SM sectors need not to be in thermal equilibrium,
otherwise the dark radiation ($\phi$ and $\chi$) would contributes the
effective number of extra neutrino species $\Delta N_{\rm eff}$ and exceed the
current observational upper bound. In fact, lowering the dark sector
temperature than SM sector alleviates the tension between $\Delta N_{\rm eff}$
and observations Marfatia:2021twj . Thus, we include the temperature ratio
$r_{T}\equiv T_{\star}/T_{\rm SM_{\star}}$ at $t_{\star}$ as one of the input
parameter for FOPT.

The dark fermion particles $\chi$’s aggregated to form macroscopic FBs due to
the FOPT. Subsequently, the attractive Yukawa force between $\chi$’s mediated
by $\phi$ destabilizes the FBs to form PBHs. In the following section, we
quantitatively discuss these processes.

###  3.1 FB properties and PBH formations

The derivation for the FB mass and radius can be found in Ref.
Marfatia:2021hcp . Including the Fermi gas kinetic energy, Yukawa potential
energy, and the temperature-dependent potential energy difference between the
false and true vacua; the total energy of a FB is approximately written as

$\displaystyle E_{\rm FB}$ $\displaystyle=$
$\displaystyle\frac{3\pi}{4}\left(\frac{3}{2\pi}\right)^{2/3}\frac{Q^{4/3}_{\rm
FB}}{R}\left[1+\frac{4\pi}{9}\left(\frac{2\pi}{3}\right)^{1/3}\frac{R^{2}T^{2}}{Q^{2/3}_{\rm
FB}}\left(1+\frac{3}{2\pi^{2}}\frac{M^{2}_{i}}{T^{2}}\right)\right]$ (13)
$\displaystyle-\frac{3g^{2}_{\chi}}{8\pi}\frac{Q^{2}_{\rm
FB}L^{2}_{\phi}}{R^{3}}+\frac{4\pi}{3}V_{0}R^{3}\left(1+\frac{T^{2}M^{2}_{i}}{12V_{0}}\right)\,,$

where $V_{0}(T)\equiv V_{\rm eff}(0,T)-V_{\rm eff}(v_{\phi}(T),T)$, which at
zero temperature returns into $B$. Since the FB is a macroscopic object, we
ignore the contribution from the surface tension. Inside a FB, each $\chi$ has
mass $M_{i}$ and couples to others via the attractive Yukawa interaction
$g_{\chi}\phi\bar{\chi}\chi$. The length of Yukawa interaction is given as

$\displaystyle
L_{\phi}(T)\equiv\left(\frac{d^{2}V_{\mathrm{eff}}}{d\phi^{2}}\bigg{|}_{\phi=0}\right)^{-1/2}=(2D(T^{2}-T^{2}_{0}))^{-1/2}\,.$
(14)

To find $R_{\rm FB}$, we require $dE_{\rm FB}/dR=0$, which yields a cubic
polynomial equation with three solutions. The largest among the three roots
gives the FB radius $R_{\rm FB}$. Evaluating $E_{\rm FB}$ at $R_{\rm FB}$
gives rise the mass of the FB, we have

$\displaystyle R_{\rm FB}$ $\displaystyle=$
$\displaystyle\left[\frac{1}{4}\left(\frac{3}{2\pi}\right)^{2/3}\frac{Q^{4/3}_{\rm
FB}}{V_{0}}\right]^{1/4}X^{1/2}$ (15) $\displaystyle M_{\rm FB}$
$\displaystyle=$ $\displaystyle\frac{3}{4}Q_{\rm
FB}\left(9\pi^{2}V_{0}\right)^{1/4}X^{-3/2}$ (16) $\displaystyle\times$
$\displaystyle\left\\{X+\frac{4}{9}\left(1+\frac{T^{2}M^{2}_{i}}{12V_{0}}\right)X^{3}+\left(\frac{2\pi
T^{2}}{9V^{1/2}_{0}}+\frac{M^{2}_{i}}{3\pi
V^{1/2}_{0}}\right)X^{2}-\frac{2g^{2}_{\chi}L^{2}_{\phi}V^{1/2}_{0}}{3\pi}\right\\}\,,$

where denote

$\displaystyle X$ $\displaystyle\equiv$
$\displaystyle\left(1+\frac{T^{2}M^{2}_{i}}{12V_{0}}\right)^{-1}$
$\displaystyle\times$
$\displaystyle\Bigg{[}\left(1+\frac{13}{108}\frac{T^{2}M^{2}_{i}}{V_{0}}+\frac{\pi^{2}}{81}\frac{T^{4}}{V_{0}}+\frac{1}{36\pi^{2}}\frac{M^{4}_{i}}{V_{0}}\right)^{1/2}\cos\theta-\frac{\pi}{18}\frac{T^{2}}{V^{1/2}_{0}}\left(1+\frac{3}{2\pi^{2}}\frac{M^{2}_{i}}{T^{2}}\right)\Bigg{]}\,.$

To make a FB stable in the false vacuum, it further requires the conditions
Kawana:2021tde

$\displaystyle\frac{dM_{\rm FB}}{dQ_{\rm FB}}\leq
g_{\chi}v_{\phi}+M_{i}\,,~{}~{}~{}~{}\frac{d^{2}M_{\rm FB}}{dQ^{2}_{\rm
FB}}\leq 0\,,$ (17)

where the second one is automatically satisfied by the surface tension.
However, in our region of interest, non-zero $M_{i}$ is always necessary to
full fill the first condition with perturbative Yukawa coupling
$|g_{\chi}|\leq\sqrt{4\pi}$.

If we neglect the Yukawa energy, $X$ reduces into

$\displaystyle X$ $\displaystyle\simeq$
$\displaystyle\frac{\sqrt{3}}{2}\left(1-\frac{1}{4\sqrt{3}\pi}\frac{M^{2}_{i}}{V^{1/2}_{0}}-\frac{\pi}{6\sqrt{3}}\frac{T^{2}}{V^{1/2}_{0}}-\frac{5}{216}\frac{M^{2}_{i}T^{2}}{V_{0}}\right)\,.$
(18)

Under the assumption $V_{0}\gg T$ and $M_{i}$, Eq. (16) can be simplified into

$\displaystyle R_{\rm FB}$ $\displaystyle=$
$\displaystyle\left[\frac{3}{16}\left(\frac{3}{2\pi}\right)^{2/3}\frac{Q^{4/3}_{\rm
FB}}{V_{0}}\right]^{1/4}\left(1-\frac{1}{8\sqrt{3}\pi}\frac{M^{2}_{i}}{V^{1/2}_{0}}-\frac{\pi}{12\sqrt{3}}\frac{T^{2}}{V^{1/2}_{0}}-\frac{5}{432}\frac{M^{2}_{i}T^{2}}{V_{0}}\right)\,,$
(19) $\displaystyle M_{\rm FB}$ $\displaystyle=$ $\displaystyle Q_{\rm
FB}\left(12\pi^{2}V_{0}\right)^{1/4}\left(1+\frac{\sqrt{3}}{8\pi}\frac{M^{2}_{i}}{V^{1/2}_{0}}+\frac{\pi}{4\sqrt{3}}\frac{T^{2}}{V^{1/2}_{0}}-\frac{1}{16}\frac{M^{2}_{i}T^{2}}{V_{0}}\right)\,.$
(20)

Due to the fact that the magnitude of Yukawa energy

$\displaystyle|E_{Y}|\simeq\frac{3g^{2}_{\chi}}{8\pi}\frac{Q^{2}_{\rm
FB}}{R}\,\left(\frac{L_{\phi}}{R}\right)^{2}$ (21)

increases as the temperature decreases, at temperature $T_{\phi}$ when
$|E_{Y}|$ equals to the Fermi-gas kinetic energy, the FB becomes unstable and
starts collapsing to a PBH. This roughly coincides when $L_{\phi}$ equals to
the mean separation distance of $\chi$’s, i.e., $L_{\phi}\simeq R_{\rm
FB}/Q^{1/3}_{\rm FB}$. Therefore, there are two scenarios of PBH formation: i)
If $T_{\phi}\leq T_{\star}$, the FB forms as an intermediate state before
collapsing into PBH. ii) If $T_{\phi}>T_{\star}$, the $\chi$’s enclosed in a
critical volume directly collapse into PBH without FB formation. In the first
scenario, we can adopt the above formulas of FB to compute the mass of PBH,
i.e. $M_{\rm PBH}(T_{\phi})=M_{\rm FB}(T_{\phi})$, and the number density
follows the adiabatic evolution of Universe $n_{\rm PBH}|_{T_{\phi}}=n_{\rm
FB}|_{T_{\star}}s(T_{\phi})/s(T_{\star})$ Kawana:2021tde . Consequently, the
PBH relic abundance and fraction at present Universe is given by

$\displaystyle\Omega_{\rm PBH}h^{2}=\frac{M_{\rm PBH}|_{0}n_{\rm
PBH}|_{0}}{3M^{2}_{\rm Pl}(H_{0}/h)^{2}}\,,~{}f_{\rm
PBH}\equiv\frac{\Omega_{\rm PBH}h^{2}}{\Omega_{\rm DM}h^{2}}\,,$ (22)

where the Hubble constant $H_{0}=2.13h\times 10^{-42}~{}{\rm GeV}$. In the
following section, we will select the benchmark points belong to the first
scenario.

##  4 Correlated signals

Figure 3: The extragalactic photon contributions of BPs in Table 1.  Figure 4:
Gravitational wave power spectra for BPs in Table 1.

Since the FOPT behaves as common origin for both PBHs and GW, it predicts the
correlated signals among the 511 keV line, extragalactic photon spectrum, and
GW spectrum. More specifically, the quartic potential dictates the FOPT and
then the Yukawa interaction determines the period of PBHs formation.
Therefore, we select and scan the input parameters: including the coefficients
of the effective potential, the asymmetry parameter $\eta_{\mathrm{DM}}$, the
temperature ratio of the dark and SM sectors $r_{T}$, the Yukawa coupling
$g_{\chi}$, and the bare DM mass $M_{i}$, in the ranges

$\displaystyle 0.1\leq B^{1/4}/{\rm MeV}\leq 10^{4},\mkern 15.0mu0.1\leq D\leq
10,\mkern 15.0mu0.05\leq\lambda\leq 0.2$ $\displaystyle 0.01\leq C/{\rm
MeV}\leq 10^{4},\mkern 15.0mu0.3\leq r_{T}\leq 1,\mkern 15.0mu0.01\leq
g_{\chi}\leq\sqrt{4\pi}$ $\displaystyle 10^{-3}\leq M_{i}/B^{1/4}\leq 10\,,$
(23)

where $A$ is fixed to be $0.1$. We pick six benchmark points (BPs) satisfying
the 511 keV excess with the corresponding PBH mass and $f_{\mathrm{PBH}}$ are
listed in Table 1 and Fig. 2.

| BP-1 | BP-2 | BP-3 | BP-4 | BP-5 | BP-6
---|---|---|---|---|---|---
$B/\rm MeV^{4}$ | $5.94\times 10^{4}$ | $4.13\times 10^{5}$ | $5.69\times 10^{2}$ | $4.80\times 10^{5}$ | $3.99\times 10^{7}$ | $5.19\times 10^{6}$
$\lambda$ | $0.1693$ | $0.0888$ | $0.0999$ | $0.0795$ | $0.0685$ | $0.1903$
$D$ | $1.9103$ | $0.5360$ | $3.5249$ | $1.6074$ | $1.3413$ | $0.8706$
$\eta_{\mathrm{DM}}$ | $4.70\times 10^{-14}$ | $3.60\times 10^{-13}$ | $6.70\times 10^{-11}$ | $7.60\times 10^{-11}$ | $4.35\times 10^{-11}$ | $1.08\times 10^{-9}$
$r_{T}$ | $0.4377$ | $0.3512$ | $0.3149$ | $0.3845$ | $0.3891$ | $0.4707$
$C/\rm MeV$ | $1.5257$ | $1.3976$ | $0.0920$ | $0.5830$ | $2.8082$ | $3.6441$
$g_{\chi}$ | $1.2728$ | $1.1556$ | $1.0275$ | $1.1958$ | $0.9141$ | $1.3278$
$M_{i}/B^{1/4}$ | $0.5570$ | $0.0516$ | $0.4973$ | $0.0047$ | $0.3048$ | $1.4341$
$M_{\mathrm{PBH}}/M_{\bigodot}$ | $9.32\times 10^{-19}$ | $4.39\times 10^{-18}$ | $1.52\times 10^{-17}$ | $3.72\times 10^{-17}$ | $4.90\times 10^{-17}$ | $7.04\times 10^{-17}$
$f_{\mathrm{PBH}}$ | $4.74\times 10^{-6}$ | $8.20\times 10^{-5}$ | $0.0013$ | $0.0100$ | $0.0235$ | $0.4350$

Table 1:  The benchmark points form PBHs after FOPT with $A=0.1$ fixed.

According to the BPs in Table 1, in order to produce the desired PBH mass, one
requires the energy scale $\mathcal{O}(1)\lesssim B^{1/4}/{\rm
MeV}\lesssim\mathcal{O}(100)$ of FOPT. By incorporating $\eta_{\rm DM}$, we
can obtain the correct values of $f_{\rm PBH}$ to generate 511 keV line. Under
our scenario, FB was formed as an intermediate state and treated as a
progenitor of PBH. In order to satisfy the stability condition of FB, i.e
Eq.(17), all BPs tend to have large values of $g_{\chi}$, close to the
perturbative limit, and nonzero values of $M_{i}$ are necessary. In addition,
the observational upper bound of $\Delta N_{\rm eff}\leq 0.5$ was imposed to
all the BPs, and thereby restricts the temperature ratio to be $0.3\leq
r_{T}\leq 0.5$. Here, we assumed the dark sector and SM sector are thermally
decoupled, such that it suppresses the dark sector contribution to the light
degree of freedom.

The extragalactic photon contributions of these BPs are shown in Fig. 3. BP-4,
BP-5 and BP-6 haven’t been rule out by the present observations
COMPTEL/EGRET/FermiLAT. This is because the associated PBH masses from BP-4,
BP-5 and BP-6 are heavier comparing with BP-1, BP-2, and BP-3. Base on the
property of Hawking radiation, the gamma-ray spectra of them peak at lower
energy window between 0.1 MeV to 1 MeV, where the present observations do not
have sufficient sensitivities. However, this window will be explored by future
AMEGO and e-ASTROGAM observations. On the other hand, the GW spectra of BPs,
mainly arising from sound waves in the plasma after bubbles collision during
the FOPT Caprini:2015zlo , are shown in Fig. 4. They cover the frequency
region from $10^{-5}$ Hz to $10^{-3}$ Hz, and their amplitudes can be
substantially detected by future $\mu$Ares telescope.

##  5 Conclusion

PBH with mass $3\times 10^{-17}\lesssim M_{\rm PBH}/m_{\odot}\lesssim 7\times
10^{-17}$ and fractional abundance $0.01\lesssim f_{\rm PBH}\lesssim 0.5$ is
favored to explain the galaxy center 511 keV line excess. PBH evaporation
emits not only positrons but also photons contributing to the extragalactic
gamma-ray flux. We investigated the novel PBH production mechanism from the
cosmological FOPT aggregating $\chi$’s into FB as intermediate state. Through
the attractive Yukawa interaction, FBs become unstable as temperature
decreasing and eventually collapse to form PBHs. To realise this scenario, we
selected the BPs satisfying several conditions: $\Delta N_{\rm eff}\leq 0.5$,
stability of FB in the false vacuum, and upper limits of gamma-ray flux from
EGERT/COMPTEL/FermiLAT. Incorporate the $\chi$ asymmetry $\eta_{\rm DM}\simeq
10^{-10}$, the phase transition with vacuum energy $\mathcal{O}(1)\lesssim
B^{1/4}/{\rm MeV}\lesssim\mathcal{O}(100)$ produces the desired PBH mass and
abundance for 511 keV excess. Consequently, correlated signals of
extragalactic $\mathcal{O}({\rm MeV})$ gamma-ray and GW spectrum with peak
frequency from $10^{-5}$ Hz to $10^{-3}$ Hz are predicted. In the future,
these signals can be either confirmed or ruled out by AMEGO/e-ASTROGAM and
$\mu$Ares.

## Acknowledgment

P.Y.Tseng is supported in part by the Ministry of Sciences and Technology with
Grant No. MoST-111-2112-M-007-012-MY3. Y.M.Yeh is supported in part by
Ministry of Education with Grant No. 111J0382I4.

## References

  * (1) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * (2) G. F. Chapline, Nature 253, no.5489, 251-252 (1975).
  * (3) M. Y. Khlopov, Res. Astron. Astrophys. 10, 495-528 (2010), [arXiv:0801.0116 [astro-ph]].
  * (4) B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no.8, 083504 (2016), [arXiv:1607.06077 [astro-ph.CO]].
  * (5) B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Rept. Prog. Phys. 84, no.11, 116902 (2021) [arXiv:2002.12778 [astro-ph.CO]].
  * (6) B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355-394 (2020), [arXiv:2006.02838 [astro-ph.CO]].
  * (7) A. M. Green and B. J. Kavanagh, J. Phys. G 48, no.4, 043001 (2021), [arXiv:2007.10722 [astro-ph.CO]].
  * (8) S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142-147 (2017), [arXiv:1603.05234 [astro-ph.CO]].
  * (9) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, 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]].
  * (10) 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]].
  * (11) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.6, 061102 (2016), [arXiv:1602.03837 [gr-qc]].
  * (12) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.24, 241103 (2016), [arXiv:1606.04855 [gr-qc]].
  * (13) B. P. Abbott et al. [LIGO Scientific and VIRGO], Phys. Rev. Lett. 118, no.22, 221101 (2017) [erratum: Phys. Rev. Lett. 121, no.12, 129901 (2018)], [arXiv:1706.01812 [gr-qc]].
  * (14) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399-415 (1974)
  * (15) M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no.6, 063001 (2018), [arXiv:1801.05235 [astro-ph.CO]].
  * (16) S. W. Hawking, I. G. Moss and J. M. Stewart, Phys. Rev. D 26, 2681 (1982).
  * (17) I. G. Moss, Phys. Rev. D 50, 676-681 (1994).
  * (18) R. V. Konoplich, S. G. Rubin, A. S. Sakharov and M. Y. Khlopov, Phys. Atom. Nucl. 62, 1593-1600 (1999).
  * (19) H. Kodama, M. Sasaki and K. Sato, Prog. Theor. Phys. 68, 1979 (1982).
  * (20) C. Gross, G. Landini, A. Strumia and D. Teresi, JHEP 09, 033 (2021) [arXiv:2105.02840 [hep-ph]].
  * (21) M. J. Baker, M. Breitbach, J. Kopp and L. Mittnacht, [arXiv:2105.07481 [astro-ph.CO]].
  * (22) K. Kawana and K. P. Xie, Phys. Lett. B 824, 136791 (2022) [arXiv:2106.00111 [astro-ph.CO]].
  * (23) J. P. Hong, S. Jung and K. P. Xie, Phys. Rev. D 102, no. 7, 075028 (2020), [arXiv:2008.04430 [hep-ph]].
  * (24) D. Marfatia and P. Y. Tseng, JHEP 11, 068 (2021), [arXiv:2107.00859 [hep-ph]].
  * (25) E. Witten, Phys. Rev. D 30, 272 (1984).
  * (26) Y. Bai, A. J. Long and S. Lu, Phys. Rev. D 99, no. 5, 055047 (2019), [arXiv:1810.04360 [hep-ph]].
  * (27) L. Bouchet, J. P. Roques and E. Jourdain, Astrophys. J. 720, 1772-1780 (2010), [arXiv:1007.4753 [astro-ph.HE]].
  * (28) D. Marfatia and P. Y. Tseng, JHEP 08, 001 (2022), [arXiv:2112.14588 [hep-ph]].
  * (29) G. Weidenspointner, V. Lonjou, J. Knodlseder, P. Jean, M. Allain, P. von Ballmoos, M. J. Harris, G. K. Skinner, G. Vedrenne and B. J. Teegarden, et al. [arXiv:astro-ph/0406178 [astro-ph]].
  * (30) E. Churazov, R. Sunyaev, S. Sazonov, M. Revnivtsev and D. Varshalovich, Mon. Not. Roy. Astron. Soc. 357, 1377-1386 (2005), [arXiv:astro-ph/0411351 [astro-ph]].
  * (31) G. Weidenspointner, J. Knoedlseder, P. Jean, G. K. Skinner, J. P. Roques, G. Vedrenne, P. Milne, B. J. Teegarden, R. Diehl and A. Strong, et al. ESA Spec. Publ. 622, 25 (2007) [arXiv:astro-ph/0702621 [astro-ph]].
  * (32) P. Jean, J. Knodlseder, W. Gillard, N. Guessoum, K. Ferriere, A. Marcowith, V. Lonjou and J. P. Roques, Astron. Astrophys. 445, 579-589 (2006), [arXiv:astro-ph/0509298 [astro-ph]].
  * (33) N. Prantzos, Astron. Astrophys. 449, 869-878 (2006), [arXiv:astro-ph/0511190 [astro-ph]].
  * (34) N. Prantzos, C. Boehm, A. M. Bykov, R. Diehl, K. Ferriere, N. Guessoum, P. Jean, J. Knoedlseder, A. Marcowith and I. V. Moskalenko, et al. Rev. Mod. Phys. 83, 1001-1056 (2011), [arXiv:1009.4620 [astro-ph.HE]].
  * (35) R. J. Wilkinson, A. C. Vincent, C. Bœhm and C. McCabe, Phys. Rev. D 94, no.10, 103525 (2016), [arXiv:1602.01114 [astro-ph.CO]].
  * (36) S. W. Hawking, Commun. Math. Phys. 43, 199-220 (1975) [erratum: Commun. Math. Phys. 46, 206 (1976)]
  * (37) S. W. Hawking, Nature 248, 30-31 (1974).
  * (38) S. W. Hawking, Commun. Math. Phys. 43, 199-220 (1975).
  * (39) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2738-2751 (1977)
  * (40) C. Keith and D. Hooper, Phys. Rev. D 104, no.6, 063033 (2021), [arXiv:2103.08611 [astro-ph.CO]].
  * (41) P. H. Frampton and T. W. Kephart, Mod. Phys. Lett. A 20, 1573-1576 (2005), [arXiv:hep-ph/0503267 [hep-ph]].
  * (42) C. Bambi, A. D. Dolgov and A. A. Petrov, Phys. Lett. B 670, 174-178 (2008) [erratum: Phys. Lett. B 681, 504-504 (2009)], [arXiv:0801.2786 [astro-ph]].
  * (43) R. G. Cai, Y. C. Ding, X. Y. Yang and Y. F. Zhou, JCAP 03, 057 (2021), [arXiv:2007.11804 [astro-ph.CO]].
  * (44) R. Laha, Phys. Rev. Lett. 123, no.25, 251101 (2019), [arXiv:1906.09994 [astro-ph.HE]].
  * (45) W. DeRocco and P. W. Graham, Phys. Rev. Lett. 123, no.25, 251102 (2019), [arXiv:1906.07740 [astro-ph.CO]].
  * (46) A. C. Robin, C. Reyle, S. Derriere and S. Picaud, Astron. Astrophys. 409, 523 (2003), [arXiv:astro-ph/0401052 [astro-ph]].
  * (47) A. Arbey and J. Auffinger, Eur. Phys. J. C 79, no.8, 693 (2019), [arXiv:1905.04268 [gr-qc]].
  * (48) A. Arbey and J. Auffinger, Eur. Phys. J. C 81, 10 (2021) [arXiv:2108.02737 [gr-qc]].
  * (49) M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde and D. A. Linde, Phys. Rev. D 46, 550-571 (1992) [arXiv:hep-ph/9203203 [hep-ph]].
  * (50) F. C. Adams, Phys. Rev. D 48, 2800 (1993), [hep-ph/9302321].
  * (51) C. Caprini et al., JCAP 1604, 001 (2016), [arXiv:1512.06239 [astro-ph.CO]].
  * (52) Y. Nakai, M. Suzuki, F. Takahashi and M. Yamada, Phys. Lett. B 816, 136238 (2021), [arXiv:2009.09754 [astro-ph.CO]].
  * (53) S. J. Huber and T. Konstandin, JCAP 0809, 022 (2008), [arXiv:0806.1828 [hep-ph]].
  * (54) J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, JCAP 1006, 028 (2010), [arXiv:1004.4187 [hep-ph]].
  * (55) J. Kehayias and S. Profumo, JCAP 1003, 003 (2010), [arXiv:0911.0687 [hep-ph]].
  * (56) D. Marfatia and P. Y. Tseng, JHEP 2102, 022 (2021), [arXiv:2006.07313 [hep-ph]].


#  Data-Driven Quantum Approximate Optimization Algorithm for Cyber-Physical
Power Systems

Hang Jing1, Ye Wang2,3*, Yan Li1†
1Department of Electrical Engineering, The Pennsylvania State University,
University Park, PA, 16802 USA.
2Duke Quantum Center, Duke University, Durham, NC, 27708 USA.
3Department of Electrical & Computer Engineering, Duke University, Durham, NC,
27708 USA.
*Corresponding author<EMAIL_ADDRESS>
†Corresponding author<EMAIL_ADDRESS>
###### Abstract

Quantum technology provides a ground-breaking methodology to tackle
challenging computational issues in power systems, especially for Distributed
Energy Resources (DERs) dominant cyber-physical systems that have been widely
developed to promote energy sustainability. The systems’ maximum power or data
sections are essential for monitoring, operation, and control, while high
computational effort is required. Quantum Approximate Optimization Algorithm
(QAOA) provides a promising means to search for these sections by leveraging
quantum resources. However, its performance highly relies on the critical
parameters, especially for weighted graphs. We present a data-driven QAOA,
which transfers quasi-optimal parameters between weighted graphs based on the
normalized graph density, and verify the strategy with $39,774$ instances.
Without parameter optimization, our data-driven QAOA is comparable with the
Goemans-Williamson algorithm. This work advances QAOA and pilots the practical
application of quantum technique to power systems in noisy intermediate-scale
quantum devices, heralding its next-generation computation in the quantum era.

##  I Introduction

Quantum technology is emerging as a new hope to address challenging
computational tasks in power systems, including quantum chemistry simulation
for new type batteries [1, 2, 3], efficient power system analysis by solving
linear systems of equations [4, 5, 6, 7, 8], forecasting highly chaotic
systems [9], scheduling and dispatching power grids [10], unit commitment
[11], optimal reconfiguration of distribution grids [12], etc. However, the
existing algorithms require substantial quantum resources, limiting their
near-term utilization on noisy intermediate-scale quantum (NISQ) devices [13].
Even though specific instances of quantum algorithms have been demonstrated on
various quantum processors with tens of qubits [14, 15, 16], practical
applications to address power system problems will still require further
advances in algorithmic design.

In power systems, one emerging quantum application is to analyze the
Distributed Energy Resources(DERs) dominant power system, which provides a
potent solution to seek an edge toward energy sustainability. We illustrate a
typical DER dominant power system in Fig. 1, where the physical system is
energized by DERs and the cyber-layer enables the communication among DERs for
coordination and control.

Figure 1: Illustration of the cyber-physical power system. It includes
physical layer and cyber layer. The physical layer is DER dominant power grid.
The cyber layer is used for the communication among DERs and the control
center for the system’s operation and control.

To improve the resiliency of the cyber-physical system, it is critical to
efficiently obtain the maximum sections of power energy in the physical layer
and data traffic in the cyber-layer [17, 18, 19]. Mathematically, finding the
maximum section of power energy or data traffic is to solve the Max-Cut
problem, an NP-hard problem [20]. Therefore people implement classical
approximation algorithms [21, 22, 23, 24] to address the Max-Cut problem in
practical applications. However, for specific instances, classical algorithms
can only guarantee an approximation ratio of $0.878$ [21, 25].

The Quantum Approximate Optimization Algorithm (QAOA), a hybrid quantum-
classical algorithm, is expected to obtain better approximate solutions than
any existing classical algorithms [26, 27]. QAOA utilizes a classical computer
trains the parameters for quantum circuit [26]. The parameterized quantum
circuit approximates the adiabatic evolution from an initial Hamiltonian,
whose ground energy state is easy to prepare, to a final Hamiltonian, whose
ground energy state encodes the solution of the Max-Cut problem. With an ideal
approximation, people expect to obtain the exact solution of the Max-Cut
problem with high probabilities [28]. Consequently, the parameters involved in
the quantum circuit play an essential role in getting high-quality
approximations [29, 30, 31]. However, how to efficiently obtain appropriate
parameters is still an open question.

This work presents a data-driven QAOA with parameter transfer strategy to
tackle the challenging issue of efficiently obtaining appropriate parameters.
Based on the normalized graph densities [32], the quasi-optimal parameters for
seed (or existing) graphs are transferred to target (or new) graphs, which can
be directly applied to obtain the Max-Cut or as an initial guess for further
optimizing iterations. The transfer strategy is extendable and designed for
generic weighted graph, which is extension studies of regular graphs [33, 34].
We numerically justify the strategy by the QAOA’s performance on $1710$ random
instances under transferred parameters and the ones after optimization.
Afterwards, we perform the data-driven QAOA in practical cyber-physical power
systems to get the maximum sections. Simulations on $996$ case studies have
validated that the data-driven QAOA can efficiently obtain comparable results
to Goemans-Williamson(GW) algorithm, a famous classical algorithm [21].
Additionally, simulation results show that the near-term achievable noise of
quantum processor is negligible to the data-driven QAOA’s performance. This
work would significantly reduce the computational effort for training QAOA
parameters, advance the development of QAOA, and highly promote its wide
applications for solving engineering problems. This is the first practical
quantum application that is feasible shortly in the NISQ era to address
problems in power systems.

##  II Maximum Sections Problem Formulation

The maximum power or data section is a segment along which the power delivery
in the power grid or data traffic in the communication network reaches the
maximum value. It is critical to efficiently obtain the maximum sections
because of three reasons. First, the maximum power energy section offers a
cost-effective way to monitor the dynamics and power delivery capability of
the physical system considering the fluctuations of DERs and/or the frequent
changes of system topology due to the join or removal of subsystems (e.g.,
microgrids) [17]. Second, the maximum power energy sections cast light on the
dynamic system’s control and operations. Dispatchable DERs can be coordinated
for reducing the electric power over the maximum section to avoid system
collapse [18]. Third, the maximum data traffic sections provide an insight
into enhancing the overall system’s resilience through strategically designing
and managing the communication network [19], e.g., packet routing and traffic
control.

Mathematically, finding the maximum section is to solve a Max-Cut problem of a
weighted graph $G=(V,E)$, where $\absolutevalue{V}=n$ is the vertex number,
$\absolutevalue{E}=m$ is the edge number, and $w_{ij}$ represents the
normalized weight of the edge $\left\langle i,j\right\rangle\in E$, where
${\rm{Max}}({w_{ij}})=1$. The Max-Cut solutions are identical before and after
normalization. The edge weight is obtained via power flow calculation for the
physical layer and means the data traffic in the cyber layer. Modeling details
are provided in the Methods Section F.

The objective is to find a subset $S\subset V$ that maximizes $\sum_{i\in
S,j\notin S}{w_{ij}}$ for cyber or physical layers, respectively. Suppose an
$n$-bit string $Z=z_{1}\cdots z_{i}\cdots z_{j}\cdots z_{n}\in\\{-1,1\\}^{n}$
can denote the status of vertices $V$, showing each bit $z_{i}$ will be equal
to $1$ if the $i^{\text{th}}$ vertex is in the subset $S$, otherwise be $-1$.
We can exhibit the partition of vertices for obtaining the maximum section.
Thus, the classical cost function of the Max-Cut problem can be defined as,

$C(Z)=\sum_{{\left\langle i,j\right\rangle}\in
E}w_{ij}\frac{1-z_{i}z_{j}}{2}=\sum_{{\left\langle i,j\right\rangle}\in
E}w_{ij}C_{ij}(Z),$ (1)

where $C_{ij}(Z)$ represents the contribution of $w_{ij}$ to the cost
function. The Max-Cut problem translates into finding the $n$-bit string $Z$
to maximize the cost function $C(Z)$. Given the $n$-bit string $Z$, we define
the approximation ratio to be $C(Z)/C(Z_{\rm Max-Cut})$, where $Z_{\rm Max-
Cut}$ is the exact Max-Cut solution. The goal of approximate algorithms is to
find the solution with a high approximation-ratio.

Figure 2: Schematic of the proposed data-driven QAOA for obtaining the maximum
sections of cyber-physical power systems. The key idea includes the following
five steps: ① By modeling the _Cyber-physical systems_ into two normalized
weighted graphs, we compute the _Normalized graph_ density from adjacency
matrix $W_{\rm{adj}}$. ② In the _Parameter transfer module_ , we first
determine the seed graph size $n_{\rm{s}}$ and layer number $p$. Then,
according to the mapping table $M_{p}{(n_{\rm{s}}\mapsto n_{\rm{t}})}$ in
Extended Fig. References, we can obtain the quasi-optimal parameters
$(\gamma,\beta)$ from seed graphs, whose normalized densities can be expressed
by an interval $[\tilde{D}_{l},\tilde{D}_{r}]$. ③ The transferred parameters
$(\gamma,\beta)$ are directly passed to the _Quantum circuit with multiple
layers for QAOA_. By measurement, it generates the probability distribution
$|\alpha_{k}|^{2}$ in (4), from which we can obtain $C(\ket{Z})$ and select a
high approximation-ratio solution. ④ If a better performance is desired, we
will further _Optimize $(\gamma,\beta)$ for $C(\ket{Z})_{\rm{max}}$_. This
step is optional. ⑤ The obtained pair $(n,D,\gamma,\beta)$ can be used to
develop an _Expandable quasi-optimal parameter database_ to provide quasi-
optimal parameter for new target graphs. If step ④ is performed, by obtaining
the optimized parameters, we can store a new pair
$(n_{\rm{t}},D,\gamma_{\rm{new}},\beta_{\rm{new}})$ in a new mapping table for
target graphs; otherwise, we can add one new entry $C(\ket{Z})$ to the
original mapping table, as exampled in Extended Fig. References. This step is
also optional.

##  III QAOA for Max-cut Problem

On quantum computers, we use $n$ quantum bits (qubits)
$\ket{Z}=\ket{z_{1}\cdots z_{i}\cdots z_{j}\cdots z_{n}}$ to represent the
status of $n$ vertexes. Each qubit $\ket{z_{i}}$ can be a superposition of
quantum states $\ket{0}$ and $\ket{1}$, denoted as
$\ket{z_{i}}=a_{i}\ket{0}+b_{i}\ket{1}$, where $\ket{0}$ and $\ket{1}$ are the
eigenstates of the Pauli-Z operator $\sigma^{z}$ with the eigenvalues of $1$
and $-1$ respectively. When we measure the qubit in the computational basis,
which is $z$ basis, according to the quantum mechanics, the qubit could
collapse to the state $\ket{0}$ with probability of $|a_{i}|^{2}$ and the
state $\ket{1}$ with probability of $|b_{i}|^{2}$. Therefore, unlike classical
computers, the measurement results could vary even though the qubit is
identical at each execution. If we consider that measuring $\ket{0}$
represents $z_{i}=1$ and measuring $\ket{1}$ represents $z_{i}=-1$, we can
obtain various $n$-bit strings $Z=z_{1}\cdots z_{i}\cdots z_{j}\cdots
z_{n}\in\\{-1,1\\}^{n}$ and calculate $C(Z)$ in (1) after every single quantum
computer execution.

On the other hand, we could obtain the deterministic $n$-bit string $Z_{k}$
out of the measurements on $2^{n}$ $n$-qubit eigenstates in the computational
basis, denoted as $\ket{Z_{k}}$ with $\ket{z_{k,i}}=\ket{0}$ or $\ket{1}$ and
$z_{k,i}=\bra{z_{k,i}}\sigma_{i}^{z}\ket{z_{k,i}}$. Therefore, we can have

$\displaystyle C(Z_{k})$ $\displaystyle=$ $\displaystyle\sum_{{\left\langle
i,j\right\rangle}\in E}w_{ij}\frac{1-z_{k,i}z_{k,j}}{2}$ (2) $\displaystyle=$
$\displaystyle\sum_{{\left\langle i,j\right\rangle}\in
E}w_{ij}\frac{1-\bra{z_{k,i}}\sigma_{i}^{z}\ket{z_{k,i}}\bra{z_{k,j}}\sigma_{i}^{z}\ket{z_{k,j}}}{2}$
$\displaystyle=$ $\displaystyle\bra{Z_{k}}H_{\rm{C}}\ket{Z_{k}}$
$\displaystyle\equiv$ $\displaystyle C(\ket{Z_{k}}),$

where

$H_{\rm{C}}=\sum_{{\left\langle i,j\right\rangle}\in
E}w_{ij}\frac{I-\sigma_{i}^{z}\sigma_{j}^{z}}{2}.$ (3)

We consider $C(\ket{Z})=\bra{Z}H_{\rm{C}}\ket{Z}$ as the quantum analog of
$C(Z)$. Then $2^{n}$ classical cost functions $C(Z_{k})$ are one-on-one mapped
to $2^{n}$ quantum cost functions $C(\ket{Z_{k}})$. The Max-Cut problem
translates into finding the quantum state $\ket{Z_{k}}$ to maximize the cost
function $C(\ket{Z_{k}})$.

The $2^{n}$ $\ket{Z_{k}}$ states form the complete basis of the $2^{n}$
Hilbert space for $n$-bit quantum states. Therefore we can decompose an
arbitrary state $\ket{Z}$ into a linear combination of $\ket{Z_{k}}$, denoted
as $\ket{Z}=\sum_{k=1}^{2^{n}}\alpha_{k}\ket{Z_{k}}$ with
$\sum_{k=1}^{2^{n}}|\alpha_{k}|^{2}=1$. The quantum cost function of an
arbitrary state $\ket{Z}$ can be written as

$\displaystyle C(\ket{Z})$ $\displaystyle=$
$\displaystyle\bra{Z}H_{\rm{C}}\ket{Z}$ (4) $\displaystyle=$
$\displaystyle(\sum_{k=1}^{2^{n}}\alpha_{k}^{*}\bra{Z_{k}})H_{\rm{C}}(\sum_{k=1}^{2^{n}}\alpha_{k}\ket{Z_{k}})$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{2^{n}}|\alpha_{k}|^{2}C(\ket{Z_{k}}).$

Since $C(\ket{Z})=\bra{Z}H_{\rm{C}}\ket{Z}\geq 0$, we have $\max
C(\ket{Z})=\max_{k=1}^{2^{n}}C(\ket{Z_{k}})$. Notably, in quantum mechanics,
$C(\ket{Z})=\bra{Z}H_{\rm{C}}\ket{Z}$ is the expectation value of system
energy for a quantum system described by Hamiltonian $H_{\rm{C}}$. The Max-Cut
problem translates into finding the maximum energy state for the quantum
system described by Hamiltonian $H_{\rm{C}}$.

QAOA utilizes a quantum circuit running on the quantum computer to approximate
an adiabatic evolution from the maximum energy state of an initial
Hamiltonian, $H_{\rm{B}}$, to the maximum energy state of the final
Hamiltonian, $H_{\rm{C}}$. For the Max-Cut problem, we particular define the
$H_{\rm{B}}$ as

$H_{\rm{B}}=\sum_{j=1}^{n}\sigma_{j}^{x}.$ (5)

According to adiabatic theorem [35], with an ideal approximation, we expect to
obtain the maximum energy state of $H_{\rm{C}}$, which leads to the exact Max-
Cut solution, with a high probability.

To implement QAOA on quantum computers, we first prepare the maximum energy
state of $H_{\rm{B}}$, ${\ket{+}}^{\otimes n}$, as the initial state for the
quantum circuit. Then we run the quantum circuit with $2p$ trainable
parameters ${{\gamma}}=\left(\gamma_{1},\gamma_{2},\ldots,\gamma_{p}\right)$
and ${{\beta}}=\left(\beta_{1},\beta_{2},\ldots,\beta_{p}\right)$ to
approximate the $p$-step Trotter expansion of the adiabatic evolution. We
measure the output state, obtain the classical $n$-bit string, and estimate
the quantum cost function using (4) with multiple executions. After that, we
use the classical-quantum hybrid optimizer iterating $2p$ parameters to
maximize the quantum cost function. Ideally, when $p$ tends to infinity, the
probability of obtaining the exact Max-Cut solution will tend to be 1. Even
with a finite $p$, measuring the final state $\ket{Z}$ of the optimized
circuit could generate high approximation-ratio solutions. More details are
discussed in the Methods Section A.

##  IV  Data-Driven QAOA

People have studied QAOA’s efficiency and accuracy in regular graphs with
constant circuit depth [26, 34, 36, 33, 37]. On the other side, the
performance of QAOA on generic weighted graphs is an open question and
challenging to estimate rigorously [38]. Heuristic strategies show potentials
to find quasi-optimal parameters with high approximation-ratio solutions, a
claim backed by numerical evidences [29, 33, 34]. However, researchers have
not exhaustively explore the heuristic strategies on generic weighted graphs
[38].

Our data-driven QAOA for generic weighted graphs can provide a high
approximation-ratio solution without parameter optimization to avoid expensive
computational effort. The data-driven QAOA is based on normalized weighted
graph density $D$ [32], which is defined as,

$D=\sum_{{\left\langle i,j\right\rangle}\in E}{\frac{2w_{ij}}{n(n-1)}},$ (6)

The data-driven QAOA includes five steps as shown in Fig. 2 and introduced as
follows, where an innovative parameter transfer strategy is the key idea.

Step 1: Formulate the search of maximum section to the Max-Cut problem of the
normalized weighted graph and calculate its $D$.

Step 2: Obtain the quasi-optimal parameters $(\gamma,\beta)$ based on $D$ via
the parameter transfer strategy, and then pass these parameters to the quantum
processors. The transfer strategy works for generic weighted graphs.

Step 3: Construct the quantum circuit with the adjacency matrix $W_{\rm{adj}}$
and parameters $(\gamma,\beta)$, and run it in a quantum processor. Then,
measure the output state of the quantum circuit to get the probability
distribution and calculate the cost function value.

Step 4: Optimize the parameters by the classical optimizer for a better result
when necessary. Step $4$ is optional.

Step 5: Expand the database by adding more pairs denoted by
$(n,D,\gamma,\beta)$ from verified cases, to provide more quasi-optimal
parameters. Step $5$ is optional.

Decent initial guesses obtained in Step 2 can also help to handle noise-free
barren plateaus, which are linked to random parameter initialization [39].
These initial guesses also significantly reduce the iterations between the
classical optimizer and the quantum processor, saving the running time of the
algorithm.

##  V Parameter Transfer Strategy

Figure 3: The mean approximation-ratio of randomly generated target graphs
under parameters of seed graphs with $n_{\rm{s1}}=10$. The mean approximation-
ratios are computed by the probability distribution, with details given in (4)
and the Methods Section B. The probability with respect to approximation-ratio
is fitted by normal distribution. The $\sigma_{\rm{window}}$ is the standard
deviation of the scatters with the same color in the $0.1$ scan window
regarding to $D$.  Figure 4: The mean approximation-ratio of randomly
generated target graphs under parameters of seed graphs with $n_{\rm{s2}}=24$.
Figure 5: The mean approximation-ratio of QAOA results after optimization.

The essential idea of the data-driven QAOA in Fig. 2 is the parameter transfer
strategy, as summarized in the Step $2$. It includes the following three
substeps, which highly improves the effectiveness of transfer.

Substep 1: Establish the initial database. Several seed graphs are randomly
generated, with their normalized graph densities spreading over $[0,1]$.
Considering the small size of these graphs, we can calculate the quasi-optimal
parameters $(\gamma,\beta)$ [33, 37] as shown in Methods Section C. These
parameters provide potentially quasi-optimal parameters for new graphs. This
feature is particularly appealing to relatively larger target graphs. The
database can then be established, based on the pairs $(n,D,\gamma,\beta)$.

Substep 2: Develop the mapping table. The mapping table is designed for
transferring quasi-optimal parameters from seed graphs to target graphs with
the same circuit layer number, $p$. For creating the mapping table, several
target graphs are also randomly generated, with normalized graph densities
spread over $[0,1]$. With the parameters obtained from the seed graphs, QAOA
calculation is performed for each target graph to get the cost function value
$C(\ket{Z})$. Assume there are $\mathcal{N}$ seed graphs with $n_{\rm{s}}$
vertices and $\mathcal{M}$ target graphs with $n_{\rm{t}}$ vertices, the
values of $C(\ket{Z})$ are then organized into a
$\mathcal{N}\times\mathcal{M}$ matrix $M_{p}{(n_{\rm{s}}\mapsto n_{\rm{t}})}$,
i.e., the mapping table. In the table, each column is corresponding to one
target graph and each row is corresponding to one seed graph. Note a mapping
table only needs to be prepared once in advance for the same $p$,
$n_{\rm{s}}$, $n_{\rm{t}}$.

Substep 3: Transfer parameters to new graphs. For a new graph with normalized
graph density $D^{\prime}$ and size $n^{\prime}_{\rm{t}}$, several appropriate
seed graphs will be selected from the mapping table $M_{p}{(n_{\rm{s}}\mapsto
n^{\prime}_{\rm{t}})}$, whose size $n_{\rm{s}}$ needs to be equal or close to
$n^{\prime}_{\rm{t}}$. Then, in the mapping table, one (or more) column, whose
corresponding $D$ is equal or close to the new graph’s $D^{\prime}$, will be
selected. Since each entry of this column is associated with a pair
$(n_{\rm{s}},\tilde{D},\gamma,\beta)$ of the seed graph, we can choose entries
that are bigger than a threshold to get quasi-optimal parameters for the new
graph. Specifically, based on the obtained entries, the parameters in the pair
corresponding to each entry will be identified and then transferred to the new
graph. In this sense, we can use an interval
$[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ to summarize the identified pairs
and denote the mapping as
$D^{\prime}\mapsto[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$, as shown in Fig.
2. To improve the result’s accuracy, the layer number $p$ can be accordingly
increased, with parameters obtained by above transfer strategy.

The presented quasi-optimal parameter database is expendable, as mentioned in
the Step $5$ of the data-driven QAOA. In the Substep $3$, we obtain the
identified quasi-optimal parameters. For one thing, these parameters can be
directly applied to the QAOA calculation of the new graph. The new result can
then be added to the current mapping table $M_{p}{(n_{\rm{s}}\mapsto
n^{\prime}_{\rm{t}})}$ as a new entry, which is associated with the pair
$(n_{\rm{s}},\tilde{D},\gamma,\beta)$. For another, these parameters can also
be decent initial guesses for further optimizing the parameters. Then based on
the optimized result, a new pair
$(n^{\prime}_{\rm{t}},D^{\prime},\gamma^{\prime},\beta^{\prime})$ can be added
to the database to provide potential parameters for new graphs, which is
equivalent to the Substep 1.

##  VI Numerical Justification of the Parameter Transfer Strategy

Since rigorously estimating the QAOA performace on generic weighted graphs is
still an open question [38], we provide numerical examples to justify the
effectiveness of the parameter transfer strategy. We verify the efficacy of
the transferred parameters from three aspects by comparing the approximation
ratios with the ones obtained by QAOA using random parameters, QAOA using
optimized parameters and GW algorithm.

For developing the mapping tables, we randomly generate $9$ unweighted graphs
with $n_{\rm{s1}}=10$ and $9$ weighted graphs with $n_{\rm{s2}}=24$ as seed
graphs, as given in Extended Table Data-Driven Quantum Approximate
Optimization Algorithm for Cyber-Physical Power Systems. The $1710$ non-planar
target graphs with $n_{\rm{t}}=24$ are also randomly generated, including
$714$ unweighted graphs and $996$ weighted ones. The justifications involve
$39,744$ QAOA expectation value calculations and at least $16,146,548,640$
shots. We apply two classical optimizers based on Newton and COBYLA methods
[40] to get the mean approximation-ratio of the seed graphs and their
corresponding quasi-optimal parameters. The values for $p=1,2,3$ are
summarized in Extended Table Data-Driven Quantum Approximate Optimization
Algorithm for Cyber-Physical Power Systems. These parameters $(\gamma,\beta)$
provide the initial data for the expandable database as introduced in Fig. 2.

According to the Substep 2, we develop $12$ mapping tables as examples, among
which the $6$ mapping tables in Extended Fig. References are for the
unweighted seed graphs under both weighted and unweighted target graphs with
$p=1,2,3$, respectively; the other $6$ mapping tables in Extended Fig.
References are for the weighted seed graphs.

###  VI-A Comparison with Using Random Parameters

We compare the QAOA results from the transferred parameters and from random
parameters to verify the parameter transfer strategy. We apply each seed
graph’s quasi-optimal parameters to the QAOA calculation for the $1,710$
target graphs, respectively. Fig. 3 summarizes the mean approximation-ratios
of QAOA in the unweighted and weighted graphs for $p=1,2,3$. In Fig. 3, each
black circle represents the mean approximation-ratio for a target graph with
the identified parameters obtained from the $9$ groups parameters in the
mapping tables developed from unweighted $n_{\rm{s1}}=10$ seed graphs. These
black circles show the high approximation-ratios, which are
$0.8501,0.8911,0.9125$ in average for $p=1,2,3$, respectively. Comparing these
black circles with the pink circles showing the approximation ratios with
random parameters, we can see that the parameter transfer strategy
significantly improves the approximation ratio, especially for low density
graphs, which verifies the effectiveness of the transfer strategy. In
addition, in Fig. 4, we also observe high approximation-ratio with parameters
transferred from weighted $n_{\rm{s2}}=24$ seed graphs.

Meanwhile, we have two insights in both Fig. 3 and 4. First, the significant
improvement in low density graph is particularly appealing to cyber-physical
systems, which usually have low densities. Second, it is challenging for the
random parameters method to efficiently handle barren plateaus, while our
method can address this issue by providing quasi-optimal initial guesses [41].
Further explanations for the effectiveness of the transfer strategy are shown
in Methods Section E.

###  VI-B Comparison with Using Optimized Parameters

We verify that the transferred parameter can provide warm starting for QAOA.

Fig. 5 compares the mean approximation-ratio of QAOA using the transferred and
unfavorable parameters in the mapping table with and without further
optimization.

On the one hand, the transferred parameters can be decent initial guesses.
Fig. 5 shows that after further optimizing the transferred parameters (blue
scatters), the result (orange scatters) has no significant improvement. On the
other hand, those transferred parameters also can be quasi-optimal parameters.
Most of blue scatters without optimization are better than the purple
scatters, which are optimized from yellow scatters with lots of computational
effort. There is a significant increment when comparing the transferred
parameters with the worst parameters in the mapping table, as shown by the
blue and yellow scatters in Fig. 5. These comparisons validate the
effectiveness of the parameter transfer strategy.

In addition, the optimized parameters associated to the orange scatters in
Fig. 5 can be adopted to expand the database.

###  VI-C Comparison with the GW Algorithm

Figure 6: The comparison of the GW algorithm and our data-driven QAOA with
different layer numbers without parameter optimization.

For further verifying our strategy can provide the promising results, we
compare the results of using the GW algorithm with the ones via the data-
driven QAOA. $p=1,2,3,10$ are adopted as examples. In Fig. 6, when $p$
increases, the overall performance of transferred parameters increases. when
$p=10$, the transferred parameters without any optimization have better mean
approximation-ratio than GW algorithm in the $113$ graphs out of total $996$
graphs. It is encouraging that, without any parameter optimization, the data-
driven QAOA is competitive with GW algorithm. We expect that proper
optimization and larger $p$ could improve the approximation ratio further.

The drop trend of approximation ratio with $p=10$ in the large graph density
regime, as shown in Fig. 6, is due to the overfitting on the seed graphs with
$D=0.9111$ and $D=1$. The $20$ parameters in the $10$-layer QAOA circuits
could be excessive to be justified for some $10$-vertex seed graphs. Notably,
$p=\mathcal{O}(\log(n))$ could be sufficient to obtain high approximation-
ratio solutions [38]. The overfitting issue could be resolved in large seed
graphs with $n$ vertexes and $p=\mathcal{O}(\log(n))$.

Table I:  The QAOA Results in Four Test Graphs with Different Normalized Graph
Densities

No. |  seed graph $\tilde{D}$ | $C(\ket{Z})(D=0.0525)$ | $C(\ket{Z})(D=0.1053)$ | $C(\ket{Z})(D=0.1143)$ | $C(\ket{Z})(D=0.3280)$
---|---|---|---|---|---
$p=1$ | $p=2$ | $p=3$ | $p=1$ | $p=2$ | $p=3$ | $p=1$ | $p=2$ | $p=3$ | $p=1$ | $p=2$ | $p=3$
1 | 0.2667 | 0.6865 | 0.7393 | 0.7724 | 0.7840 | 0.8180 | 0.8315 | 0.7828 | 0.8260 | 0.8495 | 0.8159 | 0.7581 | 0.7475
2 | 0.5333 | 0.6316 | 0.7059 | 0.7445 | 0.7577 | 0.7985 | 0.8216 | 0.7590 | 0.8150 | 0.8448 | 0.8448 | 0.8878 | 0.9107
3 | 0.6444 | 0.6135 | 0.6987 | 0.7368 | 0.7459 | 0.7941 | 0.8163 | 0.7431 | 0.8098 | 0.8384 | 0.8466 | 0.8940 | 0.9177
4 | 0.7333 | 0.6010 | 0.6991 | 0.7396 | 0.7373 | 0.7958 | 0.8176 | 0.7308 | 0.8110 | 0.8378 | 0.8450 | 0.8958 | 0.9189
5 | 0.8000 | 0.5920 | 0.7028 | 0.7460 | 0.7309 | 0.7965 | 0.8206 | 0.7216 | 0.8115 | 0.8363 | 0.8424 | 0.8963 | 0.9138
6 | 0.8667 | 0.5848 | 0.6892 | 0.7427 | 0.7258 | 0.7882 | 0.8209 | 0.7138 | 0.8027 | 0.8376 | 0.8393 | 0.8965 | 0.9146
7 | 0.9111 | 0.5793 | 0.6596 | 0.7259 | 0.7218 | 0.7769 | 0.8133 | 0.7077 | 0.7735 | 0.8306 | 0.8365 | 0.8669 | 0.9115
8 | 0.9556 | 0.5749 | 0.6448 | 0.6708 | 0.7183 | 0.7676 | 0.7876 | 0.7024 | 0.7537 | 0.7720 | 0.8337 | 0.8516 | 0.8607
9 | 1.0000 | 0.5710 | 0.4736 | 0.6386 | 0.7155 | 0.6593 | 0.7648 | 0.6979 | 0.6244 | 0.7454 | 0.8310 | 0.7617 | 0.8482

  * *

The bold expectation is the best value in the column.

##  VII Numerical Examples of Data-Driven QAOA on Cyber-Physical Power Systems

Figure 7: The comparison of the approximation-ratio distributions between the
GW algorithm and the date-driven QAOA with different layer numbers. Noise
model I is $0.1\%$ depolarizing error on single-qubit gates and $1\%$
depolarizing error on two-qubits gates; Noise model II is $0.01\%$
depolarizing error on single-qubit gates and $0.1\%$ depolarizing error on
two-qubits gates.

We test and verify the data-driven QAOA on a typical cyber-physical power
system. The physical system is a modified IEEE 24-bus system [42], as given in
Extended Fig. References. Eleven DERs are integrated into the system.
Considering the output fluctuations of DERs, the normalized graph density will
correspondingly change over time. So, two operational scenarios with
normalized graphs densities $D=0.0525$ and $D=0.1053$ are given as examples
for the test. The communication network also has $24$ vertices. Considering
the dynamic data traffic in the network, two scenarios are considered as
examples with $D=0.1143$ and $D=0.3280$, respectively. We provide the results
from the following two aspects.

###  VII-A Test without the Depolarizing Noise

We carry out the test according to the steps given in Fig. 2. Based on the
power flow calculation of the physical system or the data traffic measurement
of the cyber layer, four normalized weighted graphs can be obtained. With the
normalized graph densities, quasi-optimal parameters can be identified through
the mapping table in Extended Fig. References for the QAOA calculation. To
provide a comparison, Table I summarizes the mean approximation-ratios with
all the parameters in the mapping table when $p=1,2,3$ for the four graphs.
The highlighted results emphasize that the best results based on the mapping
table can be obtained with the transferred parameters. Thus, it justifies the
effectiveness of the parameter transfer strategy.

We also compare the data-driven QAOA’s results with the GW algorithm. Fig. 7
(a), (b), and (c) show the approximation-ratio distributions of different
normalized graph densities when $p=10$, with the following findings. First,
the results show that the approximation means are very close to those of the
GW algorithm. More importantly, Fig. 7 (b) and (c) show that the data-driven
QAOA’s results are better than the GW algorithm, as there is at least ten
times higher probability for the data-driven QAOA method than the GW algorithm
to get the highest approximation ratio, as shown in the zoom-in details. In
practice, we usually use the highest cut value as an approximate solution
instead of the mean approximation-ratio. The data-driven QAOA can be better
than the GW algorithm. Note these parameters are transferred from the mapping
tables without any further optimization. Hence, according to Section VI B,
when these parameters are used as initial guesses for further optimizing them,
the better mean approximation-ratio are $0.9569,0.9499,0.9751$ for the cases
in Fig. 7 (a), (b), and (c), respectively. Second, Table I shows that the mean
approximation-ratio will increase as $p$ increases; and thus, a relatively
larger $p$ is recommended for practical applications.

###  VII-B Test with the Depolarizing Noise

To verify the practicability of data-driven QAOA, we introduce the
depolarizing noise on quantum gates to simulate the realistic noise on quantum
simulators [43]. Two noise models are considered. The noise model I is with
$0.1\%$ depolarizing error on single-qubit gates and $1\%$ depolarizing error
on two-qubits gates, which is presently achievable. The noise model II is with
$0.01\%$ depolarizing error on single-qubit gates and $0.1\%$ depolarizing
error on two-qubits gates, which is achievable in the near term.

We carry out the numerical noise experiments on the test graphs. Fig. 7
(d)-(i) show the examples of three graphs with $D=0.1053$, $D=0.1143$, and
$D=0.3280$, under the transferred parameters when $p=3$. By comparing the
approximation-ratio distributions and means between the results with and
without noise, we can see that the mean approximation-ratios drop, with noise
model I, as given in Fig. 7 (d), (e), and (f). While with the smaller noise,
the reduction of mean approximate-ratio is negligible, as shown in Fig. 7 (g),
(h), and (i). Therefore, it is feasible to run the data-driven QAOA on a NISQ
quantum processor and address the Max-Cut problem in the practical power
system in the near term.

##  VIII Conclusions

We present a data-driven QAOA to efficiently search for the maximum power or
data sections in DER dominant cyber-physical power systems by leveraging
quantum advantage. The parameter transfer strategy is designed to provide
quasi-optimal parameters from seed graphs to target graphs. It addresses the
challenge of obtaining the critical parameters in QAOA; and thus, highly
improving the efficacy and efficiency of QAOA. In the transfer strategy,
normalized graph density is utilized to bridge the seed and target graphs for
developing an extendable mapping table. We have verified the transfer strategy
by comparing our approximation ratios with those obtained by QAOA using random
parameters, QAOA using optimized parameters and GW algorithm. The parameter
transferability has also been verified from two perspectives, namely between
unweighted and weighted graphs and between small scale and large scale graphs
as well as graphs with the same size. We simulate the presented method in a
modified IEEE 24-bus system and demonstrated its effectiveness in finding the
maximum sections with and without depolarizing noise. The presented method
showcases the new computation of power systems when meeting quantum
technology. As a promising early candidates for achieving quantum advantage on
NISQ systems, it can also be extended to address challenging issues in other
complex engineered systems and eventually evolve into a formal quantum
methodology.

##  IX Methods

###  IX-A Adiabatic Approximation with QAOA

According to the adiabatic evolution theorem [28], during the time interval
$[0,T]$, we can slowly change the system’s Hamiltonian from $H_{\rm{B}}$ to
$H_{\rm{C}}$ and obtain the maximum energy state of $H_{\rm{C}}$ with high
probability [26]. The changing process is exampled in (7).

$\displaystyle{H}\left(t\right)$ $\displaystyle=$
$\displaystyle\left[1-s\left(t\right)\right]{H}_{\rm{B}}+s\left(t\right){H}_{\rm{C}},$
(7)

where $s\left(t\right)$ is a smooth function, $s\left(0\right)=0$ and
$s\left(T\right)=1$. We then use Trotterization technique to emulate the
evolution process [44].

We discretize the total time interval $[0,T]$ into intervals $[j\Delta
t,(j+1)\Delta t]$ with small enough $\Delta t$. Over the $j^{\text{th}}$
interval, the Hamiltonian is approximately constant, i.e., $H(t)=H((j+1)\Delta
t)$. Therefore, the total time evolution operator $U(T,0)$ can be
approximately discretized into $2p$ implementable operators with constant
Hamiltonian [44], as written in (8). The approximation will improve as $\Delta
t$ gets smaller or, equivalently, as $p$ gets bigger.

$\displaystyle U\left(T,0\right)$ $\displaystyle=U(T,T-\Delta t)U(T-\Delta
t,T-2\Delta t)\cdots U(\Delta t,0)$ (8)
$\displaystyle=\prod_{j=0}^{p-1}U((j+1)\Delta t,j\Delta
t)\approx\prod_{j=1}^{p}e^{-i{H}\left({j}\Delta t\right)\Delta t},$

where $U((j+1)\Delta t,j\Delta t)$ represents the time evolution from $j\Delta
t$ to $(j+1)\Delta t$. Inserting (7) to (8) and using
$e^{i(A_{1}+A_{2})x}=e^{iA_{1}x}e^{iA_{2}x}+\mathcal{O}(x^{2})$, the time
evolution operator can be expressed as,

$\displaystyle U\left(T,0\right)$
$\displaystyle\approx\prod_{j=1}^{p}e^{-i\left[\left(1-s\left(j\Delta
t\right)\right){H}_{\rm{B}}+s\left(j\Delta t\right){H}_{\rm{C}}\right]\Delta
t}$ (9) $\displaystyle\approx\prod_{j=1}^{p}{e^{-i\left(1-s\left(j\Delta
t\right)\right){H}_{\rm{B}}\Delta t}e^{-is\left(j\Delta
t\right){H}_{\rm{C}}\Delta t}}+O(\Delta t^{2})$
$\displaystyle\approx\prod_{j=1}^{p}{U_{\rm{B}}^{\left(j\right)}U_{\rm{C}}^{\left(j\right)}},$

where $U_{\rm{B}}^{\left(j\right)}$ and $U_{\rm{C}}^{\left(j\right)}$ are the
time evolution operators, evolving the system under the Hamiltonian
$H_{\rm{B}}$ for the time period of $\beta_{j}=(1-s(j\Delta t))\Delta t$ and
the Hamiltonian $H_{\rm{C}}$ for the time period of $\gamma_{j}=s(j\Delta
t)\Delta t$, respectively, as defined in (10).

$\left\\{\begin{array}[]{lr}U_{\rm{C}}^{\left(j\right)}=e^{-is\left(j\Delta
t\right){H}_{\rm{C}}\Delta t}=e^{-i\gamma_{j}{H}_{\rm{C}}}&\\\
U_{\rm{B}}^{\left(j\right)}=e^{-i\left[1-s\left(j\Delta
t\right)\right]{H}_{\rm{B}}\Delta
t}=e^{-i{\beta_{j}}{H}_{\rm{B}}}&\end{array}\right.$ (10)

In the evolution, $\ket{\varphi}$ represents the quantum state $\ket{Z}$ in
section III. Through applying $U_{\rm{B}}^{\left(j\right)}$ and
$U_{\rm{C}}^{\left(j\right)}$ to the initial state
$\ket{\varphi(0)}={\ket{+}}^{\otimes n}$ alternately, we can compute the final
state $\ket{\varphi(T)}$ in (11), which is expected to lead a high
$C(\ket{\varphi(T)})$ and collapse to maximum energy state after measurement.

$\displaystyle{\ket{\varphi\left(T,{{\gamma}},{{\beta}}\right)}=\prod_{k=1}^{p}U_{\rm{B}}^{\left(j\right)}U_{\rm{C}}^{\left(j\right)}\ket{\varphi\left(0\right)}},$
(11)

where ${{\gamma}}=\left(\gamma_{1},\gamma_{2},\ldots,\gamma_{p}\right)$ and
${{\beta}}=\left(\beta_{1},\beta_{2},\ldots,\beta_{p}\right)$ need to be
optimized, which requires expensive computational effort.

###  IX-B Measurement Outcomes for the Cost Function Value

Quantum computers perform calculations based on the probability distribution
of quantum states. In QAOA, we obtain the cost function value in (4) by
sampling the quantum states, where $|\alpha_{k}|^{2}$ is the probability that
the final state $\ket{\varphi}$ collapses on the computational basis
$\ket{Z_{k}}$, as explained below.

First, we construct the quantum circuit of QAOA for the target graphs. In our
study, the circuit is built in the Qiskit simulator [45]. The quantum circuit
prepares the initial maximum energy state and computes the final state in (11)
by using the quantum operators $U_{\rm{B}}^{\left(j\right)}$ and
$U_{\rm{C}}^{\left(j\right)}$, whose implementations are illustrated in Fig. 2
and also shown in (12) and (13), respectively.

$\displaystyle
e^{-i{\beta_{k}}{H}_{\rm{B}}}=e^{-i\beta_{k}\sum_{j=1}^{n}\sigma_{j}^{x}}=\prod_{j=1}^{n}e^{-i\beta_{k}\sigma_{j}^{x}}=\prod_{j=1}^{n}R_{\rm{X}}^{(j)}(2\beta_{k}),$
(12) $\displaystyle e^{-i\gamma_{k}{H}_{\rm{C}}}=e^{-i\gamma_{k}\sum
w_{ij}\frac{I-\sigma_{i}^{z}\sigma_{j}^{z}}{2}}=\prod_{{\left\langle
i,j\right\rangle}\in E}R_{\rm{ZZ}}^{\left\langle
i,j\right\rangle}(-\gamma_{k}w_{ij}),$ (13)

where $R_{\rm{X}}^{(j)}$ means only applying $R_{\rm{X}}$ gate to the
$j^{\text{th}}$ qubit without changing other qubits; and
$R_{\rm{ZZ}}^{\left\langle i,j\right\rangle}$ means only applying
$R_{\rm{ZZ}}$ gate to the $i^{\text{th}}$ and $j^{\text{th}}$ qubits.

Second, we run the quantum circuit $N_{\rm{shot}}$ times and measure the final
state for the probability distribution. Suppose the final state collapses on
the $\ket{Z_{k}}$ with $N_{k}$ times, the approximation of $|\alpha_{k}|^{2}$
is $|\tilde{\alpha}_{k}|^{2}=N_{k}/N_{\rm{shot}}$. The approximation will
improve as $N_{\rm{shot}}$ gets bigger. In this study, to get an accurate
probability distribution and mean approximation-ratio, we perform
$N_{\rm{shot}}=2^{19}$ to approximate the distribution coefficients
$|\alpha_{k}|^{2}$. In practice, $2,048$ shots works well and is recommended.

Third, we calculate the cost function $C(\ket{\varphi})$ as shown in (4),
which is the weighted summation of $C(\ket{Z_{k}})$ with the non-zero
coefficients $|\tilde{\alpha}_{k}|^{2}$. Since the standard deviation of
$C(\ket{\varphi})$ is in the order of $\sqrt{m}$ [26], $N_{\rm{shot}}$ is in
the polynomial order. So, it is efficient to compute $C(\ket{Z_{k}})$ with
non-zero coefficients $|\tilde{\alpha}_{k}|^{2}$. Based on the calculation of
$C(\ket{Z_{k}})$ with non-zero coefficients $|\tilde{\alpha}_{k}|^{2}$, we
select the $\ket{Z_{\rm{opt}}}$ with the maximal cut value as the final
solution $Z_{\rm{opt}}$.

###  IX-C Optimization of Parameters

This subsection explains the parameter optimization involved in three
perspectives of Section VI and Section VII, where we need to optimize the
parameters for high cost function values in seed graphs with $n_{\rm{s}}=10$
and in target graphs with $n_{\rm{t}}=24$, when $p=1,2,3,10$. The three
perspectives are introduced below.

First, we get the optimal parameters for the seed graphs with $n_{\rm{s}}=10$
when $p=1,2,3$ by classical optimization method. In our study, the Newton
method is used to get the exact cost function values. Considering the non-
convex landscapes of the cost function, we adopt multiple initial guesses for
$(\gamma,\beta)$ within ${[0,2\pi]}^{p}\times{[0,\pi]}^{p}$. The number of
initial guesses is designed in the polynomial order of $n$ and $m$, which is
proved to be sufficient for obtaining the optimal parameters [26].

Second, we get the quasi-optimal parameters for the seed graphs with
$n_{\rm{s}}=10$ when $p=10$ by FOURIER heuristic strategy [29]. In the FOURIER
strategy, the time complexity of obtaining quasi-parameters is reduced into
$\mathcal{O}(\text{poly}(p))$ to avoid computational burden [29], thus the
parameters with high $p$ can be obtained efficiently.

Third, after getting the transferred parameters, we further optimize them for
verifying the efficacy of the transferred parameters. Specifically, we use the
COBYLA method [40] to further optimize the transferred parameters due to the
fluctuations of the cost function values. As mentioned in Methods Section B,
the quantum computer estimates the cost function values by sampling copies of
output quantum state, which results in the fluctuations of the cost function
values. The COBYLA method is used to address the optimization with this issue,
with the optimized results given in Fig. 5 and Extended Table Data-Driven
Quantum Approximate Optimization Algorithm for Cyber-Physical Power Systems.
The COBYLA method is also used to carry out the optimization of transferred
parameters for the test cases without the depolarizing noise. Note that some
certain gradient-based methods might also be able to find the quasi-optimal
parameters with fluctuating cost function value, where the inaccurate gradient
estimation may cause the escape from a local maximum and allow the converge
towards a better one [37].

###  IX-D Test Graphs Preparation

Without losing generality, the test graphs are randomly generated, as
introduced below. First, we randomly generate adjacency matrices. Second, we
set the entries to be zero with different probability to ensure the densities
of the test graphs spread over $[0,1]$. Third, since there exist a polynomial
algorithm for Max-Cut problem for the planar graph [46], we check the
generated graphs’ planar property by Kuratowski’s Theorem, such that all of
the test graphs are not planar. Finally, we generate $11,840$ graphs, sort
them by normalized graph densities, and then uniformly pick out $1710$ graphs
as test graphs.

###  IX-E Findings for Transfer Principles

Here we show two important findings in Fig. 3 to further explain the transfer
principles.

First, the mean approximation-ratio of QAOA for the target graphs are
correlated to a Lipschitz continuous curve with respect to their normalized
graph densities, although these target graphs are randomly generated. It is
justified by a scan window with the size $0.1$, as given in Fig. 3. The window
shows that the upper limit of the standard deviations of the scatters with the
same parameters is $0.057$, which further indicates the scatters approximately
follow a curve. It also indicates the normalized graph density is a effective
metric. According to these curves, we can directly estimate the mean
approximation-ratio of new graphs with parameters in the database. Thereafter
the parameters with outstanding performance can be quickly identified for the
QAOA circuit, avoiding the high computing effort.

Second, the parameters of seed graphs with low density perform better in
target graphs with low density than in the ones with high density, and vice
versa. This property is also uncovered in the mapping tables, where the yellow
area denoting the high approximation-ratio will increase as $D$ increases.
Specifically, when the sizes of the seed and target graphs are very close, the
yellow area will be around the diagonal line as shown in Extended Fig.
References; while, when the size of the target graph is much bigger than that
of the seed graphs, the area will be above the diagonal line as shown in
Extended Fig. References. With this property, the quasi-optimal parameters can
be effectively identified.

###  IX-F Modeling the Cyber-Physical Power System

We model the physical layer and then get multiple normalized weighted graphs
through the power flow calculation when disturbances from DERs are considered.
Power flow calculates the bus voltages for a given load, generation, and
network condition, based on which the line powers (weights) can be obtained.
The power flow equations are given in (14).

$\displaystyle
P_{i}-jQ_{i}=\dot{V}_{i}^{*}\sum_{k=1}^{n}{\dot{y}_{ik}\dot{V}_{k}},$ (14)

where $*$ denotes conjugate, $\dot{V}_{i}\in\mathbb{C}$ is the $i^{\text{th}}$
bus (vertex) voltage in the physical grid, $P_{i},Q_{i}\in\mathbb{R}$ is the
injection active and reactive power of the $i^{\text{th}}$ bus, and
$\dot{y}_{ik}\in\mathbb{C}$ is the admittance of the line between the
$i^{\text{th}}$ bus and $k^{\text{th}}$ bus.

After solving the power flow equations, we can obtain the complex power over
each line. In our study, the apparent power is used as the edge weight. Due to
the complex landscape of parameters [18], the edge weight is then normalized.
Thus, the modeling graph for the physical system is a normalized weighted
graph.

The modeling graph of the cyber layer is based on the communication network
data traffic that is flexible and random. In our study, we randomly generate
$n_{\rm{t}}=24$ graphs to represent the communication network. In practical
applications, we can monitor the data traffic to set up the edge weights for
the cyber graphs.

###  IX-G Depolarizing Noise Model

For demonstrating the potential of our method to be a promising candidate for
achieving quantum advantage on NISQ systems, we introduce the depolarizing
noise for the quantum gates in QAOA circuits. In our study, Qiskit [45] is
used to simulate the depolarizing noise and investigate the influences. The
simulator needs to calculate the density matrix after each quantum gate to
include the noise model, which costs exponentially more calculation resources
than the noiseless vector simulation.

###  IX-H The Approximation-ratio Distribution of GW Algorithm

In Fig. 7 (a)-(c), we obtain the approximation-ratio distributions of the GW
algorithm, as summarized below.

The GW algorithm relaxes the constraint of Max-Cut problem from discrete
variables to the vectors on a unit sphere. The relaxed problem then becomes a
semidefinite programming (SDP) problem. By solving the SDP problem, we obtain
the optimal vector distribution. By randomly cutting the unit sphere into two
parts, we correspondingly separate the vectors into two groups and obtain an
approximate solution. When we cut the sphere several times, there is a
guarantee that we have at least $0.878$ expected approximation ratio.

Similarly to the shots in QAOA, with certain cut times, the GW algorithm
outputs a probability distribution with respect to the approximation ratio,
e.g., the approximation-ratio distribution. In the theoretical research, we
usually compare the expectation of the approximation-ratio distribution to
evaluate the algorithm performance, while in the practical applications to
cyber-physical systems, we can take the maximum approximation-ratio as the
final approximation solution.

##  X Acknowledgment

Y.W. is primarily supported by the Office of the Director of National
Intelligence—Intelligence Advanced Research Projects Activity through ARO
Contract No. W911NF-16- 1-0082 and DOE BES award de-sc0019449 (quantum
algorithm analysis). The author would like to thank Dr. Rui Chao at Duke
University for the detailed discussion on the QAOA algorithm.

Recently, we became aware of a similar work by Shaydulin _et al._ about the
transferability between weighted graphs [47], which was carried out
independently.

## References

  * [1] Alán Aspuru-Guzik, Anthony D Dutoi, Peter J Love and Martin Head-Gordon  “Simulated quantum computation of molecular energies”  In _Science_ 309.5741 American Association for the Advancement of Science, 2005, pp. 1704–1707
  * [2] Matthew B Hastings, Dave Wecker, Bela Bauer and Matthias Troyer  “Improving quantum algorithms for quantum chemistry”  In _arXiv preprint arXiv:1403.1539_ , 2014
  * [3] Peter JJ O’Malley et al.  “Scalable quantum simulation of molecular energies”  In _Physical Review X_ 6.3 APS, 2016, pp. 031007
  * [4] Aram W Harrow, Avinatan Hassidim and Seth Lloyd  “Quantum algorithm for linear systems of equations”  In _Physical Review Letters_ 103.15 APS, 2009, pp. 150502
  * [5] Rozhin Eskandarpour, Kumar Ghosh, Amin Khodaei and Aleksi Paaso  “Experimental Quantum Computing to Solve Network DC Power Flow Problem”  In _arXiv preprint arXiv:2106.12032_ , 2021
  * [6] Yifan Zhou, Fei Feng and Peng Zhang  “Quantum electromagnetic transients program”  In _IEEE Transactions on Power Systems_ IEEE, 2021
  * [7] Rozhin Eskandarpour et al.  “Quantum computing for enhancing grid security”  In _IEEE Transactions on Power Systems_ 35.5 IEEE, 2020, pp. 4135–4137
  * [8] Fei Feng, Yifan Zhou and Peng Zhang  “Quantum Power Flow”  In _IEEE Transactions on Power Systems_ 36.4, 2021, pp. 3810–3812  DOI: 10.1109/TPWRS.2021.3077382
  * [9] Michael Lubasch et al.  “Variational quantum algorithms for nonlinear problems”  In _Phys. Rev. A_ 101 American Physical Society, 2020, pp. 010301  DOI: 10.1103/PhysRevA.101.010301
  * [10] Annarita Giani and Zachary Eldredge  “Quantum computing opportunities in renewable energy”  In _SN Computer Science_ 2.5 Springer, 2021, pp. 1–15
  * [11] Samantha Koretsky et al.  “Adapting Quantum Approximation Optimization Algorithm (QAOA) for Unit Commitment”  In _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , 2021, pp. 181–187  IEEE
  * [12] Filipe F.. Silva, Pedro M.. Carvalho, Luis A… Ferreira and Yasser Omar  “A QUBO Formulation for Minimum Loss Spanning Tree Reconfiguration Problems in Electric Power Networks”, 2021  arXiv:2109.09659 [quant-ph]
  * [13] John Preskill  “Quantum computing in the NISQ era and beyond”  In _Quantum_ 2 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2018, pp. 79
  * [14] Matthew P Harrigan et al.  “Quantum approximate optimization of non-planar graph problems on a planar superconducting processor”  In _Nature Physics_ 17.3 Nature Publishing Group, 2021, pp. 332–336
  * [15] Guido Pagano et al.  “Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator”  In _Proceedings of the National Academy of Sciences_ 117.41 National Acad Sciences, 2020, pp. 25396–25401
  * [16] TM Graham et al.  “Demonstration of multi-qubit entanglement and algorithms on a programmable neutral atom quantum computer”  In _arXiv preprint arXiv:2112.14589_ , 2021
  * [17] Farshid Shariatzadeh et al.  “Real-time implementation of intelligent reconfiguration algorithm for microgrid”  In _IEEE Transactions on sustainable energy_ 5.2 IEEE, 2014, pp. 598–607
  * [18] Hang Jing, Ye Wang and Yan Li  “Interoperation Analysis of Reconfigurable Networked Microgrids through Quantum Approximate Optimization Algorithm”  In _2022 IEEE Power & Energy Society General Meeting (PESGM)_, 2022, pp. 1–5  IEEE
  * [19] Mark EJ Newman  “Modularity and community structure in networks”  In _Proceedings of the national academy of sciences_ 103.23 National Acad Sciences, 2006, pp. 8577–8582
  * [20] M.R. Garey, D.S. Johnson and L. Stockmeyer  “Some simplified NP-complete graph problems”  In _Theoretical Computer Science_ 1.3, 1976, pp. 237–267  DOI: https://doi.org/10.1016/0304-3975(76)90059-1
  * [21] Michel X Goemans and David P Williamson  “Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming”  In _Journal of the ACM (JACM)_ 42.6 ACM New York, NY, USA, 1995, pp. 1115–1145
  * [22] Song Mei, Theodor Misiakiewicz, Andrea Montanari and Roberto Imbuzeiro Oliveira  “Solving SDPs for synchronization and MaxCut problems via the Grothendieck inequality”  In _Conference on Learning Theory_ , 2017, pp. 1476–1515  PMLR
  * [23] Sihong Shao, Dong Zhang and Weixi Zhang  “A simple iterative algorithm for maxcut”  In _arXiv preprint arXiv:1803.06496_ , 2018
  * [24] Weichi Yao, Afonso S Bandeira and Soledad Villar  “Experimental performance of graph neural networks on random instances of max-cut”  In _Wavelets and Sparsity XVIII_ 11138, 2019, pp. 111380S  International Society for OpticsPhotonics
  * [25] Subhash Khot, Guy Kindler, Elchanan Mossel and Ryan O’Donnell  “Optimal inapproximability results for MAX-CUT and other 2-variable CSPs?”  In _SIAM Journal on Computing_ 37.1 SIAM, 2007, pp. 319–357
  * [26] Edward Farhi, Jeffrey Goldstone and Sam Gutmann  “A quantum approximate optimization algorithm”  In _arXiv preprint arXiv:1411.4028_ , 2014
  * [27] Joao Basso et al.  “The Quantum Approximate Optimization Algorithm at High Depth for MaxCut on Large-Girth Regular Graphs and the Sherrington-Kirkpatrick Model”  In _arXiv preprint arXiv:2110.14206_ , 2021
  * [28] Edward Farhi, Jeffrey Goldstone, Sam Gutmann and Michael Sipser  “Quantum computation by adiabatic evolution”  In _arXiv preprint quant-ph/0001106_ , 2000
  * [29] Leo Zhou et al.  “Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices”  In _Phys. Rev. X_ 10 American Physical Society, 2020, pp. 021067  DOI: 10.1103/PhysRevX.10.021067
  * [30] Haibin Wang, Jiaojiao Zhao, Bosi Wang and Lian Tong  “A Quantum Approximate Optimization Algorithm with Metalearning for MaxCut Problem and Its Simulation via TensorFlow Quantum”  In _Mathematical Problems in Engineering_ 2021 Hindawi, 2021
  * [31] Jiahao Yao, Marin Bukov and Lin Lin  “Policy gradient based quantum approximate optimization algorithm”  In _Mathematical and Scientific Machine Learning_ , 2020, pp. 605–634  PMLR
  * [32] Takeshi Tokuyama  “Algorithms and computation”  Springer, 2007
  * [33] Fernando GSL Brandao et al.  “For fixed control parameters the quantum approximate optimization algorithm’s objective function value concentrates for typical instances”  In _arXiv preprint arXiv:1812.04170_ , 2018
  * [34] Alexey Galda et al.  “Transferability of optimal QAOA parameters between random graphs”  In _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , 2021, pp. 171–180  IEEE
  * [35] Max Born and Vladimir Fock  “Beweis des adiabatensatzes”  In _Zeitschrift für Physik_ 51.3 Springer, 1928, pp. 165–180
  * [36] Jonathan Wurtz and Peter Love  “MaxCut quantum approximate optimization algorithm performance guarantees for p¿ 1”  In _Physical Review A_ 103.4 APS, 2021, pp. 042612
  * [37] Gian Giacomo Guerreschi and Mikhail Smelyanskiy  “Practical optimization for hybrid quantum-classical algorithms”  In _arXiv preprint arXiv:1701.01450_ , 2017
  * [38] Boaz Barak and Kunal Marwaha  “Classical algorithms and quantum limitations for maximum cut on high-girth graphs”  In _arXiv preprint arXiv:2106.05900_ , 2021
  * [39] Marco Cerezo et al.  “Cost function dependent barren plateaus in shallow parametrized quantum circuits”  In _Nature communications_ 12.1 Nature Publishing Group, 2021, pp. 1–12
  * [40] Michael JD Powell  “A view of algorithms for optimization without derivatives”  In _Mathematics Today-Bulletin of the Institute of Mathematics and its Applications_ 43.5 Citeseer, 2007, pp. 170–174
  * [41] Samson Wang et al.  “Noise-induced barren plateaus in variational quantum algorithms”  In _Nature Communications_ 12.1 Nature Publishing Group, 2021, pp. 1–11
  * [42] Christos Ordoudis, Pierre Pinson, Juan M Morales and Marco Zugno  “An updated version of the IEEE RTS 24-bus system for electricity market and power system operation studies”  In _Technical University of Denmark_ 13, 2016
  * [43] Michael A Nielsen and Isaac L Chuang  “Quantum Computation and Quantum Information”  Cambridge University Press, 2011
  * [44] L-A Wu, MS Byrd and DA Lidar  “Polynomial-time simulation of pairing models on a quantum computer”  In _Physical Review Letters_ 89.5 APS, 2002, pp. 057904
  * [45] MD SAJID ANIS et al.  “Qiskit: An Open-source Framework for Quantum Computing”, 2021  DOI: 10.5281/zenodo.2573505
  * [46] Wei-Kuan Shih, Sun Wu and Yue-Sun Kuo  “Unifying maximum cut and minimum cut of a planar graph”  In _IEEE Transactions on Computers_ 39.5 IEEE, 1990, pp. 694–697
  * [47] Ruslan Shaydulin et al.  “Parameter Transfer for Quantum Approximate Optimization of Weighted MaxCut”  In _arXiv preprint arXiv:2201.11785_ , 2022

[]  The Parameters Performance of Seed Graphs with Different Layer Numbers $p$
and Vertex Size $n_{\rm{s}}$ seed graph $C(\ket{Z})$ when $n_{s}=10$
$C(\ket{Z})$ when $n_{s}=24$ $D$ $p=1$ $p=2$ $p=3$ $D$ $p=1$ $p=2$ $p=3$ 1
0.2667 0.7783 0.8727 0.9287 0.0831 0.7395 0.8117 0.8573 2 0.5333 0.8537 0.9193
0.9516 0.1901 0.7931 0.8526 0.8866 3 0.6444 0.8447 0.9009 0.9296 0.2869 0.8122
0.8681 0.8950 4 0.7333 0.8182 0.8601 0.8885 0.3527 0.8363 0.8872 0.9114 5
0.8000 0.8797 0.9206 0.9493 0.5163 0.8518 0.8974 0.9178 6 0.8667 0.8672 0.8943
0.9170 0.6256 0.8793 0.9168 0.9325 7 0.9111 0.9044 0.9223 0.9396 0.6756 0.8743
0.9105 0.9259 8 0.9556 0.9420 0.9553 0.9639 0.8316 0.9152 0.9261 0.9353 9
1.0000 0.9804 0.9977 0.9999 0.9608 0.9588 0.9588 0.9661

  * •

We randomly generate several seed graphs, including $9$ unweighted
$n_{\rm{s1}}=10$ graphs with ${D}$ spreading over $[0.2667,1]$, and $9$
weighted $n_{\rm{s2}}=24$ graphs with ${D}$ spreading over $[0.0831,0.9608]$.
In each seed graph, we use the classical optimizers mentioned in Methods
Section C to obtain the quasi-optimal parameters, which have better
approximation ratio as the layer number $p$ increases.

[!t]  The four scenarios of the test cyber-physical power system.

[!t]  The mapping table developed by the parameters from unweighted seed
graphs with $n_{\rm{s1}}=10$. Each sub-figure is a mapping table for
$D\mapsto[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ in the parameter transfer
module of the proposed data-driven QAOA in Fig. 2. This mapping
$D\mapsto[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ is performed based on the
entries in the mapping table. The entries are the scaling approximation
ratios, which can be obtained by applying the parameters from seed graphs to
the QAOA for the target graphs, as shown in Fig. 3. The procedure of getting
$D\mapsto[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ is that we first fix the
horizontal axis as shown in the orange arrow according to $D$ and then find
the interval $[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ as shown in purple line
according to the contour. After mapping, we identify the quasi-optimal
parameters from the seed graphs whose $\tilde{D}$ are within
$[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$.

[!t]  The mapping table developed by the parameters from weighted seed graphs
with $n_{\rm{s2}}=24$. Each sub-figure is a mapping table for
$D\mapsto[\tilde{D}_{\rm{l}},\tilde{D}_{\rm{r}}]$ in the parameter transfer
module of the proposed data-driven QAOA in Fig. 2. The developing procedure of
these mapping tables is the same with that of the mapping tables in Extended
Fig. References. Note that since we identify the parameters from
$n_{\rm{s2}}=24$ seed graphs and apply them to the same size $n_{\rm{t}}=24$
target graphs, the maximal performances are along the diagonal dash line.