# Constraints on Hadronic Contribution to LHAASO Sources with Neutrino
Observations

Tian-Qi Huang Department of Astronomy, School of Physics, Peking University,
Beijing 100871, China  Zhuo Li  Department of Astronomy, School of Physics,
Peking University, Beijing 100871, China  Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, China

###### Abstract

LHAASO detected 12 gamma-ray sources above 100 TeV which are the possible
origins of Galactic cosmic-rays. We summarize the neutrino measurements by
IceCube and ANTARES in the vicinity of LHAASO sources to constrain the
contribution of hadronic gamma-rays in these sources. We find that the current
observations constrain that the hadronic gamma-rays contribute no more than
$\sim 60\%$ of the gamma-rays from Crab Nebula. Gamma-rays from two LHAASO
sources, LHAASO J1825-1326 and LHAASO J1907+0626, are dominated by leptonic
components up to $\sim 200$ TeV, under the hypotheses in the analysis by
IceCube. The uncertainties of the constraint on the hadronic gamma-ray
emission are discussed. We also constrain the total 100 TeV gamma-ray emission
from TeV PWNe relying on the remarkable sensitivity of LHAASO at that
energies.

Gamma-ray astronomy (628), Cosmic ray sources (328), Neutrino astronomy (1100)

##  1 Introduction

Cosmic-rays extending to several PeV energies are believed to originate from
Galactic sources, called PeVatrons. It is hard to localize PeVatrons with only
cosmic-ray observations, because charged cosmic-rays below several PeV are
deflected frequently by magnetic field while traveling in the Galaxy.
Fortunately, PeV cosmic-rays produce high energy ($\rm>$100 TeV) gamma-rays
and neutrinos in the collisions with baryon background through hadronuclear
($pp$) interactions or radiation background through photohadronic ($p\gamma$)
interactions. Those gamma-rays and neutrinos are good tracers of PeVatrons, as
both of them travel in a straight line after escaping from the sources.

Gamma-rays can be produced through hadronic process ($pp$ or $p\gamma$) and/or
leptonic process (inverse Compton scattering). The production mechanisms of
observed gamma-rays usually are unclear. However, compared with TeV gamma-
rays, the production of $\gtrsim 100$TeV gamma-rays with inverse Compton
scatterings suffers more stringent Klein–Nishina suppression. The sources with
$\sim 100$ TeV emission have more chance to be cosmic ray PeVatrons (e.g. HESS
Collaboration et al., 2016; Albert et al., 2020a; Tibet AS$\gamma$
Collaboration et al., 2021).

The Large High Altitude Air Shower Observatory (LHAASO) is a hybrid extensive
air shower array for cosmic-ray and gamma-ray studies. Relying on the
unprecedented sensitivity to gamma-rays around hundreds TeV, the Kilometer
Square Array (KM2A) of LHAASO detected 12 gamma-ray sources above 100 TeV (Cao
et al., 2021). They also observed PeV gamma-rays from the direction of Cygnus
region and Crab nebula (The LHAASO Collaboration et al., 2021). These LHAASO
gamma-ray sources provide us a group of target sources to find out the cosmic
ray PeVatrons. We should look for neutrino signals from these sources to pin
down the answer.

The IceCube neutrino observatory is a cubic-kilometer array in the deep ice at
the South Pole and currently the most sensitive detector for neutrinos from
TeV to PeV. IceCube has carried out some approaches to search for neutrino
sources from known Galactic sources, but only gives upper limits to the flux.
First, in the Source List Search (SLS), they select some known astrophysical
sources to be the Neutrino Source Candidates (NSCs), and searches for neutrino
signals from them individually (e.g. Aartsen et al., 2019). The signals are
not significant enough. Second, since stacking may help to enhance signal to
noise ratio and discover weak signals, IceCube also tries to search for
neutrinos over some catalogs of known sources, i.e., Stacked Source Search
(SSS) (e.g. Aartsen et al., 2020a). Still, only upper limits are given to the
neutrino flux so far. However, as neutrinos are produced together with gamma-
rays in the hadronic interactions, one can derive the upper limits of neutrino
flux by that of hadronic gamma-ray flux, and hence constrain the contribution
of hadronic process.

With the goal of identifying the Galactic cosmic ray Pevatrons from the LHAASO
detected 100 TeV sources, in this work we use the neutrino measurements of
NSCs in the vicinity of LHAASO sources, and constrain the hadronic
contributions to gamma-ray emissions. In section 2, we introduce the neutrino
source searches used in this paper and the criteria to select the NSCs. In
section 3, we show the method to derive the upper limits on hadronic gamma-ray
flux. In section 4, we compare the gamma-ray observations with upper limits of
hadronic gamma-ray flux. In section 5, we discuss previous works and evaluate
the uncertainties in the analyses. We finally give conclusions in section 6.

##  2 Source Selection

IceCube uses track-like neutrino events to measure the astrophysical neutrino
flux from NSCs. The median location errors of these neutrinos range from 0.25
degree (above PeV) to 1 degree (sub-TeV) (Aartsen et al., 2020a), while most
LHAASO sources represent diffuse structures with angular extensions as large
as 1 degree except for the Crab Nebula and LHAASO J2108+5157 (Cao et al.,
2021). Astrophysical neutrinos from the directions of LHAASO sources probably
have the same origins with those of gamma-rays.

We select the NSCs within 1 degree from the measured center of a LHAASO source
(see Table 1). The NSCs are taken from three SLSes: the ten-year search
(Aartsen et al., 2020a), the eight-year search (Aartsen et al., 2019) and the
combined search by ANTARES & IceCube (Albert et al., 2020b). The upper limits
of neutrino flux given by SLSes depend on the hypotheses on spectral shapes
(usually assumed as a single power law) and point-source hypothesis. We also
choose the NSCs from two SSSes: the TeV Pulsar Wind Nebulae (PWNe) search
(Aartsen et al., 2020b) and the Galactic catalog searches (Aartsen et al.,
2017a) including five catalogs (one Milagro catalog, one HAWC catalog and
three Supernova Remnant (SNR) catalogs). The upper limits given by SSSes
depend on not only the spectral shapes and assumed source extensions, but also
the weighting factors of neutrino flux from individual sources in the certain
catalog. The weighting factor is the fraction of neutrino flux from individual
source and depends on the weighting scheme used. The details about the
weighting schemes are explained in Appendix A.

Nine LHAASO sources have NSCs in the neighborhood. In order to compare the
gamma-ray emission with neutrino emission in the range of TeV to PeV, we
utilize previous TeV gamma-ray measurements of the certain NSCs to obtain the
energy spectra at lower energies. The details about the TeV gamma-ray
observations of the NSCs are described in Appendix B.

##  3 Method

In hadronuclear scenarios, the differential flux of gamma-rays and neutrinos
($\Phi_{\gamma,\,\nu}=dN_{\gamma,\,\nu}/dE_{\gamma,\,\nu}$) are related as

$\frac{1}{3}\sum_{\rm\alpha}\frac{dN_{\nu_{\rm\alpha}+\bar{\nu}_{\rm\alpha}}}{dE_{\rm\nu}}=K_{\rm\pi}\frac{dN_{\rm\gamma}}{dE_{\rm\gamma}}$
(1)

where $\rm\alpha$ represents the neutrino flavor, $E_{\rm\gamma}=2E_{\rm\nu}$
and $K_{\rm\pi}\simeq 2$ is the number ratio of charged to neutral pions. For
an $dN/dE\propto E^{-\gamma}$ spectrum, the upper limit on hadronic gamma-ray
flux at 100 TeV can be expressed as

$\Phi_{\rm\gamma}^{\rm UL}(100~{}{\rm
TeV})=\frac{1}{2}\Phi_{\rm\nu_{\mu}+\bar{\nu}_{\mu}}^{\rm
UL}(E_{\rm\nu})\left(\frac{50~{}{\rm TeV}}{E_{\rm\nu}}\right)^{-\gamma}$ (2)

where $\Phi_{\rm\nu_{\mu}+\bar{\nu}_{\mu}}^{\rm UL}$ is the 90% upper limit of
neutrino flux.

Gamma-rays are partially absorbed in the propagation through interstellar
radiation field (ISRF) and cosmic microwave background (CMB) due to pair
productions ($\gamma\gamma\rightarrow e^{+}e^{-}$). Such absorption is
considered in the conversion of neutrino upper limit into gamma-ray upper
limit, if the distance to the NSC is available. The ISRF energy density is
taken from Popescu et al. (2017). Gamma-ray absorption around the source is
not considered.

##  4 Results

Figure 1:  Comparison between gamma-ray flux by LHAASO observations and the
upper limits of hadronic gamma-ray flux derived from SLSes. The grey lines
present the 90% upper limits of hadronic gamma-ray flux, where the dark parts
mark the central 90% energy ranges according to neutrino observations. The
solid, dot-dashed and dashed lines correspond to $E^{-2.0}$, $E^{-2.5}$ and
$E^{-3.0}$ neutrino spectra, respectively. Blue symbols are LHAASO
observations (Cao et al., 2021; The LHAASO Collaboration et al., 2021). Other
symbols and shaded areas are the gamma-ray measurements of (a) Crab Nebula by
Tibet AS$\gamma$ (Amenomori et al., 2019), H.E.S.S. (Aharonian et al., 2006)
and MAGIC (Albert et al., 2008); (b) HESS J1826-130 (H. E. S. S. Collaboration
et al., 2020), HESS J1825-137 (H. E. S. S. Collaboration et al., 2018a), eHWC
J1825-134 (Abeysekara et al., 2020); (c) HESS J1841-055 (Aharonian et al.,
2008); (d) HESS J1843-033 and HESS J1844-030 (H. E. S. S. Collaboration et
al., 2018b), 2HWC J1844-032 (Abeysekara et al., 2017); (e) HESS J1849-000 (H.
E. S. S. Collaboration et al., 2018b); (f) HESS J1908+063 (Aharonian et al.,
2009), eHWC J1907+063 (Abeysekara et al., 2020); (g) VER J2019+368 (Abeysekara
et al., 2018), eHWC J2019+368 (Abeysekara et al., 2020); (h) Cyg X-3
(Sinitsyna & Sinitsyna, 2020), TeV J2032+4130 (Aliu et al., 2014), MGRO
J2031+41 (Abdo et al., 2007, 2009), 2HWC J2031+415 (Abeysekara et al., 2017).
The gamma-ray absorption due to ISRF and CMB is considered if the distance is
available.  Figure 2: Comparison between gamma-ray flux by LHAASO observations
and the upper limits of hadronic gamma-ray flux derived from SSSes. Lines show
the 90% upper limits of hadronic gamma-ray flux, where the dark parts mark the
central 90% energy ranges according to neutrino observations. The solid,
dotted, dot-dashed and dashed lines correspond to $E^{-2.0}$, $E^{-2.19}$,
$E^{-2.5}$ and $E^{-3.0}$ neutrino spectra respectively. The colors of lines
represent different source catalogs: black lines correspond to the TeV PWNe
catalog (inverse age weighting); red lines to the catalog of SNRs with
associated PWN; green lines to the HAWC catalog and purple lines to the
Milagro catalog. The symbols in panel (i) are gamma-ray flux from G106+2.7
observed by VERITAS (Acciari et al., 2009) and Tibet $\rm AS\gamma$ (Tibet
AS$\gamma$ Collaboration et al., 2021), MRGO C4 (Abdo et al., 2007, 2009) and
LHAASO J2227+6057 (Cao et al., 2021). Symbols and shaded areas in other panels
are the same as those in Figure 1.  Figure 3: The cumulative gamma-ray flux
from TeV PWNe compared with the upper limits of hadronic gamma-ray flux. There
are 35 PWNe in the catalog. The cumulative gamma-ray flux from all the PWNe
(black solid line) is 26% lower than the cumulative gamma-ray flux
extrapolated from low energy before LHAASO measurement (grey dashed line). The
light grey shaded area indicates the uncertainty of the extrapolated spectrum.
Colorful lines show the upper limits on the hadronic gamma-ray emission for
four weighting schemes (Aartsen et al., 2020b).

The constraints on the hadronic gamma-ray emissions derived from SLSes are
shown in Figure 1. If there are more than one NSCs associated with the LHAASO
source, we only show the upper limit derived from the NSC closest to the
center. The upper limit for $E^{-2}$ spectrum is generally consistent with or
higher than the observed gamma-ray flux around 100 TeV in each panel. The
upper limits for $E^{-3}$ and $E^{-2.5}$ are also generally consistent with
the observed gamma-ray flux from 10 TeV to 100 TeV but much higher at the
energies below 10 TeV (see panel b, c, d, f, and g). We also show in the
figure the energy range corresponding to where 90% of signal neutrino events
will concentrate in. This energy range is determined by the assumed neutrino
spectrum, as well as the effective area of the neutrino detector. The gamma-
ray absorption due to ISRF and CMB is not significant around 100 TeV.

The constraints on the hadronic gamma-ray emissions derived from SSSes are
shown in Figure 2. The hadronic gamma-ray flux from individual source in the
catalog should be lower than the upper limit derived from the stacking
neutrino flux. The upper limit given by the HAWC catalog set strong
constraints on LHAASO J1825-1326 (panel b) and LHAASO J1908+0621 (panel f)
above 100 TeV. The leptonic component is even dominant around 200 TeV. For the
TeV PWNe, the only shown are the upper limits for inverse age weighting,
because they are lower than those for the other weighting schemes, and give
stronger constraints on the hadronic gamma-rays (Aartsen et al., 2020b). The
constraints above are reliable when the hypotheses in the SSSes are consistent
with the real case. In the HAWC catalog search, the sources are assumed to be
$0.5^{\circ}$ extended and the neutrino flux following an $E^{-2.5}$ spectrum
is assumed to be proportional to the gamma-ray flux at pivot energy for each
source. The assumption on neutrino flux ratio brings large uncertainty to the
upper limit estimate (see subsection 5.2).

We further estimate the cumulative gamma-ray flux from TeV PWNe with LHAASO
observations and compare it to the upper limits given by the TeV PWNe search
in the Figure 3. In the TeV PWNe search, Aartsen et al. (2020b) compared the
differential upper limits for different weightings with the cumulative gamma-
ray flux extrapolated from TeV observations to constrain the hadronic
component from TeV to 100 TeV. Six of these PWNe are associated with LHAASO
sources, while seven of them are not detected but in the data taken region of
LHAASO-KM2A (sky declination band $-15^{\circ}<\delta<75^{\circ}$). With the
observations at 100 TeV, we correct the extrapolated gamma-ray spectrum as the
black solid line in Figure 3 (details in Appendix C). The cumulative gamma-ray
flux at 100 TeV is 26% lower than the extrapolated value (grey dashed line)
but still consistent with the uncertainty (grey shaded area).

##  5 Discussion

###  5.1 Compare with previous works

We compare the neutrino spectrum model of Crab Nebula in Kappes et al. (2007)
with the neutrino observation results in the ten-year search. First of all,
the upper limit on neutrino flux (with a certain spectral shape)
$\Phi^{90\%}_{\nu}$ is related to the upper limit on the signal event number
$n_{90}$ as

$n_{\rm 90}=\int dt\int
dE_{\rm\nu}\Phi_{\rm\nu_{\mu}+\bar{\nu}_{\mu}}^{90\%}A_{\rm
eff}(E_{\rm\nu},\delta_{s})$ (3)

where $A_{\rm eff}$ is the effective area for neutrinos from the direction of
source ($E_{\nu}$ is the neutrino energy and $\delta_{s}$ is the decliantion
of source). The effective area is available in the IceCube muon-track data
from 2008 to 2018 (IceCube Collaboration et al., 2021).

In the seven-year search one has $n_{90}=7.0$ and $21.2$ for the spectra of
$E^{-2}$ and $E^{-3}$, respectively (Aartsen et al., 2017b). The two values
both decrease in the ten-year search (Aartsen et al., 2020a, $n_{90}=5.8$ and
$20.4$, respectively). Note that the local pre-trial p-value for Crab Nebula
is larger in the ten-year search (from 0.34 in the seven-year search to 0.49),
indicating more background like. Therefore, we assume that $n_{90}$ using
Kappes model decreases between two searches. For the Kappes model a value
$n_{90}=20.5$ is derived from the seven-year limit (Aartsen et al., 2017b)
with Equation 3, which can be considered as the upper limit for ten-year
search, and hence results in a flux limit (red line in Figure 4) by Equation 3
and Equation 1.

Figure 4:  The best-fit gamma-ray spectrum $\Phi_{\gamma}$ (black line) of the
Crab Nebula and the upper limit of hadronic gamma-ray flux $\Phi^{\rm
Kappes}_{\rm\gamma}$ (red line). The dark part of red line shows the central
90% energy range of signal events.

We fit the gamma-ray observations with a log-parabola model
$(E/10\,\rm{TeV})^{-\gamma_{1}-\gamma_{2}\,\rm{log}(E/10\,\rm{TeV})}$ (black
line in Figure 4). The number of hadronic gamma-rays in the central 90% energy
range contributes $\rm 62^{+5}_{-4}\%$ to the total gamma-rays at most, while
the energy of hadronic gamma-rays contributes $\rm 62{\pm}4\%$ at most. Thus,
the current observations suggest that hadronic contribution can not account
for the bulk of the gamma-ray emission from Crab Nebula.

Furthermore, there have been many discussion in the literatures on the
hadronic gamma-ray and neutrino emission from potential Galactic PeVatrons and
the Galactic plane, say, for PWNe (e.g. Amato et al., 2003), for SNRs (e.g.
Mandelartz & Becker Tjus, 2015), for young massive star clusters (YMCs) (e.g.
Anchordoqui et al., 2007), and for diffuse Galactic gamma-ray and neutrino
emission (Ahlers & Murase, 2014; Neronov et al., 2014). Some TeV gamma-ray
sources have been well studied to evaluate the probability of being PeVatrons,
like the Cygnus Region (e.g. Anchordoqui et al., 2007), MGRO J2019+37 (e.g.
Beacom & Kistler, 2007) and MGRO J1908+06 (e.g. Halzen et al., 2017). However,
the LHAASO detected $100$ TeV sources provide a group of sample sources for
candidate PeVatrons. The high sensitivity of LHAASO at 100 TeV ranges also
help to put stringent constraint on the PeVatron models. A wide range of
parameter space in the Crab Nebula model by Amato et al. (2003) can be
constrained by LHAASO observation at 100 TeV.

###  5.2 Uncertainties in analysis

The upper limit on the hadronic gamma-ray flux depends on the hypotheses of
source extension. As the signals are not that significant over the background,
the upper limit ratio for Extended-Source (ES) hypothesis approximates as
$\Phi^{90\%}_{\rm PS}\times(\Omega_{\rm ES}/\Omega_{\rm PSF})^{1/2}$, where
$\Phi^{90\%}_{\rm PS}$ is the upper limit for Point-Source (PS) hypothesis,
$\Omega_{\rm ES}$ is the angular size of extended source and $\Omega_{\rm
PSF}$ is the angular size of point spread function. The angular size is
defined as $\Omega_{\rm ES}=\pi(\sigma_{s}^{2}+\sigma_{\rm IC}^{2})$ and
$\Omega_{\rm PSF}=\pi\sigma_{\rm IC}^{2}$, where $\sigma_{s}$ is the source
extension and $\sigma_{\rm IC}=0.64^{\circ}$ is given by Halzen et al. (2017).
Three LHAASO sources (LHAASO J1825-1326, LHAASO J1908+0621 and LHAASO
J2226+6057) have extension measurements with Gaussian template. They are
$0.30^{\circ}$, $0.58^{\circ}$ and $0.36^{\circ}$ respectively (Cao et al.,
2021). If their neutrino counterparts have the same extensions, their upper
limits will be 10%, 35% and 15% higher in comparison with the PS hypothesis.
The upper limits given by SLSes will be harder to constrain the hadronic
component for extended sources (panel b-h of Figure 1), while source
extensions have been considered in the SSSes (e.g., HAWC sources in the
Galactic catalog searches are assumed to be $0.5^{\circ}$ extended).

The upper limit on the hadronic gamma-ray flux also relies on the hypotheses
of spectral shapes. A Bayesian analysis of the IceCube muon-track data from
2008 to 2018 is carried out to test the impact of spectral shapes on the upper
limits (Huang & Li, 2021, in prep). The hadronic component is assumed to
follow the best-fit gamma-ray spectrum with a scaling factor. Five LHAASO
sources are tested: LHAASO J0534+2202 (Crab Nebula), LHAASO J1825-1326, LHAASO
J1908+0621, LHAASO J2018+3651 and LHAASO J2226+6057. Only the hadronic gamma-
ray component of Crab Nebula is well constrained to be no more than 70% of
total gamma-ray flux. As for other sources, the upper limits on hadronic
gamma-ray flux are higher than the flux observed.

The weighting schemes used in the SSSes influence strongly the upper limit
estimate. In the TeV PWNe search, the upper limit for the equal weighting is
around 3 times of the one for the gamma-ray flux weighting (see the right
panel of Figure 3). In the stacking search for neutrinos from AGN and
starburst galaxies, the upper limits for the equal weighting vary from 0.9 to
6.1 times of the upper limits for the gamma-ray flux weighting when changing
source samples. In the HAWC catalog search, the neutrino flux from each source
is assumed to be proportional to the gamma-ray flux at the pivot energy. If
the upper limits derived from the HAWC catalog search (green lines in Figure
2) increase by a factor of six, we cannot simply use the upper limit on the
stacking neutrino flux to constrain the hadronic gamma-ray flux from LHAASO
J1908+0621 and LHAASO J1825-1326.

Kappes et al. (2007) approximate the neutrino and gamma-ray relation as
$\Phi_{\nu}(E)=(0.694-0.16\gamma)\Phi_{\gamma}(E)$, especially for primary
proton spectrum being $E^{-\gamma}$ ($1.8<\gamma<3.0$) with a high energy
exponential cutoff. If using this approximation, our upper limits on hadronic
gamma-ray flux will be increased by 34% (17%) for the $E^{-2}$ ($E^{-3}$)
spectrum.

The radiation background around LHAASO sources also influence the constraints
on hadronic gamma-ray component. The gamma-ray sources associated with young
massive stars may have a denser infrared background compared with ISRF. The
higher optical depths around 100 TeV induce the lower upper limits on the
hadronic gamma-ray flux observed. In photon-rich environment, $p\gamma$
interactions should also be considered in addition to the $pp$ interactions.

##  6 Conclusions

We adopt the neutrino observation results in the directions within one degree
from LHAASO sources, and transfer the neutrino upper limits into those of
hadronic gamma-ray flux. The upper limits derived from the SLSes are
marginally consistent with the gamma-ray flux for most LHAASO sources and
cannot constrain the hadronic contributions stringently. The upper limits
derived from the SSSes set strong constraints on the hadronic components above
100 TeV. These constraints depend on hypotheses of not only the spectral
shapes and source extensions but also the weighting schemes.

The main conclusions of this paper are summarized below:

i) The hadronic gamma-ray component contributes no more than 62% for the Crab
Nebula.

ii) LHAASO J1825-1326 and LHAASO J1908+0621 are leptonic dominant up to 200
TeV, following the assumption that each source in the HAWC catalog has an
intrinsic extension $\lesssim 0.5$ degree and follows an $E^{-2.5}$ spectrum
with neutrino flux proportional to its gamma-ray flux at pivot energy.

iii) The cumulative gamma-ray flux from TeV PWNe is lower by 26% at 100 TeV
than the extrapolation from low energy before LHAASO measurement, leaving
larger room for hadronic dominated emission.

The constraints on hadronic gamma-ray flux will be stronger as the neutrino
statistics increase in the following 10 years. For IceCube, the time evolution
of the $\rm 5\,\sigma$ discovery potential is close to $1/T$ ($T$ is the
observation time) (Aartsen et al., 2017b). We take the discovery potential in
the ten-year point source search as the baseline, corresponding to the data
from 2008 to 2018. IceCube will be able to discover the point-like neutrino
source two times fainter by 2028. If the intrinsic gamma-ray spectrum of
LHAASO J1849-0003 follows $E^{-2}$ up to PeV, the neutrino observations with
sensitivity two times better will offer more clues to answer the origins of
these high energy gamma-rays.

A stacking search for astrophysical neutrinos from the twelve LHAASO sources
is required, which will increase the significance of signals and constrain the
hadronic component more strictly. As for the LHAASO sources in the southern
sky ($\delta<-5^{\circ}$), the joint search of IceCube and ANTARES will
improve the sensitivity by a factor ${\sim 2}$ compared to both individual
analyses (Albert et al., 2020b). The instrumented volume of high energy
neutrino detectors will be three times larger when KM3NeT/ARCA (two blocks)
(Adrián-Martínez et al., 2016) and Baikal-GVD (15 clusters) (Safronov, 2020)
fully operate in the next 5 to 10 years.

## Acknowlegdments

We thank Qinrui Liu, Jon Dumm and Ruo-Yu Liu for useful discussions. This work
is supported by the Natural Science Foundation of China (No. 11773003,
U1931201).

## References

  * Aartsen et al. (2017a) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017a, ApJ, 849, 67, doi: 10.3847/1538-4357/aa8dfb
  * Aartsen et al. (2017b) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017b, ApJ, 835, 151, doi: 10.3847/1538-4357/835/2/151
  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, European Physical Journal C, 79, 234, doi: 10.1140/epjc/s10052-019-6680-0
  * Aartsen et al. (2020a) —. 2020a, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Aartsen et al. (2020b) —. 2020b, ApJ, 898, 117, doi: 10.3847/1538-4357/ab9fa0
  * Abdo et al. (2007) Abdo, A. A., Allen, B., Berley, D., et al. 2007, ApJ, 664, L91, doi: 10.1086/520717
  * Abdo et al. (2009) Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ, 700, L127, doi: 10.1088/0004-637X/700/2/L127
  * Abeysekara et al. (2016) Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2016, ApJ, 817, 3, doi: 10.3847/0004-637X/817/1/3
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, ApJ, 843, 40, doi: 10.3847/1538-4357/aa7556
  * Abeysekara et al. (2018) Abeysekara, A. U., Archer, A., Aune, T., et al. 2018, ApJ, 861, 134, doi: 10.3847/1538-4357/aac4a2
  * Abeysekara et al. (2020) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2020, Phys. Rev. Lett., 124, 021102, doi: 10.1103/PhysRevLett.124.021102
  * Acciari et al. (2009) Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, ApJ, 703, L6, doi: 10.1088/0004-637X/703/1/L6
  * Adrián-Martínez et al. (2016) Adrián-Martínez, S., Ageron, M., Aharonian, F., et al. 2016, Journal of Physics G Nuclear Physics, 43, 084001, doi: 10.1088/0954-3899/43/8/084001
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561, doi: 10.1038/s41550-019-0724-0
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 457, 899, doi: 10.1051/0004-6361:20065351
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 477, 353, doi: 10.1051/0004-6361:20078516
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 499, 723, doi: 10.1051/0004-6361/200811357
  * Ahlers & Murase (2014) Ahlers, M., & Murase, K. 2014, Phys. Rev. D, 90, 023010, doi: 10.1103/PhysRevD.90.023010
  * Albert et al. (2020a) Albert, A., Alfaro, R., Alvarez, C., et al. 2020a, ApJ, 896, L29, doi: 10.3847/2041-8213/ab96cc
  * Albert et al. (2020b) Albert, A., André, M., Anghinolfi, M., et al. 2020b, ApJ, 892, 92, doi: 10.3847/1538-4357/ab7afb
  * Albert et al. (2021) Albert, A., Alfaro, R., Alvarez, C., et al. 2021, ApJ, 907, L30, doi: 10.3847/2041-8213/abd77b
  * Albert et al. (2008) Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ, 674, 1037, doi: 10.1086/525270
  * Aliu et al. (2014) Aliu, E., Aune, T., Behera, B., et al. 2014, ApJ, 783, 16, doi: 10.1088/0004-637X/783/1/16
  * Amato et al. (2003) Amato, E., Guetta, D., & Blasi, P. 2003, A&A, 402, 827, doi: 10.1051/0004-6361:20030279
  * Amenomori et al. (2019) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Phys. Rev. Lett., 123, 051101, doi: 10.1103/PhysRevLett.123.051101
  * Anchordoqui et al. (2007) Anchordoqui, L. A., Beacom, J. F., Goldberg, H., Palomares-Ruiz, S., & Weiler, T. J. 2007, Phys. Rev. D, 75, 063001, doi: 10.1103/PhysRevD.75.063001
  * Beacom & Kistler (2007) Beacom, J. F., & Kistler, M. D. 2007, Phys. Rev. D, 75, 083001, doi: 10.1103/PhysRevD.75.083001
  * Cao et al. (2021) Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33, doi: 10.1038/s41586-021-03498-z
  * Ge et al. (2020) Ge, C., Liu, R.-Y., Niu, S., Chen, Y., & Wang, X.-Y. 2020, arXiv e-prints, arXiv:2012.11531.  https://arxiv.org/abs/2012.11531
  * H. E. S. S. Collaboration et al. (2018a) H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2018a, VizieR Online Data Catalog, J/A+A/621/A116
  * H. E. S. S. Collaboration et al. (2018b) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018b, A&A, 612, A1, doi: 10.1051/0004-6361/201732098
  * H. E. S. S. Collaboration et al. (2020) H. E. S. S. Collaboration, Abdalla, H., Adam, R., et al. 2020, A&A, 644, A112, doi: 10.1051/0004-6361/202038851
  * Halzen et al. (2017) Halzen, F., Kheirandish, A., & Niro, V. 2017, Astroparticle Physics, 86, 46, doi: 10.1016/j.astropartphys.2016.11.004
  * HESS Collaboration et al. (2016) HESS Collaboration, Abramowski, A., Aharonian, F., et al. 2016, Nature, 531, 476, doi: 10.1038/nature17147
  * Huang & Li (2021, in prep) Huang, T.-Q., & Li, Z. 2021, in prep
  * HURVICH & TSAI (1991) HURVICH, C. M., & TSAI, C.-L. 1991, Biometrika, 78, 499, doi: 10.1093/biomet/78.3.499
  * IceCube Collaboration et al. (2021) IceCube Collaboration, Abbasi, R., Ackermann, M., et al. 2021, arXiv e-prints, arXiv:2101.09836.  https://arxiv.org/abs/2101.09836
  * Kappes et al. (2009) Kappes, A., Halzen, F., & Murchadha, A. Ó. 2009, Nuclear Instruments and Methods in Physics Research A, 602, 117, doi: 10.1016/j.nima.2008.12.049
  * Kappes et al. (2007) Kappes, A., Hinton, J., Stegmann, C., & Aharonian, F. A. 2007, ApJ, 656, 870, doi: 10.1086/508936
  * Mandelartz & Becker Tjus (2015) Mandelartz, M., & Becker Tjus, J. 2015, Astroparticle Physics, 65, 80, doi: 10.1016/j.astropartphys.2014.12.002
  * McCollough et al. (2016) McCollough, M. L., Corrales, L., & Dunham, M. M. 2016, ApJ, 830, L36, doi: 10.3847/2041-8205/830/2/L36
  * Mel’Nik & Dambis (2009) Mel’Nik, A. M., & Dambis, A. K. 2009, MNRAS, 400, 518, doi: 10.1111/j.1365-2966.2009.15484.x
  * Moskalenko et al. (2006) Moskalenko, I. V., Porter, T. A., & Strong, A. W. 2006, ApJ, 640, L155, doi: 10.1086/503524
  * Neronov et al. (2014) Neronov, A., Semikoz, D., & Tchernin, C. 2014, Phys. Rev. D, 89, 103002, doi: 10.1103/PhysRevD.89.103002
  * Popescu et al. (2017) Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539, doi: 10.1093/mnras/stx1282
  * Safronov (2020) Safronov, G. 2020, arXiv e-prints, arXiv:2012.03373.  https://arxiv.org/abs/2012.03373
  * Sinitsyna & Sinitsyna (2020) Sinitsyna, V. G., & Sinitsyna, V. Y. 2020, in Journal of Physics Conference Series, Vol. 1468, Journal of Physics Conference Series, 012092, doi: 10.1088/1742-6596/1468/1/012092
  * Terrier et al. (2008) Terrier, R., Mattana, F., Djannati-Atai, A., et al. 2008, in American Institute of Physics Conference Series, Vol. 1085, American Institute of Physics Conference Series, ed. F. A. Aharonian, W. Hofmann, & F. Rieger, 312–315, doi: 10.1063/1.3076669
  * The LHAASO Collaboration et al. (2021) The LHAASO Collaboration, Cao, Z., Aharonian, F., et al. 2021, Science, 373, 425, doi: 10.1126/science.abg5137
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration, Amenomori, M., Bao, Y. W., et al. 2021, Nature Astronomy, doi: 10.1038/s41550-020-01294-9

LHAASO source | Neutrino Source | Distance | Angular Distance | Extension | $\Phi_{0,\gamma=2.0}^{90\%}$ | $\Phi_{0,\gamma=2.19}^{90\%}$ | $\Phi_{0,\gamma=2.50}^{90\%}$ | $\Phi_{0,\gamma=3.0}^{90\%}$ | Neutrino Source Search
---|---|---|---|---|---|---|---|---|---
| Candidate | [kpc] | [deg] | [deg] | $\rm[10^{-17}\,TeV^{-1}\,cm^{-2}\,s^{-1}]$ |
LHAASO J0534+2202 | Crab Nebula | 2.0a | 0.13 | PS | 3.7 | — | — | 0.76 | Ten-yeard
|  |  | 0.13 | PS | 4.74 | — | — | — | Eight-yeare
|  |  | 0.13 | PS | 10.7 | 11.8 | 10.7 | — | TeV PWNeg
|  |  | 0.10 | PS | — | — | 117 | — | SNR with PWNh
LHAASO J1825-1326 | HESS J1826-130 | 4.0a | 0.44 | PS | 13 | — | 25 | — | $\rm ANTARES\,\&\,IceCube^{f}$
| HESS J1825-137 | 3.9a | 0.39 | 0.461 | 10.7 | 11.8 | 10.7 | — | TeV PWNe
| 1HWC J1825-133 | 3.9a | 0.21 | 0.5 | — | — | 0.948 | — | HAWC Hotspotsh
LHAASO J1839-0545 | HESS J1841-055 | — | 0.52 | PS | 4.8 | — | — | 3.28 | Ten-year
| 1HWC J1838-060 | — | 0.29 | 0.5 | — | — | 0.948 | — | HAWC Hotspots
LHAASO J1843-0338 | HESS J1843-033 | — | 0.35 | PS | 2.5 | — | — | 1.09 | Ten-year
| 1HWC J1844-031c | — | 0.60 | 0.5 | — | — | 0.948 | — | HAWC Hotspots
| 1HWC J1842-046c | — | 0.98 | 0.5 | — | — | 0.948 | — | HAWC Hotspots
LHAASO J1849-0003 | HESS J1849-000 | 7a | 0.09 | PS | 2.2 | — | — | 1.01 | Ten-year
| IGR J18490-0000 | 7a | 0.11 | 0.09 | 10.7 | 11.8 | 10.7 | — | TeV PWNe
LHAASO J1908+0621 | MGRO J1908+06 | — | 0.21 | PS | 5.7 | — | — | 2.11 | Ten-year
|  |  | 0.10 | PS | 7.62 | — | — | — | Eight-year
|  |  | 0.49 | 1.3 | — | — | 3.98 | — | Milagro Sixh
| 1HWC J1907+062c | — | 0.29 | 0.5 | — | — | 0.948 | — | HAWC Hotspots
LHAASO J2018+3651 | MGRO J2019+37 | — | 0.09 | PS | 4.0 | — | — | 0.69 | Ten-year
|  |  | 0.38 | PS | 4.54 | — | — | — | Eight-year
|  |  | 0.08 | 0.75 | 10.7 | 11.8 | 10.7 | — | TeV PWNe
|  |  | 0.16 | 0.64 | — | — | 3.98 | — | Milagro Six
LHAASO J2032+4102 | 2HWC J2031+415 | 1.8a | 0.47 | PS | 9.2 | — | — | 1.42 | Ten-year
| Cyg OB2 | 1.5b | 0.18 | PS | 7.64 | — | — | — | Eight-year
| Cyg X-3 | 7.4c | 0.10 | PS | 8.20 | — | — | — | Eight-year
| MGRO J2031+41 | — | 0.39 | 1.5 | — | — | 3.98 | — | Milagro Six
| TeV J2032+4130 | 1.8a | 0.46 | 0.158 | 10.7 | 11.8 | 10.7 | — | TeV PWNe
LHAASO J2226+6057 | Boomerang | 0.8a | 0.30 | 0.22 | 10.7 | 11.8 | 10.7 | — | TeV PWNe
Table 1: Neutrino source candidates within one degree away from the measured
center of each LHAASO source. Column 3 shows the distances between the
candidates and the Earth, which are derived from a TeVCat; b Mel’Nik & Dambis
(2009); c McCollough et al. (2016). Column 4 shows the angular distances
between the candidates and LHAASO sources. Column 5 shows the extensions of
neutrino source candidates considered in the neutrino source searches. Column
6-9 show the 90% C.L. upper limits of neutrino flux parameterized as:
$dN_{\nu_{\mu}+\bar{\nu}_{\mu}}/dE_{\nu}=\Phi_{0,\gamma}^{90\%}\cdot(E_{\nu}/100\,{\rm
TeV})^{-\gamma}\times 10^{-17}\,{\rm TeV^{-1}~{}cm^{-2}~{}s^{-1}}$. The
coordinates, extensions, and neutrino flux of neutrino source candidates are
derived from SLSes: d the ten-year point source search (Aartsen et al.,
2020a), e the eight-year point source search (Aartsen et al., 2019), f the
ANTARES & IceCube combined search (Albert et al., 2020b) and SSSes: g the TeV
PWNe search (Aartsen et al., 2020b), h the Galactic catalog searches (Aartsen
et al., 2017a). The upper limits for $E^{-3}$ spectrum in the ten-year search
are taken from Figure 3 of Aartsen et al. (2020a). The data point of specific
source in Figure 3 can be identified with source declination and local pre-
trial p-value in the Table III of Aartsen et al. (2020a). The lower p-value
indicates the higher upper limit for the sources with similar declination
angles. As for TeV PWNe search, only the upper limits for the inverse age
weighting scheme are displayed.

##  Appendix A Weighting Schemes

The TeV PWNe search gave the upper limits of stacking neutrino flux from 35
identified TeV PWNe under four weighting schemes (Aartsen et al., 2020b).
Neutrino flux from each source in the catalog is assumed to be the same in the
equal weighting, proportional to gamma-ray flux at $\rm 1\,TeV$ in the gamma-
ray flux weighting, proportional to inverse age of pulsar in the inverse age
weighting and proportional to pulsar frequency in the frequency weighting.

The Galactic catalog searches gave the upper limit of stacking neutrino flux
over five Galactic catalogs respectively (Aartsen et al., 2017a). In the
Milagro catalog, neutrino flux from each source is assumed to follow the model
of Kappes et al. (2009). In the HAWC catalog, neutrino flux from each source
is assumed to be proportional to the gamma-ray flux at pivot energy. In the
other three SNR catalogs, neutrino flux from each source is assumed to be
equal.

##  Appendix B TeV Gamma-ray Observations

#### LHAASO J1825-1326

We take the gamma-ray spectrum of eHWC J1825-134 instead of that of 1HWC
J1825-133, because 1HWC J1825-133 is reported as 2HWC J1825-134 in the 2HWC
catalog (Abeysekara et al., 2017), while 2HWC J1825-134 is only 0.07 degree
from eHWC J1825-134 ($\rm>56\,TeV$) (Abeysekara et al., 2020). The flux data
points of HESS J1826-130 are raised by a factor of 2.37 from the original
values, if we consider a 2-d Gaussian intrinsic extension of $0.21^{\circ}$
with the integration region of $0.22^{\circ}$ (H. E. S. S. Collaboration et
al., 2020).

#### LHAASO J1839-0545

We only take the spectrum of HESS J1841-055, because 1HWC J1838-060 overlaps
with the extension of HESS J1841-055 and its differential flux normalization
is compatible with the previous measurements of HESS J1841-055 (Abeysekara et
al., 2016). The flux data points of HESS J1841-055 are raised by a factor of
1.14 based on its intrinsic extension of $0.41^{\circ}$ ($0.25^{\circ}$) along
the major (minor) axis and the integration region of $0.7^{\circ}$ (Aharonian
et al., 2008).

#### LHAASO J1843-0338

We take the observation of 2HWC J1844-032, because 1HWC J1844-031c is reported
as 2HWC J1844-032 in the 2HWC gamma-ray catalog (Abeysekara et al., 2017).

#### LHAASO J1849-0003

We only take the spectrum of HESS J1849-000, because HESS J1849-000 is in
spatial coincidence with IGR J18490-0000 (Terrier et al., 2008).

#### LHAASO J1908+0621

We take the spectrum of eHWC J1907+063. 1HWC J1907+062c is reported as 2HWC
J1908+063 (Abeysekara et al., 2017), while 2HWC J1908+063 is a point source
only 0.17 degree from eHWC J1907+063 ($>$56 TeV) (Abeysekara et al., 2020).
The flux data points of HESS J1908+063 are raised by a factor of 1.51 for its
intrinsic extension of $0.34^{\circ}$ and the integration region of
$1.3^{\circ}$(Aharonian et al., 2009).

#### LHAASO J2018+3651

We take the spectrum of VER J2019+368 because it constitutes the bulk of the
emission of MGRO J2019+37 (Abeysekara et al., 2018). The flux data points of
VER J2019+368 is raised by a factor of 2.53 based on its intrinsic extension
of $0.34^{\circ}$ ($0.14^{\circ}$) along the major (minor) axis and the
integration region of $0.23^{\circ}$. We also use the observation of eHWC
J2019+368 because its angular distance to LHAASO J2018+3651 is only 0.17
degree and it also has 100 TeV measurements Abeysekara et al. (2020).

#### LHAASO J2032+4102

The energy spectrum of TeV J2032+4130 is raised by different factors of 1.79,
4.64 and 2.72 for the flux data points below $\rm 1.46\,TeV$, between $\rm
1.46\,TeV$ and $\rm 2.76\,TeV$, and above $\rm 2.76\,TeV$ respectively,
according to its intrinsic extensions at different energy bands (Aliu et al.,
2014).

#### LHAASO J2226+6057

The VERITAS flux data points of G106.3+2.7 are raised by a factor of 1.59
based on its intrinsic extension of $0.27^{\circ}$ ($0.18^{\circ}$) along the
major (minor) axis with an integration region of $0.32^{\circ}$ (Acciari et
al., 2009).

For the NSCs not mentioned above, we simply take their gamma-ray measurements
without adjustments.

##  Appendix C Gamma-ray Flux from TeV PWNe

There are 6 PWNe associated with LHAASO sources. We fit the gamma-ray spectra
of these sources from TeV to sub-PeV. The best spectral fits are selected
among the single power law (PL) model
$dN/dE\propto(E/10\,\rm{TeV})^{-\gamma_{1}}$, the exponential cutoff power law
(ECPL) model $dN/dE\propto(E/10\,\rm{TeV})^{-\gamma_{1}}e^{-E/E_{c}}$ and the
log-parabola (LOGP) model
$dN/dE\propto(E/10\,\rm{TeV})^{-\gamma_{1}-\gamma_{2}\,\rm{log}(E/10\,\rm{TeV})}$
based on the corrected Akaike Information Criterion ($\rm AIC_{C}$) (HURVICH &
TSAI, 1991). Both the gamma-ray measurements of the identified PWNe and the
associating LHAASO sources are used in the fitting except for HESS J1825-137.
The observations by LHAASO are not considered for HESS J1825-137, because HESS
J1825-137 is just one of three sources composing the gamma-ray hotspot and the
gamma-ray flux at 100 TeV is not dominated by this PWN (Albert et al., 2021).
The best-fit parameters are shown in Table 2 (row 1-6).

There are 7 PWNe in the data taken region $-15^{\circ}<\delta<75^{\circ}$
($\delta$ is declination angle) but not discovered by KM2A half-array. The
upper limit on the gamma-ray flux is $1.1\times
10^{-13}\,\rm{TeV\,cm^{-2}\,s^{-1}}$ between 100 TeV and 178 TeV (Cao et al.,
2021). We fit the upper limit of energy flux at the 100-178 TeV bin with the
ECPL model, following the spectral index used in the TeV PWNe search (Aartsen
et al., 2020b). The best-fit cutoff energies are also the upper limits of
cutoff energies $E_{c}$. The best-fit $E_{c}$ is negative for SNR G054.1+00.3,
so we just follow its spectrum in the TeV PWNe search. The upper limits of
$E_{c}$ for the other 6 sources are shown in Table 2 (row 7-13).

The correction for the interstellar absorption is necessary for the 12 PWNe in
Table 2. The distances to 8 PWNe are available from TeVCat. As for the other 4
PWNe, we place them at the border of the ISRF model (a cylinder with radius 24
kpc and half-height 10 kpc) given by Popescu et al. (2017) to get the maximum
absorption to gamma-rays around 100 TeV. So we can get a conservative upper
limit on gamma-ray flux from these PWNe before absorption. The formula to
calculate gamma-ray opacity follows the equations in section 3 of Moskalenko
et al. (2006). This upper limit is also hold when the gamma-rays from LHAASO
sources are not from the PWNe, but the SNRs or YMCs associated (e.g. Aharonian
et al., 2019; Ge et al., 2020).

As for the other 22 PWNe outside the data taken region, we simply use the same
spectra in the TeV PWNe search (Aartsen et al., 2020b). The sum of gamma-ray
spectra obtained above is shown as the black solid line in Figure 3.

PWN | $\gamma_{1}$ | $\gamma_{2}$ | $E_{\rm c}$ | Model | Distance
---|---|---|---|---|---
|  |  | [TeV] |  | [kpc]
Crab Nebula | 2.86 | 0.20 | — | LOGP | 2.0
HESS J1825-137 | 2.42 | — | 31 | ECPL | 3.9
IGR J18490-0000 | 1.99 | — | — | PL | 7.0
MGRO J2019+37 | 2.2 | 0.83 | — | LOGP | 24.8∗
TeV J2032+4130 | 2.36 | 0.41 | — | LOGP | 1.8
Boomerang | 2.29 | 0.35 | — | LOGP | 0.8
CTA1 | 2.2 | — | 214 | ECPL | 1.4
Geminga | 2.23 | — | 74 | ECPL | 0.25
2HWC J0700+143 | 2.17 | — | 196 | ECPL | 17.2∗
HESS J1831-098 | 2.1 | — | 188 | ECPL | 30.9∗
HESS J1837-069 | 2.54 | — | 80 | ECPL | 6.6
MAGIC J1857.2+0263 | 2.2 | — | 119 | ECPL | 29.7∗
Table 2: The gamma-ray spectral parameters of TeV PWNe constrained by LHAASO’s
observation. Row 1-6: The best-fit spectral parameters for the PWNe associated
with LHAASO sources. Row 7-12: The upper limits of cutoff energy $E_{\rm c}$
for the PWNe in the data taken region but not discovered by LHAASO-KM2A. These
parameters are corrected for the interstellar absorption. ∗the distance to the
source at the border of ISRF given by Popescu et al. (2017).


20.00003ptV\equiv\tfrac{1}{4}\sum_{ijkl}V_{ijkl}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}.$
(82)

Additionally, three-body terms will be induced by the transformation. All
three species have the same single-particle energies:
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(0,1,20)$. We define the valence
space to consist of the lower two orbits for each flavor. The antisymmetrized
matrix elements of the perturbation $V$ are:

$\displaystyle v_{QQ}$
$\displaystyle=V_{1323}=V_{1331}=V_{1332}=V_{2331}=V_{2332}=V_{3132}$ (83)
$\displaystyle v_{PP}$ $\displaystyle=V_{1122}$ $\displaystyle v_{PQ}$
$\displaystyle=V_{2213}=V_{2231}=V_{1113}=V_{1131}.$

We take $v_{PP}$=8, $v_{QQ}$=$-8$, $v_{PQ}$=1. Essentially $v_{PP}$ and
$v_{QQ}$ mix configurations within the $P$ and $Q$ spaces, respectively,
leading to a collective $Q$ state coming down in energy as the interaction is
turned on, while a collective $P$ state is pushed up, and eventually the
states cross. The remaining term $v_{PQ}$, is initially the term we want to
suppress. It couples the $P$ and $Q$ states and makes the level crossing an
avoided crossing. To reduce somewhat the size of the problem, we exclude the
highest level for the third particle (call it the proton), reducing the three-
body Hilbert space to $3\times 3\times 2=18$ configurations. We take $v_{QQ}$
and $v_{PQ}$ to act only between neutrons, while $v_{PP}$ acts on all species.

The eigenstates of this problem may easily be found by forming the $18\times
18$ Hamiltonian matrix and diagonalizing. However, our aim here is to first
decouple the $P$ and $Q$ spaces, and then diagonalize within the decoupled
spaces.

We perform a non-perturbative decoupling in the 3-body Hilbert space using the
iterative method outlined in section 3.1.6 (this is essentially the approach
proposed by Suzuki [171] to deal with the intruder-state problem). The first
step is to construct the matrix $\mathbb{H}_{0}$, where the subscript denotes
iterations. Next, we form an anti-hermitian generator $\mathbb{G}_{n}$ which
is defined as

$\langle q|\mathbb{G}_{n}|p\rangle=\frac{\langle
q|\mathbb{H}_{n}|p\rangle}{\langle q|\mathbb{H}_{n}|q\rangle-\langle
p|\mathbb{H}_{n}|p\rangle}$ (84)

(here $p$ and $q$ label $A$-body configurations belonging to the $P$ and $Q$
spaces, respectively) and obtain the next iteration of $\mathbb{H}$ by the
Baker-Campbell-Hausdorff expansion

$\mathbb{H}_{n+1}=\mathbb{H}_{n}+\left[\mathbb{G}_{n},\mathbb{H}_{n}\right]+\tfrac{1}{2}\left[\mathbb{G}_{n},\left[\mathbb{G}_{n},\mathbb{H}_{n}\right]\right]+\ldots$
(85)

The nested commutators are evaluated until the norm of the last nested
commutator falls below $10^{-7}$. The iteration in $n$ is performed until the
norm of $\mathbb{G}_{n}$ falls below $10^{-6}$, at which point the $P$ and $Q$
spaces are decoupled.

We also perform an IMSRG decoupling171717Strictly speaking, because there is
no core, and no normal ordering is performed, there is no “medium” and so this
is really just an SRG calculation., using the flow equation formulation,
directly on the Fock space representation of the Hamiltonian (82). We use the
flow equation formulation because in the Magnus formulation for $x\gtrsim
0.5$, the Magnus operator $\Omega$ grows sufficiently large that the Baker-
Campbell-Hausdorff expansion does not converge. We perform an IMSRG(2)
calculation, discarding 3-body terms, and we also perform an IMSRG(3)
calculation, including the full three-body commutators, so the calculation is
exact for the three-body problem. The results are presented in Fig. 16 as a
function of the perturbation strength parameter $x$. The gray lines in panels
(a) and (b) are the results of the decoupling in the $A$-body space following
the iteration procedure in Eqs. (84) and (85). The purple and red lines and
symbols correspond to the IMSRG(2) and IMSRG(3) solution in Fock-space.

Figure 16: _Top:_ The eigenvalues after decoupling the $P$ and $Q$ spaces in
the IMSRG(2) approximation (the gray lines indicate the exact values),
_Middle:_ Eigenvalues after decoupling with the full IMSRG(3), _Bottom:_ The
cluster decomposition of the $P$-space component of the transformed
Hamiltonian, all as a function of the interaction strength parameter $x$.

We immediately make two observations. First, as shown in the top panel of Fig.
16, before the level crossings the IMSRG(2) eigenvalues are in agreement with
the exact ones, while after the level crossings they go astay. Second, as
shown in the bottom panel of Fig. 16, at the first level crossing near
$x\approx 0.5$ the size of the three-body contribution to the Hamiltonian
begins to grow rapidly. These two observations are clearly connected. We may
understand the error of the IMSRG(2) calculation straightforwardly as the
consequence of discarding important three-body terms. Alternatively, one can
compare the IMSRG(2) results with those presented in e.g. Fig. 1 of Ref.
[257], where a two-level system with an avoided crossing was studied in
perturbation theory. At low orders of perturbation theory, the avoided
crossing is not reproduced. With increasing orders in the perturbation series,
the levels slowly bend back to reproduce the avoided crossing. One must go to
very high orders (e.g. 50th in Ref [257]) to fully recover the avoided
crossing.

The generation of three-body terms can be understood by considering the flow
equation for a $P$-space three-body contribution:

$\frac{d}{ds}V_{111222}\sim\eta_{1123}V_{1322}-V_{1123}\eta_{1322}+\ldots$
(86)

The level crossing prevents the off-diagonal two-body matrix elements like
$V_{1322}$ from being rapidly suppressed, and keeps the door open for strength
to leak into the three-body sector via terms like (86).

This investigation of a toy problem illustrates how the effects which lead to
intruder configurations also cause problems with decoupling within the IMSRG
framework, even without small denominators, and even if the flow equation
converges. Presumably, such effects will also arise in the SMCC framework
because of the strong similarity between the two approaches. This remains an
open and important problem, and we hope that a better understanding of these
effects will lead to a solution in the near future.

##  7 Other developments

###  7.1 EFT for the shell model

As we mentioned in the introduction, the success of shell model phenomenology
strongly indicates that the shell model provides the relevant degrees of
freedom for nuclear structure. It is therefore tempting to formulate the shell
model as an effective theory with some scheme for systematic improvement. The
main difficulty is in identifying a separation of scales which one can use to
form an expansion.

20 years ago, Haxton et al. put forward the idea of formulating the shell
model as an effective theory [258, 259, 260, 261]. The method presented in
that work amounts to an effective theory for the NN interaction with the
harmonic oscillator basis serving as the regulator, and the Bloch-Horowitz
effective interaction (4) cast as an RG flow equation. Similar ideas have been
pursued in Refs. [262, 263, 264], and by the Oak Ridge group [265, 266]. This
type of approach is appealing because it is formulated in the harmonic
oscialltor basis and so yields an interaction well-suited to a number of
popular many-body methods.

Using the harmonic oscillator basis as a regulator is conceptually distinct
from formulating the standard shell model directly as an EFT. In the latter
case, one should use the shell model to define the degrees of freedom, and
write down all possible terms in the Hamiltonian consistent with the relevant
symmetries (parity, rotational invariance, charge, baryon number, etc.). Then
one should assign an importance to those terms based on some power counting. A
recent attempt more along these lines [267] employs a Weinberg chiral power
counting in a shell-model basis, modified by the Galilean invariance breaking
terms due to the presence of the core. While a rigorous basis for the use of
Weinberg’s power counting is still lacking– core excitations introduce a new
scale, for instance, that might make it more natural to treat the Fermi
momentum as a hard scale– very encouraging order-by-order convergence was
obtained.

Another possibility might be, a-la Landau-Migdal theory [268], to exploit the
similarities between the valence shell model philosophy and Landau’s Fermi
liquid theory for infinite systems. In modern parlance, the latter can be
viewed as an effective field theory for low-lying excitations (“quasi-
particles”) in the vicinity of the Fermi surface. As with any EFT, the
effective Hamiltonian of Landau’s theory incorporates the underlying
symmetries of the system, and the low-energy couplings– the Landau parameters–
can either be fixed by experiment or calculated microscopically based on the
underlying theory. Lending credence to the analogy with the valence shell
model, Shankar, Polchinski, and others have shown that Landau’s theory can be
understood as an IR fixed point of the RG as one integrates out modes away
from the Fermi surface [269, 270, 271]. Intriguingly, their analysis shows
that _(i)_ three- and higher-body quasi-particle interactions are irrelevant
in the sense of the RG, which might shed light as to why phenomenological
shell model interactions with 1- and 2-body terms are so effective, and _(ii)_
the natural small parameter is the ratio of excitation energy to the Fermi
energy, which could provide guidance for formulating an appropriate
powercounting for an EFT tailored to the shell model.

###  7.2 Uncertainty quantification

A major deficiency in essentially all shell model applications to date is the
absence of quantified theoretical uncertainties. This is no small thing; at a
fundamental level, a theoretical prediction without some confidence interval
cannot be falsified. In particular, since experimental binding energies and
excitation energies can routinely be measured at parts-per-million precision,
whether or not theory and experiment are compatible is entirely dependent on
the theoretical uncertainty.

In fact, there are various meanings one can assign to error bars. Perhaps the
most straightforward form of shell model uncertainty is the root-mean-squared
deviation from experiment (see Fig. 9). As mentioned above, the
phenomenological USDB interaction [57] has a very small rms deviation of 130
keV throughout the $sd$ shell, and this may be interpreted in the following
way. There exist a large number of states which can be interpreted as shell
model configurations mixed by the same universal residual interaction. In the
context of effective interaction theory these are the states which get mapped
to the $P$ space. We may then interpret levels where the deviation is much
larger than 130 keV as not belonging to the $P$ space. A clear illustration of
this can be found in Ref. [57], where the ground state energies of 29,30Ne,
30,31Na, and 31,32Mg have conspicuously large deviations—a clear signal that
these “island of inversion” nuclei have a different character than the others.

This notion of uncertainty has been recently made more quantitatively rigorous
by Yoshida et al. [272], who explored the various possible effective
interactions in the $p$ shell and obtained marginal distributions for each
parameter based on the deviation from experiment. These distributions were
then propagated to the calculated spectra, yielding an error bar indicating
the range of values that could be obtained by a universal $p$-shell
interaction. Comparison to experiment then yields a well-defined
quantification of how appropriate the naive shell model picture of a given
state is.

A different type of uncertainty is sought by _ab initio_ approaches employing
chiral interactions. There, one should quantify and propagate three sources of
uncertainty: _(i)_ contributions from truncated higher orders in the EFT
expansion, _(ii)_ the experimental uncertainty in data used to fit the
parameters of the EFT, and _(iii)_ uncertainty due to approximations made in
solving the many-body problem. In this case, because one has begun with the
most general Lagrangian compatible with the symmetries of the Standard Model,
the resulting theoretical error bar would then indicate compatibility with the
Standard Model. While we certainly expect such compatibility from the known
nuclear energy levels, this becomes very important for testing extensions to
the Standard Model through e.g. searches for neutinoless double beta decay
[273] or dark matter [274]. While uncertainty quantification is feasible for
quantum Monte Carlo or no-core shell model calculations, there is as yet no
rigorous means for uncertainty quantification of _ab initio_ shell model
effective interactions. This is an important avenue of future work.

###  7.3 Coupling to the Continuum

Throughout the previous sections we have briefly touched upon the need to
account for continuum effects in theoretical calculations. Of course, this
will be especially relevant as we seek to understand the structure for
increasingly neutron-rich nuclei. While a variety of approaches for coupling
the shell model and other many-body methods to the continuum exist (see, e.g.,
[275, 276, 277], as well as the reviews [278, 279]), methods based on the
Berggren basis [280] appear to be most suitable in the context of VS-IMSRG and
SMCC.

The Berggren basis adds resonant and scattering states to the single-particle
basis from which many-body states are constructed. In valence-space
configuration interaction calculations, one obtains what is colloquially known
as the Gamow Shell Model (GSM), which entails the large-scale diagonalization
of a complex symmetric Hamiltonian [278]. Applications of this method to
weakly bound nuclei have been quite successful [281, 282, 283, 284, 285, 286,
287, 288, 289], and there is a push to move from the commonly used
phenomenological interactions to fundamental ones [289]. The Berggren basis
has been used successfully in ground- and excited-state CC calculations (see,
e.g., [210, 128]), hence the inclusion in VS-IMSRG and SMCC is technically
straightforward. However, the proliferation of states due to the inclusion of
the continuum aggravates the intruder-state problem discussed in Sec. 6. If
this issue can be solved, it would allow us to properly account for the
continuum coupling in the derivation of effective interactions. For the time
being, work is underway to at least treat the impact of the continuum on the
dynamics of the valence particles via the GSM.

##  8 Concluding remarks

In this work, we have reviewed the current state of efforts to derive
effective interactions for the shell model from modern nuclear forces, with an
emphasis on the impact of RG and EFT ideas on our understanding of the shell
model itself. We have summarized the popular approaches, and discussed their
relations at a formal level. We emphasized the importance of three-body forces
in eliminating the need for phenomenological adjustments and presented the
ensemble normal ordering (ENO) approach to efficiently including three-body
effects. We presented highlights from recent applications of _ab initio_ shell
model calculations, and discussed the current challenges of low-lying
collective excitations and intruder states.

Before concluding, we would like to make some remarks and clarify some common
misconceptions about _ab initio_ valence space methods. [SUMMARY POINTS]

  * •

The shell model picture is inherently a low-momentum description of nuclear
structure. The basic assumption of the shell model is that nucleons are
(almost) independent particles moving in a mean field potential, and that
nuclear spectra can be explained by the mixing of a few valence configurations
above an inert core via a residual interaction. The bound mean-field solution
and weak (possibly even perturbative) residual interaction that are the
foundation of this intuitive picture can only be obtained if low and high
momenta are decoupled in the Hamiltonian [102, 116, 115].

Of course, nuclear observables — energies, radii, transition rates — must be
independent of the resolution scale at which a theory operates. In principle,
there is nothing that prevents one from microscopically constructing a valence
shell model Hamiltonian starting from a high-resolution description, e.g.,
using an input interaction with a highly repulsive core. However, not only
does such a choice make computations more difficult, but it complicates
interpretations as the resulting shell model wave functions bear little
resemblance to the exact ones, which contain sizable contributions from a vast
number of configurations and defy a simple interpretation. In contrast, the
exact wave functions of a low-resolution Hamiltonian at least qualitatively
resemble those that come out of the shell model diagonalization, providing a
simple and intuitive picture.

  * •

Approaches such as the VS-IMSRG, SMCC, or the $\hat{Q}$-box resummation, are
methods for solving the nuclear many-body problem, not new shell model
interactions. The approaches described in this article combine the derivation
of effective interactions with a shell model diagonalization. This should be
understood as an efficient alternative to a large-scale, full no-core
configuration interaction (e.g., NCSM) calculation that would yield exact
results for nuclear spectra, but is infeasible in most cases.

  * •

A careful comparison with experimental data or theoretical results _requires_
that the method and underlying nuclear interaction be specified. When
comparing two phenomenological shell model calculations to experimental data,
the interpretation is generally straightforward: the interaction that better
reproduces the data is the better interaction. With _ab initio_ approaches,
such an interpretation is no longer appropriate.

Disagreements between theory and experiment must be caused either by
deficiencies in the underlying nuclear interactions, or the approximations
employed in deriving the effective interaction (provided the shell model
calculation is done without further approximations of its own). It is
therefore _crucially important_ to specify both the method and the input
interaction when comparisons with experimental data or other theoretical
results are presented.

It should also be kept in mind that the major advantage of _ab initio_
approaches is the ability to systematically improve the precision of the
theoretical result by lifting approximations, or improving the input nuclear
interactions.

  * •

Three-body forces are inevitable and non-negligible in nuclear structure. So
long as we choose to use protons and neutrons as our active degrees of
freedom—excluding explicit Deltas, anti-nucleons, etc.—there will be “genuine”
(in the traditional language) many-body forces accounting for these integrated
out degrees of freedom. So long as we wish to work in a valence space, there
will be “effective” many-body forces accounting for excitations outside the
valence space. The relative importance of these many-body forces will depend
on the details of the implementation (scheme and scale). Indeed, an excellent
description of a local region of the chart can be obtained with a purely two-
body interaction, like USDB for the $sd$-shell [57]. But such an interaction
will only work locally, and will need modification (e.g., scaling of matrix
elements with mass — again think of USDB) in order to be used over a wider
range of nuclei. The theoretical evidence supports the expectation that three-
body forces are the underlying source of such ad hoc modifications (see Sec.
4).

  * •

The mass-dependence of modern effective interactions does not imply a loss of
predictive power. Effective interactions for different target nuclei are
derived from the _same_ two- plus three-nucleon force, and there are _no
parameter refits or phenomenological modifications_. Again, the derivation of
the effective interaction and subsequent shell model diagonalization are
merely an efficient alternative to a full no-core configuration interaction
calculation.

  * •

All observables, not just the Hamiltonian, must be treated consistently to
produce a true _ab initio_ result. As discussed in section 3, the effective
interaction corresponds to a similarity transformation of the original
Hamiltonian, and in order to perform a consistent calculation, _all operators_
must also be transformed.

As a simple example, consider computing the deuteron ground state by
generating an effective interaction for the $0s$ shell. By construction, the
energy obtained by a (trivial) diagonalization in the $0s$ shell would be
identical to the result from a diagonalization in the full space with the bare
Hamiltonian. Now, if one were to calculate the deuteron quadrupole moment
using the bare $E2$ operator with the valence space wave functions (which
would be pure $s$-wave), the result would be zero. Using a consistently
transformed $E2$ operator will give the full-space result (for an
illustration, see Ref. [123]).

One could argue that the bare operator gives the leading-order contribution to
the effective operator (cf. Eq. (77)). Trouble arises when the contribution of
the leading term is suppressed by a symmetry, or missing degrees of freedom.
Examples are the case for the $E2$ operator in the deuteron or in $sd$-shell
oxygen isotopes, which only have neutrons in the valence space. Along these
same lines, the use of phenomenological effective charges in conjunction with
an _ab initio_ shell model interaction should be considered inappropriate. In
general, it is difficult to make meaningful conclusions based on inconsistent
calculations.

  * •

The use of an inert core does not constitute an ad hoc approximation.
Calculations based on effective interction theory do not formally rely on an
assumption that excitations out of the core are “negligible”. Such excitations
are accounted for by the effective interaction. Certainly, there will be
states in the experimental spectrum that are not generated in the valence
space calculation even with a perfect effective interaction—these belong to
the excluded $Q$-space. However, those states that are generated will not be
improved by, e.g., allowing core excitations described by a schematic
interaction, as this would amount to double-counting. Instead, to include core
excitations explicitly, one should re-define the $P$ and $Q$ spaces and derive
a new effective interaction.

The shell model has been the primary intellectual and computational framework
for low energy nuclear structure for the past 70 years. While the
computational work has been largely phenomenological over that time, an
enormous amount of knowledge and intuition has been developed. At the same
time, perhaps no problem in nuclear structure has so stubbornly resisted a
satisfactory solution as the microscopic derivation of shell model
interactions. The general path has been more or less known for over half a
century, but it is only recently that the available computational power,
combined with a more systematic way of thinking about nuclear forces and the
many-body problem, has allowed a direct connection between the shell model,
the forces applicable to few-body scattering, and the underlying physics of
the Standard Model.

We are not quite yet in the promised land. While there are certainly many
details remaining to be worked out (including those mentioned here), and
several clear extensions to be made (continuum effects, reactions), there are
still two major hills to climb: a fully consistent and satisfactory power-
counting for the interaction, and a rigorous uncertainty quantification for
our many-body methods. We hope that progress can be made on these fronts in
the near future, enabling a broadly applicable, quantitatively predictive
theory of nuclear structure.

## ACKNOWLEDGMENTS

We would like to thank B. Alex Brown, Takayuki Miyagi, Titus Morris, Petr
Navrátil and Zhonghao Sun for helpful discussions, and Gaute Hagen for
providing coupled cluster results. S.R.S is supported by the U.S. DOE under
contract DE-FG02-97ER41014. H.H. acknowledges support by the National Science
Foundation under Grant No. PHY-1614130, as well as the U.S. Department of
Energy, Office of Science, Office of Nuclear Physics under Grants No. DE-
SC0017887 and DE-SC0018083 (NUCLEI SciDAC Collaboration). S.K.B. acknowledges
support by the National Science Foundation under Grant No. PHY-1713901, as
well as as the U.S. Department of Energy, Office of Science, Office of Nuclear
Physics under Grant No. DE-SC0018083 (NUCLEI SciDAC Collaboration).

## APPENDIX

##  Appendix A IMSRG flow equations

For reference, we now present the IMSRG(2)/VS-IMSRG(2) flow equations [167,
145, 146]. Ground-state and valence-space decoupling only differ by the choice
of the generator $\eta$ (see Refs. [168, 146]).

The system of flow equations for the zero-, one-, and two-body parts of $H(s)$
result from evaluating

$\frac{dH}{ds}=[\eta(s),H(s)]$ (87)

with normal-ordered Fock-space operators that are truncated at the two-body
level181818These expressions can be easily adapted to evaluate the nested
commutators appearing in the Magnus formulation of the IMSRG.:

$H=E_{0}+\sum_{ij}f_{ij}\\{a^{\dagger}_{i}a_{j}\\}+\frac{1}{4}\sum_{ijkl}\Gamma_{ijkl}\\{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}\\}$
(88)
$\eta=\sum_{ij}\eta_{ij}\\{a^{\dagger}_{i}a_{j}\\}+\frac{1}{4}\sum_{ijkl}\eta_{ijkl}\\{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}\\}.$
(89)

The flow equations are then

$\displaystyle\frac{dE_{0}}{ds}$
$\displaystyle=\sum_{ab}n_{a}\bar{n}_{b}(\eta_{ab}f_{ba}-f_{ab}\eta_{ba})+\frac{1}{4}\sum_{abcd}n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}(\eta_{abcd}\Gamma_{cdab}-\Gamma_{abcd}\eta_{cdab})$
(90) $\displaystyle\frac{df_{ij}}{ds}$
$\displaystyle=\sum_{a}(\eta_{ia}f_{aj}-f_{ia}\eta_{aj})+\sum_{ab}(n_{a}-n_{b})(\eta_{ab}\Gamma_{biaj}-f_{ab}\eta_{biaj})$
$\displaystyle~{}~{}+\frac{1}{2}\sum_{abc}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})\left(\eta_{ciab}\Gamma_{abcj}-\Gamma_{ciab}\eta_{abcj}\right)$
(91) $\displaystyle\frac{d\Gamma_{ijkl}}{ds}$
$\displaystyle=\sum_{a}(1-P_{ij})(\eta_{ia}\Gamma_{ajkl}-f_{ia}\eta_{ajkl})-(1-P_{kl})(\eta_{ak}\Gamma_{ijal}-f_{ak}\eta_{ijal})$
$\displaystyle~{}~{}+\frac{1}{2}\sum_{ab}(n_{a}n_{b}-\bar{n}_{a}\bar{n}_{b})(\eta_{ijab}\Gamma_{abkl}-\Gamma_{ijab}\eta_{abkl})$
$\displaystyle~{}~{}-\sum_{ab}(n_{a}-n_{b})(1-P_{ij})(1-P_{kl})\eta_{bjal}\Gamma_{aibk}$
(92)

where $P_{ij}$ exchanges indices $i$ and $j$, $n_{a}$ is the occupation of
orbit $a$ and ${\bar{n}_{a}\equiv 1-n_{a}}$.

##  Appendix B Canonical perturbation theory to second order

We partition the Hamiltonian into a zeroth-order piece and a perturbation,

$H=H_{0}+xV\,,$ (93)

and we consider a perturbative expansion in powers of the dimensionless order
parameter $x$, where in the end we will take $x=1$. We further distinguish
between “diagonal” and “off-diagonal” components, $V=V^{d}+V^{od}$, where
“off-diagonal” generically means the terms we wish to suppress by the
transformation

$H_{\mathrm{eff}}=e^{\mathcal{G}}He^{-\mathcal{G}}\,.$ (94)

As in the discussion leading to Eq. (39), we use the superoperator notation to
express a commutator with $H_{0}$ in terms of an energy denominator $\Delta$.
Through second order in $x$, we obtain for $\mathcal{G}$

$\displaystyle\mathcal{G}^{[]1]}$ $\displaystyle=\frac{V^{od}}{\Delta}$ (95)
$\displaystyle\mathcal{G}^{[2]}$
$\displaystyle=\left[\mathcal{G}^{[1]},V^{d}\right]^{od}/\Delta+\tfrac{1}{2}\left[\mathcal{G}^{[1]},V^{od}\right]^{od}/\Delta$

The second term in $\mathcal{G}^{[2]}$ will vanish for $A$-body Hilbert space
formulations, but not in general for a Fock space formulation. This is related
to the different meanings of “off-diagonal” in the two formulations. The
transformed Hamiltonian through second order is

$\displaystyle H_{\mathrm{eff}}^{[0]}$ $\displaystyle=H_{0}$ (96)
$\displaystyle H_{\mathrm{eff}}^{[1]}$ $\displaystyle=V^{d}$ $\displaystyle
H_{\mathrm{eff}}^{[2]}$
$\displaystyle=\left[\mathcal{G}^{[1]},V^{d}\right]^{d}+\tfrac{1}{2}\left[\mathcal{G}^{[1]},V^{od}\right]^{d}$

##  Appendix C Integration of Magnus flow equation to second order

Here, we integrate the flow equation (28) explicitly to second order. As in
appendix B, we partition the Hamiltonian into a zero order piece $H_{0}$, and
a perturbation $V$, and we split up the perturbation into “diagonal” and “off-
diagonal” pieces. We use the White generator, which we write as

$\eta(s)\equiv\frac{H^{od}(s)}{\Delta}$ (97)

using the super-operator notation introduced in Sec. 2.3. To first order in
$x$, the flow equation for $\Omega$ is

$\frac{d\Omega^{[1]}}{ds}=\eta^{[1]}(s)=\frac{H^{[1]od}(s)}{\Delta}=\frac{V^{od}}{\Delta}+[\Omega^{[1]}(s),H_{0}]^{od}/\Delta=\frac{V^{od}}{\Delta}-\Omega^{[1]}(s).$
(98)

A differential equation for $\Omega^{[2]}(s)$ may be obtained in a similar
manner. The solutions given the initial condition $\Omega(0)=0$ are

$\displaystyle\Omega^{[1]}(s)$ $\displaystyle=(1-e^{-s})\frac{V^{od}}{\Delta}$
(99) $\displaystyle\Omega^{[2]}(s)$
$\displaystyle=(1-e^{-s}-se^{-s})\left[\frac{V^{od}}{\Delta},V^{d}\right]^{od}/\Delta+\tfrac{1}{2}(1-e^{-2s})\left[\frac{V^{od}}{\Delta},V^{od}\right]^{od}/\Delta.$

The flowing Hamiltonian through second order is

$\displaystyle H^{[0]}$ $\displaystyle=H_{0}$ (100) $\displaystyle H^{[1]}$
$\displaystyle=V^{d}+e^{-s}V^{od}$ $\displaystyle H^{[2]}$
$\displaystyle=(1-e^{-s})\left[\frac{V^{od}}{\Delta},V^{d}\right]^{d}+\tfrac{1}{2}(1-e^{-2s})\left[\frac{V^{od}}{\Delta},V^{od}\right]^{d}$
$\displaystyle~{}~{}+se^{-s}\left[\frac{V^{od}}{\Delta},V^{d}\right]^{od}+e^{-s}(1-e^{-s})\left[\frac{V^{od}}{\Delta},V^{od}\right]^{od}.$

We see that at first order, the off-diagonal part of the perturbation is
exponentially suppressed. At $s=0$, the second order piece is by definition
zero. As $s$ increases, we initially induce both diagonal and off-diagonal
second-order terms, and eventually, the induced second-order terms are
suppressed exponentially, leaving a purely diagonal the second-order
correction.

Taking the limit $s\rightarrow\infty$, we obtain the same generator and
effective Hamiltonian as the canonical perturbation theory in appendix B. Note
that this equivalence requires the same definition of “off-diagonal”, and that
the results will differ at higher orders if we include the commutator terms on
the right hand side of the Magnus flow equation in (28).

## References

  * [1] Mayer MG.  Phys. Rev. 74:235 (1948)
  * [2] Haxel O, Jensen JHD, Suess HE.  Phys. Rev. :1766 (1949)
  * [3] Goeppert Mayer M, Jensen JHD.  Wiley (1955)
  * [4] Dawson JF, Talmi I, Walecka JD.  Ann. Phys. (N. Y). 18:339 (1962)
  * [5] Bertsch GF.  Nucl. Phys. 74:234 (1965)
  * [6] Kuo TTS, Brown GE.  Nucl. Phys. 85:40 (1966)
  * [7] Kuo TT.  Nucl. Phys. A 90:199 (1967)
  * [8] Kuo TTS, Brown GE.  Nuclear Physics A 114:241 (1968)
  * [9] Cohen S, Kurath D.  Nuclear Physics 73:1 (1965)
  * [10] Wildenthal BH.  Prog. Part. Nucl. Phys. 11:5 (1984)
  * [11] Caurier E, et al.  Phys. Rev. C 59:2033 (1999)
  * [12] DeShalit A, Talmi I, De-Shalit A, Talmi I.  New York: Academic Press (1963)
  * [13] Ring P, Schuck P.  Berlin, Heidelberg: Springer-Verlag, 1st ed. (1980)
  * [14] Heyde K.  Berlin, Heidelberg: Springer (1990)
  * [15] Talmi I.  Contemp. Concepts Phys. Vol 7. Chur, Switzerland: Harwood Academic Publishers (1993)
  * [16] Brown BA.  Prog. Part. Nucl. Phys. 47:517 (2001)
  * [17] Talmi I.  In Advances in Nuclear Physics, Volume 27, eds. JW Negele, EW Vogt, vol. 27. Springer US, 1–275 (2003)
  * [18] Caurier E, et al.  Rev. Mod. Phys. 77:427 (2005)
  * [19] Johnson CW, Ormand WE, Krastev PG.  Comput. Phys. Commun. 184:2761 (2013)
  * [20] Brown B, Rae W.  Nuclear Data Sheets 120:115 (2014)
  * [21] Sorlin O, Porquet MG.  Prog. Part. Nucl. Phys. 61:602 (2008)
  * [22] Otsuka T, et al. arXiv:1805.06501 [nucl-th] (2018)
  * [23] Lindgren I, Morrison J.  New York: Springer Verlag (1986)
  * [24] Kuo TTS, Osnes E.  vol. 364 of Lecture Notes in Physics.  Springer Berlin Heidelberg (1990)
  * [25] Shavitt I, Bartlett RJ.  New York: Cambridge University Press (2009)
  * [26] Kuo TTS.  Annu. Rev. Nucl. Sci. 24:101 (1974)
  * [27] Ellis PJ, Osnes E.  Rev. Mod. Phys. 49:777 (1977)
  * [28] Hjorth-Jensen M, Kuo TTS, Osnes E.  Phys. Rep. 261:125 (1995)
  * [29] Dean DJ, et al.  Prog. Part. Nucl. Phys. 53:419 (2004)
  * [30] Coraggio L, et al.  Prog. Part. Nucl. Phys. 62:135 (2009)
  * [31] Talmi I.  Helv. Phys. Acta 25:185 (1952)
  * [32] Lawson RD, Uretsky JL.  Phys. Rev. 106:1369 (1957)
  * [33] Talmi L, Unna I.  Annual Review of Nuclear Science 10:353 (1960)
  * [34] Talmi I, Unna I.  Nuclear Physics 30:280 (1962)
  * [35] Hamada T, Johnston ID.  Nuclear Physics 34:382 (1962)
  * [36] Reid RV.  Annals of Physics 50:411 (1968)
  * [37] Brueckner KA, Levinson CA, Mahmoud HM.  Phys. Rev. 95:217 (1954)
  * [38] Brueckner KA, Levinson CA.  Phys. Rev. 97:1344 (1955)
  * [39] Day BD.  Rev. Mod. Phys. 39:719 (1967)
  * [40] Barrett BR, Kirson MW.  Nuclear Physics A 148:145 (1970)
  * [41] Goode P, Kirson MW.  Physics Letters B 51:221 (1974)
  * [42] Kirson MW.  Annals of Physics 66:624 (1971)
  * [43] Barrett BR.  Physics Letters B 38:371 (1972)
  * [44] Kirson MW.  Annals of Physics 82:345 (1974)
  * [45] Vary JP, Sauer PU, Wong CW.  Phys. Rev. C 7:1776 (1973)
  * [46] Schucan T, Weidenmüller H.  Ann. Phys. (N. Y). 73:108 (1972)
  * [47] Schucan TH, Weidenmüller HA.  Ann. Phys. (N. Y). 76:483 (1973)
  * [48] Barrett BR.  Lect. Notes Phys. Vol. 40. Berlin: Springer Verlag (1975)
  * [49] Kuo TTS, Lee SY, Ratcliff KF.  Nucl. Physics, Sect. A 176:65 (1971)
  * [50] Krenciglowa E, Kuo T.  Nucl. Phys. A 240:195 (1975)
  * [51] Whitehead RR, Watt A, Cole BJ, Morrison I.  In Adv. Nucl. Phys. Boston, MA: Springer US, 123–176 (1977)
  * [52] Chung W. 1976.  Empirical renormalization of shell-model hamiltonians and magnetic dipole moments of sd-shell nuclei.  Phd, Michigan State University
  * [53] Poves A, Zuker A.  Phys. Rep. 70:235 (1981)
  * [54] Brown BA, Wildenthal BH.  Annu. Rev. Nucl. Part. Sci. 38:29 (1988)
  * [55] Caurier E, Zuker AP, Poves A, Martinez-Pinedo G.  Phys. Rev. C 50:225 (1994)
  * [56] Brown BA, Richter WA, Julies RE, Wildenthal BH.  Ann. Phys. (N. Y). 182:191 (1988)
  * [57] Brown BA, Richter WA.  Phys. Rev. C 74:034315 (2006)
  * [58] Pieper SC, Wiringa RB.  Ann. Rev. Nucl. Part. Sci. 51:53 (2001)
  * [59] Carlson J, et al.  Rev. Mod. Phys. 87:1067 (2015)
  * [60] Lynn JE, Tews I, Gandolfi S, Lovato A arXiv:1901.04868 [nucl-th] (2019)
  * [61] Navrátil P, Barrett BR.  Phys. Rev. C 54:2986 (1996)
  * [62] Navrátil P, Kamuntavicius GP, Barrett BR.  Phys. Rev. C 61:044001 (2000)
  * [63] Barrett BR, Navrátil P, Vary JP.  Prog. Part. Nucl. Phys. 69:131 (2013)
  * [64] Kümmel H, Lührmann KH, Zabolitzky JG.  Physics Reports 36:1 (1978)
  * [65] Zabolitzky JG, Ey W.  Nuclear Physics A 328:507 (1979)
  * [66] Emrich K, Zabolitzky JG.  Nuclear Physics A 351:439 (1981)
  * [67] Zheng DC, et al.  Phys. Rev. C 48:1083 (1993)
  * [68] Pudliner BS, Pandharipande VR, Carlson J, Wiringa RB.  Phys. Rev. Lett. 74:4396 (1995)
  * [69] Coester F, Cohen S, Day B, Vincent CM.  Phys. Rev. C 1:769 (1970)
  * [70] Bethe HA.  Ann. Rev. Nucl. Sci. 21:93 (1971)
  * [71] Weinberg S.  Phys. Lett. B 251:288 (1990)
  * [72] Weinberg S.  Nucl. Phys. B 363:3 (1991)
  * [73] Epelbaum E, Hammer HW, Meißner UG.  Rev. Mod. Phys. 81:1773 (2009)
  * [74] Machleidt R, Entem D.  Phys. Rept. 503:1 (2011)
  * [75] Machleidt R, Sammarruca F.  Physica Scripta 91:083007 (2016)
  * [76] Pastore S, et al.  Phys. Rev. C 80:034004 (2009)
  * [77] Pastore S, Girlanda L, Schiavilla R, Viviani M.  Phys. Rev. C 84:024001 (2011)
  * [78] Piarulli M, et al.  Phys. Rev. C 87:014006 (2013)
  * [79] Kölling S, Epelbaum E, Krebs H, Meißner UG.  Phys. Rev. C 80:045502 (2009)
  * [80] Kölling S, Epelbaum E, Krebs H, Meißner UG.  Phys. Rev. C 84:054008 (2011)
  * [81] Entem DR, Machleidt R.  Phys. Lett. B 524:93 (2002)
  * [82] Entem DR, Kaiser N, Machleidt R, Nosyk Y.  Phys. Rev. C 91:014002 (2015)
  * [83] Entem DR, Kaiser N, Machleidt R, Nosyk Y.  Phys. Rev. C 92:064001 (2015)
  * [84] Epelbaum E, Krebs H, Meißner UG.  Phys. Rev. Lett. 115:122301 (2015)
  * [85] Reinert P, Krebs H, Epelbaum E.  The European Physical Journal A 54:86 (2018)
  * [86] Ekström A, et al.  Phys. Rev. C 91:051301 (2015)
  * [87] Gezerlis A, et al.  Phys. Rev. C 90:054323 (2014)
  * [88] Lynn JE, et al.  Phys. Rev. Lett. 116:062501 (2016)
  * [89] Piarulli M, et al.  Phys. Rev. C 91:024003 (2015)
  * [90] Piarulli M, et al.  Phys. Rev. Lett. 120:052503 (2018)
  * [91] Navrátil P.  Few-Body Syst. 41:117 (2007)
  * [92] Gazit D, Quaglioni S, Navrátil P.  Phys. Rev. Lett. 103:102502 (2009)
  * [93] Wilson K.  Phys. Rep. 12:75 (1974)
  * [94] Wilson KG.  Rev. Mod. Phys. 47:773 (1975)
  * [95] Polchinski J.  Nucl. Phys. B 231:269 (1984)
  * [96] Birse MC, Mcgovern JA, Richardson KG.  Phys. Lett. B 464:169 (1999)
  * [97] Bogner SK, Schwenk A, Kuo TTS, Brown GE.  arxiv nucl-th/0111042 (2001)
  * [98] Barford T, Birse MC.  Phys. Rev. C 67:064006 (2002)
  * [99] Bogner SK, Kuo TTS, Schwenk a.  Phys. Rep. 386:1 (2003)
  * [100] Bogner SK, Furnstahl RJ, Perry RJ.  Phys. Rev. C 75:061001 (2007)
  * [101] Bogner S, Furnstahl R, Perry R, Schwenk A.  Phys. Lett. B 649:488 (2007)
  * [102] Bogner S, Furnstahl R, Schwenk A.  Prog. Part. Nucl. Phys. 65:94 (2010)
  * [103] Furnstahl RJ, Hebeler K.  Rept. Prog. Phys. 76:126301 (2013)
  * [104] Polyzou W, Glöckle W.  Few-Body Systems 9:97 (1990)
  * [105] Bethe HA.  Phys. Rev. 103:1353 (1956)
  * [106] Goldstone J.  Proc. R. Soc. A Math. Phys. Eng. Sci. 239:267 (1957)
  * [107] Bethe HA, Brandow BH, Petschek AG.  Phys. Rev. 129:225 (1963)
  * [108] Phillips AC, Barton G.  Phys. Lett. B 28:378 (1969)
  * [109] Tjon JA.  Phys. Lett. 56:217 (1975)
  * [110] Rustgi ML, et al.  Phys. Rev. C 4:854 (1971)
  * [111] Abzouzi A, Caurier E, Zuker AP.  Phys. Rev. Lett. 66:1134 (1991)
  * [112] Glazek SD, Wilson KG.  Phys. Rev. D 48:5863 (1993)
  * [113] Wegner F.  Ann. Phys. (Leipzig) 3:77 (1994)
  * [114] Bogner SK, Furnstahl RJ, Ramanan S, Schwenk A.  Nucl. Phys. A 773:203 (2006)
  * [115] Hoppe J, et al.  Phys. Rev. C 96:054002 (2017)
  * [116] Tichai A, Langhammer J, Binder S, Roth R.  Physics Letters B 756:283 (2016)
  * [117] Jurgenson E, Navrátil P, Furnstahl RJ.  Phys. Rev. Lett. 103:082501 (2009)
  * [118] Hebeler K.  Phys. Rev. C 85:021002 (2012)
  * [119] Wendt KA.  Phys. Rev. C 87:061001 (2013)
  * [120] Roth R, Calci A, Langhammer J, Binder S.  Phys. Rev. C 90:024325 (2014)
  * [121] Anderson ER, Bogner SK, Furnstahl RJ, Perry RJ.  Phys. Rev. C 82:054001 (2010)
  * [122] Schuster MD, et al.  Phys. Rev. C 90:011301 (2014)
  * [123] Parzuchowski NM, et al.  Phys. Rev. C 96:034324 (2017)
  * [124] Jurgenson ED, et al.  Phys. Rev. C 87:054312 (2013)
  * [125] Hergert H, et al.  Phys. Rev. Lett. 110:242501 (2013)
  * [126] Binder S, Langhammer J, Calci A, Roth R.  Phys. Lett. B 736:119 (2014)
  * [127] Hagen G, Papenbrock T, Hjorth-Jensen M, Dean DJ.  Rept. Prog. Phys. 77:096302 (2014)
  * [128] Hagen G, Jansen GR, Papenbrock T.  Phys. Rev. Lett. 117:172501 (2016)
  * [129] Morris TD, et al.  Phys. Rev. Lett. 120:152503 (2018)
  * [130] Coraggio L, et al.  Phys. Rev. C 68:034320 (2003)
  * [131] Holt JD, et al.  Phys. Rev. C 72:041304 (2005)
  * [132] Otsuka T, et al.  Phys. Rev. Lett. 105:032501 (2010)
  * [133] Holt JD, Otsuka T, Schwenk A, Suzuki T.  J. Phys. G 39:085111 (2012)
  * [134] Holt JD, Menéndez J, Schwenk A.  J. Phys. G 40:075105 (2013)
  * [135] Holt J, Menéndez J, Schwenk A.  Eur. Phys. J. A 49:1 (2013)
  * [136] Holt JD, Menéndez J, Schwenk A.  Phys. Rev. Lett. 110:022502 (2013)
  * [137] Holt JD, Menéndez J, Simonis J, Schwenk A.  Phys. Rev. C 90:024312 (2014)
  * [138] Tsunoda N, Takayanagi K, Hjorth-Jensen M, Otsuka T.  Phys. Rev. C 89:024313 (2014)
  * [139] Bogner SK, et al.  Phys. Rev. Lett. 113:142501 (2014)
  * [140] Jansen GR, et al.  Phys. Rev. Lett. 113:142502 (2014)
  * [141] Simonis J, et al.  Phys. Rev. C 93:011302 (2016)
  * [142] Stroberg SR, et al.  Phys. Rev. C 93:051301(R) (2016)
  * [143] Stroberg SR, et al.  Phys. Rev. Lett. 118:032502 (2017)
  * [144] Tsunoda N, et al.  Phys. Rev. C 95:021304 (2017)
  * [145] Hergert H, et al.  Phys. Rep. 621:165 (2016)
  * [146] Hergert H.  Phys. Scr. 92:023002 (2017)
  * [147] Tichai A, Müller J, Vobig K, Roth R :1 arXiv:1809.07571 (2018)
  * [148] Forssén C, et al.  Phys. Rev. C 97:1 (2018)
  * [149] Lanczos C.  J. Res. Nat’l Bur. Std. 45:255 (1950)
  * [150] Arnoldi WE.  Quarterly of Applied Mathematics 9:17 (1951)
  * [151] Davidson ER.  Comp. Phys. Comm. 53:49 (1989)
  * [152] Friman B, Schwenk A arXiv:1101.4858 (2011).  Three-Body Interactions in Fermi Systems, chap. 6.  World Scientific, 141–156
  * [153] Bogner SK, Furnstahl RJ, Perry RJ.  Ann. Phys. (N. Y). 323:1478 (2008)
  * [154] Bloch C, Horowitz J.  Nucl. Phys. 8:91 (1958)
  * [155] Feshbach H.  Ann. Phys. (N. Y). 5:357 (1958)
  * [156] Brandow BH.  Rev. Mod. Phys. 39:771 (1967)
  * [157] Suzuki K, Lee S.  Prog. Theor. Phys. 64:2091 (1980)
  * [158] Shurpin J, Müther H, Kuo TTS, Faessler A.  Nucl. Physics, Sect. A 293:61 (1977)
  * [159] Shurpin J, et al.  Phys. Lett. B 69:395 (1977)
  * [160] Takayanagi K.  Nucl. Phys. A 864:91 (2011)
  * [161] Takayanagi K.  Nucl. Phys. A 852:61 (2011)
  * [162] Navrátil P, Thoresen M, Barrett BR.  Phys. Rev. C 55:R573 (1997)
  * [163] Lisetskiy AF, et al.  Phys. Rev. C 78:044302 (2008)
  * [164] Lisetskiy AF, et al.  Phys. Rev. C 80:024315 (2009)
  * [165] Dikmen E, et al.  Phys. Rev. C 91:064301 (2015)
  * [166] Jansen GR, et al.  Phys. Rev. C 94:011301(R) (2016)
  * [167] Tsukiyama K, Bogner SK, Schwenk A.  Phys. Rev. Lett. 106:222502 (2011)
  * [168] Tsukiyama K, Bogner SK, Schwenk A.  Phys. Rev. C 85:061304 (2012)
  * [169] White SR.  J. Chem. Phys. 117:7472 (2002)
  * [170] Klein DJ.  J. Chem. Phys. 61:786 (1974)
  * [171] Suzuki K.  Prog. Theor. Phys. 58:1064 (1977)
  * [172] Wick GC.  Phys. Rev. 80:268 (1950)
  * [173] Kutzelnigg W, Mukherjee D.  J. Chem. Phys. 107:432 (1997)
  * [174] Morris TD, Parzuchowski NM, Bogner SK.  Phys. Rev. C 92:034331 (2015)
  * [175] Morris TD. 2016.  A Novel Approach to Solving the In-Medium SRG.  Ph.D. thesis, Michigan State University
  * [176] Yanai T, Chan GKL.  J. Chem. Phys. 127:104107 (2007)
  * [177] Taube AG, Bartlett RJ.  Int. J. Quantum Chem. 106:3393 (2006)
  * [178] Shavitt I, Redmon LT.  J. Chem. Phys. 73:5711 (1980)
  * [179] Sun ZH, et al.  Phys. Rev. C 98:054320 (2018)
  * [180] Kutzelnigg W.  J. Chem. Phys. 77:3081 (1982)
  * [181] Kutzelnigg W, Koch S.  J. Chem. Phys. 79:4315 (1983)
  * [182] Suzuki K, Okamoto R.  Prog. Theor. Phys. 70:439 (1983)
  * [183] Suzuki K, Okamoto R.  Prog. Theor. Phys. 71:1221 (1984)
  * [184] Lindgren I.  J. Phys. B At. Mol. Phys. 7:2441 (1974)
  * [185] Suzuki K.  Prog. Theor. Phys. 68:1627 (1982)
  * [186] Miyagi T, Abe T, Okamoto R, Otsuka T.  Phys. Rev. C 96:1 (2017)
  * [187] Primas H.  Rev. Mod. Phys. 35:710 (1963)
  * [188] Watson TJ, Chan GKL.  Journal of Chemical Theory and Computation 12:512 (2016)
  * [189] Brandow BH.  In Effective Interactions and Operators in Nuclei. Lecture Notes in Physics vol 40, ed. BR Barrett. Springer Berlin Heidelberg, 1–24 (1975)
  * [190] Hergert H, et al.  Phys. Rev. C 87:034307 (2013)
  * [191] Somà V, et al.  Phys. Rev. C 89:061301 (2014)
  * [192] Jansen GR.  Phys. Rev. C 88:024305 (2013)
  * [193] Primakoff H, Holstein T.  Phys. Rev. 55:1218 (1939)
  * [194] Wentzel G.  Helv. Phys. Acta 15:111 (1942)
  * [195] Drell SD, Huang K.  Phys. Rev. 91:1527 (1953)
  * [196] Skyrme THR.  Nucl. Phys. 9:615 (1959)
  * [197] Vautherin D, Brink DM.  Phys. Rev. C 5:626 (1972)
  * [198] Dechargé J, Gogny D.  Physical Review C 21:1568 (1980)
  * [199] Polls A, et al.  Nucl. Phys. A 401:124 (1983)
  * [200] Rath PK, et al.  Nucl. Phys. A 427:511 (1984)
  * [201] Osnes E.  Phys. Lett. B 26:274 (1968)
  * [202] Goldhammer P, Hill JR, Nachamkin J.  Nucl. Phys. A 106:62 (1968)
  * [203] Singh B.  Nucl. Phys. A 219:621 (1973)
  * [204] Caesar C, Simonis J, Others.  Phys. Rev. C 88:034313 (2013)
  * [205] Zuker A.  Phys. Rev. Lett. 90:042502 (2003)
  * [206] Bogner SK, et al.  Phys. Lett. Sect. B Nucl. Elem. Part. High-Energy Phys. 576:265 (2003)
  * [207] Pieper SC, Varga K, Wiringa RB.  Phys. Rev. C 66:044310 (2002)
  * [208] Navrátil P, et al.  Phys. Rev. Lett. 99:042501 (2007)
  * [209] Cipollone A, Barbieri C, Navrátil P.  Phys. Rev. Lett. 111:062501 (2013)
  * [210] Hagen G, et al.  Phys. Rev. Lett. 108:242501 (2012)
  * [211] Nomura M.  Phys. Lett. B 85:187 (1979)
  * [212] Hergert H, et al.  Phys. Rev. C 90:041302 (2014)
  * [213] Cipollone A, Barbieri C, Navrátil P.  Phys. Rev. C 92:014306 (2015)
  * [214] Perez-Martin S, Robledo LM.  Phys. Rev. C 76:064314 (2007)
  * [215] Perez-Martin S, Robledo LM.  Phys. Rev. C 78:014304 (2008)
  * [216] Matsubara T.  Prog. Theor. Phys. 14:351 (1955)
  * [217] Gaudin M.  Nucl. Phys. 15:89 (1960)
  * [218] Fetter AL, Walecka JD.  Mineola: Dover Publications, Inc (2003)
  * [219] Cole BJ, Watt A, Whitehead RR.  J. Phys. G-Nuclear Part. Phys. 1:213 (1975)
  * [220] Hebeler K, et al.  Phys. Rev. C 83:031301 (2011)
  * [221] Brown BA, Signoracci A, Hjorth-Jensen M.  Physics Letters, Section B: Nuclear, Elementary Particle and High-Energy Physics 695:507 (2011)
  * [222] Dirim H, Elliott JP, Evans JA.  Nucl. Phys. A 244:301 (1975)
  * [223] Roth R, et al.  Phys. Rev. Lett. 107:072501 (2011)
  * [224] Fukui T, et al.  Phys. Rev. C 98:044305 (2018)
  * [225] Leach KG, et al.  Phys. Rev. C 94:011304 (2016)
  * [226] Brodeur M, et al.  Phys. Rev. C 96:034316 (2017)
  * [227] Crawford HL, et al.  Phys. Rev. C 95:064317 (2017)
  * [228] Garnsworthy AB, et al.  Phys. Rev. C 96:044329 (2017)
  * [229] Reiter MP, et al.  Phys. Rev. C 96:052501 (2017)
  * [230] Steppenbeck D, et al.  Phys. Rev. C 96:064310 (2017)
  * [231] Henderson J, et al.  Phys. Lett. Sect. B Nucl. Elem. Part. High-Energy Phys. 782:468 (2018)
  * [232] Izzo C, et al.  Phys. Rev. C 97:014309 (2018)
  * [233] Leistenschneider E, et al.  Phys. Rev. Lett. 120:062503 (2018)
  * [234] Liu HN, et al. arXiv:1811.08451 [nucl-ex] (2018)
  * [235] Michimasa S, et al.  Phys. Rev. Lett. 121:022506 (2018)
  * [236] Mougeot M, et al.  Phys. Rev. Lett. 120:232501 (2018)
  * [237] Reiter MP, et al.  Phys. Rev. C 98:024310 (2018)
  * [238] Randhawa JS, et al.  Phys. Rev. C 99:021301 (2019)
  * [239] Wienholtz F, et al.  Nature 498:346 (2013)
  * [240] Lapoux V, et al.  Phys. Rev. Lett. 117:1 (2016)
  * [241] Otsuka T, et al.  Prog. Part. Nucl. Phys. 47:319 (2001)
  * [242] Stumpf C, Braun J, Roth R.  Phys. Rev. C 93:021301 (2016)
  * [243] Nndc
  * [244] Warburton E, Becker J, Brown BA.  Phys. Rev. C 41 (1990)
  * [245] Caurier E, Nowacki F, Poves A.  Phys. Rev. C 90:014302 (2014)
  * [246] Caurier E, Langanke K, Martínez-Pinedo G.  Phys. Lett. B 522:240 (2001)
  * [247] Garcia Ruiz RF, et al.  Nat. Phys. 12:594 (2016)
  * [248] Horie H, Arima A.  Phys. Rev. 623:778 (1955)
  * [249] Siegel S, Zamick L.  Phys. Lett. B 28:450 (1969)
  * [250] Siegel S, Zamick L.  Phys. Lett. B 28:453 (1969)
  * [251] Siegel S, Zamick L.  Nucl. Phys. A 145:89 (1970)
  * [252] Angeli I, Marinova K.  At. Data Nucl. Data Tables 99:69 (2013)
  * [253] Stone NJ.  At. Data Nucl. Data Tables 111-112:1 (2016)
  * [254] Wildenthal BH, Chung W.  Phys. Rev. C 22:2260 (1980)
  * [255] Watt A, Singhal RP, Storm MH, Whitehead RR.  Journal of Physics G: Nuclear Physics 7:L145 (1981)
  * [256] Hu B.  private communication
  * [257] Ellis PJ, et al.  Nucl. Phys. A 573:216 (1994)
  * [258] Haxton WC, Song CL.  Phys. Rev. Lett. 84:5484 (2000)
  * [259] Haxton WW, Luu T.  Nucl. Phys. A 690:1 (2001)
  * [260] McElvain KS, Haxton WC arXiv:1607.06863 [nucl-th] (2016)
  * [261] McElvain KS, Haxton WC arXiv:1902.03543 [nucl-th] (2019)
  * [262] Stetcu I, Barrett B, van Kolck U.  Phys. Lett. B 653:358 (2007)
  * [263] Stetcu I, Rotureau J, Barrett BR, van Kolck U.  J. Phys. G 37:064033 (2010)
  * [264] Yang CJ.  Phys. Rev. C 94:064004 (2016)
  * [265] Binder S, et al.  Phys. Rev. C 93:044332 (2015)
  * [266] Bansal A, et al.  Phys. Rev. C 98:054301 (2018)
  * [267] Huth L, Durant V, Simonis J, Schwenk A.  Phys. Rev. C 98:044301 (2018)
  * [268] Migdal AB.  Interscience (Wiley) (1967)
  * [269] Shankar R.  Rev. Mod. Phys. 66:129 (1994)
  * [270] Polchinski J. 1992.  In Proceedings, Theoretical Advanced Study Institute (TASI 92): From Black Holes and Strings to Particles: Boulder, USA, June 1-26, 1992
  * [271] Chen T, Frohlich J, Seifert M. 1995.  In NATO Advanced Study Institute: Les Houches Summer School, Session 62: Fluctuating Geometries in Statistical Mechanics and Field Theory Les Houches, France, August 2-September 9, 1994
  * [272] Yoshida S, Shimizu N, Togashi T, Otsuka T.  Phys. Rev. C 98:061301 (2018)
  * [273] Engel J, Menéndez J.  Reports on Progress in Physics 80:046301 (2017)
  * [274] Menéndez J, Gazit D, Schwenk A.  Physical Review D 86:103511 (2012)
  * [275] Volya A, Zelevinsky V.  Phys. Rev. C 74:064314 (2006)
  * [276] Papadimitriou G, Vary JP.  Physics Letters B 746:121 (2015)
  * [277] Baroni S, Navrátil P, Quaglioni S.  Phys. Rev. C 87:034326 (2013)
  * [278] Michel N, Nazarewicz W, Płoszajczak M, Vertse T.  J. Phys. G 36:013101 (2009)
  * [279] Myo T, Kikuchi Y, Masui H, Katō K.  Prog. Part. Nucl. Phys. 79:1 (2014)
  * [280] Berggren T.  Nucl. Phys. A 109:265 (1968)
  * [281] Betan RI, Nazarewicz W.  Phys. Rev. C 86:034338 (2012)
  * [282] Jaganathen Y, Michel N, Płoszajczak M.  Phys. Rev. C 89:034624 (2014)
  * [283] Fossez K, et al.  Phys. Rev. C 91:034609 (2015)
  * [284] Fossez K, et al.  Phys. Rev. C 94:054302 (2016)
  * [285] Fossez K, et al.  Phys. Rev. C 93:011305 (2016)
  * [286] Fossez K, Rotureau J, Michel N, Nazarewicz W.  Phys. Rev. C 96:024308 (2017)
  * [287] Fossez K, Rotureau J, Michel N, Płoszajczak M.  Phys. Rev. Lett. 119:032501 (2017)
  * [288] Jones MD, et al.  Phys. Rev. C 96:054322 (2017)
  * [289] Fossez K, Rotureau J, Nazarewicz W.  Phys. Rev. C 98:061302 (2018)