# Vacancy-Engineered Flat-Band Superconductivity in Holey Graphene

Matheus S. M. de Sousa  Department of Physics, PUC-Rio, 22451-900 Rio de
Janeiro, Brazil  Fujun Liu  Instituto de Física, Universidade de Brasília,
Brasília-DF, Brazil  Fanyao Qu  Instituto de Física, Universidade de Brasília,
Brasília-DF, Brazil  Wei Chen  Department of Physics, PUC-Rio, 22451-900 Rio
de Janeiro, Brazil

###### Abstract

A bipartite lattice with chiral symmetry is known to host zero energy flat
bands if the numbers of the two sublattices are different. We demonstrate that
this mechanism of producing flat bands can be realized on graphene by
introducing periodic vacancies. Using first-principle calculations, we
elaborate that even though the pristine graphene does not exactly preserve
chiral symmetry, this mechanism applied to holey graphene still produces
single or multiple bands as narrow as $\sim 0.5$eV near the Fermi surface
throughout the entire Brillouin zone. Moreover, this mechanism can combine
with vacancy-engineered nonsymmorphic symmetry to produce band structures with
coexisting flat bands and nodal lines. A weak coupling mean-field treatment
suggests the stabilization of superconductivity by these vacancy-engineered
narrow bands. In addition, superconductivity occurs predominantly on the
majority sublattices, with an amplitude that increases with the number of
narrow bands.

##  I Introduction

The superconductivity (SC) discovered recently in twisted bilayer graphene
(TBLG) revises the interest on the search for SC on graphene based
materialsCao  _et al._ (2018a). While the origin of the observed SC is still
under intense debate, the small flat band due to splitting and anticrossing of
the Dirac cones from the two graphene sheets is generally considered to play
an important roleSuárez Morell  _et al._ (2010); Bistritzer and MacDonald
(2011); Moon and Koshino (2012); Trambly de Laissardière  _et al._ (2012);
Lopes dos Santos  _et al._ (2012); Fang and Kaxiras (2016); Cao  _et al._
(2018b). In fact, the possibility of flat band induced SC has long been of
great interest, since the CuO2 plane of high temperature superconductor
materialsBednorz and Müller (1986); Lee  _et al._ (2006); Keimer  _et al._
(2015), at which the SC occursAnderson (1987); Emery (1987); Zhang and Rice
(1988), has the same structure (if Cu and O are not distinguished) as the Lieb
lattice that is known to host a flat band and has been realized in various
other systemsShen  _et al._ (2010); Guzmán-Silva  _et al._ (2014); Mukherjee
_et al._ (2015); Taie  _et al._ (2015); Vicencio  _et al._ (2015); Julku  _et
al._ (2016); Xia  _et al._ (2016); Diebel  _et al._ (2016); Slot  _et al._
(2017); Klembt  _et al._ (2017); Cui  _et al._ (2020). On the other hand, from
the density of states (DOS) point of view, it is intriguing to ask whether
there exists some generic mechanisms that can generate flat bands in a large
area of the Brillouin zone (BZ) of single-layer graphene at low energy, such
that SC may be stabilized.

In this paper, we elaborate that the a generic mechanism of producing zero-
energy flat bands (ZEFBs) on any bipartite lattice, first put forward by
LiebLieb (1989), can be realized on graphene by introducing periodic
vacancies. Lieb’s theorem states that on a bipartite lattice with chiral
(sublattice) symmetry, ZEFBs occur if the numbers of the two sublattices per
unit cell are different $N_{A}\neq N_{B}$, and the ZEFBs are at least
$|N_{A}-N_{B}|$-fold degenerate. We demonstrate that this situation can be
created on the honeycomb lattice by removing different numbers of the two
sublattices in an enlarged unit cell, and ZEFBs throughout the entire BZ occur
provided the tight-binding model of the lattice preserves the chiral symmetry.
This mechanism is then tested in a realistic single-layer graphene by first-
principle calculations. Our results indicate that despite graphene in reality
does not exactly preserve chiral symmetry, this mechanism can still produce
bands as narrow as $\sim 0.5$eV near the Fermi surface throughout the entire
BZ, and moreover can be used to engineer an exotic band structure that
contains both flat bands and nodal lines.

Our theoretical investigation is largely motivated by the fact that graphene
with vacancies, often called holey graphene or graphene nanomesh, has been
realized by various experimental techniques, such as nitrogenationPawlak  _et
al._ (2020); Mahmood  _et al._ (2015), self-aligned anisotropic etchingShi
_et al._ (2011), nano-network maskingJung  _et al._ (2014), and lithography
using copolymerBai  _et al._ (2010), nanosphereWang  _et al._ (2013), or He
ion beamArchanjo  _et al._ (2014), suggesting the feasibility of vacancy
engineering in reality. Note that our proposal considers complete removal of
carbon atoms, which is in contrast to the flat bands produced by the periodic
potentials from adatomsSkurativska  _et al._ (2021), requires much less
removal of atoms compared to the cyclicgraphyneYou  _et al._ (2019) or
aziteShima and Aoki (1993) proposals, and may also be realizable in
superlattices nanolithographed in semiconductor thin filmsTadjine  _et al._
(2016). In addition, since these vacancy-engineered narrow bands dramatically
enlarge the DOS at the Fermi surface compared to the pristine graphene, we
examine the possibility of phonon-mediated SC by means of a weak coupling
mean-field theory using the realistic phonon band widthSaito  _et al._ (2001);
Maultzsch  _et al._ (2004); Mohr  _et al._ (2007); Jorio  _et al._ (2011). The
results indeed point to the occurrence of SC in these vacancy configurations,
whose spatial pattern is highly influenced by the chiral symmetric wave
function of the ZEFBs, with a magnitude that increases with the number of
ZEFBs.

The structure of the paper is organized in the following manner. In Sec. II,
we first revisit Lieb’s theorem with an emphasis on the chiral symmetry of the
ZEFB wave functions. Two vacancy configurations are then used to demonstrate
perfect ZEFBs on a honeycomb lattice that preserves chiral symmetry. The
narrow bands of these two configurations realized in graphene are then
elaborated by first-principle calculations, and additionally another vacancy
configuration that yields coexisting narrow bands and nodal lines. In Sec.
III, we lay out the weak-coupling mean field theory to investigate SC on the
first two vacancy configurations, especially to detail their spatial pattern
and dependence on the pairing interaction and degeneracy of the ZEFBs. These
results are summarized in Sec. IV.

##  II Vacancy-engineered flat bands on graphene

###  II.1 Chiral symmetry and rank-nullity theorem

We first revisit Lieb’s theorem that is based on the rank-nullity theoremLieb
(1989), with a special emphasis on the nonspatial symmetries, localization of
the wave functions, and applications to periodic vacancies. We consider any
two- (2D) or three-dimensional (3D) bipartite-lattices described by single-
particle Hamiltonian $H({\bf k})$ in momentum space that preserves time-
reversal (TR), particle-hole (PH) chiral symmetries

$\displaystyle TH({\bf k})T^{-1}=H(-{\bf k}).$ $\displaystyle CH({\bf
k})C^{-1}=-H(-{\bf k}).$ $\displaystyle SH({\bf k})S^{-1}=-H({\bf k})$ (1)

which are nonspatial symmetries particularly relevant to topological
orderSchnyder  _et al._ (2008); Ryu  _et al._ (2010); Chiu  _et al._ (2016);
von Gersdorff  _et al._ (2021) and topological phase transitionsChen and
Schnyder (2019). In these bipartite lattices, the Hamiltonian matrix arranged
in the basis of the electron operators of the two sublattices $(c_{B{\bf
k}},c_{A{\bf k}})$ is a block-off-diagonal $2\times 2$ matrix, and the
symmetry operators are implemented by $T=K$, $C=\sigma_{3}K$ and
$S=\sigma_{3}$, where $K$ is the complex conjugation operator.

Now suppose we enlarge the unit cell to contain not 2 but $N=$ even number of
sites with the same amount of two sublattices, then the Hamiltonian matrix
arranged in the basis $(c_{B_{1}{\bf k}},...,c_{B_{N/2}{\bf k}},c_{A_{1}{\bf
k}},...,c_{A_{N/2}{\bf k}})$ remains block-off-diagonal

$\displaystyle H({\bf k})=\left(\begin{array}[]{cc}&Q({\bf k})\\\
Q^{{\dagger}}({\bf k})&\\\ \end{array}\right),$ (4)

where $Q({\bf k})$ is an $(N/2)\times(N/2)$ square matrix. The symmetry
operators are implemented by $C=\sigma_{3}\otimes I_{N/2}K$ and
$S=\sigma_{3}\otimes I_{N/2}$ in this case, with $I_{N/2}$ the
$(N/2)\times(N/2)$ identity matrix. If we now introduce periodic vacancies
into the lattice, the columns and rows in the $H({\bf k})$ in Eq. (4) that
correspond to the vacancy sites will be removed. It is then obvious that if
the number of vacancies on the $A$ and $B$ sublattices are different, then the
unit cell will contain different number of sublattices $N_{A}\neq N_{B}$. As a
result, the $Q({\bf k})$ in Eq. (4) as an $N_{A}\times N_{B}$ matrix is not a
square matrix anymore, as elaborated in Fig. 1. Nevertheless, the PH and
chiral symmetries of the systems still hold since removing the columns and
rows in $C=\sigma_{3}\otimes I_{N/2}K$ and $S=\sigma_{3}\otimes I_{N/2}$ that
correspond to the vacancy sites preserve Eq. (1). As a result, the band
structure at any vacancy configuration is PH symmetric, and the ZEFB wave
functions must be localized on one of the two sublattices since they are
eigenstates of the chiral operator $S$. In fact, we will prove below that the
ZEFB wave functions must localize in the majority sublattices.

Figure 1: Schematics of the proposed generic mechanism for ZEFBs. Starting
from a chiral symmetric bipartite lattice, the Hamiltonian matrix $H({\bf k})$
describing an enlarged unit cell containing $N_{A}=N_{B}=N/2$ sites of each
sublattice is block-off-diagonal (blue blocks). After periodic vacancies are
introduced, which correspond to removing some columns and rows in the
Hamiltonian, the off-diagonal block becomes not square if the unit cell
contains different numbers of the two sublattices $N_{B}>N_{A}$. In this case,
the rank-nullity theorem ensures that ZEFBs must occur.

Denoting $r(M)$ as the rank and $\eta(M)$ as the nullity of a matrix $M$, our
interest is how the nullity of the Hamiltonian $\eta(H)$, which counts the
number of ZEFBs, can be nonzero. Without loss of generality, we assume that
the vacancy configuration is such that the $B$ sublattices are the majority
$N_{B}>N_{A}$. In this case, the rank-nullity theorem states that

$\displaystyle r(Q)+\eta(Q)=N_{A},\;\;\;r(H)+\eta(H)=N_{A}+N_{B}.$ (5)

For $H$ of the form of Eq. (4), the rank satisfies

$\displaystyle r(Q)=r(Q^{{\dagger}})=r(QQ^{{\dagger}})=r(Q^{{\dagger}}Q),$
$\displaystyle r(H)=r(Q)+r(Q^{{\dagger}})=2r(Q).$ (6)

Using these simple identities in linear algebra, we now prove the following
propositions.

Proposition 1: $\eta(H)>0$ if $Q$ is not square. This can be proved easily
from Eqs. (5) and (6)

$\displaystyle\eta(H)=N_{A}+N_{B}-2r(Q)=N_{B}-N_{A}+2\eta(Q).$ (7)

Hence $\eta(H)>0$ if $N_{B}>N_{A}$, meaning that ZEFB must emerge if we remove
the two sublattices in different quantities.

Proposition 2: $\eta(H)=N_{B}-N_{A}$ if $\eta(Q)=0$. To prove this, we start
from

$\displaystyle
r(QQ^{{\dagger}})+\eta(QQ^{{\dagger}})=N_{B},\;\;\;r(Q^{{\dagger}}Q)+\eta(Q^{{\dagger}}Q)=N_{A},\;\;\;\;\;$
(8)

which implies

$\displaystyle r(QQ^{{\dagger}})-r(Q^{{\dagger}}Q)=N_{B}-N_{A}.$ (9)

If $Q$ itself is not singular $\eta(Q)=0$, which is true in many practical
examples, then Eq. (5) implies
$r(Q)=N_{A}=r(QQ^{{\dagger}})=r(Q^{{\dagger}}Q)$. Then from Eq. (8) one sees
that $\eta(QQ^{{\dagger}})=N_{B}-N_{A}$ and $\eta(Q^{{\dagger}}Q)=0$. Because
the square of the Hamiltonian $H^{2}={\rm
diag}(QQ^{{\dagger}},Q^{{\dagger}}Q)$ has the same nullity and ZEFBs wave
functions as $H$, we immediately see that
$\eta(H)=\eta(H^{2})=\eta(QQ^{{\dagger}})=N_{B}-N_{A}$. The propositions 1 and
2 constitute the original version of Lieb’s theorem.

Proposition 3: ZEFB wave functions are localized in the majority sublattices
if $\eta(Q)=0$. This is simply because the proof for proposition 2 shows that
the ZEFB wave functions are given by diagonalizing the first block
$QQ^{{\dagger}}$ in $H^{2}$ that belongs to the majority $B$ sublattices,
whereas wave functions in $A$ sublattices are zero.

Figure 2: Numerical results for the two vacancy configurations C15 (left
column) and C14 (right column), whose lattice structures are shown in (a) with
the two sublattices $A$ and $B$ indicated. (b) The band structures obtained
from the nearest-neighbor tight-binding model that preserves chiral symmetry,
which contains a single ZEFB for C15 and doubly degenerate ZEFBs for C14. (c)
DFT band structures and (d) the corresponding DOS for these two vacancy
configurations realized by graphene that does not exactly preserve the chiral
symmetry, which still yield very narrow bands of mainly $p_{z}$ at low energy.

We aim to examine these propositions in the spinless honeycomb lattice with
nearest-neighbor hopping

$\displaystyle H=\sum_{\langle ij\rangle}t\,c_{i}^{{\dagger}}c_{j}+U\sum_{i\in
v}c_{i}^{{\dagger}}c_{i}$ (10)

owing to its relevance to the $p_{z}$ orbital of single-layer graphene that
will be addressed later, where $c_{i}$ denotes the electron annihilation
operator at site $i$, and $U$ is the large on-site potential that can be used
to conveniently project out the vacancy sites $i\in v$. Increasing from $U=0$
to $U>100t$ can continuously change the band structure from Dirac points to
ZEFBs, similar to that proposed recently in the Lieb-kagome latticesLim  _et
al._ (2020), but we will focus on the large on-site potential regime $U>100t$
that completely removes the vacancy sites. We denote these vacancy
configurations by C${}_{N_{A}+N_{B}}$. In the left panel of Fig. 2 (a), we
show a C15 example of removing a single $A$ sublattice from a $N=16$
rectangular unit cell such that $N_{A}=7$ and $N_{B}=8$. In this case, the BZ
is rectangular, and the PH symmetric tight-binding band structure plotted
along a high-symmetry line (HSL) clearly indicates a single ZEFB throughout
the BZ, as shown in the left panel of Fig. 2 (b). In contrast, the C14
configuration shown in the right panel of Fig. 2 (a) and (b) that removes two
$A$ sublattices on the same 16-site unit cell, such that $N_{A}=6$ and
$N_{B}=8$, has doubly degenerate ZEFBs $\eta(H)=N_{B}-N_{A}=2$, consistent
with proposition 2. In Appendix A, we also elaborate proposition 3 by
presenting the wave functions of the ZEFBs for C15 and C14 and show that
indeed both are localized on the majority $B$ sublattices. Finally, we remark
that although we focus on periodic vacancies on an infinite graphene, the
propositions 1 to 3 also explain the number of zero eigenenergies and their
wave functions of a finite size graphene with random vacanciesBouzerar and
Mayou (2021).

###  II.2 Application to realistic graphene

To examine the proposed mechanism on the realistic single-layer graphene, it
should be first noted that graphene in reality does not exactly preserve the
PH and chiral symmetries in Eq. (1) owing to the complications such as longer
range hopping and hybridization between different orbitals, although
magnitudes of these factors are much smaller than the nearest-neighbor
hoppingCastro Neto  _et al._ (2009). To investigate the effect of these
symmetry breaking factors, we turn to density functional theory (DFT) to
obtain the band structure of graphene with periodic vacancies. In Fig. 2 (c),
we show the DFT band structures and DOS for C15 and C14, which indicate that
very narrow bands do occur near the Fermi surface. Although not perfectly
flat, these bands are as narrow as $\sim 0.5$eV, reminisce the ZEFBs.
Moreover, the DOS at the Fermi surface is dramatically enhanced by these
narrow bands compared to the pristine graphene, as indicated by Fig. 2 (d).

###  II.3 Coexistence of ZEFBs and nodal lines

We proceed to demonstrate that the proposed mechanism can coexist with another
vacancy engineering principle proposed recently, namely the nodal-line
semimetals caused by 2D nonsymmorphic vacancy configurationsLiu  _et al._
(2021, ); de Sousa  _et al._ . In these nonsymmorphic configurations, every
two carbon atoms map to each other under glide plane operation, resulting in
nodal lines at the BZ boundary regardless the details of the Hamiltonian,
which serves as a vacancy-engineering principle to obtain symmetry-enforced
nodal linesYoung and Kane (2015); Yamakage  _et al._ (2016); Zhao and Schnyder
(2016); Wieder  _et al._ (2018). As an example, in Fig. 3 we show a C22
configuration that contains two vacancies on the $A$ sublattice, and
additionally a glide plane along ${\hat{\bf x}}$ direction. The resulting DFT
band structure indicates that both mechanisms prevail in this situation,
yielding a band structure that contains two low energy narrow bands that
reminisce doubly degenerate ZEFBs, and in addition every two pairs of bands
(each pair is spin degenerate) stick together at the BZ boundary $k_{x}=0$
(the $X-V$ section of Fig. 3 (b)) to form symmetry-enforced nodal lines as
predicted. This example indicates that vacancy engineering can combine
different crystalline and nonspatial symmetries to produce very exotic band
structures that may not be easily found in nature.

Figure 3: (a) A vacancy configuration C22 that contains two missing $A$
sublattices and a glide plane $G_{x}$. (b) The resulting band structure
contains two low energy narrow bands of mainly $p_{z}$ orbital origin, and in
addition every two pairs of spin degenerate bands stick together to form four-
fold degenerate nodal lines at the $X-V$ section of the BZ boundary.

##  III Mean-field treatment of SC stabilized by ZEFBs

###  III.1 Bogoliubov-de Gennes equations for graphene with periodic vacancy

The enlarged DOS near the Fermi surface due to the narrow bands is expected to
dramatically change the electronic and magnetic properties of the holey
graphene compared to the pristine graphene. However, these narrow bands are
neither perfectly flat nor entirely located at zero energy, and hence it is
unclear at present whether strong correlations will play an important role on
the electronic properties in the same way as in TBLGCao  _et al._ (2018a). On
the other hand, the phonon band width in graphene is about $\omega_{D}\approx
0.25$eVSaito  _et al._ (2001); Maultzsch  _et al._ (2004); Mohr  _et al._
(2007); Jorio  _et al._ (2011), which means that the narrow bands with band
width $\sim 0.5$eV in a large area of the BZ are within the Debye frequency,
suggesting a large phase space for phonon-mediated Cooper pairingKopnin  _et
al._ (2011); Heikkilä and Volovik (2016). These features motivate us to
examine whether a conventional phonon-mediated $s$-wave SC phase emerges in
the holey graphene with low energy narrow bands. For this purpose, we resort
to the following spinful mean-field model of $s$-wave SC

$\displaystyle H$ $\displaystyle=$ $\displaystyle\sum_{\langle
ij\rangle\sigma}t\,c_{i\sigma}^{{\dagger}}c_{j\sigma}+\sum_{\langle\langle
ij\rangle\rangle\sigma}t^{\prime}\,c_{i\sigma}^{{\dagger}}c_{j\sigma}-\sum_{i\sigma}\mu\,c_{i\sigma}^{{\dagger}}c_{i\sigma}$
(11)
$\displaystyle+\sum_{i}\left(\Delta_{i}c_{i\uparrow}^{{\dagger}}c_{i\downarrow}^{{\dagger}}+\Delta_{i}^{\ast}c_{i\downarrow}c_{i\uparrow}\right)+U\sum_{i\in
v}c_{i\sigma}^{{\dagger}}c_{i\sigma},$

where $t=2.8$eV is the nearest-neighbor hopping on the honeycomb lattice, and
the on-site potential $U>100t$ is used to project out the vacancy sites $i\in
v$. We use the next-nearest-neighbor hopping $t^{\prime}$ and chemical
potential $\mu$ to simulate the breaking of chiral symmetry in realistic
graphene, and find that the values $t^{\prime}=-0.2$eV and $\mu=0.2$eV can
give a reasonable fit to the narrow bands obtained by DFT in both C15 and C14,
as demonstrated in the Appendix A. The mean field Hamiltonian is diagonalized
into $H={\rm const.}+\sum_{\bf k\alpha}E_{\bf k}\gamma_{\bf
k\alpha}^{{\dagger}}\gamma_{\bf k\alpha}$ by a Bogoliubov transformation

$\displaystyle c_{i\uparrow}=\sum_{\bf k}\gamma_{\bf k\uparrow}u_{\bf
k}(i)-\gamma_{\bf k\downarrow}^{{\dagger}}v_{\bf k}^{\ast}(i),$ $\displaystyle
c_{i\downarrow}=\sum_{\bf k}\gamma_{\bf k\downarrow}u_{\bf k}(i)+\gamma_{\bf
k\uparrow}^{{\dagger}}v_{\bf k}^{\ast}(i),$ (12)

where $i=1,2...N_{A}+N_{B}$ denotes the site inside a unit cell, and
$\gamma_{\bf k\sigma}$ is the annihilation operator of the Bogoliubov
quasiparticles. The wave functions $\left\\{u_{\bf k}(i),v_{\bf
k}(i)\right\\}$ and eigenenergy $E_{\bf k}$ satisfy

$\displaystyle E_{\bf k}u_{\bf k}(i)$ $\displaystyle=$
$\displaystyle\sum_{\langle ij\rangle}t\,u_{\bf k}(j)+\sum_{\langle\langle
ij\rangle\rangle}t^{\prime}\,u_{\bf k}(j)$ $\displaystyle+U\delta_{i\in
v}u_{\bf k}(i)+\Delta_{i}v_{\bf k}(i),$ $\displaystyle E_{\bf k}v_{\bf k}(i)$
$\displaystyle=$ $\displaystyle-\sum_{\langle ij\rangle}t\,v_{\bf
k}(j)-\sum_{\langle\langle ij\rangle\rangle}t^{\prime}\,v_{\bf k}(j)$ (13)
$\displaystyle-U\delta_{i\in v}v_{\bf k}(i)+\Delta_{i}^{\ast}u_{\bf k}(i).$

After the Hamiltonian is diagonalized, the pairing amplitude at site $i$ is
calculated by

$\displaystyle\Delta_{i}=\sum_{\bf k}V\theta(\omega_{D}-E_{\bf
k})\left[2f(E_{\bf k})-1\right]u_{\bf k}(i)v_{\bf k}^{\ast}(i),$ (14)

where $f(E_{\bf k})=(e^{E_{\bf k}/k_{B}T}+1)^{-1}$ is the Fermi distribution
and $V<0$ is the pairing interaction acting within Debye frequency
$\omega_{D}=0.25$eV, as ensured by the step function $\theta(\omega_{D}-E_{\bf
k})$. Equations (13) and (14) are solved self-consistently until the local
pairing amplitude $\Delta_{i}$ converges at a given pairing interaction and
temperature $\left\\{V,T\right\\}$.

Figure 4: (a) Local pairing amplitude $\Delta_{i}$ represented by circle size
for the C15 configuration in Fig. 2, calculated at zero temperature and
pairing interaction $V=-2.2$eV by our weak coupling mean-field model. The
largest circles correspond to $\Delta_{i}=0.036$eV. (b) $\Delta_{i}$ for the
C14 configuration in Fig. 2 at zero temperature and pairing interaction
$V=-0.5$eV, where the largest circles correspond to $\Delta_{i}=0.074$eV. In
these two figures, one sees that the $B$ sublattices that the vacancy sites do
not belong to have much larger pairing amplitude. (c) Spatially averaged
pairing amplitude at zero temperature $\Delta(0)$ versus the pairing
interaction $|V|$. (d) The spatially average pairing amplitude $\Delta(T)$ as
a function of temperature.

###  III.2 Numerical results for the local pairing

Numerical calculation of our mean-field theory applied to the pristine
graphene yields $\Delta_{i}=0$ everywhere, consistent with the experimental
observation that the single-layer pristine graphene has no SC. In contrast,
Fig. 4 (a) and (b) show the local gap $\Delta_{i}$ at zero temperature for the
two vacancy configurations C15 and C14 illustrated in Fig. 2. Our result
indicates a finite $\Delta_{i}$ in the holey graphene, and moreover
$\Delta_{i}$ on the $B$ sublattices is about two orders of magnitude larger
than that on the $A$ sublattices. This feature stems from the fact that the
chiral symmetry breaking terms are relatively small
$\left\\{t^{\prime},\mu\right\\}\ll t$, so the narrow band wave functions
still roughly satisfy the proposition 3 in Sec. II.1 and hence are mainly
localized on the $B$ sublattices, as discussed in Appendix A. For C15, the
spatially averaged gap $\Delta(T)=\sum_{i}\Delta_{i}/(N_{A}+N_{B})$ at zero
temperature $\Delta(0)$ is vanishingly small at small pairing potential $V$.
Only when the pairing potential has the same order of magnitude as the hopping
$|V|\sim t$ does a sizable gap emerge, implying that a sufficiently strong
electron-phonon interaction is needed to support SC in C15. On the other hand,
C14 requires much smaller $|V|$ to trigger SC in comparison with C15,
suggesting that increasing the number of narrow bands does help to create the
SC phase. Concerning the temperature dependence, the spatially averaged gap
shows a trend similar to the usual weak coupling superconductors, which
vanishes at a critical temperature $T_{c}$ that is higher at larger pairing
potential $V$. In addition, $T_{c}$ is generally raised in the configurations
with more narrow bands, consistent with that expected from an enlarged DOS.
Since the proposed engineering mechanism in principle has no restriction on
the number of narrow bands $N_{B}-N_{A}$ (times $2$ if including spin), we
anticipate that the vacancy configurations with very different numbers of the
two sublattices may yield a very high $T_{c}$, which is presumably more likely
to engineer in configurations with a larger unit cell.

##  IV Conclusions

In summary, we demonstrate that Lieb’s theorem of ZEFBs can be realized on
graphene by introducing periodic vacancies, which can be experimentally
relevant to holey graphene or graphene nanomesh. Although graphene in reality
does not preserve the chiral symmetry required by the theorem, periodic
vacancies can still induce bands as narrow as $\sim 0.5$eV near the Fermi
surface because the symmetry breaking factors are relatively weak compared to
the nearest-neighbor hopping. Moreover, our results suggest that periodic
vacancies can be used to combine various nonspatial and crystalline symmetries
to produce very exotic band structures, such as the coexisting ZEFBs and nodal
lines revealed in the present work, paving a way to engineer band structures
of 2D materials beyond the limitation set by the underlying crystalline
structures.

The vacancy-engineered narrow bands dramatically enlarge the DOS near the
Fermi surface, and hence are expected to significantly change the electronic
and magnetic properties of the holey graphene compared to the pristine one,
which await further investigations. In particular, our mean-field theory
survey reveals that a phonon-mediated conventional SC is feasible due to the
enlarged DOS. The local pairing amplitude is highly localized on the majority
sublattices, a feature originated from the chiral symmetry of the flat band
wave functions. The minimal electron-phonon interaction required to create the
SC phase varies significantly with the number of narrow bands, but increasing
the numbers of narrow bands generally reduces the minimal electron-phonon
interaction, as expected from the DOS point of view. As a result, we
anticipate that a certain experimental effort to search for the appropriate
vacancy configuration that has a sufficient number of narrow bands is needed
to observe the SC.

##  Appendix A Flat band wave functions with and without chiral symmetry

The proposition 3 in Sec. II.1 states that the ZEFB wave functions for a
chiral symmetric system in any vacancy configuration must localize on the
majority sublattices. As an example, in Fig. 5 we show the single flat band
wave function for C15 and the two degenerate flat band wave functions for C15
at momentum ${\bf k}=(0.15,0.37)$, both described by a spinless nearest-
neighbor hopping Hamiltonian $H=\sum_{\langle
ij\rangle}t\,c_{i}^{{\dagger}}c_{j}$ that preserves chiral symmetry. We find
that all these wave functions are localized in the majority $B$ sublattices
where the vacancies do not belong to, satisfying the proposition 3 in Sec.
II.1.

Figure 5: Wave functions $|\psi_{i}|^{2}$ at momentum ${\bf k}=(0.15,0.37)$ of
(a) the single ZEFB of C15 and (b) the two degenerate ZEFBs of C14 with only
nearest-neighbor hopping, whose band structure is that shown in Fig. 2 (b).
The largest circles correspond to $|\psi_{i}|^{2}=0.313$. All these wave
functions preserve chiral cymmetry and hence are localized only on the
majority $B$ sublattices.  Figure 6: Tight-binding band structure of (a) C15
and (b) C14 with $t=2.8$eV, $t^{\prime}=-0.2$eV, and $\mu=0.2$eV, which yields
a reasonable fit to the $p_{z}$ orbital bands obtained from DFT shown in Fig.
2 (c). The wave functions $|\psi_{i}|^{2}$ of the narrow bands at momentum
${\bf k}=(0.15,0.37)$ are shown in (c) and (d), which are still highly
localized on the $B$ sublattices since the chiral symmetry breaking terms are
relatively weak $\left\\{t^{\prime},\mu\right\\}\ll t$.

For the realistic graphene that does not preserve chiral symmetry, we rely on
numerical calculation to investigate the wave functions. Using a tight-binding
model with nearest $t$ and next-nearest-neighbor hopping $t^{\prime}$, and
additionally a chemical potential $\mu$ (described by the Hamiltonian in Eq.
11 without pairing and spin), we find that the parameters $t=2.8$eV,
$t^{\prime}=-0.2$eV and $\mu=0.2$eV can fit the $p_{z}$ orbital bands of the
DFT band structure shown in Fig. 2 (c) reasonably well, as shown in Fig. 6 (a)
for C15 and (b) for C14. The wave functions of the narrow bands close to Fermi
surface are shown in Fig. 6 (c) and (d) at the same momentum ${\bf
k}=(0.15,0.37)$ as that shown in Fig. 5 (a) and (b). We find that because
these chiral symmetry breaking terms $\left\\{t^{\prime},\mu\right\\}$ are
relatively small compared to $t$, the wave functions on the majority $B$
sublattices are about two orders of magnitude larger than that on the minority
$A$ sublattices. In other words, the wave functions still roughly preserves
the chiral symmetry and approximately satisfies the proposition 3 in Sec.
II.1. As a result, the local pairing amplitude in the SC phase is much larger
on $B$ sublattices, as shown in Fig. 4.

## References

  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Suárez Morell  _et al._ (2010) E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 (2011).
  * Moon and Koshino (2012) P. Moon and M. Koshino, Phys. Rev. B 85, 195458 (2012).
  * Trambly de Laissardière  _et al._ (2012) G. Trambly de Laissardière, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012).
  * Lopes dos Santos  _et al._ (2012) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * Fang and Kaxiras (2016) S. Fang and E. Kaxiras, Phys. Rev. B 93, 235153 (2016).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018b).
  * Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Zeitschrift fÃ¼r Physik B Condensed Matter 64, 189 (1986).
  * Lee  _et al._ (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
  * Emery (1987) V. J. Emery, Phys. Rev. Lett. 58, 2794 (1987).
  * Zhang and Rice (1988) F. C. Zhang and T. M. Rice, Phys. Rev. B 37, 3759 (1988).
  * Shen  _et al._ (2010) R. Shen, L. B. Shao, B. Wang, and D. Y. Xing, Phys. Rev. B 81, 041410 (2010).
  * Guzmán-Silva  _et al._ (2014) D. Guzmán-Silva, C. Mejía-Cortés, M. A. Bandres, M. C. Rechtsman, S. Weimann, S. Nolte, M. Segev, A. Szameit, and R. A. Vicencio, New J. Phys. 16, 063061 (2014).
  * Mukherjee  _et al._ (2015) S. Mukherjee, A. Spracklen, D. Choudhury, N. Goldman, P. Öhberg, E. Andersson, and R. R. Thomson, Phys. Rev. Lett. 114, 245504 (2015).
  * Taie  _et al._ (2015) S. Taie, H. Ozawa, T. Ichinose, T. Nishio, S. Nakajima, and Y. Takahashi, Sci. Adv. 1, e1500854 (2015).
  * Vicencio  _et al._ (2015) R. A. Vicencio, C. Cantillano, L. Morales-Inostroza, B. Real, C. Mejía-Cortés, S. Weimann, A. Szameit, and M. I. Molina, Phys. Rev. Lett. 114, 245503 (2015).
  * Julku  _et al._ (2016) A. Julku, S. Peotta, T. I. Vanhala, D.-H. Kim, and P. Törmä, Phys. Rev. Lett. 117, 045303 (2016).
  * Xia  _et al._ (2016) S. Xia, Y. Hu, D. Song, Y. Zong, L. Tang, and Z. Chen, Opt. Lett. 41, 1435 (2016).
  * Diebel  _et al._ (2016) F. Diebel, D. Leykam, S. Kroesen, C. Denz, and A. S. Desyatnikov, Phys. Rev. Lett. 116, 183902 (2016).
  * Slot  _et al._ (2017) M. R. Slot, T. S. Gardenier, P. H. Jacobse, G. C. P. van Miert, S. N. Kempkes, S. J. M. Zevenhuizen, C. M. Smith, D. Vanmaekelbergh, and I. Swart, Nat. Phys. 13, 672 (2017).
  * Klembt  _et al._ (2017) S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, H. Suchomel, J. Beierlein, M. Emmerling, C. Schneider, and S. HÃ¶fling, Applied Physics Letters 111, 231102 (2017).
  * Cui  _et al._ (2020) B. Cui, X. Zheng, J. Wang, D. Liu, S. Xie, and B. Huang, Nature Communications 11, 66 (2020).
  * Lieb (1989) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * Pawlak  _et al._ (2020) R. Pawlak, X. Liu, S. Ninova, P. D’Astolfo, C. Drechsel, S. Sangtarash, R. Häner, S. Decurtins, H. Sadeghi, C. J. Lambert, U. Aschauer, S.-X. Liu, and E. Meyer, J. Am. Chem. Soc. 142, 12568 (2020).
  * Mahmood  _et al._ (2015) J. Mahmood, E. K. Lee, M. Jung, D. Shin, I.-Y. Jeon, S.-M. Jung, H.-J. Choi, J.-M. Seo, S.-Y. Bae, S.-D. Sohn, N. Park, J. H. Oh, H.-J. Shin, and J.-B. Baek, Nat. Commun. 6, 6486 (2015).
  * Shi  _et al._ (2011) Z. Shi, R. Yang, L. Zhang, Y. Wang, D. Liu, D. Shi, E. Wang, and G. Zhang, Adv. Mater. 23, 3061 (2011).
  * Jung  _et al._ (2014) I. Jung, H. Y. Jang, J. Moon, and S. Park, Nanoscale 6, 6482 (2014).
  * Bai  _et al._ (2010) J. Bai, X. Zhong, S. Jiang, Y. Huang, and X. Duan, Nat. Nanotechnol. 5, 190 (2010).
  * Wang  _et al._ (2013) M. Wang, L. Fu, L. Gan, C. Zhang, M. RÃ¼mmeli, A. Bachmatiuk, K. Huang, Y. Fang, and Z. Liu, Sci. Rep. 3, 1238 (2013).
  * Archanjo  _et al._ (2014) B. S. Archanjo, B. Fragneaud, L. Gustavo CanÃ§ado, D. Winston, F. Miao, C. Alberto Achete, and G. Medeiros-Ribeiro, Appl. Phys. Lett. 104, 193114 (2014).
  * Skurativska  _et al._ (2021) A. Skurativska, S. S. Tsirkin, F. D. Natterer, T. Neupert, and M. H. Fischer, Phys. Rev. Research 3, L032003 (2021).
  * You  _et al._ (2019) J.-Y. You, B. Gu, and G. Su, Sci. Rep. 9, 20116 (2019).
  * Shima and Aoki (1993) N. Shima and H. Aoki, Phys. Rev. Lett. 71, 4389 (1993).
  * Tadjine  _et al._ (2016) A. Tadjine, G. Allan, and C. Delerue, Phys. Rev. B 94, 075441 (2016).
  * Saito  _et al._ (2001) R. Saito, A. Jorio, A. G. Souza Filho, G. Dresselhaus, M. S. Dresselhaus, and M. A. Pimenta, Phys. Rev. Lett. 88, 027401 (2001).
  * Maultzsch  _et al._ (2004) J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phys. Rev. Lett. 92, 075501 (2004).
  * Mohr  _et al._ (2007) M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Phys. Rev. B 76, 035439 (2007).
  * Jorio  _et al._ (2011) A. Jorio, M. S. Dresselhaus, R. Saito, and G. Dresselhaus,  _Raman Spectroscopy in Graphene Related Systems_ (Wiley, 2011).
  * Schnyder  _et al._ (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * Ryu  _et al._ (2010) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New J. Phys. 12, 065010 (2010).
  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * von Gersdorff  _et al._ (2021) G. von Gersdorff, S. Panahiyan, and W. Chen, Phys. Rev. B 103, 245146 (2021).
  * Chen and Schnyder (2019) W. Chen and A. P. Schnyder, New J. Phys. 21, 073003 (2019).
  * Lim  _et al._ (2020) L.-K. Lim, J.-N. Fuchs, F. Piéchon, and G. Montambaux, Phys. Rev. B 101, 045131 (2020).
  * Bouzerar and Mayou (2021) G. Bouzerar and D. Mayou, Phys. Rev. B 103, 075415 (2021).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Liu  _et al._ (2021) F. Liu, F. Qu, I. Žutić, S. Xie, D. Liu, A. L. A. Fonseca, and M. Malard, J. Phys. Chem. Lett. 12, 5710 (2021).
  * (51) F. Liu, F. Qu, I. Žutić, and M. Malard, arXiv:2108.05502 .
  * (52) M. S. M. de Sousa, F. Liu, M. Malard, F. Qu, and W. Chen, arXiv:2109.12468 .
  * Young and Kane (2015) S. M. Young and C. L. Kane, Phys. Rev. Lett. 115, 126803 (2015).
  * Yamakage  _et al._ (2016) A. Yamakage, Y. Yamakawa, Y. Tanaka, and Y. Okamoto, J. Phys. Soc. Jpn. 85, 013708 (2016).
  * Zhao and Schnyder (2016) Y. X. Zhao and A. P. Schnyder, Phys. Rev. B 94, 195109 (2016).
  * Wieder  _et al._ (2018) B. J. Wieder, B. Bradlyn, Z. Wang, J. Cano, Y. Kim, H.-S. D. Kim, A. M. Rappe, C. L. Kane, and B. A. Bernevig, Science 361, 246 (2018).
  * Kopnin  _et al._ (2011) N. B. Kopnin, T. T. Heikkilä, and G. E. Volovik, Phys. Rev. B 83, 220503 (2011).
  * Heikkilä and Volovik (2016) T. T. Heikkilä and G. Volovik,  _Basic Physics of Functionalized Graphite_ (Springer, 2016).


# ANN-assisted CoSaMP Algorithm for Linear Electromagnetic Imaging of
Spatially Sparse Domains

Ali I. Sandhu  Division of Computer, Electrical, and Mathematical Science and
Engineering (CEMSE)
King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi
Arabia, 23955-6900
{aliimran.sandhu, salman.shaukat<EMAIL_ADDRESS>Salman A. Shaukat
Division of Computer, Electrical, and Mathematical Science and Engineering
(CEMSE)
King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi
Arabia, 23955-6900
{aliimran.sandhu, salman.shaukat<EMAIL_ADDRESS>Abdulla Desmal
Department of Electrical Engineering, Higher Colleges of Technology, (HCT)
Ras Al-Khaimah, United Arab Emirates<EMAIL_ADDRESS>Hakan Bagci  Division
of Computer, Electrical, and Mathematical Science and Engineering (CEMSE)
King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi
Arabia, 23955-6900
{aliimran.sandhu, salman.shaukat<EMAIL_ADDRESS>
###### Abstract

Greedy pursuit algorithms (GPAs) are widely used to reconstruct sparse
signals. Even though many electromagnetic (EM) inverse scattering problems are
solved on sparse investigation domains, GPAs have rarely been used for this
purpose. This is because (i) they require a priori knowledge of the sparsity
level in the investigation domain, which is often not available in EM imaging
applications, and (ii) the EM scattering matrix does not satisfy the
restricted isometric property. In this work, these challenges are respectively
addressed by (i) using an artificial neural network (ANN) to estimate the
sparsity level, and (ii) adding a Tikhonov regularization term to the diagonal
elements of the scattering matrix. These enhancements permit the compressive
sampling matching pursuit (CoSaMP) algorithm to be efficiently used to solve
the two-dimensional EM inverse scattering problem, which is linearized using
the Born approximation, on spatially sparse investigation domains. Numerical
results, which demonstrate the efficiency and applicability of the proposed
ANN-enhanced CoSaMP algorithm, are provided.

Keywords: Compressed sensing, CoSaMP, electromagnetic imaging, inverse
problems, neural network, Born approximation, sparse reconstruction, machine
learning

##  1 Introduction

In the last two decades, greedy pursuit algorithms (GPAs) making use of
compressed sensing (CS) have been successfully used in reconstruction of
sparse signals[1, 2, 3, 4, 5]. In general, a CS-based scheme seek for the
sparsest solution by constraining the optimization problem with the
$\ell_{0}$-norm of the solution. It is well known that a direct solution of
this optimization problem is NP-hard, however, GPAs, under certain conditions,
provide a good approximation to its solution in an efficient way [3, 4, 5].
GPAs that are widely used in the image and signal processing communities rely
on the orthogonal matching pursuit (OMP) [3] and its different variants
including the regularized OMP [5] and the compressive sampling matching
pursuit (CoSaMP) [4].

Many electromagnetic (EM) imaging applications involve naturally sparse
investigation domains. Here, sparsity means that the scatterer, which is
defined by a contrast function (difference between the dielectric permittivity
of the scatterer and the background medium), occupies only a small fraction of
the investigation domain, i.e., contrast function is spatially sparse[6]. Note
that sparsity can also be defined with respect to a different basis expansion
or in a transformed domain [7, 8, 9]. Having said that, application of GPAs to
the solution of the EM inverse scattering problem has been rather limited [10,
11, 12]. This can be explained by two fundamental reasons: (i) For accurate
and efficient reconstruction, GPAs require the exact number of non-zero
elements (samples of the contrast function) in the investigation domain
(denoted by $k$ in this work) to be known a priori and this information is not
available in many EM imaging applications. (ii) The GPAs require the
measurement matrix to satisfy the restricted isometric property (RIP) [3, 4,
5]. For EM imaging applications, the measurement matrix is the scattering
matrix and is obtained by sampling the Green function of the background medium
between the measurement/receiver locations and the investigation domain. This
matrix does not satisfy the RIP [10, 11, 12].

In [10], a simultaneous OMP algorithm (S-OMP) is developed to reconstruct
sparse investigation domains, however, the algorithm is provided with the
value of $k$, which would not be possible in many EM imaging applications. In
[11], a phase-less reconstruction scheme is developed to reconstruct low-
contrast point like scatterers. The optimization problem is linearized under
the Born approximation [13] and made convex using $\ell_{1}$-norm
regularization. In [12], a flexible tree search based OMP algorithm (FTB-OMP)
is developed to reconstruct closely-spaced point-like scatterers and $k$ is
predicted using the data misfit at each stage of the tree search. This scheme
suffers from a tradeoff between the reconstruction accuracy and the
computational cost, both of which increase with the search tree size. Another
limitation of the FTB-OMP is that if the very first estimate of the solution
component is incorrect, the algorithm would converge to a local minimum or an
entirely incorrect solution.

It should be mentioned here that the EM inversion schemes making use of
Bayesian compressive sensing (BCS) frameworks (some of which also operate
under first-order linearization schemes such as Born and Rytov approximations)
have been developed [14, 15, 16, 17, 18]. Since these algorithms belong to the
family of probabilistic inversion schemes, they are not discussed here. The
readers are referred to [14, 15, 16, 17, 18] and the references cited in these
papers for details.

In this work, CoSaMP is used to solve the EM inverse scattering problem, which
is linearized using the Born approximation, on two-dimensional (2D) spatially
sparse domains. Unlike OMP, CoSaMP works by estimating multiple elements in
the support set of the scatterer (i.e., indices of non-zero contrast samples)
instead of one at a time and refines this support set iteratively.
Consequently, it offers a faster rate of convergence and can reconstruct large
connected scatterers, which is a limitation associated with the OMP-based
algorithms. The two main challenges associated with the application of GPAs to
EM imaging (as described above) are respectively addressed by (i) using an
artificial neural network (ANN) to estimate $k$ and (ii) applying Tikhonov
regularization to ensure that the RIP is satisfied. The ANN is trained
beforehand using synthetic but noisy (scattered field) measurements on many
investigation domains (with known $k$) which include various scatterers with
different contrast levels and are discretized using different number of
elements. This ANN is then used to predict $k$ given the fields scattered from
an unknown domain. The Tikhonov ($\ell_{2}$-norm) regularization term is added
to the diagonal elements of the scattering matrix to ensure that the RIP is
satisfied. This enables the use of CoSaMP and reduces the sensitivity of the
solution to the noise in the measurements.

The advantages of this Tikhonov and ANN-enhanced CoSaMP algorithm are
threefold: (i) It is computationally efficient since the least squares
problem, which is solved for unknown values of the contrast samples, is
executed on a reduced set rather than the full investigation domain, (ii) it
does not require tuning of a thresholding parameter [19], and (iii) the
reconstructed images are more accurate and sharper than those produced by
smoothness promoting inverse algorithms.

It should be emphasized here that the ANN used in this work predicts only the
number of non-zero elements ($k$) but neither their locations nor values. From
this perspective, it can be considered as a “qualitative” imaging approach and
its structure is much simpler compared to the ANNs used for “quantitative”
imaging where the contrast is fully reconstructed [20, 21, 22, 23, 24]. Having
said that, it still comes with a computational overhead due to training (even
though it is done only once) and one can extend other “qualitative” methods
[25, 26, 27, 28, 29] to predict $k$ while avoiding this computational cost.

It should also be noted here that a preliminary version of the CoSaMP
framework proposed here has been described in a conference contribution [30].
This inversion scheme work on current samples and is provided the exact value
of $k$.

##  2 Formulation

###  2.1 EM Inverse Scattering Problem

Let $S$ represent the support of a 2D inhomogeneous investigation domain that
resides in an unbounded background medium. The permittivity and permeability
in $S$ and in the background medium are
$\\{\varepsilon({\mathbf{r}}),{\mu_{0}}\\}$ and
$\\{{\varepsilon_{0}},{\mu_{0}}\\}$, respectively. $S$ is individually
illuminated by $N^{\mathrm{T}}$ line-source transmitters that generate
transverse magnetic (to $z$) $\mathrm{TM}^{z}$ incident fields
$E_{i}^{\mathrm{inc}}({\mathbf{r}})$, $i=1,\ldots,N^{\mathrm{T}}$ at frequency
$\omega=2\pi f$. Upon excitation by $E_{i}^{{\mathrm{inc}}}({\mathbf{r}})$, an
electric current density is induced on $S$ and this current density generates
the scattered electric field $E_{i}^{{\mathrm{sca}}}({\mathbf{r}})$. Let
$E_{i}^{{\mathrm{tot}}}({\mathbf{r}})=E_{i}^{{\mathrm{inc}}}({\mathbf{r}})+E_{i}^{{\mathrm{sca}}}({\mathbf{r}})$
represent the total electric field, then
$E_{i}^{{\mathrm{sca}}}({\mathbf{r}})$ can be expressed as [13]:

$E_{i}^{{\mathrm{sca}}}({\mathbf{r}})={k_{0}^{2}}\int\nolimits_{S}{\tau({\mathbf{r^{\prime}}})E_{i}^{{\mathrm{tot}}}({\mathbf{r^{\prime}}})G({\mathbf{r}}_{m}^{\mathrm{R}},{\mathbf{r^{\prime}}})ds^{\prime}}.$
(1)

Here, ${k_{0}}=\omega\sqrt{{\varepsilon_{0}}{\mu_{0}}}$ is the wavenumber in
the background medium,
$\tau(\mathbf{r})=\varepsilon(\mathbf{r})/\varepsilon_{0}-1$ is the dielectric
contrast, and
$G({\mathbf{r}},{\mathbf{r^{\prime}}})=H_{0}^{2}(k_{0}\left|{\mathbf{r}}-{\mathbf{r^{\prime}}}\right|)/(4j)$
is the 2D Green function.

The EM inverse scattering problem (Fig. 1) is defined as finding
$\tau(\mathbf{r})$ given measured
$E_{i}^{{\mathrm{sca}}}({\mathbf{r}_{m}^{\mathrm{R}}})$, where
${\mathbf{r}_{m}^{\mathrm{R}}}$, $m=1,\ldots,N^{\mathrm{R}}$ are receiver
locations away from $S$. To solve this problem numerically, first, $S$ is
discretized into $N$ number of square elements. Let $\mathbf{r}_{n}$,
$n=1,\ldots,N$ represent the centers of these elements. Then,
$E_{i}^{{\mathrm{tot}}}({\mathbf{r}})$ and $\tau({\mathbf{r}})$ are
approximated as

$\displaystyle E_{i}^{{\mathrm{tot}}}({\mathbf{r}})$
$\displaystyle=\sum\nolimits_{n=1}^{N}{{\\{\bar{E}_{i}^{\mathrm{tot}}\\}_{n}}{p_{n}}({\mathbf{r}})}$
(2) $\displaystyle\tau({\mathbf{r}})$
$\displaystyle=\sum\nolimits_{n=1}^{N}{{{\\{\bar{\tau}\\}}_{n}}{p_{n}}({\mathbf{r}})}$

where
${\\{\bar{E}_{i}^{\mathrm{tot}}\\}_{n}}=E_{i}^{{\mathrm{tot}}}({{\mathbf{r}}_{n}})$,
${\\{\bar{\tau}\\}_{n}}=\tau({{\mathbf{r}}_{n}})$, and ${p_{n}}({\mathbf{r}})$
is the pulse basis function on element $n$ with support ${S_{n}}$. Note that
${p_{n}}({\mathbf{r}})$ is non-zero only for ${\mathbf{r}}\in{S_{n}}$ with
unit amplitude. Substituting (2) into (1) and evaluating the expression at
${\mathbf{r}_{m}^{\mathrm{R}}}$, $m=1,\ldots,N^{\mathrm{R}}$ yields a matrix
equation

$\bar{E}_{i}^{\mathrm{sca}}=\bar{G}\bar{D}\\{\bar{E}_{i}^{\mathrm{tot}}\\}\bar{\tau},\;i=1,\ldots,N^{\mathrm{T}}$
(3)

where
${\\{\bar{E}_{i}^{\mathrm{sca}}\\}_{m}}=E_{i}^{{\mathrm{sca}}}({\mathbf{r}}_{m}^{\mathrm{R}})$,
$\bar{\bar{D}}\\{\bar{E}_{i}^{\mathrm{tot}}\\}$ is a diagonal matrix with
entries ${\\{\bar{E}_{i}^{\mathrm{tot}}\\}_{n}}$, and the entries of the
matrix ${\bar{G}}$ are
${\\{\bar{G}\\}_{m,n}}=k_{0}^{2}\int_{{S_{n}}}{G({\mathbf{r}}_{m}^{\mathrm{R}},{\mathbf{r^{\prime}}})ds^{\prime}}$.

Matrix equation (3) is nonlinear in $\bar{\tau}$ since
$\bar{E}_{i}^{\mathrm{tot}}$ is a function of $\bar{\tau}$. Under the
assumption of weak scattering, (3) can be linearized using the Born
approximation $\bar{E}_{i}^{\mathrm{tot}}\approx\bar{E}_{i}^{\mathrm{inc}}$,
where
${\\{\bar{E}_{i}^{\mathrm{inc}}\\}_{n}}=E_{i}^{{\mathrm{inc}}}({{\mathbf{r}}_{n}})$
[13]. Inserting this approximation into (3) yields a linear matrix equation

$\bar{E}_{i}^{\mathrm{sca}}=\bar{H_{i}}\bar{\tau}=\bar{G}\bar{D}\\{\bar{E}_{i}^{\mathrm{inc}}\\}\bar{\tau},\;i=1,\ldots,N^{\mathrm{T}}.$
(4)

This equation is “inverted” for $\bar{\tau}$ using the CoSaMP algorithm [4] as
described in the next section.

###  2.2 CoSaMP for EM Inverse Scattering Problem

In this work, the CoSaMP algorithm [4] is used to solve (4) for sparse
$\bar{\tau}$ as described next. First, (4) is expressed as an optimization
problem

$\bar{\tau}=\mathop{\min}\limits_{\bar{\tau}}{\left\|{\bar{\tau}}\right\|_{0}}\qquad\mathrm{s.t.}\left\|{\bar{E}^{\mathrm{meas}}-{\bar{H}}\bar{\tau}}\right\|_{2}^{2}\leq\epsilon.$
(5)

Here, $\bar{E}^{\mathrm{meas}}\approx\bar{E}^{\mathrm{sca}}+\bar{\eta}$, where
$\bar{\eta}$ represents the additive white Gaussian noise, and
$\bar{E}^{\mathrm{sca}}$ and $\bar{H}$ are obtained by respectively cascading
$\bar{E}_{i}^{\mathrm{sca}}$ and $\bar{H}_{i}$ for all
$i=1,\ldots,N^{\mathrm{T}}$.

For the CoSaMP algorithm to converge, the scattering matrix ${\bar{H}}$ has to
satisfy the RIP[4]. This condition states that for any vector $\bar{y}$ there
should be $\delta\in(0,1)$, such that
$(1-\delta)\|\bar{y}\|^{2}_{2}\leq\|{\bar{H}}\bar{y}\|^{2}_{2}\leq(1+\delta)\|\bar{y}\|^{2}_{2}$
holds. The largest value of $\delta$ that satisfies the RIP is known as the
restricted isometric constant (RIC). For any system that satisfies the RIP
with RIC value $\delta_{r}$, the following condition holds [4]:

$\sqrt{1-\delta_{r}}\leq\sigma({\bar{H}})_{\mathrm{min}}\leq\sigma({\bar{H}})_{\mathrm{max}}\leq\sqrt{1+\delta_{r}}.$
(6)

Here, $\sigma({\bar{H}})_{\mathrm{min}}$ and
$\sigma({\bar{H}})_{\mathrm{min}}$ are the largest and smallest singular
values of $\bar{H}$. For the EM inverse scattering problem, the lower bound in
(6) is not satisfied since $\bar{H}$ is ill-conditioned with
$\sigma({\bar{H}})_{\mathrm{min}}\to 0$, which enforces $\delta_{r}\to 1$, and
consequently results in violation of the RIP. To alleviate this problem and
enable the use of CoSaMP in solving the optimization problem (5), the data
misfit $||{\bar{E}^{\mathrm{meas}}-{\bar{H}}\bar{\tau}}||$ is replaced with
its Tikhonov-regularized version as

$\bar{\tau}=\mathop{\min}\limits_{\bar{\tau}}{\left\|{\bar{\tau}}\right\|_{0}}\quad\mathrm{s.t.}\left\|{\tilde{\bar{E}}^{\mathrm{meas}}-{{\bar{H}}}^{\lambda}\bar{\tau}}\right\|_{2}^{2}\leq\epsilon.$
(7)

Here,
$\tilde{\bar{E}}^{\mathrm{meas}}={\bar{H}}^{\dagger}\bar{E}^{\mathrm{meas}}$,
${\bar{H}}^{\lambda}={\bar{H}}^{\dagger}{\bar{H}}+\lambda\bar{\bar{I}}$,
${\bar{H}}^{\dagger}$ is the complex conjugate of ${\bar{H}}$, and $\lambda$
is a positive real constant. Note that (7) can also be written as

$\bar{\tau}=\mathop{\min}\limits_{\bar{\tau}}{\left\|{\bar{\tau}}\right\|_{0}}\quad\mathrm{s.t.}\left\|{\tilde{\bar{E}}^{\mathrm{meas}}-{{\bar{H}}}^{\dagger}{\bar{H}}\bar{\tau}}\right\|_{2}^{2}+\lambda^{2}{\left\|{\bar{\tau}}\right\|_{2}^{2}}\leq\epsilon.$
(8)

The form of the optimization problem provided in (8) clearly shows that it is
not only constrained by the sparsity requirement but also regularized using
the $\ell_{2}$-norm of the solution weighted with $\lambda^{2}$. Note that
$\sigma({\bar{H}}^{\lambda})=\sigma({{\bar{H}}}^{\dagger}{\bar{H}})+\lambda$,
and since ${{\bar{H}}}^{\dagger}{\bar{H}}$ is ill-conditioned with
$\sigma({{\bar{H}}}^{\dagger}{\bar{H}})_{\mathrm{\mathrm{min}}}\to 0$,
$\sigma({{{\bar{H}}}^{\lambda}})_{\mathrm{\mathrm{min}}}\to\lambda$.
Consequently, since $\lambda$ is a positive real constant, the lower bound of
the RIP can now be satisfied with $\delta_{r}\in(0,1)$.

The CoSaMP algorithm is applied to the optimization problem (8) to yield a
sparse solution. The steps of this algorithm read:

$\displaystyle{\rm Step}\;1:\;{\rm
initialize}\;\bar{r}^{(0)}={\bar{E}}^{\mathrm{meas}},n=0,k,\lambda$
$\displaystyle{\rm Step}\;2:\;\rm repeat$ $\displaystyle{\rm
Step}\;2.1:\;n=n+1$ $\displaystyle{\rm
Step}\;2.2:\;\bar{y}^{(n)}={\bar{H}}^{\dagger}\bar{r}^{(n-1)}$
$\displaystyle{\rm
Step}\;2.3:\;\Omega^{(n)}=\mathrm{supp}_{2k}(|\bar{y}^{(n)}|)$
$\displaystyle{\rm Step}\;2.4:\;F^{(n)}=\Omega^{(n)}\cup\Gamma^{(n)}$
$\displaystyle{\rm Step}\;2.5:\;{\rm
solve}\;({\bar{H}}_{:,F^{(n)}}^{\dagger}{\bar{H}}_{:,F^{(n)}}+\lambda\bar{I})\bar{s}^{(n)}=\tilde{\bar{E}}^{\mathrm{meas}}$
$\displaystyle{\rm
Step}\;2.6:\Gamma^{(n)}=\mathrm{supp}_{k}(|\bar{s}^{(n)}|),\;\bar{\tau}^{(n)}_{\Gamma^{(n)}}=\bar{s}^{(n)},\;\bar{\tau}^{(n)}_{\Gamma_{c}^{(n)}}=\bar{0}$
$\displaystyle{\rm
Step}\;2.7:\;\bar{r}^{(n)}={\bar{E}}^{\mathrm{meas}}-{\bar{H}}\bar{\tau}^{(n)}$

At Step $1$, several parameters are initialized. The number of non-zero
entries in $\bar{\tau}$ is denoted as $k$ and predicted by providing
${\bar{E}}^{\mathrm{meas}}$ to an ANN trained as described in the next
section. As explained above, the parameter $\lambda$ is primarily used to
ensure that the RIP is satisfied. But it also smoothens the high frequency
components since it acts as a Tikhonov-type regularizer. Therefore, it is
selected heuristically to increase the robustness of the solution to noise
while maintaining its sharpness within the support set of the scatterer (see
below). At Step $2.2$ the residual $\bar{r}^{(n-1)}$ from the last iteration
is projected onto the model subspace to determine which components of the
unknown model are yet to be determined. At Step $2.3$, function
$\mathrm{supp}_{L}(.)$ function choses indices of the largest $L$ entries of
its input vector. As a result, index/support set $\Omega^{(n)}$ stores $2k$
column-indices from ${\bar{H}}$ which contribute maximally towards the
projection at Step $2.2$. At Step $2.4$, the newly identified support set
$\Omega^{(n)}$ is unified with the final support set $\Gamma^{(n)}$ from Step
$2.6$ (of the previous iteration) to eliminate any repetitions in the support
elements. At Step $2.5$, $\bar{s}^{(n)}$ is computed by solving a least
squares problem over the merged support set $F^{(n)}$. Note that
${\bar{H}}_{:,F^{(n)}}^{\lambda}={\bar{H}}_{:,F^{(n)}}^{\dagger}{\bar{H}}_{:,F^{(n)}}+\lambda\bar{I}$
contain only the columns of ${\bar{H}}$ whose indices are in the merged
support set $F^{(n)}$. This significantly reduces the computational cost in
contrast to solving the least squares problem involving the “whole” matrix
${\bar{H}}$. At Step $2.6$, support set $\Gamma^{(n)}$ stores the indices of
the largest $k$ entries of $|\bar{s}^{(n)}|$. These are stored in the correct
entries of $\bar{\tau}^{(n)}$ while its remaining entries are set to zero.
Finally, at Step $2.7$, $\bar{r}^{(n)}$ is updated. The algorithm is set to
terminate if the residual between successive iterations does not change
significantly, i.e. $||\bar{r}^{(n)}-\bar{r}^{(n-1)}||/||\bar{r}^{(n)}||\leq
10^{-6}$.

###  2.3 Sparsity Estimation

Greedy algorithms (including CoSaMP) are very efficient in reconstructing
sparse signals however their reconstruction accuracy degrades if the number of
non-zero elements (denoted as $k$ in the previous section) is not known a
priori. Unfortunately, this is often the case for a typical EM inverse
scattering problem. In this work, an ANN, which is constructed, trained, and
tested as described below, is used to predict $k$. The value of this
prediction is denoted by $\hat{k}$.

The ANN used here is a simple two-layer feedforward perceptron network [31].
As shown in Fig. 2(a), the two hidden layers consist of 2048 and 64 neurons.
“ReLU” is the activation function applied to the outputs of each layer [31,
32]. The training step minimizes a loss function that is defined as the mean
squared error (MSE) between (i) $\hat{k}$ and $k$ (both of which are
normalized with $N$) and (ii) the true and estimated contrast level. Only
homogenous scatterers are used to train the ANN, i.e., contrast profile have
the same value for all elements that represent the scatterer. This value is
referred to as “contrast value” in the rest of the text. Note that the
contrast level is not the reconstructed image. It is estimated at the ANN
output and is used only to train the ANN better, which helps in producing a
more accurate $\hat{k}$. This information is not used by the CoSaMP algorithm
to reconstruct the contrast profile. The optimization algorithm is “RMSprop”
with a learning rate of 0.001 and default momentum of 0.0 [32].

This ANN directly accepts the scattered electric field measured/sampled at the
receiver locations (together with added noise) at its input layer. In all
training scenarios, the frequency of the transmitters $f=100$ MHz and A 2D
volume integral equation solver (for $\mathrm{TM}^{z}$ fields) [13], which
uses the discretization scheme described by (2) in the previous section, is
used to compute the scattered field samples at the receiver locations. Then,
$25\mathrm{dB}$ of Gaussian noise is added to these samples to generate a set
of measurements. Four different ANNs are generated. For the first ANN, ANN-1,
a total of $14450$ scattering scenarios including uniformly distributed
scatterers (i) of different shapes, in particular, circular rings with random
radii, single and double cylinders with random radii, and with varied
separation distance, (ii) with contrast level selected from set
$\\{0.2,0.4,0.6,0.8\\}$, and (iii) with varying total number of discretization
elements $N\in\\{784,3136,12544\\}$, are created. A randomly selected 70% of
this set is used to train the ANN, and the remaining 30% is used to test it.
These training and test sets are denoted by $T_{1}$ and $S_{1}$, respectively.
For ANN-1, $N^{\rm T}=32$ and $N^{\rm R}=32$.

Fig. 2(b) plots the convergence of the loss function for ANN-1 on $T_{1}$.
Even though the contrast level is not used in the reconstruction by the CoSaMP
algorithm, for the sake of completeness, the mean and variance of the error in
the contrast level estimated by ANN-1 in $S_{1}$, are provided here: $0.088$
and $0.007$. Fig. 2(c) shows the normalized probability histogram of the
difference $|k-\hat{k}|$ obtained over $S_{1}$. It is observed that for $90\%$
of the test cases $|k-\hat{k}|=0$ and the maximum value of $|k-\hat{k}|$ is
$3$, which occurs for less than $1\%$ of the test cases in $S_{1}$. The
histogram in Fig. 2(c) presents only $|k-\hat{k}|$ without providing any
information about the actual value of $k$. This figure is complemented by Fig.
2(d)-(e) that plots $k_{\rm min}$ vs $|k-\hat{k}|$. Here, $k_{\rm min}$ is the
smallest number of non-zero elements in the investigation domains of the test
cases that produce the given value of $|k-\hat{k}|$. For example, in Fig.
2(d), for $|k-\hat{k}|=1$, $k_{\rm min}=11$, which means that the largest
relative error is $|k-\hat{k}|/k_{\rm min}=9\%$, and for $|k-\hat{k}|=3$,
$k_{\rm min}=116$ and the largest relative error is $2.5\%$.

Two more test sets are generated to demonstrate the accuracy of ANN-1. Test
set $S_{2}$ is created by selecting the scattering scenarios with $N=784$ from
$S_{1}$ and changing their contrast level in $\\{0.3,0.7\\}$. Note that these
contrast levels are not used in $T_{1}$ which is used to train ANN-1. Test set
$S_{3}$ is created using Austria-like profiles [33] (also see Section III-C),
which do not exist in $T_{1}$. $S_{3}$ includes a total of $400$ scattering
scenarios where the scatterers include (i) a circular ring with constant
radius and the two outer cylinders rotated around for $360^{0}$ in steps of
$1^{0}$ (resulting in $360$ examples), and (ii) circular rings with varied
radii and outer cylinders rotated in steps of $90^{0}$ (resulting in $40$
examples). The contrast level of all scatterers is set to $0.2$ and
$N\in\\{3136,12544\\}$. The maximum value of $|k-\hat{k}|$ is 3 over all test
examples in $S_{2}$ and $S_{3}$. Fig. 2(d) plots $k_{\rm min}$ versus
$|k-\hat{k}|$ for both sets demonstrating the performance of ANN-1 for
predicting the sparsity level of Austria-like profiles and when the test
scatterers have contrast levels different from those of the training
scatterers.

For ANN-2, ANN-3, and ANN-4, $\\{N^{\rm T},N^{\rm R}\\}=$ $\\{16,16\\}$,
$\\{8,8\\}$, and $\\{4,4\\}$, respectively. All three ANNs are trained on set
$T_{2}$ which is created by selecting the scattering scenarios with contrast
level $0.4$ and $N=3136$ from $T_{1}$, and they are tested on set $S_{4}$
which is created from $S_{1}$ in the same way. Fig. 2(e) shows the three
normalized probability histograms of the difference $|k-\hat{k}|$ obtained
over $S_{4}$ using ANN-2, ANN-3, and ANN-4. The figure shows that the ANN is
accurate in predicting $\hat{k}$ even when $N^{\rm T}$ and $N^{\rm R}$ are
reduced to 4.

One last test is carried out to demonstrate the dependence of the prediction
accuracy on the frequency. A test case with $N=12544$, $k=288$, and contrast
level of $0.4$ is selected from $S_{1}$. ANN-1 (which is trained at $100$ MHz)
is used to compute $\hat{k}$ while the frequency of the test scattering
scenario is changed from $80$ MHz to $120$ MHz. Fig. 2(f) plots $|k-\hat{k}|$
versus the frequency. The figure shows that the prediction accuracy
deteriorates at higher frequencies but ANNs can still be used to provide $k$
for the CoSAMP algorithm at the frequencies that are closer to the frequency
at which they are trained at [see Fig. 4(d) and 5(d)].

##  3 Numerical results

This section demonstrates the accuracy and efficiency of the proposed scheme
via numerical experiments. In all experiments,
${{\tau}^{\mathrm{ref}}}(\mathbf{r})$ and
${{\\{{{\bar{\tau}}^{\mathrm{ref}}}\\}}_{n}}={{\tau}^{\mathrm{ref}}}({{\mathbf{r}}_{n}})$,
$n=1,\ldots,N$ are the actual contrast and its samples, respectively. A 2D
volume integral equation solver (for $\mathrm{TM}^{z}$ fields) [13], which the
discretization described in Section II-A is used to compute
$\bar{E}^{\mathrm{sca}}$ from ${{\bar{\tau}}^{\mathrm{ref}}}$, then
$\bar{E}^{\mathrm{meas}}$ is generated by adding Gaussian noise to
$\bar{E}^{\mathrm{sca}}$. Unless otherwise stated, the SNR for this noise is
$25$dB and $N^{\mathrm{T}}=N^{\mathrm{R}}=32$ and $f=100$ MHz in all
experiments. The quality of reconstruction is measured using

$err={{{\|{{{\bar{\tau}}}}-{{{\bar{\tau}}}^{\mathrm{ref}}}\|}_{2}}}/{{{\|{{{\bar{\tau}}}^{\mathrm{ref}}}\|}_{2}}}$
(9)

where ${{\bar{\tau}}}$ stores the samples of the reconstructed contrast.

###  3.1 Closely Spaced Point Like Targets

The investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ with $N=1369$ and
$k=4$ and the transmitter-receiver configuration with $N^{\mathrm{T}}=30$,
$N^{\mathrm{R}}=50$, and $f=300$ MHz are shown in Fig. 3(a). Note that this
example is taken from [12], where the FTB-OMP algorithm is introduced. Twelve
reconstructions are carried out using CoSaMP and FTB-OMP for six different SNR
values of the noise in $\bar{E}^{\mathrm{meas}}$. Fig. 3(c) plots $err$ versus
the SNR value. It is clear from the figure that CoSaMP is more accurate than
FTB-OMP over the whole range of SNR values considered for this example. Fig.
3(b) shows the reconstruction obtained using CoSaMP with $30$dB noise in
$\bar{E}^{\mathrm{meas}}$. Note that, for this example, $k$ is not predicted
using an ANN, instead its actual value is provided to CoSaMP.

###  3.2 Closely Spaced Cylinders

The investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ with $N=3136$ and
$k=72$ is shown in Fig. 4(a). ANN-1 estimates $k$ exactly, i.e., $\hat{k}=k$.
Fig. 4(b) and (c) show reconstructions obtained using CoSaMP and the Born
approximation with soft thresholding [19], respectively. Clearly the image
produced by CoSaMP is sharper and more accurate with $err=28\%$. It is
discussed in Section II-C that the maximum value of $|\hat{k}-k|$ for all test
examples is $3$ (even though for this specific example $\hat{k}=k$).
Therefore, CoSaMP is executed with $\bar{k}\in\\{k-3,k-2\ldots k+2,k+3\\}$
just to demonstrate the effect of $\bar{k}$ on the accuracy of reconstruction.
As expected, as shown in Fig. 4(d), the accuracy is highest when $\hat{k}=k$.

Next, the effect of using smaller $N^{\rm T}$ and $N^{\rm R}$ on the accuracy
of the solution is demonstrated. The reconstruction of the investigation
domain described above is carried out for $\\{N^{\rm T},N^{\rm R}\\}$ =
$\\{16,16\\}$, $\\{8,8\\}$, and $\\{4,4\\}$ (in addition to $\\{N^{\rm
T},N^{\rm R}\\}=\\{32,32\\}$ case above). The transmitters and receivers
surround the investigation domain. Just like ANN-1, ANN-2, ANN-3, and ANN-4
estimate $k$ exactly for these three configurations, respectively. Fig. 4(e)
plots the reconstruction error versus $N^{\rm T}\times N^{\rm R}$. The figure
shows that the accuracy of reconstruction degrades significantly for $N^{\rm
T}=4$ and $N^{\rm R}=4$. This demonstrates that the CoSaMP algorithm might
fail to achieve a good solution when the number of measurements is too small
(even though the ANN predicts the exact $k$).

###  3.3 Austria

The third example is the well-known Austria profile [33]. The investigation
domain ${{\bar{\tau}}^{\mathrm{ref}}}$ with $N=784$ and $k=66$ is shown in
Fig. 5(a). ANN-1 estimates $k$ exactly, i.e., $\hat{k}=k$. Reconstructions
obtained using CoSaMP and the Born approximation with soft thresholding [19],
are shown in Fig. 5(b) and (c) , respectively. CoSaMP is more accurate with
$err=43\%$. Similar to the previous example, just to demonstrate the effect of
$\hat{k}$ on the accuracy of reconstruction, CoSaMP is executed with
$\bar{k}\in\\{k-4,k-3\ldots k+3,k+4\\}$ (even though the ANN predicts
$\hat{k}=k$) and Fig. 5(d) shows that reconstruction accuracy degrades as
$|\hat{k}-k|$ increases.

###  3.4 L-Shaped Object

The last example is an L-shaped object. The investigation domain
${{\bar{\tau}}^{\mathrm{ref}}}$ with $N=961$ and $k=65$ is shown in Fig. 6(a).
CoSaMP reconstructs this object with $err=31.6\%$ as shown in Fig. 6(b).

##  4 Conclusion

A greedy algorithm is used together with a simple ANN for efficient and
accurate EM imaging of 2D sparse investigation domains. To enable the
application of CoSaMP to solving the EM inverse scattering problem (i) a
simple ANN is used to predict the number of non-zero contrast samples in the
investigation domain and (ii) a constant is added to the diagonal entries of
the scattering matrix (resulting in a Tikhonov-type regularization), which
ensures that the RIP condition is satisfied. The resulting EM inversion scheme
is computationally efficient since it calls for solution of a smaller least
squares problem on a reduced set determined by the number of non-zero contrast
samples. Numerical results, which demonstrate that the proposed scheme
produces more accurate and sharper images than Born approximation with soft
thresholding, are provided.

## Acknowledgment

The authors would like to thank the King Abdullah University of Science and
Technology (KAUST) Supercomputing Laboratory (KSL) for providing the required
computational resources.

## References

  * [1] E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” _IEEE Signal Process. Mag._ , vol. 25, no. 2, pp. 21–30, 2008.
  * [2] R. G. Baraniuk, “More is less: signal processing and the data deluge,” _Science_ , vol. 331, no. 6018, pp. 717–719, 2011.
  * [3] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 12, pp. 4655–4666, 2007.
  * [4] D. Needell and J. A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” _Appl. Comput. Harmon. A._ , vol. 26, no. 3, pp. 301–321, 2009.
  * [5] D. Needell and R. Vershynin, “Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit,” _IEEE J. Sel. Topics Signal Process._ , vol. 4, no. 2, pp. 310–316, 2010.
  * [6] A. Massa, P. Rocca, and G. Oliveri, “Compressive sensing in electromagnetics- a review,” _IEEE Antennas Propag. Mag._ , vol. 57, no. 1, pp. 224–238, 2015\.
  * [7] N. Anselmi, M. Salucci, G. Oliveri, and A. Massa, “Wavelet-based compressive imaging of sparse targets,” _IEEE Trans. Antennas Propag._ , vol. 63, no. 11, pp. 4889–4900, 2015.
  * [8] M. Li, O. Semerci, and A. Abubakar, “A contrast source inversion method in the wavelet domain,” _Inverse Probl._ , vol. 29, no. 2, p. 025015, 2013.
  * [9] A. I. Sandhu, A. Desmal, and H. Bagci, “A sparsity-regularized born iterative method for reconstruction of two-dimensional piecewise continuous inhomogeneous domains,” in _European Conference on Antennas and Propagation (EuCAP)_. IEEE, 2016, pp. 1–3.
  * [10] T. C. Ye and S. Y. Lee, “Non-iterative exact inverse scattering using simultaneous orthogonal matching pursuit (S-OMP),” in _Proc. IEEE Int. Conf. Acoustics, Speech and Signal Process._ , 2008, pp. 2457–2460.
  * [11] X. C. L. Pan and S. P. Yeo, “A compressive-sensing-based phaseless imaging method for point-like dielectric objects,” _IEEE Trans. Antennas Propag._ , vol. 60, no. 11, pp. 5472–5475, 2012.
  * [12] R. V. Senyuva, O. Ozdemir, G. K. Kurt, and E. Anarim, “Electromagnetic imaging of closely spaced objects using matching pursuit based approaches,” _IEEE Antennas Wireless Propag. Lett._ , vol. 15, pp. 1179–1182, 2015.
  * [13] M. Pastorino, _Microwave Imaging_. Hoboken, New Jersey, US: Wiley, 2010.
  * [14] G. Oliveri, L. Poli, N. Anselmi, M. Salucci, and A. Massa, “Compressive sensing-based Born iterative method for tomographic imaging,” _IEEE Trans. Microw. Theory Techn._ , vol. 67, no. 5, pp. 1753–1765, May 2019.
  * [15] A. E. Fouda and F. L. Teixeira, “Bayesian compressive sensing for ultrawideband inverse scattering in random media,” _Inverse Probl._ , vol. 30, no. 11, p. 114017, Oct 2014.
  * [16] L. Poli, G. Oliveri, P. P. Ding, T. Moriyama, and A. Massa, “Multifrequency Bayesian compressive sensing methods for microwave imaging,” _J. Opt. Soc. Am. A_ , vol. 31, no. 11, pp. 2415–2428, 2014.
  * [17] G. O. L. Poli and A. Massa, “Microwave imaging within the first-order Born approximation by means of the contrast-field Bayesian compressive sensing,” _IEEE Trans. Antennas Propag._ , vol. 60, no. 6, pp. 2865–2879, 2012.
  * [18] G. Oliveri, L. Poli, P. Rocca, and A. Massa, “Bayesian compressive optical imaging within the Rytov approximation,” _Opt. Lett._ , vol. 37, no. 10, pp. 1760–1762, 2012.
  * [19] A. Desmal and H. Bagci, “Shrinkage-thresholding enhanced Born iterative method for solving 2D inverse electromagnetic scattering problem,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 7, pp. 3878–3884, 2014.
  * [20] L. Li, L. G. Wang, F. L. Teixeira, C. Liu, A. Nehorai, and T. J. Cui, “Deepnis: Deep neural network for nonlinear electromagnetic inverse scattering,” _IEEE Trans. Antennas Propag._ , vol. 67, no. 3, pp. 1819–1825, 2019.
  * [21] Z. X. Y. Sun and U. Kamilov, “Efficient and accurate inversion of multiple scattering with deep learning,” _Opt. Express_ , vol. 26, no. 11, pp. 14 678–14 688, 2018.
  * [22] E. Bermani, S. Caorsi, and M. Raffetto, “An inverse scattering approach based on a neural network technique for the detection of dielectric cylinders buried in a lossy half-space,” _Prog. Electromagn. Res._ , vol. 26, pp. 67–87, 2000.
  * [23] M. T. McCann, K. H. Jin, and M. Unser, “Convolutional neural networks for inverse problems in imaging: A review,” _IEEE Signal Process Mag_ , vol. 34, no. 6, pp. 85–95, 2017.
  * [24] A. Massa, G. Oliveri, M. Salucci, N. Anselmi, and P. Rocca, “Learning-by-examples techniques as applied to electromagnetics,” _J. Electromagn. Waves Appl._ , vol. 32, no. 4, pp. 516–541, 2018.
  * [25] A. Litman, D. Lesselier, and F. Santosa, “Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set,” _Inverse Probl._ , vol. 14, no. 3, p. 685, 1998.
  * [26] O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” _Inverse Probl._ , vol. 22, no. 4, p. R67, 2006.
  * [27] M. Benedetti, D. Lesselier, M. Lambert, and A. Massa, “Multiple-shape reconstruction by means of multiregion level sets,” _IEEE Trans. Geosci. Remote Sens._ , vol. 48, no. 5, pp. 2330–2342, 2010.
  * [28] F. Cakoni, M. Fares, and H. Haddar, “Analysis of two linear sampling methods applied to electromagnetic imaging of buried objects,” _Inverse Probl._ , vol. 22, no. 3, p. 845, 2006.
  * [29] D. Colton, H. Haddar, and M. Piana, “The linear sampling method in inverse electromagnetic scattering theory,” _Inverse Probl._ , vol. 19, no. 6, p. S105, 2003.
  * [30] A. I. Sandhu and H. Bagci, “A modified CoSaMP algorithm for electromagnetic imaging of two dimensional domains,” in _Proc. Int. Appl. Comput. Electromagn. Soc. (ACES) Symp._ , 2017.
  * [31] N. Purkait, _Hands-On Neural Networks with Keras: Design and Create Neural Networks using Deep Learning and Artificial Intelligence Principles_. Packt Publishing Ltd, 2019.
  * [32] M. J. Kochenderfer and T. A. Wheeler, _Algorithms for optimization_. Mit Press, 2019.
  * [33] P. Van Den Berg and A. Abubakar, “Contrast source inversion method: state of art,” _J. Electromagn. Waves and Appl._ , vol. 15, no. 11, pp. 1503–1505, 2001.

Figure 1: Description of the 2D EM inverse scattering problem.

(a)  (b)

(c)  (d)

(e)  (f)

(g)

Figure 2:  (a) ANN used to estimate the sparsity level $k$. (b) Convergence of
the loss function for ANN-1 on $T_{1}$. (c) Normalized probability histogram
of $|k-\hat{k}|$ obtained by ANN-1 over $S_{1}$. (d) $k_{\rm min}$ versus
$|k-\hat{k}|$ over $S_{1}$ for ANN-1. (e) $k_{\rm min}$ versus $|k-\hat{k}|$
over $S_{2}$ and $S_{3}$ for ANN-1. (f) Normalized probability histogram of
$|k-\hat{k}|$ obtained by ANN-2, ANN-3 and ANN-4 (trained on $T_{2}$) over
$S_{4}$. (f) $|k-\hat{k}|$ versus the frequency of a test scenario for ANN-1.

(a)  (b)

(c)

Figure 3: (a) Investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ and
transmitter-receiver configuration for two closely spaced point like targets.
(b) $\bar{\tau}$ reconstructed using CoSaMP with $30$dB noise in
$\bar{E}^{\mathrm{meas}}$. (c) $err$ in $\bar{\tau}$ reconstructed using
CoSaMP and FTB-OMP versus six different levels of noise in
$\bar{E}^{\mathrm{meas}}$.

(a)  (b)  (c)

(d)  (e)

Figure 4: (a) Investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ and
transmitter-receiver configuration for two closely spaced cylinders.
$\bar{\tau}$ reconstructed using (b) CoSaMP with $\hat{k}=k$ (as estimated by
ANN-1) and (c) Born approximation with soft-thresholding. (d) $err$ in
$\bar{\tau}$ reconstructed using CoSaMP versus different values of $\hat{k}$.
(e) $err$ in $\bar{\tau}$ reconstructed using CoSaMP with $\hat{k}=k$ (as
estimated by ANN-4, ANN-3, ANN-2 and ANN-1) versus $N^{\rm T}\times N^{\rm
R}$.

(a)  (b)  (c)

(d)

Figure 5:  (a) Investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ and
transmitter-receiver configuration for Austria. $\bar{\tau}$ reconstructed
using (b) CoSaMP with $\hat{k}=k$ (as estimated by ANN-1) and (c) Born
approximation with soft-thresholding. (d) $err$ in $\bar{\tau}$ reconstructed
using CoSaMP versus different values of $\hat{k}$.

(a)  (b)

Figure 6: (a) Investigation domain ${{\bar{\tau}}^{\mathrm{ref}}}$ and
transmitter-receiver configuration for L-shaped object. (b) $\bar{\tau}$
reconstructed using CoSaMP.


# Effect of $\Sigma$-beam Asymmetry Data on Fits to Single Pion
Photoproduction off Neutron

I.I. Strakovsky111Electronic address<EMAIL_ADDRESS>R.A. Arndt222Deceased,
W.J. Briscoe, M.W. Paris, R.L. Workman  The George Washington University,
Washington, DC 20052, USA

###### Abstract

We investigate the influence of new GRAAL $\Sigma$-beam asymmetry measurements
on the neutron in multipole fits to the single-pion photoproduction database.
Results are compared to those found with the addition of a double-polarization
quantity associated with the sum rule.

###### pacs:

13.60.Le, 25.20.Lj, 13.88.+e, 11.80.Et

##  I Introduction

Only with good data on both proton and neutron targets one can hope to
disentangle the isoscalar and isovector EM couplings of the various baryon
resonances, as well as the isospin properties of the non-resonant background
amplitudes. In particular, the simple quark model predicts several resonances
that couple much stronger to the neutron than to the proton. The lack of
$\gamma n\to\pi^{-}p$ and $\gamma n\to\pi^{0}n$ data does not allow us to be
as confident about the determination of neutron couplings relative to those of
the proton.

Some of the N∗ baryons, N(1675)5/2-, for instance, have stronger
electromagnetic couplings to the neutron than to the proton but parameters are
very uncertain. PDG PDG  estimate for the A1/2 and A3/2 decay amplitudes of
the N(1720)3/2+ state are consistent with zero, while the recent SAID
determination Dugger  gives small but non-vanishing values. The reason for the
disagreement between the PDG estimate for the A1/2 decay amplitude and the
recent SAID determination Dugger  is also unclear. Other unresolved issues
relate to the N(1700)$3/2^{-}$ and second P11, N(1710)1/2+, that are not seen
in the recent GW $\pi$N partial-wave analysis (PWA) Ar06 , contrary to other
PWAs used by the Particle Data Group PDG .

New, high quality data on $\gamma n\to\pi^{-}p$ and $\gamma n\to\pi^{0}n$ are
needed to shed light on these issues, and the tagged-photon hall at GRAAL
offered a state-of-the-art facility to obtain such data. Here we report on an
analysis included novel $\Sigma$-measurements, covering incident photon
energies from threshold ($E_{\gamma}=707$ MeV) up to $E_{\gamma}=1500$ MeV.
The present measurement of $\Sigma$s for $\vec{\gamma}n\to\pi^{-}p$ Ma10  and
for $\vec{\gamma}n\to\pi^{0}n$ Sa09  is part of an extensive program at the
GRAAL to provide data of unrivaled quality on charged and neutral meson
photoproduction on the neutron, which includes polarized beam observable in
addition to the cross sections.

##  II New GRAAL Measurements for $\Sigma$ on the Neutron

Figure 1: $\Sigma$-beam asymmetry for $\vec{\gamma}n\to\pi^{0}n$. Data from
GRAAL Collaboration Sa09 . Solid lines correspond to the SAID-MA09 solution
(GRAAL data included in the database) Ma10 . Dash-dotted (dashed) lines show
the SAID-SP09 Dugger  (MAID2007 maid ) (GRAAL data excluded in the database).

To gauge the influence of new GRAAL data and their compatibility with previous
measurements, the GRAAL $\Sigma$s have been included in a number of fits using
the full SAID database for $\gamma N\to\pi N$ up to $E_{\gamma}=2.7$ GeV said
. The impact of new data on the SAID PWA can be understood from the comparison
of the new SAID fit MA09 Ma10 , which involves new GRAAL data, with the
previous SAID fit SP09 Dugger  and MAID2007 results maid .

216 $\Sigma$s for $\pi^{0}n$ final state at Eγ=703– 1475 MeV and
$\theta$=53–164∘ with 99 $\Sigma$s for $\pi^{-}p$ final state at Eγ=753–1439
MeV and $\theta$=33–163∘ GRALL data have been added to the GW SAID database
said . We have to notice that this GRAAL $\pi^{0}n$ contribution doubled the
World database for this reaction. Our best fit MA09 Ma10  for $\pi^{0}n$ and
$\pi^{-}p$, reduced initial $\chi^{2}$/dp=223 and 89 (SP09 Dugger ) to 3.1 and
4.9, respectively. It shows, in particular, that previous $\pi^{-}p$
measurements provided a better constraint vs. $\pi^{0}n$ case.

In Figs. 1 and 2, we show the excitation functions for several production
angles. The number of the distributions shown is enough to illustrate the
quality of new GRAAL data, the main features of the $\gamma n\to\pi N$
dynamics at the measured energy range, and the impact of the present data on
PWAs. The most noticeable effect of the present data on the new MA09 is due to
very good measurements of the medium-angle (65–140∘) $\Sigma$s for W in the
range above 1650 MeV. Earlier, this angular region either had been measured
with worse accuracy or could only be reached by extrapolation.

Figure 2: $\Sigma$-beam asymmetry for $\vec{\gamma}n\to\pi^{-}p$. Data from
GRAAL Collaboration Ma10 . The notation of the PWA solutions is the same as in
Fig. 1

The difference between our MA09 and SP09 results for the neutron target is
visible specifically for $S_{11}$nE (Fig. 3). It is observed above
E${}_{\gamma}\sim$400 MeV while modified MAID2007 shown a significant changes
vs. MAID2007 Sa09  above 1 GeV (see Fig. 7 at Ref. Sa09 ).

Figure 3: The multipole amplitude $S_{11}$nE (${}_{n}E_{0+}^{1/2}$). (a) Re
and (b) Im parts show. The vertical arrows indicate $W_{R}$ (Breit-Wigner
mass) and the horizontal bars show the full and partial width $\Gamma$ and
$\Gamma_{\pi N}$ associated with the SAID solution SP06 for $\pi N$ Ar06 .

The difference between previous pion photoproduction and new GRAAL
measurements may result in significant changes in the neutron couplings.

##  III Helicity-Dependent Photoabsorption Cross Sections on the Neutron

The amplitudes obtained in our analyses can be used to evaluate the single-
pion production component of several sum rules, in particular GDH, Baldin, and
forward spin polarizability Ar02 . In Table 1, we summarized our results for
the neutron target.

The running integrals are shown in Fig. 4. The evaluation of sum rules (GDH,
Baldin, and forward spin polarizability) for the neutron target and for a
single pion contribution exhibits convergence by 1 GeV. Agreement with Mainz
is good. Clearly, calculations above 450 MeV have to take into account
contributions beyond single-pion photoproduction.

Figure 4: Running (a) GDH, (b) Baldin, and (c) forward spin polarizability $\gamma_{0}$ integrals for the neutron target. The solid (dash-dotted) lines represent the SAID-MA09 (SAID-SP09) solution. Dashed lines show the MAID2007 predictions.  Table 1: Comparison of the recent SAID-MA09, SAID-SP09, and MAID2007 calculations for the GDH, Baldin and the forward spin polarizability from threshold up to 2.5 GeV in W (for MAID up to 2 GeV) and displayed as MA09/SP09/MAID2007. Reaction | GDH | Baldin | $\gamma_{0}$
---|---|---|---
| ($\mu b$) | ($10^{-4}fm^{3}$) | ($10^{-4}fm^{4}$)
$\gamma n\to\pi^{-}p$ | 21/ 21/ 20 | 8.4/8.4/8.7 | 1.4/ 1.5/ 1.4
$\gamma n\to\pi^{0}n$ |  $-$159/$-$157/$-$151 | 4.8/4.8/4.6 |  $-$1.5/$-$1.5/$-$1.5

###### Acknowledgements.

This work was supported in part by the U. S. Department of Energy under Grant
DE-FG02-99ER41110.

## References

  * (1) K. Nakamura, et al. (Particle Data Group), _J. Phys. G_ 37, 075021 (2010).
  * (2) M. Dugger, et al. (CLAS Collaboration), _Phys. Rev. C_ 79, 065206 (2009).
  * (3) R.A. Arndt, W.J. Briscoe, I.I. Strakovsky, and R.L. Workman, _Phys. Rev. C_ 74, 0605082 (2006).
  * (4) G. Mandaglio, et al. (GRAAL Collaboration), _Phys. Rev. C_ 82, 045209 (2010).
  * (5) R. Di Salvo, et al. (GRAAL Collaboration), _Eur. Phys. J. A_ 42, 151 (2009).
  * (6) D. Drechsel, S.S. Kamalov, and L. Tiator, _Eur. Phys. J. A_ 34, 69 (2007).
  * (7) The SAID website contains data and fits for this and a number of other medium-energy reactions: http://gwdac.phys.gwu.edu.
  * (8) R.A. Arndt, W.J. Briscoe, I.I. Strakovsky, and R.L. Workman, _Phys. Rev. C_ 66, 055213 (2002).