# Nuclear Many-Body Effect on Particle Emissions Following Muon Capture on
${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$

Futoshi Minato<EMAIL_ADDRESS>Nuclear Data Center, Japan Atomic
Energy Agency, Tokai, Ibaraki 319-1195 Japan  RIKEN Nishina Center, Wako,
Saitama, 351-0198 Japan  Tomoya Naito  RIKEN Interdisciplinary Theoretical and
Mathematical Sciences (iTHEMS), Wako, Saitama, 351-0198 Japan  Department of
Physics, Graduate School of Science, The University of Tokyo, Tokyo, 113-0033
Japan  Osamu Iwamoto  Nuclear Data Center, Japan Atomic Energy Agency, Tokai,
Ibaraki 319-1195 Japan

###### Abstract

Muon captures on nuclei have provided us with plenty of knowledge of nuclear
properties. Recently, this reaction attracts attention in electronics, because
it is argued that charged particle emissions following muon capture on silicon
trigger non-negligible soft errors in memory devices. To investigate the
particle emissions from a nuclear physics point of view, we develop a new
approach using a microscopic model of muon capture and up-to-date particle
emission models. We paid attention to the muon capture rates, the particle
emission spectra, and the multiplicities that have a close interrelation with
each other, and found that the nuclear many-body correlation including two-
particle two-hole excitations is a key to explaining them simultaneously.

††preprint: RIKEN-iTHEMS-Report-22

Introduction—Negative muon capture on nuclei recently attracts high attention
due to its various feature. Muon is firstly captured into outer atomic
orbitals generated by the nuclear Coulomb potential and then transits to lower
orbitals emitting characteristic muonic X-rays and Auger electrons. These
X-rays are utilized to accurately determine nuclear charge radii [1, 2, 3] and
are also applied to non-destructive analyses inside vessels [4]. The muon
eventually settles in the lowest $s$ orbital and decays via
$\mu^{-}\rightarrow e^{-}+\nu_{\mu}+\bar{\nu}_{e}$ or is captured by a
nucleus. The latter case is analogue to electron capture of neutron-deficient
nuclei, but can occur even for stable nuclei due to the large muon mass. Thus,
the muon capture is utilized to analyze theoretically calculated nuclear
matrix elements of double $\beta$ decay [5].

Another importance of the muon capture is to bring nuclei that captured muon
to highly excited states. Those nuclei deexcite by emitting various particles,
and some of which ionize surrounding materials. Recently, great attention is
paid to this process in electronics field because non-negligible soft errors
are evoked in memory devices by charged particle and recoiled nucleus after
muon capture on silicon [6]. This issue becomes more serious with reducing the
scale of memory devices (the so-called die shrink) and operating them at low
voltage [7]. Currently, the soft errors in memory devices are studied with
Monte-Carlo transport simulations [8, 9] in which, however, the muon capture
process is greatly simplified by omitting to solve the nuclear many-body
problem. To date, there is no theoretical framework based on the nuclear
structure that describes the muon captures to the particle emissions
comprehensively, so developing a more practical model is highly demanded.

Particle emissions following the muon capture were extensively studied in the
70–80’s (see review article of Ref. [10]), and it was discussed that high
energy spectra are mainly attributed to instant particle emissions from the
preequilibrium state where excitation energy is shared only with a few
nucleons in nuclei. At this stage, it was also pointed out that two-body
meson-exchange current (MEC) becomes essential to explain high energy neutron
and proton emissions [11, 12]. In contrast, low energy spectra are attributed
to particle evaporation from a compound state where excitation energy is
shared with many nucleons in nuclei. Experimental data of particle emission
spectra were explained qualitatively by phenomenological models of
preequilibrium and compound states. However, they have not yet been reproduced
accurately enough to be used for a practical application Moreover, the muon
capture process was not discussed carefully from the nuclear structure point
of view.

Recent theoretical studies that consider the nuclear many-body systems more
appropriately pointed out that correlations resulting from the interaction
between nucleons play a significant role in the muon capture [13, 14, 15].
This fact motivated us to revisit the problem of particle emission following
muon capture. To understand the mechanism, we develop a new method that
considers the muon capture with a microscopic nuclear model and the particle
emission with an up-to-date model of preequilibrium and compound states. We
will demonstrate that many-body correlations resulting from the residual two-
body interactions are essential to describe particle emission spectra and
multiplicities, as well as muon capture rates. In particular, the effect of
two-particle two-hole ($2p$-$2h$) states, that is to say, the door-way state,
is significant for the particle emission spectra. The target nucleus of this
study is ${}^{28}\mathrm{Si}$, the main material of semiconductors. In
addition, we study ${}^{40}\mathrm{Ca}$ that have experimental data of
emission spectra from low to high energies.

Model—We assume that one-neutron particle one-proton hole ($1p_{\nu}1h_{\pi}$)
states are produced in nuclei by the muon capture at first. We describe this
process with the Tamm-Dancoff approximation (TDA) and second TDA (STDA) which
is the extension of TDA to $2p$-$2h$ model spaces [16, 17]. They are
extensively applied to study nuclear states at low to high excitation energies
including charge-exchange reactions. Then, we assume that the
$1p_{\nu}$-$1h_{\pi}$ state generated by the one-body weak interaction evolves
to more complicated multi-particle multi-hole ($mp$-$mh$) states leading to
the compound state. This process, the so-called preequilibrium state, is still
one of the challenging subjects in the nuclear physics. There are mainly two
approaches that have been investigated to describe the preequilibrium states,
which are semi-classical [18] and quantum models [19] (see also Ref. [20] for
the current status). The latter model has a good predictive power on several
experimental data without phenomenological parameters; however, it is still
difficult to describe various kinds of multiple particle emissions in low to
high energies comprehensively. In contrast, the semi-classical model shows a
good performance of reproducing various particle emissions, but with some
phenomenological ingredients of partial level densities, collision matrices,
and so on. We choose in this work the two-component exciton model [21, 22],
one of the semi-classical approaches, because various charged particle
emissions are the present scope and the targets in interest are stable nuclei
for which the model is well established through the nuclear data evaluations
[23, 24]. When nuclei reach the compound state, we describe the particle
emission with the Hauser-Feshbach statistical model [25] that is also applied
extensively to study nuclear reactions and particle emission after
$\beta$-decay [26, 27].

In STDA, excited states of daughter nuclei with spin $J$ are created by
operating a phonon creation operator of the vibrational states $Q_{\lambda
J}^{\dagger}$ to the target nuclear ground state $\left|0\right\rangle$ as

$\left|\lambda J\right\rangle=Q_{\lambda
J}^{\dagger}\left|0\right\rangle,\qquad Q_{\lambda J}\left|0\right\rangle=0.$
(1)

We calculate the ground state $\left|0\right\rangle$ with the Skyrme-Hartree-
Fock (SHF) method [28] in coordinate space assuming spherical symmetry. Note
that ${}^{28}\mathrm{Si}$ is considered to form a oblate shape, and the
deformation effect and the pairing correlation need to be taken into account
practically. However, it has been studied that those influences are not
striking because of the conservation of closed shell structure for oblate
${}^{28}\mathrm{Si}$ [29]. We studied with using two effective forces, which
are SGII [30] and SkO’ [31]. Those forces provide a reasonable strength
distribution of time-odd channels and the order of low-lying states of
${}^{28}\mathrm{Al}$ and ${}^{40}\mathrm{K}$. Continuum states are discretized
by a box size $R=14\,\mathrm{fm}$. The STDA phonon operator of the vibrational
states is given by

$Q_{\lambda J}^{\dagger}=\sum_{mi}X_{mi}^{\lambda
J}\mathcal{O}_{mi}^{J\dagger}+\sum_{m\leq n,i\leq
j}\mathcal{X}_{mnij}^{\lambda J}\mathcal{O}_{mnij}^{J\dagger},$ (2)

where $m$ and $n$ denote particle states, while $i$ and $j$ denote hole
states. The operators $\mathcal{O}_{mi}^{JM\dagger}$ and
$\mathcal{O}_{mnij}^{JM\dagger}$ create $1p$-$1h$ and $2p$-$2h$ states,
respectively. Omitting the second term of Eq. (2) corresponds to the TDA. The
coefficients $X$ and $\mathcal{X}$ in Eq. (2) that effectively reflect the
effects of the residual two-body interaction are obtained by solving the TDA
and STDA equation [16, 17]. The model space for $1p$-$1h$ state is set to be
$\varepsilon_{m}-\varepsilon_{i}\leq 100\,\mathrm{MeV}$, where $\varepsilon$
is the single-particle energy, and that for $2p$-$2h$ state is restricted to
$2n_{q}+l_{q}\leq 4$ except the $0g$ states, where $n_{q}$ and $l_{q}$
($q=\pi$,$\nu$) are the principle number and orbital angular momentum of a
single-particle state.

Muon capture rates are written as

$\displaystyle\omega\left(E\right)=\sum_{\lambda,\,J}\frac{2G^{2}\nu^{2}}{1+\nu/M}$
$\displaystyle\times\left[\left|\sum_{\nu\pi}X_{\nu\pi}^{\lambda
J}\left\langle
j_{\nu}l_{\nu}\middle\|\phi\left(\hat{\mathcal{M}}_{J}-\hat{\mathcal{L}}_{J}\right)\middle\|j_{\pi}l_{\pi}\right\rangle\right|^{2}\right.$
$\displaystyle\left.+\left|\sum_{\nu\pi}X_{\nu\pi}^{\lambda J}\left\langle
j_{\nu}l_{\nu}\middle\|\phi\left(\hat{\mathcal{T}}_{J}^{\text{el}}-\hat{\mathcal{T}}_{J}^{\text{mag}}\right)\middle\|j_{\pi}l_{\pi}\right\rangle\right|^{2}\right]\delta\left(E-E_{\text{x}}^{\text{($\lambda$,
$J$)}}\right)$ (3)

where $G=1.166\times 10^{-11}\,\mathrm{MeV}^{-2}$ is the Fermi coupling
constant, $M$ is the mass of the target nucleus, $\nu=m_{\mu}-\Delta
M+\varepsilon_{\mu}-E_{\text{x}}^{\text{($\lambda$, $J$})}$ is the muon
neutrino energy, $\Delta M$ is the mass difference between parent and
descendant nuclei taken from AME2020 [32], $j_{\nu,\,\pi}$ are the total
angular momentum, and $\varepsilon_{\mu}$ is the muon binding energy. The
excitation energy $E_{\text{x}}^{\text{($\lambda$,$J$)}}$ with respect to the
ground states of the daughter nuclei (${}^{28}\mathrm{Al}$ and
${}^{40}\mathrm{K}$) resulting from the muon capture on the parent nuclei
(${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$), respectively, are
approximated by $E_{\text{x}}^{\text{($\lambda$,
$J$)}}=E_{\text{TDA}}^{\text{($\lambda$,
$J$)}}-E_{\text{TDA}}^{\text{($0$)}}$, where
$E_{\text{TDA}}^{\text{($\lambda$, $J$)}}$ and $E_{\text{TDA}}^{\text{($0$)}}$
are the STDA or TDA phonon energies and the lowest energies, respectively. We
consider spin-parity up to $J^{\pi}\leq 5^{\pm}$. The lowest states of TDA and
STDA are, respectively, $3^{+}$ and $4^{-}$ for ${}^{28}\mathrm{Al}$ and
${}^{40}\mathrm{K}$, which are consistent with the experimental data. The one-
body operators of charge $\hat{\mathcal{M}}_{J}$, longitudinal
$\hat{\mathcal{L}}_{J}$, transverse electric
$\hat{\mathcal{T}}^{\text{el}}_{J}$, and transverse magnetic
$\hat{\mathcal{T}}^{\text{mag}}_{J}$ fields are found in Ref. [33]. For the
axial-vector coupling, the free nucleon value is quenched to some extent in
nuclei; however, the exact value is still unknown exactly. Hence, we study
with $g_{A}=-1.26$ (free nucleon) and the quenched value of $g_{A}=-1$. This
quenching is partly explained by the coupling with $2p$-$2h$ states, and we
will discuss this point later. For numerical purpose, $\delta$-function in Eq.
(3) is replaced by the Lorentzian function with a width of $1\,\mathrm{MeV}$.

The muon binding energy of $1s_{1/2}$ orbital, $\varepsilon_{\mu}$, and the
wave function, $\phi\equiv\phi\left(r\right)$, are calculated by solving the
Dirac equation under the Coulomb potential formed by the atomic nucleus,
$V_{\text{$N$-$\mu$}}^{\text{Coul}}$, and the $\left(Z-1\right)$ electrons,
$V_{\text{$e$-$\mu$}}^{\text{Coul}}$, and the vacuum polarization of the
Coulomb potential formed by the atomic nucleus, $V_{\text{$e$-$\mu$}}$. Here,
$Z$ denotes the atomic number of the atom. In this work, the electron density
distribution is calculated by the density functional theory [34, 35] with the
local density approximation (LDA), where the PZ81 LDA correlation functional
[36] is used. It should be noted that the electron density and the muon wave
function are solved simultaneously and self-consistently; thus, the effect of
the muon is also considered in the electron density distribution as well. The
vacuum polarization between the nucleus and the muon
$V_{\text{$N$-$\mu$}}^{\text{VP}}$ is considered by using the Uehling
effective potential [37, 38]. Consequently, the muon wave function reads

$\left[T+V_{\text{{\hbox{ N }}-{\hbox{
\mu}}}}^{\text{Coul}}\left(\bm{r}\right)+V_{\text{{\hbox{ N }}-{\hbox{
\mu}}}}^{\text{VP}}\left(\bm{r}\right)+V_{\text{{\hbox{ e }}-{\hbox{
\mu}}}}^{\text{Coul}}\left(\bm{r}\right)\right]\phi\left(\bm{r}\right)=\varepsilon\phi\left(\bm{r}\right),$
(4)

where $T$ is the Dirac kinetic operator. For more detail, see Supplemental
Material.

Assuming that the $1p_{\nu}$-$1h_{\pi}$ state generated by the one-body
operators of the muon capture as an initial state, we carry out the two-
component exciton model calculation. The master equation of the two-component
exciton model is [22]

$\displaystyle\frac{dP\left(p_{\pi},p_{\nu},t\right)}{dt}$
$\displaystyle=P\left(p_{\pi}-1,p_{\nu},t\right)\lambda_{\pi+}\left(p_{\pi}-1,p_{\nu}\right)$
$\displaystyle+P\left(p_{\pi},p_{\nu}-1,t\right)\lambda_{\nu+}\left(p_{\pi},p_{\nu}-1\right)$
$\displaystyle+P\left(p_{\pi}-1,p_{\nu}+1,t\right)\lambda_{\nu\pi}\left(p_{\pi}-1,p_{\nu}+1\right)$
$\displaystyle+P\left(p_{\pi}+1,p_{\nu}-1,t\right)\lambda_{\pi\nu}\left(p_{\pi}+1,p_{\nu}-1\right)$
$\displaystyle-P\left(p_{\pi},p_{\nu},t\right)\left[\lambda\left(p_{\pi},p_{\nu}\right)+W\left(p_{\pi},p_{\nu}\right)\right],$
(5)

where
$\lambda\left(p_{\pi},p_{\nu}\right)=\sum_{x}\lambda_{x}\left(p_{\pi},p_{\nu}\right)$
is the total transition rate, $\lambda_{q+}$ is the creation rate of a
particle-hole pair, $\lambda_{\pi\nu}$ is the exchange rate of proton and
neutron particle-hole pairs, and $P\left(p_{\pi},p_{\nu},t\right)$ is the
occupation probability of the exciton state having proton particle number
$p_{\pi}$ and neutron particle number $n_{\nu}$ at time $t$.
$W\left(p_{\pi},p_{\nu}\right)=\sum_{b}W_{b}\left(p_{\pi},p_{\nu}\right)$ is
the total particle emission rate, and $W_{b}$ is the particle emission rate
for the particle $b$ calculated with the inverse reaction cross section and
the partial level density for residual nuclei. Here, notation of hole state is
omitted for simplicity. By solving Eq. (5) with initial condition of
$1p_{\nu}$-$1h_{\pi}$, the probability for emitting particle $b$ can be
calculated by

$\sigma\left(E\right)=R\left(E\right)\sum_{p_{\pi},p_{\nu}}Q\left(p_{\pi},p_{\nu}\right)W_{b}\left(p_{\pi},p_{\nu}\right)$
(6)

where $Q$ is the cumulative occupation probability defined as
$Q\left(p_{\pi},p_{\nu}\right)=\int_{t_{0}}^{t_{1}}P\left(p_{\pi},p_{\nu}\right)\,dt$,
$R$ the normalized capture rate as
$R\left(E\right)=\omega\left(E\right)/\int\omega\left(E^{\prime}\right)\,dE^{\prime}$
and $W_{b,p}\left(p_{\pi},p_{\nu}\right)$ is the particle emission rate from
the exciton state $\left(p_{\pi},p_{\nu}\right)$. As the number of exciton is
greater than $12$, the calculation of preequilibrium state terminates and that
of compound state initiates. We use a standard parameter set for
preequilibrium and compound states that is globally used in evaluating nuclear
data (see Ref. [22] for more detail). Only for the proton and neutron single-
particle state densities $g$, we adjust to $g_{\pi}=Z/19$ and $g_{\nu}=N/19$
from the standard value of the Fermi gas, $g_{\pi}=Z/15$ and $g_{\nu}=N/15$.
This adjustment is reasonable considering the semi-magic structure of
${}^{28}\mathrm{Si}$ and the magic structure of ${}^{40}\mathrm{Ca}$ (see
Supplemental Material).

Result and discussion—Table 1 lists the calculated muon capture rates of
${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$. In addition to TDA and STDA, we
also show the result of “FREE” that is obtained by assuming that nucleons move
independently in the nuclear potential. The experimental data for natural
silicon and calcium are also listed. Note that the natural abundances of 28Si
and 40Ca are about $92$% and $97$%, respectively, and their muon capture rates
are expected to be close to those for the natural elements. The range of
calculated muon capture rates is estimated with the axial-vector coupling of
$g_{A}=-1$ and $-1.26$, which correspond to the lower and upper values of the
calculated muon capture rates, respectively. We find that the muon capture
rates of FREE overestimate experimental data of ${}^{28}\mathrm{Si}$ and
${}^{40}\mathrm{Ca}$ both for SGII and SkO’, while TDA and STDA reproduces
reasonably well within the uncertainties of $g_{A}$. One may notice that the
lower value of the capture rate, that is the result of $g_{A}=-1.0$, is
favorable for TDA, while the result of weakly quenched axial-vector coupling
is supportive for STDA. This result is consistent to a picture that the
quenching of the axial-vector coupling is partly explained by the coupling
with $2p$-$2h$ states [39].

Table 1: Calculated muon capture rates for ${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$ (in the unit of $10^{6}\,\mathrm{s}^{-1}$). The lower and upper values of model calculations are obtained by using $g_{A}=-1$ and $-1.26$, respectively. The experimental data for natural silicon and calcium are also listed [40]. The numbers in parentheses are the uncertainties of the corresponding last digits. Nucl. | Force | FREE | TDA | STDA | Expt.
---|---|---|---|---|---
${}^{28}\mathrm{Si}$ | SGII |  $1.02$–$1.46$ |  $0.87$–$1.26$ |  $0.81$–$1.18$ | $0.8712\left(18\right)$
| SkO’ |  $1.04$–$1.49$ |  $0.87$–$1.26$ |  $0.72$–$1.04$ |
${}^{40}\mathrm{Ca}$ | SGII |  $3.12$–$4.35$ |  $2.54$–$3.58$ |  $2.44$–$3.32$ | $2.557\left(14\right)$
| SkO’ |  $2.58$–$3.63$ |  $2.07$–$2.95$ |  $1.90$–$2.58$ |

The agreement with experimental data is related to the many-body correlations
resulting from the residual two-body interaction. To explain it, we show in
Fig. 1 the normalized capture rates $R\left(E\right)$ of ${}^{28}\mathrm{Al}$
and ${}^{40}\mathrm{K}$ calculated by SkO’ (see Supplemental Material for that
by SGII). The functions of $R\left(E\right)$ for FREE, TDA, and STDA are
similar to each other in excitation energies less than $20$ MeV, where are the
major part of the capture rates having more than $1\,\%$. However, those of
TDA and STDA distribute in higher excitation energies than FREE by a few MeV,
and moreover, significant enhancements above $E=25\,\mathrm{MeV}$ are
observed. This is because the residual two-body interaction works repulsively
for most of the $J^{\pi}$ channels. Such a feature substantially decreases the
momentum of outgoing muon neutrino, reducing the muon capture rate through the
factor $\nu^{2}$ in Eq. (3).

STDA gives additional enhancements above $E=30\,\mathrm{MeV}$ as compared with
TDA as seen in Fig. 1. This is because the $1p$-$1h$ states couple with the
$2p$-$2h$ ones, and some of them at high energies receive substantial
strengths from those at low energies. The enhancement of transition
probabilities at high energies induces the further reduction of muon capture
rates as found in Table 1, and the weak quenching of the axial-vector coupling
shows a favorable agreement with the experimental data accordingly. Although
the capture rates above $30\,\mathrm{MeV}$ are about $1\,\%$ at most and the
enhancement due to the $2p$-$2h$ states is a little, they affect particle
emission spectra appreciably as discussed next. In general, the effect of the
residual interaction becomes small with increasing energy, and the capture
rates of FREE, TDA, and STDA get closer to each other.

Figure 1: Normalized capture rates $R\left(E\right)$ for ${}^{28}\mathrm{Si}$
(top) and ${}^{40}\mathrm{Ca}$ (bottom) obtained by FREE, TDA, and STDA. The
results of SkO’ force with $g_{A}=-1.0$ are shown. The horizontal axis
represents the excitation energy of ${}^{28}\mathrm{Al}$ and
${}^{40}\mathrm{K}$. Contributions of MEC approximated by a Gaussian form are
also shown by the dashed line (see text).

Figure 2 shows the particle emission spectra of muon capture on
${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$. The results of SkO’ with the
axial-vector coupling $g_{A}=-1$ are illustrated together with the
experimental data. The peaks formed around $E=2$–$4\,\mathrm{MeV}$ are
resulting from the particle emission from the compound state, while high
energy tails of spectra are due to the particle emission from the
preequilibrium state. The result of FREE underestimates largely the
experimental data. This shortcoming is improved by TDA because of the
enhancement of the capture rates at high energies as seen in Fig. 1. STDA
further raises the calculated spectra and the results get more closer to
experimental data. We would like to stress here that the effect of the
coupling of $2p$-$2h$ states was only less than $1\,\%$ for muon capture rates
as seen in Fig. 1; however, its influence on particle emission spectra (Fig.
2) is non-negligible.

Figure 2: Particle yields after the muon capture on ${}^{28}\mathrm{Si}$ (top) and ${}^{40}\mathrm{Ca}$ (bottom). The result of SkO’ for $g_{A}=-1$ is shown. Experimental data of neutron (filled symbols) and proton (open symbols) are taken from Refs. [41, 42, 43] and [44, 45, 46], respectively. Note that the unit is not given in the original paper of Budyashov [44], so that we normalized the second point from the low energy to STDA+MEC. Table 2: Calculated multiplicities per $10^{3}$ muon captures for ${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$, compared with the experimental data for the natural silicon and calcium. The results of SGII and SkO’ forces with $g_{A}=-1$ are shown. Experimental data of neutron are taken from Ref. [47], while those of charged particles are from Ref. [45]. The numbers in parentheses are the uncertainties of the corresponding last digits. Nucl. | Particle | Energy | FREE | TDA | STDA | STDA+MEC | Expt.
---|---|---|---|---|---|---|---
|  | Range ($\mathrm{MeV}$) | SGII | SkO’ | SGII | SkO’ | SGII | SkO’ | SGII | SkO’ |
${}^{28}\mathrm{Si}$ | $n$ | Entire | $639$ | $573$ | $631$ | $812$ | $676$ | $938$ | $729$ | $983$ | $864\left(72\right)$
| $p$ |  $4$–$20$ | $14.9$ | $13.5$ | $23.9$ | $23.7$ | $31.8$ | $41.4$ | $54.4$ | $63.5$ | $52.46\left(192\right)$
| $d$ |  $5$–$17$ | $1.99$ | $2.01$ | $3.98$ | $4.13$ | $5.89$ | $7.13$ | $7.96$ | $9.15$ | $9.80\left(46\right)$
| $t$ |  $6$–$17$ | $0.303$ | $0.309$ | $0.714$ | $0.720$ | $1.18$ | $1.54$ | $1.63$ | $1.97$ | $1.70\left(13\right)$
| $\alpha$ |  $15$–$20$ | $0.288$ | $0.285$ | $0.598$ | $0.607$ | $0.944$ | $1.14$ | $1.28$ | $1.47$ | $0.57\left(10\right)$
${}^{40}\mathrm{Ca}$ | $n$ | Entire | $352$ | $473$ | $583$ | $659$ | $559$ | $631$ | $613$ | $681$ | $764\left(32\right)$

In Fig. 2, STDA still underestimates the proton spectra largely. The one-body
operators of muon capture in Eq. (3) create only $1p_{\nu}$-$1h_{\pi}$ state.
To enhance the proton emission spectra from this initial state, neutron must
give its energy instantly to proton during the preequilibrium state; however,
it is difficult to simulate this energy transfer enough to reproduce
experimental data. In this respect, Lifshitz and Singer have discussed that
MEC plays an important role to explain the experimental proton spectra [12].
We also consider the effect phenomenologically in the present framework as
follows. The capture rate due to MEC is almost universal for nuclei according
to Lifshitz and Singer, and we approximate it by a Gaussian function, where
the mean excitation energy and the width of which are given as $56$ and
$15\,\mathrm{MeV}$, respectively. The contribution of MEC to the capture rate
is set to be $5\,\%$ as estimated by Lifshitz and Singer [12], while that of
STDA is reduced by $95\,\%$. The capture rate of MEC is illustrated in Fig. 1.
Since MEC is two-body current and has a form of two-body isospin operator
$\left[\tau_{1}\otimes\tau_{2}\right]^{1,-1}$, we assume that two
configurations that are $1p_{\pi}$-$2h_{\pi}$-$1p_{\nu}$-$0h_{\nu}$ and
$0p_{\pi}$-$1h_{\pi}$-$2p_{\nu}$-$1h_{\nu}$ states are created with the equal
probability after the muon capture, and set them as the initial preequilibrium
state. The result of considering MEC (STDA+MEC) is shown in Fig. 2. Particle
emission spectra at high energies are enhanced further, in particular, a
remarkable improvement is obtained for the proton spectra.

Table 2 lists the calculated and experimental data of multiplicities of
emitted particles for ${}^{28}\mathrm{Si}$ and ${}^{40}\mathrm{Ca}$. Here, we
used $g_{A}=-1$ in the calculation; however, we confirmed that the difference
from $g_{A}=-1.26$ is less than $3\,\%$. The results of FREE underestimate the
experimental data, while we obtain improvements with increasing the many-body
correlations of $1p$-$1h$ mixture (TDA), and the coupling with $2p$-$2h$
states (STDA). The effect of MEC is comparable to TDA and STDA, making the
calculated results even closer to the experimental data. On the contrary, the
calculated $\alpha$ particle multiplicity of ${}^{28}\mathrm{Si}$ and neutron
multiplicity of ${}^{40}\mathrm{Ca}$ deviate from the experimental data. The
measured energy range of $\alpha$ particle multiplicity for
${}^{28}\mathrm{Si}$ is limited only to $15$–$20\,\mathrm{MeV}$ and the
available experimental data of neutron multiplicity for ${}^{40}\mathrm{Ca}$
is only one, which was measured more than $50$ years ago [47]. The
multiplicity strongly depends on the distribution of capture rates that
modulates particle emissions from the preequilibrium state which emits only a
few particles and the compound state which emits multiple particles. For
further understanding of particle emissions and validations of the nuclear
model, more experimental studies that cover the spectra from low to high
energies are required. The calculated multiplicities of charged particles for
entire energy range are provided in Supplemental Material.

Summary—We demonstrated that the many-body correlation resulting from the
residual two-body interaction operating between nucleons played an important
role in the muon capture rates of ${}^{28}\mathrm{Si}$ and
${}^{40}\mathrm{Ca}$. In particular, the coupling with $2p$-$2h$ states was
essential to describe the particle emission spectra and the multiplicities. In
other words, this indicates that particle emissions following muon captures
have information on the nuclear structure at high energies. The present study
showed that the combination of the microscopic approach of muon capture and
the two-component exciton model of particle emission is an effective tool to
describe particle emission following the muon captures, giving the nuclear
structure information additionally. Only the effect of MEC was considered
within a simple manner so that it is demanded to take into account in STDA in
future. The present outcomes are expected to contribute to the development of
Monte-Carlo transport simulations of muon captures [8, 9] and further
understanding of the nuclear structure.

The $\alpha$ particle multiplicity of 28Si and neutron multiplicity of
${}^{40}\mathrm{Ca}$ are the remaining questions of this work. Currently,
three facilities providing negative muon beams (TRIUMF, RAL, and J-PARC) are
running in the world, and new experiments are planned there. In addition, a
new negative muon facility is going to launch in RCNP at Osaka Univ. Those
activities will increase information on the interaction between muon and
nuclei, and help us not only to verify the present framework but also to
develop theoretical models.

###### Acknowledgements.

The authors thanks Prof. S. Kawase, Dr. A. Abe, and Prof. G. Colò for fruitful
discussion. FM thanks Dr. T. Fukui for supporting this work. TN thanks Dr. S.
A. Sato for discussion on the code for electron wave functions in the density
functional theory. TN acknowledges the financial support from the RIKEN
Special Postdoctoral Researchers Program. The numerical calculations were
partly performed on cluster computers at the RIKEN iTHEMS program.

## References

  * Fricke  _et al._ [1995] G. Fricke, C. Bernhardt, K. Heilig, L. A. Schaller, L. Schellenberg, E. B. Shera, and C. W. De Jager, At. Data Nucl. Data Tables 60, 177 (1995).
  * Antognini  _et al._ [2013] A. Antognini _et al._ , Science 339, 417 (2013).
  * Saito  _et al._ [2022] T. Y. Saito _et al._ , “Muonic X-Ray Measurement for the Nuclear Charge Distribution: the Case of Stable Palladium Isotopes,” (2022), arXiv:2204.03233 [nucl-ex] .
  * Shimada-Takaura  _et al._ [2021] K. Shimada-Takaura, K. Ninomiya, A. Sato, N. Ueda, M. Tampo, S. Takeshita, I. Umegaki, Y. Miyake, and K. Takahashi, J. Nat. Med. 75, 532 (2021).
  * Hashim and Ejiri [2021] I. H. Hashim and H. Ejiri, Front. Astron. Space Sci. 8, 82 (2021).
  * Serre  _et al._ [2012] S. Serre, S. Semikh, S. Uznanski, J. L. Autran, D. Munteanu, G. Gasiot, and P. Roche, IEEE Trans. Nucl. Sci. 59, 714 (2012).
  * Liao  _et al._ [2019] W. Liao, M. Hashimoto, S. Manabe, S.-i. Abe, and Y. Watanabe, IEEE Trans. Nucl. Sci. 66, 1390 (2019).
  * Abe and Sato [2017] S. Abe and T. Sato, J. Nucl. Sci. Tech. 54, 101 (2017).
  * Kossov [2007] M. V. Kossov, Eur. Phys. J. A 33, 7 (2007).
  * Measday [2001] D. F. Measday, Phys. Rep. 354, 243 (2001).
  * Dautry  _et al._ [1976] F. Dautry, M. Rho, and D. O. Riska, Nucl. Phys. A 264, 507 (1976).
  * Lifshitz and Singer [1988] M. Lifshitz and P. Singer, Nucl. Phys. A 476, 684 (1988).
  * Lovato  _et al._ [2019] A. Lovato, N. Rocco, and R. Schiavilla, Phys. Rev. C 100, 035502 (2019).
  * Jokiniemi  _et al._ [2019] L. Jokiniemi, J. Suhonen, H. Ejiri, and I. H. Hashim, Phys. Lett. B 794, 143 (2019).
  * Ciccarelli  _et al._ [2020] M. Ciccarelli, F. Minato, and T. Naito, Phys. Rev. C 102, 034306 (2020).
  * Da Providência [1965] J. Da Providência, Nucl. Phys. 61, 87 (1965).
  * Minato [2016] F. Minato, Phys. Rev. C 93, 044319 (2016).
  * Griffin [1966] J. J. Griffin, Phys. Rev. Lett. 17, 478 (1966).
  * Feshbach  _et al._ [1980] H. Feshbach, A. Kerman, and S. Koonin, Ann. Phys. 125, 429 (1980).
  * Carlson  _et al._ [2014] B. V. Carlson, J. E. Escher, and M. S. Hussein, J. Phys. G 41, 094003 (2014).
  * Koning and Duijvestijn [2004] A. J. Koning and M. C. Duijvestijn, Nucl. Phys. A 744, 15 (2004).
  * Iwamoto  _et al._ [2016] O. Iwamoto, N. Iwamoto, S. Kunieda, F. Minato, and K. Shibata, Nucl. Data Sheets 131, 259 (2016).
  * Hetrick  _et al._ [1997] D. M. Hetrick, D. C. Larson, N. M. Larson, L. C. Leal, and S. J. Epperson,  _Evaluation of 28,29,30Si neutron induced cross sections for ENDF/B-VI_, Tech. Rep. (Oak Lidge National Laboratory, 1997) TM-11825.
  * Shibata and Kunieda [2008] K. Shibata and S. Kunieda, J. Nucl. Sci. Tech. 45, 123 (2008).
  * Hauser and Feshbach [1952] W. Hauser and H. Feshbach, Phys. Rev. 87, 366 (1952).
  * Mumpower  _et al._ [2018] M. R. Mumpower, T. Kawano, T. M. Sprouse, N. Vassh, E. M. Holmbeck, R. Surman, and P. Möller, Astrophys. J. 869, 14 (2018).
  * Minato  _et al._ [2021] F. Minato, T. Marketin, and N. Paar, Phys. Rev. C 104, 044321 (2021).
  * Vautherin and Brink [1972] D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 (1972).
  * Ha and Cheoun [2017] E. Ha and M.-K. Cheoun, Eur. Phys. J. A 53, 26 (2017).
  * Van Giai and Sagawa [1981] N. Van Giai and H. Sagawa, Phys. Lett. 106B, 379 (1981).
  * Reinhard  _et al._ [1999] P.-G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys. Rev. C 60, 014316 (1999).
  * Wang  _et al._ [2021] M. Wang, W. J. Huang, F. G. Kondev, G. Audi, and S. Naimi, Chin. Phys. C 45, 030003 (2021).
  * O’Connell  _et al._ [1972] J. S. O’Connell, T. W. Donnelly, and J. D. Walecka, Phys. Rev. C 6, 719 (1972).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Uehling [1935] E. A. Uehling, Phys. Rev. 48, 55 (1935).
  * Wayne Fullerton and Rinker [1976] L. Wayne Fullerton and G. A. Rinker, Phys. Rev. A 13, 1283 (1976).
  * Ichimura  _et al._ [2006] M. Ichimura, H. Sakai, and T. Wakasa, Prog. Part. Nucl. Phys. 56, 446 (2006).
  * Suzuki  _et al._ [1987] T. Suzuki, D. F. Measday, and J. P. Roalsvig, Phys. Rev. C 35, 2212 (1987).
  * Sundelin and Edelstein [1973] R. M. Sundelin and R. M. Edelstein, Phys. Rev. C 7, 1037 (1973).
  * Kozlowski  _et al._ [1985] T. Kozlowski, W. Bertl, H. P. Povel, U. Sennhauser, H. K. Walter, A. Zglinski, R. Engfer, C. Grab, E. A. Hermes, H. P. Isaak, A. Van Der Schaaf, J. Van Der Pluym, and W. H. A. Hesselink, Nucl. Phys. A 436, 717 (1985).
  * Van Der Pluym  _et al._ [1986] J. Van Der Pluym, T. Kozlowski, W. H. A. Hesselink, A. Van Der Schaaf, C. Grab, E. A. Hermes, and W. Bertl, Phys. Lett. B 177, 21 (1986).
  * Budyashov  _et al._ [1971] Y. G. Budyashov, V. G. Zinov, A. D. Konin, A. I. Mukhin, and A. M. Chatrchyan, Sov. Phys. JETP 33, 11 (1971).
  * Edmonds  _et al._ [2022] A. Edmonds _et al._ (AlCap Collaboration), Phys. Rev. C 105, 035501 (2022).
  * Manabe [2020] S. Manabe,  _Study of Muon-induced Single Event Upsets in SRAMs and their Fundamental Physical Process_ , Ph.D. thesis, Kyushu Univ. (2020).
  * Macdonald  _et al._ [1965] B. Macdonald, J. A. Diaz, S. N. Kaplan, and R. V. Pyle, Phys. Rev. 139, B1253 (1965).


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Efficient Classical Computation of Quantum Mean Values for Shallow QAOA
Circuits

Wei-Feng Zhuang  Beijing Academy of Quantum Information Sciences, Beijing
100193, China  Ya-Nan Pu  Beijing Academy of Quantum Information Sciences,
Beijing 100193, China  Hong-Ze Xu  Beijing Academy of Quantum Information
Sciences, Beijing 100193, China  Xudan Chai  Beijing Academy of Quantum
Information Sciences, Beijing 100193, China  Yanwu Gu  Beijing Academy of
Quantum Information Sciences, Beijing 100193, China  Yunheng Ma  Beijing
Academy of Quantum Information Sciences, Beijing 100193, China  Shahid Qamar
Beijing Academy of Quantum Information Sciences, Beijing 100193, China  Chen
Qian  Beijing Academy of Quantum Information Sciences, Beijing 100193, China
Peng Qian  Beijing Academy of Quantum Information Sciences, Beijing 100193,
China  Xiao Xiao  Beijing Academy of Quantum Information Sciences, Beijing
100193, China  Meng-Jun Hu<EMAIL_ADDRESS>Beijing Academy of Quantum
Information Sciences, Beijing 100193, China  Dong E. Liu
<EMAIL_ADDRESS>State Key Laboratory of Low Dimensional Quantum
Physics, Department of Physics, Tsinghua University, Beijing, 100084, China
Beijing Academy of Quantum Information Sciences, Beijing 100193, China
Frontier Science Center for Quantum Information, Beijing 100184, China

###### Abstract

The Quantum Approximate Optimization Algorithm (QAOA), which is a variational
quantum algorithm, aims to give sub-optimal solutions of combinatorial
optimization problems. It is widely believed that QAOA has the potential to
demonstrate application-level quantum advantages in the noisy intermediate-
scale quantum(NISQ) processors with shallow circuit depth. Since the core of
QAOA is the computation of expectation values of the problem Hamiltonian, an
important practical question is whether we can find an efficient classical
algorithm to solve quantum mean value in the case of general shallow quantum
circuits. Here, we present a novel graph decomposition based classical
algorithm that scales linearly with the number of qubits for the shallow QAOA
circuits in most optimization problems except for complete graph case.
Numerical tests in Max-cut, graph coloring and Sherrington-Kirkpatrick model
problems, compared to the state-of-the-art method, shows orders of magnitude
performance improvement. Our results are not only important for the
exploration of quantum advantages with QAOA, but also useful for the
benchmarking of NISQ processors.

##  I Introduction

The rapid development of quantum computing technologies in past decades has
attracted a lot of interests from both academia and industrial community. In
2019, Google demonstrated the so-called “quantum supremacy” by using the 53-
qubit superconducting quantum processor Sycamore sycamore . Although the
Google’s claim seems controversial now claim1 ; claim2 ; claim3 , stronger
quantum advantage on the same problem has been verified by the USTC team with
a higher quality 60-qubit superconducting quantum processor Zuchongzhi zu1 ;
zu2 . These works imply that we have arrived the so-called noisy intermediate-
scale quantum (NISQ) era NISQ , where the quantum processors contain about
fifty to a few hundreds noisy qubits. Quantum error correction and fault
tolerance could in principle support the quantum computation with arbitrary
accuracy, but on the other hand, they require a much larger scale processor
with an error rate in each physical element mitigated to an extremely low
level. Unfortunately, due to high error rate of qubit operations, those NISQ
processors are limited to shallow depth circuit without error corrections. In
that case, a naturally important question arises: whether or not application-
level quantum advantage can be demonstrated in NISQ processors.

Variational quantum algorithms (VQA), which includes variational quantum
eigensolver (VQE), quantum machine learning (QML) and quantum approximate
optimization algorithm (QAOA), has been widely believed as a promising
approach for demonstrating NISQ applications, such as quantum chemistry chem1
; chem2 ; chem3 , machine learning ml1 ; ml2 ; ml3  and combinatorial
optimization op1 ; op2 ; op3 . The basic idea of VQA is to estimate the cost
function by sampling from parameterized quantum circuit and calls for
classical optimizer to find iterative parameters until convergence condition
is reached VQA . The bitstring as solution is obtained in the last stage by
sampling from quantum circuit with the optimalized parameters. The cost
function is in general expectation value of a Hamiltonian that can be written
as a linear combination of poly$(n)$ Pauli operators, e.g., Ising Hamiltonian
in the QAOA. Since the quantum circuit sampling has been shown classically
intractable as the number of qubits increase sampling1 ; sampling2 ; sampling3
, it is natural to ask whether or not quantum mean value of shallow depth
circuit can be calculated efficiently by a classical computer. If the answer
is positive, then NISQ processor is only necessary in the quantum sampling
part of QAOA.

Figure 1: Example illustration of basic procedure of graph decomposition
algorithm. a: The cost operator $\hat{C}$ is mapped to a weighed graph $G$, in
which $w_{ij}$ and $h_{k}$ are the weight of edge $(ij)$ and bias of node $k$
respectively. According to the value of $p$, we only show $p=1$ case here, the
G can be decomposed into independent subgraphs that correspond to terms
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ and
$\langle\hat{Z}_{k}\rangle_{p}$. The number of subgraphs equals to the number
of individual terms in $\hat{C}$. b: There is an one-to-one correspondence
between the subgraphs and the quantum circuits for computation of
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ and
$\langle\hat{Z}_{k}\rangle_{p}$. Detail realization is given in Algorithm 1.

In a recent paper, Bravyi, Gosset and Movassagh show that there exist the
classical algorithms scales linearly with the number of qubits for the special
case of geometrically local two-dimensional quantum circuits mean . Whether
the quantum mean value problems can be solved efficiently on a classical
computer in the case of general shallow circuits, however, remains a central
open question. Here, we show that there does exist a classical algorithm
capable of efficiently calculating quantum mean values in the case of general
shallow QAOA circuits. The core of our algorithm is to divide a large QAOA
instance into many independent small instances by using methods of graph
decomposition and graph-to-circuit mapping. Since small instances can be
handled independently, our algorithm is a natural parallel algorithm. Analysis
and numerical tests show that the run time of our algorithm scales linearly
with qubits number for most of optimization problems. For the special QAOA
instances that its connections in problem Hamiltonian scales with size of
instance, e.g. Sherrington-Kirkpatrick (S-K) model, though the run time scales
exponentially, our algorithm still behaves much better than the state-of-the-
art tensor network algorithm. Our results definitely show that current NISQ
processors has no advantages in quantum mean value problem in QAOA.

##  II The Quantum Approximate Optimization Algorithm

In the seminal paper of Farhi, Goldstone, and Gutmann, the QAOA is proposed as
a variational quantum algorithm to produce approximate solutions for
combinatorial optimization (CO) problems op1 . Since then numerous research
works on QAOA have been shown both theoretically theo1 ; theo2 ; theo3 ; theo4
; theo5 ; theo6 ; theo7 ; theo8 ; theo9 ; theo10 and experimentally ex1 ; ex2
; ex3 ; ex4 ; ex5 ; op3 . Similar to quantum annealing (QA), in which CO
problems are modeled as the form of Ising Hamiltonian, the QAOA also starts
with Ising form as cost function. The cost function is quadratic and its
general operator form reads

$\hat{C}=\sum_{i,j}w_{ij}\hat{Z}_{i}\hat{Z}_{j}+\sum_{k}h_{k}\hat{Z}_{k},$ (1)

where $\hat{Z}_{i}$ are Pauli-Z operators with eigenvalues $\pm 1$, $w_{ij}$
and $h_{k}$ represent weights and bias respectively. If we denote
computational basis vector as $|z\rangle$ with $z\equiv z_{1}z_{2}...z_{n}$
are $\\{1,-1\\}$ bitstrings, the $\hat{C}$ is diagonal in the computational
basis vector $|z\rangle$ such that $\langle z|\hat{C}|z\rangle=C(z)$.

The bitstring $z$ that minimize (or maximize) $\langle z|\hat{C}|z\rangle$ is
the optimal solution of optimization problems. Since most optimization
problems are NP hard, there is still a lack of efficient algorithms. The QAOA,
however, approximates the optimal solution by generating quantum state in
quantum circuit that as close as possible to the optimal basis vector. This is
done by introducing another operator $\hat{B}=\sum_{i}\hat{X}_{i}$ and
generating quantum state with $2p$ parameters as

$|\vec{\gamma},\vec{\beta}\rangle=\hat{U}_{B}(\beta_{p})\hat{U}_{C}(\gamma_{p})\cdots\hat{U}_{B}(\beta_{1})\hat{U}_{C}(\gamma_{1})|+\rangle,$
(2)

where $\hat{U}_{C}(\gamma)=e^{-i\gamma\hat{C}}$,
$\hat{U}_{B}(\beta)=e^{-i\beta\hat{B}}$, and the initial state
$|+\rangle=\sum_{z}|z\rangle/\sqrt{2^{n}}$ be the uniform superposition of
computational basis vectors. The expectation value of cost operator $\hat{C}$

$E_{p}(\vec{\gamma},\vec{\beta})=\langle\vec{\gamma},\vec{\beta}|\hat{C}|\vec{\gamma},\vec{\beta}\rangle$
(3)

is then optimized via outer-loop classical optimizer to find optimal
parameters $\\{\vec{\gamma}_{opt},\vec{\beta}_{opt}\\}$. Once the optimized
state $|\vec{\gamma}_{opt},\vec{\beta}_{opt}\rangle$ is obtained, we can
perform sampling from the optimized quantum circuit and output the bitstring
that minimized the $C(z)$ as approximate solution.

The core of QAOA is the computation of quantum mean value
$E_{p}(\vec{\gamma},\vec{\beta})$, and there are two different
implementations. The first one, which estimates the
$E_{p}(\vec{\gamma},\vec{\beta})$ from sampling results, is realized by making
repeated calls to the quantum processor. The alternative one is to calculate
the $E_{p}(\vec{\gamma},\vec{\beta})$ directly on the classical computers. In
both cases, constrained by the current capabilities of quantum processors and
classical simulations, the computation is limited to the case of shallow depth
circuit such that in general $p\leq 4$. In the following we focus on the
efficient classical computation of quantum mean values for shallow depth
circuit.

##  III Graph decomposition based classical algorithm

Since the cost operator $\hat{C}$ is consist of independent terms
$\hat{Z}_{i}\hat{Z}_{j}$ and $\hat{Z}_{k}$, the
$E_{p}(\vec{\gamma},\vec{\beta})$ can be recast as

$E_{p}(\vec{\gamma},\vec{\beta})=\sum_{ij}w_{ij}\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}+\sum_{k}h_{k}\langle\hat{Z}_{k}\rangle_{p},$
(4)

where
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}\equiv\langle\vec{\gamma},\vec{\beta}|\hat{Z}_{i}\hat{Z}_{j}|\vec{\gamma},\vec{\beta}\rangle$,
and the same form for $\langle\hat{Z}_{k}\rangle_{p}$. The computation of
$E_{p}(\vec{\gamma},\vec{\beta})$ now depends on the computation of individual
terms $\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ and
$\langle\hat{Z}_{k}\rangle_{p}$. The sum form of Eq. (4) indicates that an
efficient parallel algorithm is feasible if we can come up with a procedure to
calculate individual terms. This critical procedure is the graph decomposition
networkx  that be introduced below.

It should be noted here that the initial idea of graph decomposition has
already been mentioned in the seminal paper of QAOA when dealing with Max-cut
problem op2 . Those authors realized that there exists only finite subgraphs
for each edge $(ij)$ in Max-cut graph for finite $p$ and the expectation value
$E_{p}(\vec{\gamma},\vec{\beta})$ is thus determined by expectation values
related to subgraphs and the number of occurrences of the subgraphs op2 .
Surprisingly, no specific algorithm implementation is given in the follow-up
study. Besides, there are no further discussions about extending graph
decomposition of Max-cut to general combinatorial optimization problems. Our
contribution here is connecting Ising cost operator $\hat{C}$ of optimization
problems to weighted graph representation $G$ and develop specific graph
decomposition based algorithm to realize efficient calculations for shallow
depth QAOA circuits.

Algorithm 1 Graph Decomposition Algorithm

1:$\hat{C},\quad p,\quad(\vec{\gamma},\vec{\beta})$

2:$E_{p}(\vec{\gamma},\vec{\beta})$

3:weighted graph $\mathcal{G}$ $\leftarrow$ $\hat{C}$

4:subgraph_set $\leftarrow$ $\varnothing$, expectation_set $\leftarrow$
$\varnothing$

5:$E_{p}(\vec{\gamma},\vec{\beta})$ $\leftarrow$ $0$

6:for nodes and edges in $\mathcal{G}$ do

7: elements $\leftarrow$ { nodes or edges }

8: for i $\leftarrow$ $1$ to $p$ do

9: elements $\leftarrow$ elements $\cup$ adjacency nodes and edges

10: subgraph $\mathcal{G}(i,j,p)$ or $\mathcal{G}(k,p)\leftarrow$ elements

11: subgraph_set $\leftarrow$ subgraph $\cup$ subgraph_set

12:for subgraph in subgraph_set do

13: qubits number $N\leftarrow$ node numbers of subgraph

14: prepare $|+\rangle^{\otimes N}$ initial state

15: for i $\leftarrow$ 1 to $p$ do

16: for edge in subgraph do

17: $(q_{m},q_{n})$ $\leftarrow$ nodes of edge

18: $w_{mn}\leftarrow$ weight of edge

19: apply $RZZ(\gamma_{i}\times\omega_{mn})$ gate to $(q_{m},q_{n})$

20: for node in subgraph do

21: $q_{m},h_{k}\leftarrow$ $node$, weight of node

22: apply $RZ(\gamma_{i}\times h_{k})$ gate to $q_{m}$

23: apply $RX(\beta_{i}\times h_{k}$) gate to $q_{m}$

24: expectation $\leftarrow$ calculating
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ or $\langle\hat{Z}_{k}\rangle_{p}$

25: expectation_set $\leftarrow$ expectation $\cup$ expectation_set

26:for expectation in expectation_set do

27: $E_{p}(\vec{\gamma},\vec{\beta})$ $\leftarrow$
$\sum_{ij}w_{ij}\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}+\sum_{k}h_{k}\langle\hat{Z}_{k}\rangle_{p}$

Figure 2: Benchmarking results for random 3-regular Max-cut, graph coloring
and S-K model. (a) and (b) show the results of Max-cut, in which $p=1$ in (a)
and $p=2$ in (b). For $p=3,4$, we give the result in the case of fixed $1000$
qubits number. (c) show the result of coloring random generated 3-degree graph
with three colors in case of $p=1$. In the case of Max-cut and graph coloring,
Qcover scales linearly with qubits number while other methods scale
exponentially. Since the graph of S-K model is complete graph, both Qcover and
other methods scale exponentially but Qcover behaves much better as can be
seen in (d) with $p=1$.

The idea of our graph decomposition is based on two observations. The first
one is that any factors in the operators
$\hat{U}_{c}(\vec{\gamma}),\hat{U}_{B}(\vec{\beta})$ that do not involve
qubits $i$ or $j$ commute through $\hat{Z}_{i}\hat{Z}_{j},\hat{Z}_{i}$, and
can thus be cancelled out without affecting the calculation of
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p},\langle\hat{Z}_{i}\rangle_{p}$. The
other one is that the cost operator $\hat{C}$ in Eq. (1) has one-to-one
correspondence with the weighted graph $G$. The $\hat{Z}_{i}$ corresponds to
the node $i$ in $G$, while $w_{ij}$ and $h_{k}$ are the weight of edge $(ij)$
and the bias of node $k$ respectively. From the above two observations, it can
be concluded that individual term like
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ depends only on the subgraph
$G(i,j,p)$, which involves nodes $i,j$ and other nodes connecting to them with
no more than $p$ distance. Also the term like $\langle\hat{Z}_{k}\rangle$
depends on the subgraph $G(k,p)$. For each subgraph $G(i,j,p)$ or $G(k,p)$,
there exists an one-to-one correspondence with the quantum circuits
representation for the computation of
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ or $\langle\hat{Z}_{k}\rangle_{p}$,
in which the qubits number equals to nodes number in subgraphs. The edge
$(ij)$ with weight $w_{ij}$ in the subgraph corresponds to two-qubit gate
$e^{-i\gamma w_{ij}\hat{Z}_{i}\hat{Z}_{j}}$ operating on qubits $i,j$. The
node $k$ with bias $h_{k}$ in the subgraph corresponds to one-qubit operation
$e^{-i\gamma h_{k}\hat{Z}_{k}}$. We give a simple example in Fig. 1 to
illustrate the basic procedure of graph decomposition algorithm. The size of
subgraphs is determined by $p$ and degree of connection of $G$, which is
usually independent of $n$ unless the $G$ is a complete graph. The weighted
graphs of optimization problems are sparse or medium density in most cases,
making the size of subgraphs not so big for small $p$. For example, in the
case of 3-regular Max-cut, the subgraphs for $p=1$ contains only $6$ qubits
which can be easily simulated classically, regardless of the size of the
problem. Since the computational complexity of graph decomposition is $o(n)$,
our algorithm scales linearly with sparse density optimization problems in the
case of shallow circuits. The size of subgraphs grows as $p$ increases, which
increase the run time for calculation of each individual term in Eq. (4).
Suppose that the weighted graph $G$ has an average degree $d=l/n$ with $l,n$
are total number of edges and nodes respectively. If $d$ keeps constant as $n$
grows, e.g., 3-regular Max-cut, then the size of subgraphs is constant for
arbitrary finite $p$, which implies the time complexity of our algorithm be
$C(d,p)o(n)$ with $C(d,p)$ is a constant depending on the value of $(d,p)$.
For large $d$ or $p$, the value of $C(d,p)$ may large enough that beyond
current classical simulation but the resource requirement is independent of
$n$. For special case that $d$ depends on $n$, e.g., $d=n-1$ in S-K model,
$C(d,p)$ now depends on $n$ with exponential complexity in general for current
classical algorithms. For most optimization problems, $d$ is constant or
$o(\log n)$, the run time of graph decomposition algorithm scales linearly or
polynomially with nodes number $n$.

##  IV Numerical results

The pseudo code implementation of graph decomposition algorithm is given in
Algorithm 1. The input is the cost operator $\hat{C}$ and $p$, while the
output is the expectation value $E_{p}(\vec{\gamma},\vec{\beta})$. The first
step is convert the $\hat{C}$ to corresponding weighted graph $G$, and then
decompose the $G$ starting from nodes and edges of $G$ according to the value
of $p$. The total number of subgraphs is the sum of numbers of nodes and
edges. For each subgraph $G(i,j,p)$ or $G(k,p)$, it can be mapped to a
corresponding QAOA quantum circuit for the calculation of expectation value
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ or $\langle\hat{Z}_{k}\rangle_{p}$.
The simulation of quantum circuits can be implemented by kinds of method, e.g.
statevector or tensor contraction. At last the
$E_{p}(\vec{\gamma},\vec{\beta})$ is obtained by adding every terms
$\langle\hat{Z}_{i}\hat{Z}_{j}\rangle_{p}$ or $\langle\hat{Z}_{k}\rangle_{p}$.
Based on Algorithm 1, we have developed a quantum software package named
QCover Qcover  to help finding optimal parameters of shallow QAOA circuits
more faster than existing software .

To test our algorithm, we perform numerical experiments on three classical
optimization problems, i.e. 3-regular Max-cut, graph coloring and S-K model by
comparing the QCover with the IBM Qiskit qiskit  and the Quimb quimb . The IBM
Qiskit uses statevector method to realize QAOA simulation, while the Quimb
takes the method of tensor network contraction. The most important part of
tensor method is to search the optimal contraction path. The Quimb provides
many heuristic algorithms to determine contraction path and we choose the best
one in our instance tests. The engine of the QCover to compute subgraphs could
be statevector or tensor, depending on the characteristics of subgraphs. For
example, statevector engine is more efficient than tensor one in the case of
3-regular Max-cut with $p=1$, while it is better to use tensor engine for
$p\geq 2$. The QCover can automatically choose the best engine to complete the
computation. The cost operator of Max-cut and S-K are given with the same form

$\hat{C}=\sum_{i\neq j}w_{ij}\hat{Z}_{i}\hat{Z}_{j},$ (5)

where $w_{ij}=1$ for Max-cut and $w_{ij}$ is randomly chosen to be $\pm 1$ for
S-K model. For graph coloring problem, we first model it with quadratic
unconstrained binary optimization (QUBO) form and then transformed it into
Isng cost form $\hat{C}$ QUBO . The QCover has an instance library to
automatically generate the weighted graph $G$ corresponding to $\hat{C}$ for
above three problems.

The weighted graph $G$ of S-K model is a complete graph SK , while it is
constant density graph for Max-cut or graph coloring problem QUBO . The size
of subgraphs for complete graph equals to $G$ even for $p=1$ case, which
constrains the calculations for small size. For constant density graph, the
size of subgraphs depend only on the density and the value of $p$, which makes
us capable of dealing with large size instances given that the size of
subgraphs is within the scope of computing power. The Fig. 2 shows the test
results in which horizontal axis represents the node numbers (or qubit
numbers) of weighted graph $G$ and longitudinal axis represents the time to
complete calculation of $E_{p}(\vec{\gamma},\vec{\beta})$ with
$(\vec{\gamma},\vec{\beta})$ is given randomly. For each data, we show the
average time value of 5 independent tests. The QCover scales linearly with the
node numbers in Max-cut and graph coloring problems as shown in sub-figures
(a), (b) and (c) of Fig. 2, while the qiskit and Quimb behave exponentially.
For S-K model result shown in the sub-figure (d) of Fig. 2, the QCover also
has better performance but scales exponentially due to the size of subgraphs
equal to the node numbers. It should be noted here that all date from the
QCover are obtained from the personal computer run without using parallel
computing. Since the subgraphs obtained from weighted graph can be handled
individually with diffient computing cores, the QCover supports high
performance parallel computing and better speed-up can be achieved.

##  V Discussion and conclusion

Although we have shown quantum mean values for shallow QAOA circuits can be
efficiently solved by classical computers, it remains an open question for
other VQAs, e.g., VQE and QML. Besides, the sampling in the last stage of
QAOA, which has been shown classically intractable, still requires NISQ
processor. Our algorithm and software, however, can be used as a powerful
classical-assisted tool to help finding and realizing the possible
optimization problems that is suitable to demonstrate application-level
quantum advantage in NISQ processor. It is also helpful of our algorithm for
the verification and benchmarking of NISQ computers. Since QAOA deals with
Ising Hamiltonian, our algorithm has the potential to assist realizing fast
approximate ground state preparation of arbitrary long-range Ising-type
Hamiltonian ex3 .

In summary, we have present a graph decomposition algorithm to realize
efficient classical computation of quantum mean values for shallow QAOA
circuits. For most optimization problems, the run time of our algorithm scales
linearly with the instance size in the case of shallow circuits. When combined
with high-performance parallel computing in the next stage, the performance of
our algorithm can be further significantly improved. Our algorithm and the
related software Qcover, when used with NISQ processor, may accelerate the
demonstration of application-level quantum advantage.

###### Acknowledgements.

Acknowledgments: The authors would like to thank Bai-Ting Liu, Xingyao Wu,
Diqing Chen and Yuxuan Wang for valuable discussions. This work is supported
by the Beijing Academy of Quantum Information Sciences. Data availability: All
codes and data can be found in GitHub https://github.com/BAQIS-Quantum/Qcover.

## References

  * (1) F. Arute et al., Quantum supremacy using a programmable superconducting processor, Nature 574, 505-510 (2019).
  * (2) Cupjin Huang et al., Classical simulation of quantum supremacy circuits, arXiv:2005.06787 (2020).
  * (3) Feng Pan and Pan Zhang. Simulating the sycamore quantum supremacy circuits, arXiv:2103.03074 (2021).
  * (4) Yong Liu et al., Closing the “Quantum Supremacy” Gap: Achieving Real-Time Simulation of a Random Quantum Circuit Using a New Sunway Supercomputer, arXiv:2110.14502 (2021).
  * (5) Yulin Wu et al., Strong Quantum Computational Advantage Using a Superconducting Quantum Processor, Phys. Rev. Lett. 127, 180501 (2021).
  * (6) Qingling Zhu et al., Quantum computational advantage via 60-qubit 24-cycle random circuit sampling, Science Bulletin (2021).
  * (7) J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * (8) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Simulated Quantum Computation of Molecular Energies, Science 309, 1704-1707 (2005).
  * (9) P. J. J. O’Malley et al., Scalable Quantum Simulation of Molecular Energies, Phys. Rev. X 6, 031007 (2016).
  * (10) F. Arute et al., Hartree-Fock on a superconducting qubit quantum computer, Science 365, 1084-1089 (2020).
  * (11) J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195-202 (2016).
  * (12) M. Benedetti, E, Lloyd, S. Sack, and M. Fiorentini, Parameterized quantum circuits as machine learning models, Quantum Sci. Technol. 4, 043001 (2019).
  * (13) He-Liang Huang et al., Experimental Quantum Generative Adversarial Networks for Image Generation, Phys. Rev. Applied 16, 024051 (2021).
  * (14) E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv:1411.4028 (2014).
  * (15) G. E. Crooks, Performance of the Quantum Approximate Optimization Algorithm on the Maximum Cut Problem, arXiv:1811.08419 (2018).
  * (16) M. P. Harrigan et al., Quantum approximate optimization of non-planar graph problems on a planer superconducting processor, Nat. Phys. 17, 332-336 (2021).
  * (17) M. Cerezo et al., Variational quantum algorithms, Nat. Rev. Phys. 3, 625-644 (2021).
  * (18) E. Farhi and A. W. Harrrow, Quantum Supremacy through the Quantum Approximate Optimization Algorithm, arXiv:1602.07674 (2016).
  * (19) A. P. Lund, M. J. Bremner, and T. C. Ralph, Quantum sampling problems, BosonSampling and quantum supremacy, npj Quantum Information 3, 15 (2017).
  * (20) Sergio Boixo et al., Characterizing quantum supremacy in near-term devices, Nat. Phys. 14, 595-600 (2018).
  * (21) S. Bravyi, D. Gosset, and R. Movassagh, Classical algorithms for quantum mean values, Nat. Phys. 17, 337-341 (2021).
  * (22) S. Lloyd, Quantum approximate optimization is computationally universal, arXiv:1812.11075 (2018).
  * (23) M. Streif and M. Leib, Training the Quantum Approximate Optimization Algorithm without access to a Quantum Processing Unit, arXiv:1908.08862 (2019).
  * (24) C. Huang, M. Szegedy, F, Zhang, X. Gao, J. Chen, and Y. Shi, Alibaba Cloud Quantum Development Platform: Applications to Quantum Algorithm Design, arXiv:1909.02559 (2019).
  * (25) V. Akshay, H. Philathong, M. E. S. Morales, and J. D. Biamonte, Reachability Deficits in Quantum Approximate Optimization, Phys. Rev. Lett. 124, 090504 (2020).
  * (26) P. Vikstål, M. Grönkvist, M. Svensson, M. Andersson, G. Johansson, and G. Ferrini, Phys. Rev. Applied 14, 034009 (2020).
  * (27) L. Zhu et al., An Adaptive Quantum Approximate Optimization Algorithm for Solving Combinatorial Problems on a Quantum Computer, arXiv:2005.10258 (2020).
  * (28) L. Zhou, S. T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Phys. Rev. X 10, 021067 (2020).
  * (29) A. Ozaeta, W. V. Dam, and P. L. McMahon, Expectation Values from the Single-Layer Quantum Approximate Optimization Algorithm on Ising Problems, arXiv:2012.03421 (2021).
  * (30) J. Basso, E. Farhi, K. Marwaha, B. Villalonga, and L. Zhou, The Quantum Approximate Optimization Algorithm at High Depth for MaxCut on Large-Girth Regular Graphs and the Sherrington-Kirkpatrick Model, arXiv:2110.14206 (2021).
  * (31) M. Medvidović and G. Carleo, Classical variational simulation of the Quantum Approximate Optimization Algorithm, npj Quantum Inf 7, 101 (2021).
  * (32) J. S. Otterbach et al., Unsupervised Machine Learning on a Hybrid Quantum Computer, arXiv:1712.05771 (2017).
  * (33) X. Qiang et al., Large-scale silicon quantum photonics implementing arbitrary two-qubit processing, Nat. Photon. 12, 534-539 (2018).
  * (34) G. Pagano et al., Quantum Approximate Optimization of the Long-Range Ising Model with a Trapped-Ion Quantum Simulator, PNAS 117, 25396-25401 (2020).
  * (35) A. Bengtsson et al., Improved Success Probability with Greater Circuit Depth for the Quantum Approximate Optimization Algorithm, Phys. Rev. Applied 14, 034010 (2020).
  * (36) M. Willsch, D. Willsch, F. Jin, H. D. Raedt, and K. Michielsen, Benchmarking the Quantum Approximate Optimization Algorithm, Quantum Inf. Process. 19, 197 (2020).
  * (37) https://networkx.org.
  * (38) https://github.com/BAQIS-Quantum/Qcover.
  * (39) M. S. Anis et al., Qiskit: An Open-source Framework for Quantum Computing, DOI: 10.5281/zenodo.2573505 (2021) https://github.com/Qiskit.
  * (40) J. Gray, quimb: A python package for quantum information and many-body calculations, Journal Of Open Source Software, 3, 819 (2018) https://github.com/jcmgray/quimb.
  * (41) F. Glover, G. Kochenberger, and Y. Du, A Tutorial on Formulating and Using QUBO Models, arXiv:1811.11538 (2018).
  * (42) E. Farhi, J. Goldstone, S. Gutmann, and L. Zhou, The Quantum Approximate Optimization Algorithm and the Sherrington-Kirkpatrick Model at Infinite Size, arXiv:1910.08187 (2019).