# Projective Quasiparticle Interference of a Single Scatterer
to Analyze the Electronic Band Structure of ZrSiS

Wenhao Zhang  Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Kunliang Bu  Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Fangzhou Ai  Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Zongxiu Wu  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Ying Fei
Zhejiang Province Key Laboratory of Quantum Technology and Device, Department
of Physics, Zhejiang University, Hangzhou, 310027, China  Yuan Zheng  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Jianhua Du  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Minghu Fang  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Collaborative
Innovation Center of Advanced Microstructures, Nanjing University, Nanjing
210093, China  Yi Yin<EMAIL_ADDRESS>Zhejiang Province Key Laboratory of
Quantum Technology and Device, Department of Physics, Zhejiang University,
Hangzhou, 310027, China  Collaborative Innovation Center of Advanced
Microstructures, Nanjing University, Nanjing 210093, China

###### Abstract

Quasiparticle interference (QPI) of the electronic states has been widely
applied in scanning tunneling microscopy (STM) to analyze electronic band
structure of materials. Single-defect-induced QPI reveals defect-dependent
interaction between a single atomic defect and electronic states, which
deserves special attention. Due to the weak signal of single-defect-induced
QPI, the signal-to-noise ratio (SNR) is relatively low in a standard two-
dimensional QPI measurement. In this paper, we introduce a projective
quasiparticle interference (PQPI) method, in which a one-dimensional
measurement is taken along high-symmetry directions centered on a specified
defect. We apply the PQPI method to a topological nodal-line semimetal ZrSiS.
We focus on two special types of atomic defects that scatter the surface and
bulk electronic bands. With enhanced SNR in PQPI, the energy dispersions are
clearly resolved along with high symmetry directions. We discuss the defect-
dependent scattering of bulk bands with the non-symmorphic symmetry-enforced
selection rules. Furthermore, an energy shift of the surface floating band is
observed and a new branch of energy dispersion (${\bm{q}}_{6}$) is resolved.
This PQPI method can be applied to other complex materials to explore defect-
dependent interactions in the future.

##  I Introduction

In scanning tunneling microscopy (STM), quasiparticle interference (QPI) of
electronic states has been a powerful tool to analyze the electronic band
structure of condensed matter materials Crommie1993 ; Hasegawa93 ; LCDavis91 ;
Avouris94 ; Petersen98 ; Sprunger1764 ; PhHofmann97 ; Pascual04 ; Weismann1190
; Hoffman1148 ; Wang03 ; kostin2018imaging ; Aynajian2012 ; Allan2013 ;
Roushan2009 ; ZhangTong09 ; Zeljkovic2014 ; Chang2016 ; Inoue1184 ; Zheng2016
; Batabyale1600709 . QPI arises when the electronic state with initial
momentum ${\bm{k}}_{i}$ is elastically scattered to a final momentum
${\bm{k}}_{f}$. The potential barrier of scattering is often induced by point
defects, steps or other local perturbations in materials. The scattering
process leads to a spatial oscillation of electronic state with wave vector
${\bm{q}}={\bm{k}}_{f}-{\bm{k}}_{i}$. The wave vector can be extracted from
Fourier transform of QPI oscillations. As a function of energy E, the
${\bm{q}}(E)$ dispersion reflects the electronic band structure in
${\bm{k}}$-space.

The QPI study initially focused on electronic surface state, whose QPI
oscillation (or Friedel oscillation) decays slowly with increasing distance
from the scattering center Crommie1993 ; Hasegawa93 ; LCDavis91 ; Avouris94 .
QPI has been thereafter applied in both surface states of materials, and
electronic structure of two-dimensional (2D) materials Petersen98 ;
Sprunger1764 ; PhHofmann97 ; Pascual04 . On the other hand, parallel features
in Fermi surface structure may cause anisotropic propagation of a three-
dimensional (3D) band, which can also result in a strong QPI oscillation on
the sample surface Weismann1190 . The standard QPI measurement requires a 2D
grid measurement, while in some special case it can be reduced to $1D$
measurement due to a quasi $1D$ electronic structure near an edge or
homogeneous electronic structure induced by a single defect Haim2018am ;
Liu_2017 ; Drozdov_2014 . The development of QPI technique enables extensive
analysis of band structure of complex materials, including
high-$T_{\mathrm{c}}$ superconductors Hoffman1148 ; Wang03 ; kostin2018imaging
, heavy fermion systems Aynajian2012 ; Allan2013 , and topological materials
Roushan2009 ; ZhangTong09 ; Zeljkovic2014 ; Chang2016 ; Inoue1184 ; Zheng2016
; Batabyale1600709 .

Although all QPI oscillations are related to the underlying electronic band
structure, QPI induced by a single scatterer deserves special attention
Inoue1184 ; Zheng2016 ; Simon2009 ; Derry2015 . Different types of point
defects trigger defect-dependent interaction between the defect and electronic
states. The QPI analysis around specified point defects can reveal a selective
scattering. For example, in topological nodal-line materials ZrSiS and ZrSiSe
xu2015two ; hu2016evidence ; schoop2016dirac , two different types of point
defects are found to scatter electronic states of the floating surface band
ZZ2018NatCom ; Surface_floating  and bulk band, respectively
lodge2017observation ; butler2017quasiparticle ; Bu2018PRB ; ZZ2018NatCom ;
surface_termination . In ZrSiSe, both the surface and bulk bands were observed
to be scattered by a single defect Bu2018PRB ; ZZ2018NatCom , which has not
been detected yet in ZrSiS.

However, the QPI signal around a single scatterer is relatively weak,
resulting in a low signal-to-noise ratio (SNR) in the Fourier-transformed QPI
pattern. The vague band structure in this analysis limits the data-based
discussion of physical properties. Here in this paper, a new type of point
defect is discovered in ZrSiS, and a projective QPI (PQPI) method is
introduced to analyze the scattered electronic bands in two different point
defects, with a much higher SNR. The first new type of Zr-site defect scatters
both the surface and bulk band. The second type S-site defect only scatters
the bulk band, which has been observed before lodge2017observation  with a
different interpretation. A preliminary 2D QPI measurement shows that the QPI
pattern induced by a single defect is anisotropic and highly concentrated
along high symmetry directions. With the PQPI method, we could clearly resolve
the dispersion branches and compare them with the density functional theory
(DFT) calculation. We discuss the selective scattering with non-symmorphic
symmetry-related selection rules. We also observe a possible defect-induced
energy shift of the floating surface band, and an extra bulk band dispersion
of ${\bm{q}}_{6}$ branch. PQPI is a simple and intuitive method that can be
applied in general single scatterer induced QPI studies of different
materials.

##  II Experimental Method

Single crystals of ZrSiS were grown by a two-step chemical vapor transport
method using iodine as a transport agent sankar2017crystal . In the first
step, a stoichiometric amount of 99.9 % purity precursors of Zr:Si:S = 1:1:1
molar ratio was pressed into a tablet, and put in an alumina crucible. After
sealed in an evacuated quartz ampoule, the sample was treated at 1100 ∘C for
two days and then furnace-cooled to room temperature. In the second step, the
tablet of ZrSiS was ground into a fine powder and then sealed in an evacuated
quartz ampoule with 5mg/cm3 iodine. The quartz ampoule was treated in a two-
zone tube furnace with a thermal gradient of about 1100 ∘C - 950 ∘C. After a
period of 8 days, single crystals of ZrSiS were obtained.

STM measurements were carried out in a commercial ultra-high vacuum system
Bu2018PRB ; zheng2017study ; FeiY_O_vacancy ; FeiY_Charge_order . An
electrochemically etched tungsten tip was treated with field emission on a
single crystalline of the Au (111) surface. The samples were cleaved in situ
at liquid nitrogen temperature and immediately inserted into the STM head. A
bias voltage $V_{\mathrm{b}}$ was applied to the sample, and the tunneling
current collected from the tip was maintained at a fixed setpoint
$I_{\mathrm{s}}$ by a feedback loop control. All data were acquired at liquid
helium temperature ($\sim$ 4.5 K). The differential conductance ($dI/dV$)
spectrum was acquired with a standard lock-in technique with modulation of 10
mV at 1213.7 Hz. The integration time is 3 ms for a single spectrum in 2D
measurement. In a grid measurement (2D or 1D), the tip was moved to a
different grid point in constant current mode. At each grid point, the
feedback was turned off while taking the corresponding $dI/dV$ spectrum.
Afterward the feedback was turned on again, and the tip was moved to the next
point for data collection. The DFT calculations were carried out using the
Vienna $\it{ab}$ $\it{initio}$ simulation package (VASP) Kresse96efficient ;
Kresse99paw ; Perdew96GGA ; Marzari97prb . A $1\times 1\times 5$ supercell
with a vacuum layer larger than $2$ nm was applied in the slab model.

##  III Results and Discussion

Figure 1: (a) Crystal structure of ZrSiS, with a cleavage plane between S
layers. Top view of the cleaved surface is shown on right. The yellow, green
and blue dots represent S, Zr and Si atoms, respectively. (b) A $25\times 25$
nm2 topography of ZrSiS taking under $I_{s}=1$ nA and $V_{\mathrm{b}}=600$ mV.
The two perpendicular white arrows represent lattice directions. (c) The
average $dI/dV$ spectra at the top (red), bridge (orange) and hollow (black)
sites in the supercell image (5.2 $\times$ 5.2 Å2), which is shown in the
inset. (d) Four different point defects in topography under the same bias
voltage $V_{\mathrm{b}}=500$ mV.

The family of ZrSiX (X $=$ S, Se, Te) shares a layered crystal structure. In
Fig. 1(a), the crystal structure of ZrSiS shows that a square lattice of Si
atoms is sandwiched between two sets of ZrS bilayers with glide mirror
symmetry xu2015two ; Wang1995Main . Then the crystal structure of ZrSiS is
non-symmorphic with Si lattice as the mirror plane. After an inversion towards
the mirror plane and a glide of the $ab$ plane with a vector of
$(1/2,1/2)a_{0}$ (where $a_{0}$ is the lattice constant of the ZrS bilayer),
the crystal structure becomes the same as the original one. This non-
symmorphic symmetry is critical to the topological properties of ZrSiS. For
STM experiment, the single crystal sample is cleaved between two ZrS bilayers,
with a S layer naturally exposed to be the surface plane. Figure 1(b) displays
a typical topography taken on the exposed S surface, with a tunneling junction
of $V_{\mathrm{b}}=600$ mV and $I_{\mathrm{s}}=1$ nA. In the topography, a
clear square lattice can be observed with an interatomic spacing of
$a_{0}\approx 3.6$ Å. For the family of ZrSiX, the density of states (DOS)
around the Fermi level is mainly contributed by $d$ electrons of Zr atoms
xu2015two . Top sites of the square lattice are determined to be at locations
of Zr atoms, even though Zr atoms are beneath the cleaved surface plane of S
atoms.

Different from that in ZrSiSe Bu2018PRB , our ZrSiS experiment shows a bias-
independent topography, without a shift of square lattice for different bias-
voltage polarities. In a clean area of the sample (12 $\times$ 12 nm2), we
performed a 2D $dI/dV$ spectrum measurement, with the topography acquired
simultaneously. A supercell image was created by overlaying portions of the
topography, following the algorithm in Ref. Lawler10nature ; Zeljkovic12NatMat
. The supercell image is shown in the inset of Fig. 1(c), based on which the
measured $dI/dV$ spectra are separately extracted over the top, hollow, and
bridge sites. As shown in Fig. 1(c), the spectra at different sites are almost
indistinguishable. They all exhibit a nonzero DOS around the Fermi level (zero
bias), while the DOS of occupied state is smaller than that of the empty
state. The spatially homogeneous spectrum is consistent with the bias-
independent topography.

In Fig. 1(b), sparse point defects of different types can be observed within
the square lattice. A typical diamond- and X-shaped defects are enlarged in
the top row in Fig. 1(d), whose centers are at Zr and S sites, respectively.
In previous STM studies of both ZrSiS and ZrSiSe lodge2017observation ;
butler2017quasiparticle ; Bu2018PRB ; ZZ2018NatCom , the diamond (X-shaped)
defects are found to selectively scatter the electronic surface band (bulk
band). For ZrSiSe, a strong scatterer is found to scatter both the surface and
bulk bands Bu2018PRB , which is hitherto not reported in ZrSrS. In the bottom
row of Fig. 1(d), we show two different types of atomic point defects in
ZrSiS, QPI around which is the main focus of this paper.

Figure 2: (a) A $16\times 16$ nm2 topography under $V_{\rm b}=300$ mV with
defect $\\#$1 at the center. The orange and red line across the defect are
along the lattice and diagonal directions, respectively. (b) The
$dI/dV(V=300~{}\mathrm{mV})$ map simultaneously taken with (a). (c) Fourier
transform of the $dI/dV$ map in $\bm{q}$-space. (d) A CCE model in
$\bm{k}$-space, and (e) the calculated $\bm{q}$-space map following this CCE
model. (f-g) QPI energy dispersions along red and orange lines shown in (c).

The bottom left defect in Fig. 1(d) is centered at the Zr site, around which a
larger topography ($16\times 16$ nm2 ) is shown in Fig. 2(a). To study QPI
around a single atomic defect, a standard method is to obtain the $dI/dV$ maps
from spectroscopy measurement. For a single $dI/dV$ map, a quick procedure is
to collect the $dI/dV$ signal at the fixed bias voltage while scanning the tip
in constant current mode. Along with the constant-current topography in Fig.
2(a), a $dI/dV$ map was taken simultaneously at $V_{\mathrm{b}}=300$ mV. As
shown in Fig. 2(b), this $dI/dV$ map exhibits an obvious pattern of standing
wave centered around the point defect, referred as a QPI image later. The
standing wave originates from the point-defect-induced scattering between
electronic states of different wave vectors (initial ${\bm{k}}_{i}$ and final
${\bm{k}}_{f}$) but the same energy. In Fig. 2(b), the QPI image is not
azimuthally symmetric but shows strong oscillations along the lattice
direction and the diagonal ($45^{\circ}$) direction. Fourier transform of the
$dI/dV$ map is calculated and drawn in Fig. 2(c). Similar to the previous
report for ZrSiSe Bu2018PRB , the Fourier-transformed QPI pattern can be
mainly partitioned into three groups: the central diamond, the concentric
square, and four sets of parallel lines around Bragg peaks.

The QPI pattern is described in the momentum $\bm{q}$ space with
${\bm{q}}={\bm{k}}_{f}-{\bm{k}}_{i}$. Figure 2(d) shows a contour of constant
energy (CCE) model similar to that in Ref. Bu2018PRB . For the floating
surface band ZZ2018NatCom ; Surface_floating , there are four pairs of short
parallel arcs around four X points. The scattering between parallel arcs in
the same pair (${\bm{q}}_{1}$) corresponds to the central diamond in the QPI
pattern. Scattering between short arcs in diagonal pairs (${\bm{q}}_{2}$)
corresponds to the parallel lines around Bragg peaks in the QPI pattern shown
in Fig. 2(e). For the bulk band, two concentric squares in the CCE model may
contribute to concentric squares in the QPI pattern, with possible wave
vectors ${\bm{q}}_{3}$, ${\bm{q}}_{4}$ and ${\bm{q}}_{5}$. In the QPI pattern
in Fig. 2(c), the scattering of both the surface band and bulk band can be
identified, confirming the discovery of a new type of scatterer in ZrSiS. For
the concentric squares in the QPI pattern, either a single square or two
squares have been found for different point defects in ZrSiX
lodge2017observation ; butler2017quasiparticle ; Bu2018PRB ; ZZ2018NatCom .
The concentric square indicated by wave vector ${\bm{q}}_{5}$ has never been
found in literature and our experimental results. Without a high SNR in the
QPI pattern, it is hard to judge whether the results are intrinsic
characteristics of the point defect or just vague and indistinct signals with
limited SNR. For simplicity, we intentionally forbid the scattering between
outer to outer squares when calculating the $\bm{q}$-space map in Fig. 2(e).

To study the energy-dependent QPI pattern, a three-dimensional (3D) dataset
has to be taken. For each spatial point in the scan area, the feedback loop is
temporarily interrupted and a $dI/dV$ spectrum is taken in a selected voltage
range, with the energy $E=eV$. After the measurement, the energy-dependent
$dI/dV$ maps can be extracted from the 3D dataset for further analysis. A long
time measurement is necessary for this process (e.g. 12-24 hours), in which
the system instability and thermal-drift affect the SNR. To display the
energy-dependent QPI result, the Fourier-transformed result is often shown
along a high symmetry direction in $\bm{q}$ space and plotted as a function of
the energy. As shown in Fig. 2(f), the Fourier-transformed result is displayed
along the high symmetry direction, from $(1,-1){\pi}/{a_{0}}$ to
$(-1,1){\pi}/{a_{0}}$ in $\bm{q}$-space [red line in Fig. 2(c)]. The
concentric square in the QPI pattern intersects with this line at the wave
vector ${\bm{q}}$, later confirmed as ${\bm{q}}_{3}$. In Fig. 2(f), the
energy-dependent dispersion of ${\bm{q}}_{3}$ can be observed, as guided by
the red solid line. Similarly, figure 2(g) shows the Fourier-transformed
result along the orange line in Fig. 2(c), from one Bragg peak of
$(1,0){2\pi}/{a_{0}}$ to the diagonal Bragg peak of $(-1,0){2\pi}/{a_{0}}$ in
$\bm{q}$-space. The diamond in the QPI pattern intersects with this line at
${\bm{q}}_{1}$, and the dispersion of ${\bm{q}}_{1}$ can be observed in Fig.
2(g). We can observe a limited SNR in the energy-dependent results, which
hinders a precise extraction of dispersions of ${\bm{q}}_{3}$ and
${\bm{q}}_{1}$ branches.

Figure 3: (a) Linecut measurement of $dI/dV$ spectrums along the red line in
Fig. 2(a) for defect $\\#$1\. (b) Fourier transform of the linecut measurement
in (a). Maximum of energy dispersion is guided with a red solid line, and
labeled by ${\bm{q}}_{3}$. A red dashed line represents the absence of
possible scattering of ${\bm{q}}_{4}$ branch. (c) Linecut measurement of
$dI/dV$ spectrums along the orange line in Fig. 2(a) for defect $\\#$1\. (d)
Fourier transform of the linecut measurement in (c). Maximum of two energy
dispersions are guided with orange (${\bm{q}}1$) and black (${\bm{q}}6$)
lines, respectively.

Here in this work, we introduce a simple and intuitive method, a projective
quasiparticle interference (PQPI) on a single defect, to study the same
energy-dependent scattering and extraction of the electronic band structure.
In the two-dimensional QPI image [Fig. 2(a)], the standing wave propagates
strongly along the lattice direction and the diagonal direction. In this PQPI
method, the $dI/dV$ spectra were measured at dense spatial points along with
two corresponding linecuts as labeled in Fig. 2(a). By decreasing dimension
from 2D to 1D in real space measurement, we could increase the average number
in the spectroscopy measurement. In the following 1D linecut measurement, each
spectrum is acquired with the same parameters as in 2D measurement but
averaged 5 times. The effect of system instability and thermal-drift is also
lessened within the short measurement time (e.g. half an hour for a single
linecut).

In Fig. 3(a), the measured $dI/dV$ spectrum is shown as a function of energy
(each vertical line), along the linecut of diagonal direction. For each
energy, the oscillating standing wave can be observed along the linecut in the
real space. The real-space signal can be further Fourier-transformed, leading
to the $\bm{q}$-space QPI pattern along the high symmetry direction from
$(1,-1){\pi}/{a_{0}}$ to $(-1,1){\pi}/{a_{0}}$. As shown in Fig. 3(b), the
${\bm{q}}_{3}$ branch is more clearly identified, from the strongly enhanced
SNR in the PQPI measurement. In the meantime, there is no clear dispersion
signal of scatter wave vector ${\bm{q}}_{4}$ (guided by a red dashed line) in
Fig. 3(b). Because of the short measurement time in PQPI, the energy range is
enlarged to $[-400,1000]$ meV with an energy resolution of $4$ meV. A similar
spectroscopy measurement was taken along the lattice direction, with the real-
space data shown in Fig. 3(c). The Fourier-transformed result in
$\bm{q}$-space is shown in Fig. 3(d), in which the ${\bm{q}}_{1}$ branch
exhibits a clearly resolved dispersion (later confirmed by DFT calculations).
With the relative high SNR, another ${\bm{q}}_{6}$ branch is also identified,
which will be discussed later. Putting Fig. 3(b) and 3(d) together, we
conclude that this special impurity scatters both electronic surface and bulk
band, similar to the special defect found in ZrSiSe Bu2018PRB . Only one
${\bm{q}}_{3}$ branch is identified for the scattering within concentric
squares.

Figure 4:  (a) A $16\times 16$ nm2 topography under $V_{\rm b}=600$ mV with
defect $\\#$2 at the center. (b) The $dI/dV(V=300~{}\mathrm{mV})$ map
simultaneously taken with (a). (c) Fourier transform of $dI/dV$ map in
$\bm{q}$-space. (d-e) The energy dispersions along the M-$\Gamma$-M direction.
The result in (d) is extracted from a standard 2D $dI/dV$ map along the red
line in (c). The result in (e) is extracted from a linecut measurement along
the red line in (a). Two red guide lines highlight the scattering of
${\bm{q}}_{3}$ and ${\bm{q}}_{4}$ branches. (f-g) The energy dispersions along
the X-$\Gamma$-X direction. The results are extracted from two different
datasets similar to that described in (d-e). An orange dashed line represents
the absence of possible scattering of ${\bm{q}}_{1}$ branch.

Now we turn to the second type point defect ($\\#2$). As shown in the bottom
right image in Fig. 1(d), this type of defect is centered at the S site. A
larger topography around this impurity was taken [Fig. 4(a)]. For the same
field of view, the $dI/dV$ map at the bias voltage $300$ mV was also
simultaneously taken. As shown in Fig. 4(b), the standing wave around this
defect is observed to mainly propagate along the $45^{\circ}$ direction with
respect to the lattice direction. With a nearby X-shaped defect, the QPI image
of our interest is partially affected by the X-shaped-defect-induced standing
wave. The $dI/dV$ map in Fig. 4(b) is Fourier transformed, leading to the QPI
pattern in $\bm{q}$ space in Fig. 4(c). The concentric squares appear in the
center of the QPI pattern, while the diamond and four sets of parallel lines
around Bragg peaks are absent. This defect seems only scatter the bulk band
which is similar to the X-shaped defect. From the 3D spectroscopy dataset,
figure 4(d) shows the extracted Fourier-transformed result along the
$\bm{q}$-space red line in Fig. 4(c), as a function of energy. We could
roughly distinguish two dispersed lines, labeled as ${\bm{q}}_{3}$ and
${\bm{q}}_{4}$ branches, respectively.

For the PQPI measurement, two real-space lines are chosen in Fig. 4(a) to be
away from the extra standing wave from the X-shaped defect. The $dI/dV$
spectrum was measured along the line of diagonal direction, whose Fourier
transform is performed and shown in Fig. 4(e). Within a large range of energy,
two dispersed branches (${\bm{q}}_{3}$ and ${\bm{q}}_{4}$) can be clearly
identified, confirming the two vague dispersions in Fig. 4(e). Similarly, a
series of $dI/dV$ spectra were measured along the line of lattice direction,
whose Fourier-transformed result is presented in Fig. 4(g). Different from the
QPI pattern in Fig. 3, the ${\bm{q}}_{1}$ branch of dispersion is obviously
absent, consistent with that in Fig. 4(c) and 4(f). The high SNR result in
Fig. 4(g) confirms that this defect does not scatter the electronic surface
band.

Figure 5: (a-c) DFT calculation of the slab band structure along the
M-$\Gamma$-M, M-X-M and $\Gamma$-X-$\Gamma$ directions in $\bm{k}$-space. The
blue, green, red and orange dots represent $d_{xy}$, $d_{z^{2}}$,
$d_{xz}/d_{yz}$ and $d_{x^{2}-y^{2}}$ orbitals of Zr atoms, respectively.
(d-e) Comparison of energy dispersions between experimental result and the DFT
calculation along M-$\Gamma$-M direction. (f-g) Comparison of two energy
dispersions between experimental result and the DFT calculation along
X-$\Gamma$-X direction for defect $\\#$1\. The surface band has to be shifted
$150$ mV up to match with the experimental result in (f).

With the PQPI method, electronic branches from scattering can be clearly
identified for different defects, which enables a precise extraction of
dispersions and a quantitative analysis of the electronic band structure. For
comparison, the electronic band structure of ZrSiS was calculated with a DFT
of a slab model. Figure 5(a) shows the calculated band structure along the
M-$\Gamma$-M direction in $\bm{k}$-space. The orbital projection has been
considered in the band calculation, as presented with different colored dots
in Fig. 5(a). From the DFT result, the bands near the Fermi surface are mainly
contributed by different orbital components of Zr atoms. The outer band above
the Fermi level with orange dots is mainly composed of $d_{x^{2}-y^{2}}$
components, meanwhile, the inner branch with red dots is composed of
degenerated $d_{xz}$ and $d_{yz}$ components. The QPI branch with wave vector
${\bm{q}}_{3}$ corresponds to the scattering between two internal bands
originating from $d_{xz}/d_{yz}$ orbital of Zr atoms, or two sides of the
internal square in the CCE model. The ${\bm{q}}_{4}$ branch corresponds to the
scattering between one band with $d_{xz}/d_{yz}$ orbital and another band with
$d_{x^{2}-y^{2}}$ orbital of Zr atoms. In the CCE model, it is equivalent to
the scattering from one side of the internal square to the opposite side of
the outer square. The ${\bm{q}}_{5}$ branch corresponds to the scattering
between two bands with $d_{x^{2}-y^{2}}$ orbital of Zr atoms, which is
indicated by a scattering between two sides of the outer square in the CCE
model.

For the defect $\\#1$, we extract ${\bm{q}}(E)$ from the dispersed line with
high intensity in Fig. 3(b). We made a constant energy shift of 100 meV for
all the bands to present our DFT calculation results. With this constant
shift, the extracted ${\bm{q}}(E)$ is well consistent with the
${\bm{q}}_{3}(E)$ calculated from the electronic band structure [Fig. 5(d)],
which proves that the single branch in Fig. 3(b) matches with ${\bm{q}}_{3}$,
instead of ${\bm{q}}_{4}$ or ${\bm{q}}_{5}$. For the defect $\\#2$, we extract
two dispersed branches in Fig. 4(e). In Fig. 5(e), the two extracted branches
are well consistent with ${\bm{q}}_{3}(E)$ and ${\bm{q}}_{4}(E)$ calculated
from the electronic band structure. With the decreasing energy, the amplitudes
of ${\bm{q}}_{3}(E)$ and ${\bm{q}}_{4}(E)$ increase, but at different speeds.
The dispersion of ${\bm{q}}_{3}(E)$ and ${\bm{q}}_{4}(E)$ merge around the
energy of $300$ meV above the Fermi level, related with the nodal line in this
nodal-line semimetal. The determination of the nodal line is consistent with
the result for ZrSiSe in the previous work Bu2018PRB .

For these two defects, we discover either a single branch of scattering with
wave vector ${\bm{q}}_{3}$ or two branches of scattering with wave vectors
${\bm{q}}_{3}$ and ${\bm{q}}_{4}$. This phenomenon is similar to that in
previous reports for ZrSiS and ZrSiSe butler2017quasiparticle ;
lodge2017observation ; surface_termination . The scattering of ${\bm{q}}_{5}$
is never discovered heretofore. The high SNR in our results prove that the
absence of scattering of ${\bm{q}}_{5}$ is not due to a limited SNR in the QPI
measurement, but an intrinsic property of the scattering. The appearance of a
single branch (${\bm{q}}_{3}$) or two branches (${\bm{q}}_{3}$ and
${\bm{q}}_{4}$) are also clearly distinguished for two different point
defects.

We next extract the scattering between electronic states in the surface band.
Figure 5(b) shows the calculated band structure along the M-$X$-M direction,
perpendicular to the two parallel arcs in the CCE model in $\bm{k}$-space.
From the calculation, the ${\bm{q}}_{1}$ branch happens between bands with
$d_{z^{2}}$ orbital. This surface band is a floating band, originating from
the surface-induced symmetry breaking from nonsymmorphic group $P4/nmm$ to
symmorphic wallpaper group $P4mm$. With the broken symmetry, the high band
degeneracies are not protected anymore and can be lifted, resulting in
floating or unpinned two-dimensional surface band Surface_floating . In Fig.
5(f), we extract ${\bm{q}}_{1}(E)$ from the dispersed line with high intensity
in Fig. 3(d). The main dispersion of ${\bm{q}}_{1}$ branch is linear from
$300$ meV up to $1$ V. However, the calculated surface band has to be shifted
$150$ meV up to match the experimental results, in addition to the constant
energy shift for all bands. This energy shift of 150 meV may show the
sensitivity of the floating band position with respect to the impurity
surface_potential ; vacancy-induced . The deviation of the calculated surface
band above $700$ meV may be related with a band bending effect in slab model
calculation.

In Fig. 3(d), with the high SNR in our PQPI method, a new branch of
${\bm{q}}_{6}(E)$ can be observed, which has never been reported in previous
STM experiments. The branch of ${\bm{q}}_{6}(E)$ is extracted and shown in
Fig. 5(g). After careful comparison, this branch is found to be consistent
with the scattering between bulk bands along $\Gamma X\Gamma$ direction [Fig.
5(c)]. Here the involved bulk bands bulk_band  correspond to the corners of
the inner concentric square in Fig. 2(d). Different from other ${\bm{q}}$
branches, ${\bm{q}}_{6}$ branch is related to a scattering process between
adjacent Brillouin zone (BZ), as illustrated in the Supplementary Material spm
. Normally the inter-BZ scattering is not detectable in QPI patterns. The
glide mirror symmetry however effectively reduces the unit cell by half and
expands the first BZ by two folds. Then the inter-BZ scattering of
${\bm{q}}_{6}$ becomes an intra-BZ scattering in the non-symmorphic reshaped
1st BZ, which makes this scattering detectable. A similar picture has been
applied to explain the anomalous half-missing Umklapp feature ZZ2018NatCom .

To explain the complicated defect-dependent scattering is very difficult. The
clear scattering signal, however, enables analysis of symmetry-enforced
selection rules stern2018prl . With preserved non-symmorphic symmetry for the
bulk band, the band with $d_{xz}/d_{yz}$ and $d_{x^{2}-y^{2}}$ orbital can be
characterized by an integer number of symmetry flavor $\nu=1$ and $\nu=0$,
respectively stern2018prl . The bands with different $\nu$ induce a direct
band crossing and protected nodal-line (or Dirac ring) in ZrSiX. The bulk band
scattering of ${\bm{q}}_{3}(E)$ branch happens between bands with
$d_{xz}/d_{yz}$ orbital, with $\Delta\nu=0$. The scattering of
${\bm{q}}_{4}(E)$ branch happens between bands with $d_{x^{2}-y^{2}}$ and
$d_{xz}/d_{yz}$ orbital, with $\Delta\nu=1$. Theoretically, the two branches
must be distinguished by matrix elements because ${\bm{q}}_{3}$ corresponds to
scattering on the same $\nu$ and ${\bm{q}}_{4}$ on different $\nu$
stern2018prl . For defect $\\#1$, only one ${\bm{q}}_{3}$ branch is induced,
which means that the orbital character of defect allows the scattering with
$\Delta\nu=0$, while forbids the scattering with $\Delta\nu=1$. For defect
$\\#2$, we see both ${\bm{q}}_{3}$ and ${\bm{q}}_{4}$ branches imply that the
defect does not impose a selection rule, and the defect should have a mixed
orbital character.

In the STM experiment, the tunneling current depends on the overlap between
the tip and sample wavefunctions STM_book . The tip-related effect should also
be discussed. The coexistence of ${\bm{q}}_{3}$ and ${\bm{q}}_{4}$ implies
that the tip also does not impose a selection rule. The partial overlap
between tip and sample wavefunctions is related to a nonuniversal value of
$\beta$, which is defined for the tip in Ref. stern2018prl . With a vertical
$z$-component, $d_{z^{2}}$ and $d_{xz}/d_{yz}$ orbitals of Zr atoms are prone
to be overlapped with the typical s-wave symmetric tip state
Spectral_decomposition . The overlap between $d_{x^{2}-y^{2}}$ orbital (in the
$xy$ plane) and the tip state is comparably smaller. Although with the same
$\Delta\nu=0$ as the ${\bm{q}}_{3}$ branch, the ${\bm{q}}_{5}$ branch has
never been observed. The minimum overlap between $d_{x^{2}-y^{2}}$ orbital and
the tip state may lead to a negligible signal of ${\bm{q}}_{5}$ in the QPI
result. We emphasize that although the tip-related tunneling is involved in
this detection of QPI, the observed phenomenon is robust against different
samples and tips. For example, there are always two branches of scattering
(${\bm{q}}_{3}$ and ${\bm{q}}_{4}$) for defect $\\#2$, detected by different
tips on different samples. The type of defect is the key to induce selective
scattering of electronic bands.

Although we cannot determine the different defect types yet, further
exploration of the impurity spectrum may provide extra information for a later
determination Zhou_impurity . The most frequently found defects are the
diamond-shaped Zr-site defect and the X-shaped S-site defect, which are
expected to be located within the top ZrS bilayer. When we measure the
impurity spectrum for both defects, no obvious different features can be
discerned when compared with the clean-area spectrum [see Fig. S2(a-b) in the
Supplementary Material]. For the Zr-site defect $\\#1$, in contrast, the
central impurity spectrum shows a peak feature around $-40$ meV while the
impurity spectrum at neighboring Zr site shows a peak feature around $350$ meV
[see Fig. S2(c) in the Supplementary Material]. For the S-site defect $\\#2$,
instead, the impurity spectrum at neighboring S site shows a peak feature
around $300$ meV [see Fig. S2(d) in the Supplementary Material]. Although
defect $\\#2$ looks like a S-site defect, we cannot avoid the possibility that
it is located at the Zr-site within the bottom ZrS bilayer. Please note that
the neighboring dark S sites and the four radiating dark S lines share some
similarities with that of bright Zr atoms in the topography of defect $\\#1$.
In this possible situation, the Zr-site defects within the top bilayer are
prone to scatter surface band, and the Zr-site defects within the bottom
bilayer only scatter the bulk band. In the future, scanning transmission
electron microscopy (STEM) may be applied to determine the defect types
Bu_FeSe . A comprehensive theoretical model and first-principle calculations
are also required to analyze the orbital character of different defects and
explain the defect-dependent scattering and selection rules.

##  IV Summary

PQPI is a complementary tool to analyze the single-defect-induced QPI pattern
and its energy dependence. With a single defect as the scattering center, the
QPI oscillation decays with increasing distance away from the defect. A
standard $dI/dV$ mapping of the QPI image is a $2D$ measurement within a small
area around the defect. Although a long time is required in the data-taking
procedure, the SNR of the $2D$ measurement is still relatively low. With an
anisotropic propagation of the scattering oscillations, some high symmetry
directions can be chosen in the PQPI method, along which a $1D$ $dI/dV$
measurement can be finished within a short time. Changing from a $2D$
measurement to a $1D$ measurement, we can increase the data-taking time of
every single spectrum and enhance the SNR of measured results.

In summary, we investigate single-defect-induced QPI oscillations in the
nodal-line semimetal ZrSiS. A new type of Zr-site defect is found to scatter
both the bulk band and surface floating band. With the PQPI method, clear QPI
dispersions along high symmetry directions have been clearly resolved. The
clear scattering signal enables a discussion about the non-symmorphic-
symmetry-enforced selection rules. An extra energy shift of the surface
floating band is determined and a new branch of ${\bm{q}}_{6}$ scattering is
discovered. The PQPI method can be generally applied in other complex
materials to explore the distinct interaction between a single atomic defect
and electronic states.

###### Acknowledgements.

We acknowledge and thank R. Queiroz, Z. Fang, X. Dai and H. Weng for
discussions and communications. This work was supported by the National Basic
Research Program of China (2019YFA0308602, 2016YFA0300402), the National
Natural Science Foundation of China (NSFC-11374260, 11374261), and the
Fundamental Research Funds for the Central Universities in China.

## References

  * (1) M. F. Crommie, C. P. Lutz, and D. M. Eigler,  Nature (London) 363, 524 (1993).
  * (2) Y. Hasegawa and P. Avouris,  Phys. Rev. Lett. 71, 1071 (1993).
  * (3) L. C. Davis, M. P. Everson, R. C. Jaklevic, and W. Shen,  Phys. Rev. B 43, 3821 (1991).
  * (4) P. Avouris, I. Lyo, R. E. Walkup, and Y. Hasegawa,  J. Vac. Sci. Technol. B: Microelectron. Nanometer Struct. Process. Meas. Phenom. 12, 1447 (1994).
  * (5) L. Petersen, P. T. Sprunger, P. Hofmann, E. Lægsgaard, B. G. Briner, M. Doering, H. P. Rust, A. M. Bradshaw, F. Besenbacher, and E. W. Plummer,  Phys. Rev. B 57, R6858 (1998).
  * (6) P. T. Sprunger, L. Petersen, E. W. Plummer, E. Lægsgaard, and F. Besenbacher,  Science 275, 1764 (1997).
  * (7) P. Hofmann, B. G. Briner, M. Doering, H. P. Rust, E. W. Plummer, and A. M. Bradshaw,  Phys. Rev. Lett. 79, 265 (1997).
  * (8) J. I. Pascual, G. Bihlmayer, Y. M. Koroteev, H. P. Rust, G. Ceballos, M. Hansmann, K. Horn, E. V. Chulkov, S. Blügel, P. M. Echenique, and P. Hofmann,  Phys. Rev. Lett. 93, 196802 (2004).
  * (9) A. Weismann, M. Wenderoth, S. Lounis, P. Zahn, N. Quaas, R. G. Ulbrich, P. H. Dederichs, and S. Blügel,  Science 323, 1190 (2009).
  * (10) J. E. Hoffman, K. McElroy, D. H. Lee, K. M. Lang, H. Eisaki, S. Uchida, and J. C. Davis,  Science 297, 1148 (2002).
  * (11) Q. H. Wang and D. H. Lee,  Phys. Rev. B 67, 020511 (2003).
  * (12) A. Kostin, P. O. Sprau, A. Kreisel, Y. X. Chong, A. E. Böhmer, P. C. Caneld, P. J. Hirschfeld, B. M. Andersen, and J. C. Davis,  Nat. Mater. 17, 869 (2018).
  * (13) P. Aynajian, E. H. da Silva Neto, A. Gyenis, R. E. Baumbach, J. D. Thompson, Z. Fisk, E. D. Bauer, and A. Yazdani,  Nature (London) 486, 201 (2012).
  * (14) M. P. Allan, F. Massee, D. K. Morr, J. Van Dyke, A. W. Rost, A. P. Mackenzie, C. Petrovic, and J. C. Davis,  Nat. Phys. 9, 468 (2013).
  * (15) P. Roushan, J. Seo, C. V. Parker, Y. S. Hor, D. Hsieh, D. Qian, A. Richardella, M. Z. Hasan, R. J. Cava, and A. Yazdani,  Nature (London) 460, 1106 (2009).
  * (16) T. Zhang, P. Cheng, X. Chen, J. F. Jia, X. Ma, K. He, L. Wang, H. Zhang, X. Dai, Z. Fang, X. Xie, and Q. K. Xue,  Phys. Rev. Lett. 103, 266803 (2009).
  * (17) I. Zeljkovic, Y. Okada, C. Y. Huang, R. Sankar, D. Walkup, W. Zhou, M. Serbyn, F. Chou, W. F. Tsai, H. Lin, A. Bansil, L. Fu, M. Z. Hasan, and V. Madhavan,  Nat. Phys. 10, 572 (2014).
  * (18) G. Chang, S. Y. Xu, H. Zheng, C. C. Lee, S. M. Huang, I. Belopolski, D. S. Sanchez, G. Bian, N. Alidoust, T. R. Chang, C. H. Hsu, H. T. Jeng, A. Bansil, H. Lin, and M. Z. Hasan,  Phys. Rev. Lett. 116, 066601 (2016).
  * (19) H. Inoue, A. Gyenis, Z.Wang, J. Li, S. W. Oh, S. Jiang, N. Ni, B. A. Bernevig, and A. Yazdani,  Science 351, 1184 (2016).
  * (20) H. Zheng, S. Y. Xu, G. Bian, C. Guo, G. Chang, D. S. Sanchez, I. Belopolski, C. C. Lee, S. M. Huang, X. Zhang, R. Sankar, N. Alidoust, T. R. Chang, F. Wu, T. Neupert, F. Chou, H. T. Jeng, N. Yao, A. Bansil, S. Jia, H. Lin, and M. Z. Hasan,  ACS Nano 10, 1378 (2016).
  * (21) R. Batabyal, N. Morali, N. Avraham, Y. Sun, M. Schmidt, C. Felser, A. Stern, B. Yan, and H. Beidenkopf,  Sci. Adv. 2, e1600709 (2016).
  * (22) N. Avraham, J. Reiner, A. Kumar-Nayak, N. Morali, R. Batabyal, B. Yan, and H. Beidenkopf  Adv. Mater. 30, 1707628 (2018).
  * (23) X. G. Liu, H. J. Du, J. F. Wang, M. Y. Tian, X. Sun, and B. Wang,  J. Phys.: Condens. Matter, 29, 185002 (2017).
  * (24) Ilya K. Drozdov, A. Alexandradinata, Sangjun Jeon, S. Nadj-Perge, H. W. Ji, R. J. Cava, B. A. Bernevig, and A. Yazdani,  Nat. Phys. 10, 664 (2014).
  * (25) L. Simon, C. Bena, F. Vonau, D. Aubel, H. Nasrallah, M. Habar, and J. C. Peruchetti,  Eur. Phys. J. B 69, 351 (2009).
  * (26) P. G. Derry, A. K. Mitchell, and D. E. Logan,  Phys. Rev. B 92, 035126 (2015).
  * (27) Q. Xu, Z. Song, S. Nie, H. Weng, Z. Fang, and X. Dai,  Phys. Rev. B 92, 205310 (2015).
  * (28) J. Hu, Z. Tang, J. Liu, X. Liu, Y. Zhu, D. Graf, K. Myhro, S. Tran, C. N. Lau, J. Wei, and Z. Mao,  Phys. Rev. Lett. 117, 016602 (2016).
  * (29) L. M. Schoop, M. N. Ali, C. Straßer, A. Topp, A. Varykhalov, D. Marchenko, V. Duppel, S. S. P. Parkin, B. V. Lotsch, and C. R. Ast,  Nat. Commun. 7, 11696 (2016).
  * (30) A. Topp, R. Queiroz, A. Grüneis, L. Müchler, A. W. Rost, A. Varykhalov, D. Marchenko, M. Krivenkov, F. Rodolakis, J. L. McChesney, B. V. Lotsch, L. M. Schoop, and C. R. Ast,  Phys. Rev. X 7, 041073 (2017).
  * (31) Z. Zhu, T. R. Chang, C. Y. Huang, H. Pan, X. A. Nie, X. Z. Wang, Z. T. Jin, S. Y. Xu, S. M. Huang, D. D. Guan, S. Wang, Y. Y. Li, C. Liu, D. Qian, W. Ku, F. Song, H. Lin, H. Zheng, and J. F. Jia,  Nat. Commun. 9, 4153 (2018).
  * (32) K. Bu, Y. Fei, W. Zhang, Y. Zheng, J. Wu, F. Chen, X. Luo, Y. Sun, Q. Xu, X. Dai, and Y. Yin,  Phys. Rev. B 98, 115127 (2018).
  * (33) M. S. Lodge, G. Chang, C. Y. Huang, B. Singh, J. Hellerstedt, M. T. Edmonds, D. Kaczorowski, M. M. Hosen, M. Neupane, H. Lin, M. S. Fuhrer, B. Weber, and M. Ishigami,  Nano. Lett. 17, 7213 (2017).
  * (34) C. J. Butler, Y. M. Wu, C. R. Hsing, Y. Tseng, R. Sankar, C. M. Wei, F. C. Chou, and M. T. Lin,  Phys. Rev. B 96, 195125 (2017).
  * (35) C. C. Su, C. S. Li, T. C. Wang, S. Y. Guan, R. Sankar, F. Chou, C. S. Chang, W. L. Lee, G. Y. Guo, and T. M. Chuang,  New J. Phys. 20, 103025 (2018).
  * (36) R. Sankar, G. Peramaiyan, I. P. Muthuselvam, C. J. Butler, K. Dimitri, M. Neupane, G. N. Rao, M. T. Lin, and F. C. Chou,  Sci. Rep. 7, 40603 (2017).
  * (37) Y. Zheng, Y. Fei, K. Bu, W. Zhang, Y. Ding, X. Zhou, J. E. Hoffman, and Y. Yin,  Sci. Rep. 7, 8059 (2017).
  * (38) Y. Fei, K. L. Bu, W. H. Zhang, Y. Zheng, X. Sun, Y. Ding, X. J. Zhou, and Y. Yin,  Sci. China Phys. Mech. Astron. 61, 127404 (2018).
  * (39) Y. Fei, Y. Zheng, K. L. Bu, W. H. Zhang, Y. Ding, X. J. Zhou, and Y. Yin,  Sci. China Phys. Mech. Astron. 63, 227411 (2020).
  * (40) G. Kresse and J. Furthmüller,  Phys. Rev. B 54, 11169 (1996).
  * (41) G. Kresse and D. Joubert,  Phys. Rev. B 59, 1758 (1999).
  * (42) J. P. Perdew, K. Burke, and M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996).
  * (43) N. Marzari and D. Vanderbilt,  Phys. Rev. B 56, 12847 (1997).
  * (44) C. Wang and T. Hughbanks,  Inorg. Chem. 34, 5524 (1995).
  * (45) M. J. Lawler, K. Fujita, Jhinhwan Lee, A. R. Schmidt, Y. Kohsaka, Chung Koo Kim, H. Eisaki, S. Uchida, J. C. Davis, J. P. Sethna, and Eun-Ah Kim  Nature (London) 466, 347 (2010).
  * (46) I. Zeljkovic, E. J. Main, T. L. Williams, M. C. Boyer, K. Chatterjee, W. D. Wise, Y. Yin, M. Zech, A. Pivonka, T. Kondo, T. Takeuchi, H. Ikuta, J. Wen, Z. Xu, G. D. Gu, E. W. Hudson, and J. E. Hoffman,  Nat. Mater. 11, 585 (2012).
  * (47) C. Mann, D. West, I. Miotkowski, Y. P. Chen, S. Zhang, and C. K. Shih,  Nat. Commun. 4, 2277 (2013).
  * (48) T. Förster, P. Krüger, and M. Rohlfing,  Phys. Rev. B 91, 035313 (2015).
  * (49) B. B. Fu, C. J. Yi, T. T. Zhang, M. Caputo, J. Z. Ma, X. Gao, B. Q. Lv, L. Y. Kong, Y. B. Huang, P. Richard, M. Shi, V. N. Strocov, C. Fang, H. M. Weng, Y. G. Shi, T. Qian, and H. Ding,  Sci. Adv. 5, eaau6459 (2019).
  * (50) See Supplementary Material at http://
  * (51) R. Queiroz, and A. Stern,  Phys. Rev. Lett. 121, 176401 (2018).
  * (52) C. J. Chen, Introduction to Scanning Tunneling Microscopy (Oxford University Press, 1993).
  * (53) J. Nieminen, I. Suominen, R. S. Markiewicz, H. Lin, and A. Bansil,  Phys. Rev. B 80, 134509 (2009).
  * (54) T. Zhou, W. Chen, Y. Gao and Z. D. Wang  Phys. Rev. B 100, 205119 (2019).
  * (55) K. L. Bu, B. Wang, W. H. Zhang, Y. Fei, Y. Zheng, F. Z. Ai, Z. X. Wu, Q. S. Wang, H. L. Wo, J. Zhao, C. H. Jin and Y. Yin,  Phys. Rev. B 100, 155127 (2019).


Creating Teams of Simple Agents for Specified Tasks:
A Computational Complexity Perspective

Todd Wareham

Department of Computer Science

Memorial University of Newfoundland

St. John’s, NL Canada

(Email<EMAIL_ADDRESS>
> Abstract: Teams of interacting and co-operating agents have been proposed as
> an efficient and robust alternative to monolithic centralized control for
> carrying out specified tasks in a variety of applications. A number of
> different team and agent architectures have been investigated, e.g., teams
> based on single vs multiple behaviorally-distinct types of agents
> (homogeneous vs heterogeneous teams), simple vs complex agents, direct vs
> indirect agent-to-agent communication. A consensus is emerging that (1)
> heterogeneous teams composed of simple agents that communicate indirectly
> are preferable and (2) automated methods for verifying and designing such
> teams are necessary. In this paper, we use computational complexity analysis
> to assess viable algorithmic options for such automated methods for various
> types of teams. Building on recent complexity analyses addressing related
> questions in swarm robotics, we prove that automated team verification and
> design are by large both exact and approximate polynomial-time intractable
> in general for the most basic types of homogeneous and heterogeneous teams
> consisting of simple agents that communicate indirectly. Our results suggest
> that tractability for these problems must be sought relative to additional
> restrictions on teams, agents, operating environments, and tasks.

##  1 Introduction

Teams of interacting and co-operating agents have been proposed as an
efficient and robust alternative to monolithic centralized control for
carrying out specified tasks in a variety of applications. A number of
different team and agent architectures have been investigated. Three
dimensions of these architectures are of particular importance:

  1. 1.

Should teams consist of a single type of agent or multiple types of agents,
i.e., should teams be homogeneous or heterogeneous?

  2. 2.

Should individual agents have simple or complex control mechanisms, i.e.,
simple reflex or model- / goal- / utility-based agents [1, Section 2.4]?

  3. 3.

Should individual agents communicate directly with each other by agent-to-
agent messages or indirectly via their sensed presences and/or environmental
modifications, i.e., via stigmergy [2]?

Based on the experience gained with various proof-of-concept experiments and
implementations, the consensus is emerging that (1) heterogeneous teams
composed of relatively simple agents that communicate indirectly are
preferable [3, 4, 5] and (2) automated methods for verifying and designing
such teams are necessary [6, 7].

A natural question at this point is what algorithmic options are and are not
available for the efficient verification and design of teams relative to the
three dimensions listed above. In this paper, we give some initial answers to
this question, building on recent work [8, 9, 10, 11, 12, 13] addressing
related questions in swarm robotics for distributed construction.111 Though
there has been other complexity-theoretic work on individual agent
verification and design [14, 15, 16], the models of agents and environments
used were too abstract to allow examination of agent teams, simple reflex
agents, or stigmergy.  In particular, we give proofs (modified from several
given previously in [9, 10, 11]222 These modifications are described along
with the proofs of our results in the appendix. ) which demonstrate that the
problems of team verification and design are by large both exact and
approximate polynomial-time intractable in general relative to the most basic
types of homogeneous and heterogeneous teams consisting of simple reflex
agents that do not use stigmergy. This in turn suggests that tractability must
be sought relative to additional restrictions on teams, agents, operating
environments, and tasks.

##  2 Methods

In this section, we first review the basic entities in our model of task
performance by robot teams — namely, environments, individual robots, robot
teams, and tasks (with the last of these being new to this paper). Though this
is a basic model in which robots sense and move without uncertainty in a
discrete and synchronous manner in a 2D grid-based environment, it is flexible
enough to allow investigations along the three team and agent architectural
dimensions listed in the introduction. In the interests of concision, most of
this review is the “short form” given in [8, 10]; readers wishing more details
should consult [9, 11]. This will be followed by formalization of
computational problems within this model corresponding to various types of
robot team verification and design.

The basic entities in our model are as follows:

  * •

Environments: Our robots operate in a finite 2D square-based environment $E$
in which each square is either a freespace (which a robot can occupy or move
through) or an obstacle (which a robot cannot occupy or move through), and has
a square-type, e.g., grass, gravel, wall, drawn from a set $E_{T}$. Let
$E_{i,j}$ denote the square that is in the $i$th column and $j$th row of $E$
such that $E_{1,1}$ is the square in the southwest-most corner of $E$.

  * •

Robots: Each robot occupies a square in $E$ and in a basic movement-action can
either move exactly one square to the north, south, east or west of its
current position or elect to stay at its current position. Each robot has a
sensing-distance bound $r$ such that the robot can sense the type of the
square at any position within Manhattan distance $r\geq 0$ of the robot’s
current position (with $r=0$ corresponding to the square on which the robot is
standing). These square-types are accessible via predicates of the form
$enval(e,pos)$ which returns $True$ if the square at position $pos$ has type
$e\in E_{T}\cup\\{e_{robot}\\}$ (with the sensor returning $e_{robot}$ if a
robot is occupying square $pos$) and $False$ otherwise, where a position $pos$
is specified in terms of a pair $(x,y)$ specifying an environment-square
$E_{i+x,j+y}$ if the robot is currently occupying $E_{i,j}$. Each robot can
change the type of the square at any position within Manhattan distance one of
the robot’s current position to type $e$ via predicates of the form
$enmod(e,pos)$ where $pos$ is specified as for $enval()$.

Each robot has a finite-state controller and is hence known as a Finite-State
Robot (FSR). Each such controller consists of a set $Q$ of states linked by
transitions, where each transition $(q,f,x,dir,q^{\prime})$ between states $q$
and $q^{\prime}$ has a propositional logic trigger-formula $f$, an environment
modification specification $x$, and a movement-direction
$dir\in\\{goNorth,goSouth,goEast,goWest,stay\\}$. Trigger-formulas and
modification specifications are typically stated in terms of predicates
$enval()$ and $enmod()$, respectively. Both of these specifications can also
be stated as a special symbol $*$, which is interpreted as follows: If $f\neq
x\neq*$ and the transition’s trigger-formula evaluates to $True$, i.e., the
transition is enabled, this causes the environment-modification specified by
$x$ to occur, the robot to move one square in direction $dir$, and the robot’s
state to change from $q$ to $q^{\prime}$. If $f=*$, the transition executes if
no other transition executes (making this in effect the default transition);
if $x=*$, no environment-modification is made.

The transitions in the FSR described above can be viewed as condition-action
rules within the agent framework given in [1, Section 2.4] such that the $*$
trigger formula on a transition leaving state $q$ can be viewed as the
negation of the disjunction of the trigger formulas of all other transitions
leaving $q$. Given this, single- and multi-state FSR correspond to simple
reflex and model-based reflex agents agents, where the model and UPDATE-
STATE() function in [1, Section 2.4.3] are implicit in the transitions between
states. Note that actions correspond to FSR movements and environmental
modifications as required, and are flexible enough to allow situations in
which FSR movements or environment modifications do not occur.

  * •

Robot teams: A team $T$ consists of a set of the robots described above, where
there may be more than one robot with the same controller on a team; as such,
we allow both homogeneous and heterogeneous teams. Let $T_{i}$ denote the
$i$th robot on the team. Each square in $E$ can hold at most one member of
$T$; if at any point in the execution of a task two robots in a team attempt
to occupy or modify the same free space or a robot attempts to occupy the same
space as an obstacle, the execution terminates and is considered unsuccessful.
A positioning of $T$ in $E$ is an assignment of the robots in $T$ to a set of
$|T|$ squares in $E$. For simplicity, team members move synchronously, and
once movement is triggered, it is atomic in the sense that the specified
movement is completed.

Note that robots in our teams do not communicate with each other directly —
rather, they can communicate with each other indirectly through their sensed
presences in and changes they make to the environment, i.e., via stigmergy
[2]. In the remainder of the paper, we will find it useful to distinguish
these two types of communication, which will be denoted as agent- and
environment-based stigmergy, respectively.

  * •

Tasks: Tasks are specified in terms of a desired set of environment square-
values, robot positions, and/or robot internal states, e.g., a $3\times 3$
square of square-type $e_{X}$ has been created at a particular location in the
environment, all robots in a team are in state $q_{F}$ and located on the
eastmost edge of the environment. Such is a specification will be denoted as a
task’s target configuration. We will assume that for each task $Tsk$ and an
environment $E$ in which a robot team $T$ is operating, $E$ can be checked for
the target configuration associated with $Tsk$ in time polynomial in the size
of $E$.

We use the notions of deterministic robot and team operation proposed in [8,
11] (i.e., requiring that at any time as the team operates in an environment,
all transitions enabled in a robot relative to the current state of that robot
perform the same environment modifications and progress to the same next
state). Given this, an individual FSR is not itself deterministic but rather
the operation of that FSR is deterministic in the context of a particular FSR
team operating in a particular environment.

Let us now consider the team verification design problems that we will analyze
in this paper, starting with verification.

Team / environment verification for Task $Tsk$ (TeamEnvVer)
Input: An environment $E$ based on square-type set $E_{T}$, a an FSR team $T$,
an initial positioning $p_{I}$ of $T$ in $E$, and an integer $\\#_{ec}\geq 0$.
Question: Does $T$ started at $p_{I}$ perform task $Tsk$ while making at most
$\\#_{ec}$ square-type changes in $E$?

We will consider two types of robot team design. Both types of design are done
relative to a given design library. In the first case, $L$ consists of
transition templates of the form $(q,f,x,move,q^{\prime})$ which are used to
construct FSR controllers from a specified set of states by instantiating
transition templates relative to those states. Note it may be the case that
$q=q^{\prime}$ in such a construction, i.e., a transition may loop back on the
same state.

Controller design by library selection for task $Tsk$ (ContDesLS)
Input: An environment $E$ based on square-type set $E_{T}$, a requested team-
size $|T|$, an initial positioning $p_{I}$ of $T$ in $E$, a transition
template library $L$, and integers $r\geq 0$, $|Q|\geq 1$, $d\geq 1$, $h\geq
1$, and $\\#_{ec}\geq 0$.
Output: A controller-set $C$ in which each controller has sensory radius $r$,
at most $|Q|$ states, and at most $d$ transitions chosen from $L$ out of any
state such that an FSR team $T$ based on $h$ controllers from $C$ started at
$p_{I}$ performs $Tsk$ while making at most $\\#_{ec}$ square-type changes in
$E$, if such a $C$ and $T$ exists, and special symbol $\bot$ otherwise.

In the second case, $L$ consists of complete FSR controllers.

Team design by library selection for task $Tsk$ (TeamDesLS)
Input: An environment $E$ based on square-type set $E_{T}$, a requested team-
size $|T|$, an FSR library $L$, an initial region $E_{I}$ of size $T$ in $E$,
and integers $h\geq 1$ and $\\#_{ec}\geq 0$.
$E$.
Output: An FSR team $T$ based on $h$ robots selected from $L$ such that $T$
started in $E_{I}$ performs $Tsk$ while making at most $\\#_{ec}$ square-type
changes in $E$, if such a a $T$ exists, and special symbol $\bot$ otherwise.

We will subsequently analyze these problems relative to three parameters:

  1. 1.

The number of different types of FSR controllers in a team ($h$);

  2. 2.

The maximum number of states in the robots in a team ($|Q|$); and

  3. 3.

The maximum allowable number of environmental changes made by a team in
performing specified tasks ($\\#_{ec}$).

Different values of these parameters will allow us to investigate the question
posed in the introduction — namely, the effects on the computational
difficulty of team verification and design when using (1) teams based on
single ($h=1$) and multiple ($h>1$) types of FSR controllers (i.e.,
homogeneous and heterogeneous teams), (2) simple reflex ($|Q|=1$) and model-
based reflex ($|Q|>1$) agents, and (3) agents that do ($\\#_{ec}>0$) and do
not ($\\#_{ec}=0$) use environment-based stigmergy in the performance of
tasks.

##  3 Results

In this section, we analyze the computational difficulty of our team
verification and design problems relative to the three restrictions proposed
at the end of the previous section. We evaluate this difficulty relative to
several criteria of efficient algorithm operation using two standard
techniques – namely, proving tractability by giving algorithms and
intractability by giving reductions from known intractable problems (see [17,
18] for details of these techniques). All proofs of results are relegated to
the appendix.

Let us first consider exact polynomial-time solvability. An exact polynomial-
time algorithm is an algorithm which always produces the correct output for a
given input and whose runtime is asymptotically upper-bounded, i.e., upper-
bounded when $|x|$ goes to infinity, by $c|x|^{c^{\prime}}$, where $|x|$ is
the size of the input $x$ and $c$ and $c^{\prime}$ are constants. A problem
that has a polynomial-time algorithm is said to be polynomial-time tractable.
Polynomial-time tractability is desirable because runtimes increase slowly as
input size increases, and hence allow the solution of larger inputs. We start
with team verification.333 Results A, B, D, and E hold relative to some
combination of the conjectures $P\neq NP$ and $P=BPP$, which though unproven
are widely believed to be true within computer science [19, 20].

Result A (Modified from Lemma 4 in the Supplementary Materials of [11]):
TeamEnvVer is not exact polynomial-time tractable when $|Q|=1$, $h=1$, and
$\\#_{ec}=0$. Moreover, this intractability holds for any version of
TeamEnvVer when $h>1$, $|Q|>1$, and $\\#_{ec}>0$.

This result demonstrates that verification is polynomial-time intractable in
general for homogeneous teams of simple reflex agents that do not use
environment-based stigmergy. This in turn motivates the notion (introduced in
[11]) of $(c_{1},c_{2})$-completability, which requires that each robot team
complete its task within $c_{1}(|E|+|Q|)^{c_{2}}$ timesteps for constants
$c_{1}$ and $c_{2}$.444 For technical reasons that are described in detail in
[10], $c_{1}$ and $c_{2}$ are not part of the problem input but are specified
beforehand To ensure generous but still low-order polynomial team runtime
bounds, we will assume that $c_{1}=10$ and $c_{2}=3$.  Let the versions of
ContDesLS and TeamDesLS with this completability restriction be denoted by
ContDesLSres and TeamDesLSres. It turns out that this restriction does not
always help in general.

Result B (Modified from Lemma 5 in [11]): ContDesLSres is not exact
polynomial-time tractable when $|Q|=1$, $h=1$, and $\\#_{ec}=0$. Moreover,
this intractability holds for any version of ContDesLSres when $h>1$, $|Q|>1$,
and $\\#_{ec}>0$.

Result C (Modified from Result A in [9]): TeamDesLSres is exact polynomial-
time tractable when $|Q|\geq 1$, $h=1$, and $\\#_{ec}\geq 0$.

Result D (Modified from Lemma A.7 in [10]): TeamDesLSres is not exact
polynomial-time tractable when $|Q|=1$, $h=2$, and $\\#_{ec}=0$. Moreover,
this intractability holds for any version of TeamDesLSres when $h>2$, $|Q|\geq
1$, and $\\#_{ec}\geq 0$.

The above demonstrates that (1) restricted controller design (like
verification) is polynomial-time intractable in general for homogeneous teams
of simple reflex agents that do not use environment-based stigmergy (Result B)
and (2) though restricted team design is polynomial-time tractable for any
type of homogeneous team (Result C), it is polynomial-time intractable in
general for the simplest heterogeneous teams based on simple reflex agents
that do not use environment-based stigmergy (Result D).

Let us now consider polynomial-time approximate solvability. This type of
solvability may be acceptable in situations where always getting the correct
output for an input is not required. Three popular types of polynomial-time
approximation algorithms are:

  1. 1.

algorithms that always run in polynomial time but are frequently correct in
that they produce the correct output for a given input in all but a small
number of cases (i.e., the number of errors for input size $n$ is bounded by
function $err(n)$) [21];

  2. 2.

algorithms that always run in polynomial time but are frequently correct in
that they produce the correct output for a given input with high probability
[22]; and

  3. 3.

algorithms that run in polynomial time with high probability but are always
correct [23].

Algorithms of type (2) are of particular interest as they include evolutionary
algorithms. Unfortunately, none of these options are in general open to us
either courtesy of the following result.

Result E (Modified from Results A.4 and A.5 in [10]): None of the intractable
versions of TeamEnvVer, ContDesLSres, or TeamDesLSres described in Results A,
B, and D are polynomial-time approximable in senses (1–3).

Note that all intractability results above hold relative to the simplest types
of tasks (do some subset of the robots in team $T$ reach particular positions
in $E$?) and (in the case of ContDesLSres and TeamDesLSres) the most
restrictive type of completability, i.e., $(1,1)$-completability.

##  4 Discussion

In the previous section, we demonstrated that agent team verification and two
types of design by library selection (agent controller and team) are by large
both exact (Results A, B, and D) and approximate (Result E) polynomial-time
intractable in general relative to the simplest possible types of teams and
agents, i.e., homogeneous teams consisting of simple reflex agents that do not
use environment-based stigmergy. Even in the one case where we have
polynomial-time tractability (Result C; homogeneous team design by library
selection), intractability asserts itself when as few as two types of agents
co-exist on a team (Result D).555 It is tempting to think that polynomial-time
intractability of team verification and design follows from the well-known
combinatorial explosion in the number of possible team states and designs in
multi-agent systems. However, there are many examples of problems with such
exponential-size search spaces that are nonetheless solvable in polynomial
time by algorithms that exploit structure in those spaces, e.g., Minimum
spanning tree [24, Chapter 23]. Hence, definitive proof of polynomial-time
intractability requires proofs such as those given here.  That the various
intractabilities we have demonstrated cannot be vanquished by invoking
verification and design relative to teams with higher values of $h$, $|Q|$,
and $\\#_{ec}$ suggests that we are seeing a tight frontier of tractability
[18, Section 4.1] relative to the lowest possible values of these
parameters.666 It is also intriguing that agent-based stigmergy is critical to
some (Results A and D) but not all (Result B) of our intractability results.
It is all too often assumed in discussions about stigmergy that environmental
modifications (including “smart” materials [25]) are key. Our proofs suggest
that in certain situations, environmental modifications are unnecessary given
sufficiently large and appropriately-structured groups of mobile agents. As
large groups of agents are of increasing interest in certain applications,
this warrants further investigation.

As all of our results are derived relative to a simplified team operation
model in which deterministic agents operate in a synchronous and discrete
manner within a 2D grid-based environment, these results are not immediately
applicable to probabilistic agents that operate in an asynchronous and
continuous manner in the real world.777 Results like ours do nonetheless have
a surprising generality; the interested reader is referred to Sections 5 and
6.2 of [9] and Section 5 of [11] for details.  That being said, our results do
for now offer some reasons for real-world roboticists to be cautious. In
particular, the fact that team verification and design are polynomial-time
intractable even when agent motion and sensing are error-free and occur in
completely-observable environments hints that there may be additional sources
of computational difficulty in these problems that are not associated with
agent motion and sensing under partial observability and uncertainty [26].
These sources should be acknowledged and investigated, particularly if team
verification and design must behave both efficiently and correctly when fully
automated without human oversight.

This last point highlights a crucial caveat when interpreting our
intractability results – namely, these results hold relative to a simplified
model of agent team operation and, perhaps more importantly, restrictions on
the values of $h$, $|Q|$, and $\\#_{ec}$. Our frequent proviso that
intractability results apply “in general” was not mere rhetoric — tractability
may well hold for the cases we considered when additional restrictions are in
place; then again, it may not. In either case, this must be determined by
future complexity analyses. Hence, our results should be seen not as final
statements on the intractability of team verification and design but rather as
interim guidelines suggesting where in the universe of restriction
possibilities tractability does and does not hold.

Given the above, future research into team verification and design should
perhaps more closely incorporate computational complexity analyses like those
given here. Such research could initially focus on more fully characterizing
those combinations of restrictions that do and do not render team verification
and design tractable. Such work has already been started [8, 9, 10, 12, 11,
13] for team verification and design in swarm robotics relative a variety of
restrictions using more advanced analysis techniques (e.g., parameterized
complexity analysis [27]). Additional restrictions of particular interest here
are those that “break” the reductions underlying our intractability results,
e.g., restrictions on the degree and type of structure encountered by agents
in their environments (including the presences of other agents). Once these
initial intractability maps have been derived for simplified team operation
models, they should be extended to more realistic models incorporating
stochasticity and uncertainty. Part of this can be done by using complexity
analysis techniques that explicitly incorporate stochasticity [22, 28, 29].
Complexity-based frameworks that incrementally build on simplified operation
models in a systematic and principled manner to create results applicable to
more realistic models (analogous to those developed in linguistics [30] and
cognitive science [31]) may also be of use in this endeavour.

## Acknowledgments

Funding for this work was provided by a National Science and Engineering
Research Council (NSERC) Discovery Grant to TW (grant # 228105-2015).

## References

  * [1] S. Russell and P. Norvig, _Artificial Intelligence: A Modern Approach_ , 3rd ed. Prentice Hall, 2010.
  * [2] E. Bonabeau, M. Dorigo, and G. Theraulaz, _Swarm intelligence: From natural to artificial systems_. Oxford University Press, 1999.
  * [3] M. Dorigo, D. Floreano, L. M. Gambardella, F. Mondada, S. Nolfi, T. Baaboura, M. Birattari, M. Bonani, M. Brambilla, A. Brutschy _et al._ , “Swarmanoid: A novel concept for the study of heterogeneous robotic swarms,” _IEEE Robotics & Automation Magazine_, vol. 20, no. 4, pp. 60–71, 2013.
  * [4] D. Kengyel, H. Hamann, P. Zahadat, G. Radspieler, F. Wotawa, and T. Schmickl, “Potential of heterogeneity in collective behaviors: A case study on heterogeneous swarms,” in _International Conference on Principles and Practice of Multi-agent Systems_. Springer, 2015, pp. 201–217.
  * [5] S. Schmid, D. Schraudner, and A. Harth, “Performance comparison of simple reflex agents using stigmergy with model-based agents in self-organizing transportation,” in _2021 IEEE International Conference on Autonomic Computing and Self-Organizing Systems Companion (ACSOS-C)_. Los Alamitos, CA: IEEE Press, 2021, pp. 93–98.
  * [6] A. Bagnato, R. K. Bíró, D. Bonino, C. Pastrone, W. Elmenreich, R. Reiners, M. Schranz, and E. Arnautovic, “Designing swarms of cyber-physical systems: The H2020 CPSwarm project,” in _Proceedings of the Computing Frontiers Conference_ , 2017, pp. 305–312.
  * [7] M. Dorigo, G. Theraulaz, and V. Trianni, “Reflections on the future of swarm robotics,” _Science Robotics_ , vol. 5, no. 49, p. eabe4385, 2020.
  * [8] M. Timmar and T. Wareham, “The computational complexity of controller-environment co-design using library selection for distributed construction,” in _Distributed Autonomous Robotic Systems: The 14th International Symposium_ , ser. Springer Proceedings in Advanced Robotics, N. Correll, M. Schwager, and M. Otte, Eds., vol. 9. Springer Nature Switzerland AG, 2019, pp. 51–63.
  * [9] T. Wareham, “Designing robot teams for distributed construction, repair, and maintenance,” _ACM Transactions on Autonomous and Adaptive Systems_ , vol. 14(1), pp. 2:1–2:29, 2019.
  * [10] T. Wareham, R. de Haan, A. Vardy, and I. van Rooij, “Swarm control for distributed construction: A computational complexity perspective,” _ACM Transactions on Human-Robot Interaction_ , Submitted.
  * [11] T. Wareham and A. Vardy, “Putting it together: The computational complexity of designing robot controllers and environments for distributed construction,” _Swarm Intelligence_ , vol. 12, no. 2, pp. 111–128, 2018.
  * [12] ——, “Viable algorithmic options for designing reactive robot swarms,” _ACM Transactions on Autonomous Adaptive Systems_ , vol. 13, no. 1, pp. 5:1–5:23, 2018.
  * [13] ——, “The computational complexity of designing scalar-field sensing robot teams and environments for distributed construction,” in _2021 IEEE International Conference on Autonomic Computing and Self-Organizing Systems Companion (ACSOS-C)_. Los Alamitos, CA: IEEE Press, 2021, pp. 232–237.
  * [14] P. E. Dunne, M. Laurence, and M. Wooldridge, “Complexity results for agent design,” _Annals of Mathematics, Computing & Teleinformatics_, vol. 1, no. 1, pp. 19–36, 2003.
  * [15] I. A. Stewart, “The complexity of achievement and maintenance problems in agent-based systems,” _Artificial Intelligence_ , vol. 2, no. 146, pp. 175–191, 2003.
  * [16] M. Wooldridge and P. E. Dunne, “The computational complexity of agent verification,” in _Intelligent Agents VIII_. Springer, 2002, pp. 115–127.
  * [17] S. Arora and B. Barak, _Computational complexity: A modern approach_. Cambridge University Press, 2009\.
  * [18] M. R. Garey and D. S. Johnson, _Computers and Intractability_. W.H. Freeman, 1979.
  * [19] L. Fortnow, “The Status of the P Versus NP Problem,” _Communications of the ACM_ , vol. 52, no. 9, pp. 78–86, 2009.
  * [20] A. Wigderson, “P, NP and mathematics — A computational complexity perspective,” in _Proceedings of ICM 2006: Volume I_. Zurich: EMS Publishing House, 2007, pp. 665–712.
  * [21] L. A. Hemaspaandra and R. Williams, “Complexity Theory Column 76: An atypical survey of typical-case heuristic algorithms,” _ACM SIGACT News_ , vol. 43, no. 4, pp. 70–89, 2012.
  * [22] R. Motwani and P. Raghavan, _Randomized Algorithms_. Chapman & Hall/CRC, 2010.
  * [23] J. Gill, “Computational complexity of probabilistic Turing machines,” _SIAM Journal on Computing_ , vol. 6, no. 4, pp. 675–695, 1977.
  * [24] T. Cormen, C. Leiserson, R. Rivest, and C. Stein, _Introduction to Algorithms_ , 3rd ed. Boston, MA: The MIT Press, 2009.
  * [25] J. Werfel and R. Nagpal, “Extended stigmergy in collective construction,” _IEEE Intelligent Systems_ , vol. 21, no. 2, pp. 20–28, 2006.
  * [26] S. M. Hezavehi, D. Weyns, P. Avgeriou, R. Calinescu, R. Mirandola, and D. Perez-Palacin, “Uncertainty in self-adaptive systems: A research community perspective,” _ACM Transactions on Autonomous and Adaptive Systems_ , vol. 15, no. 4, pp. 1–36, 2021.
  * [27] R. G. Downey and M. R. Fellows, _Fundamentals of Parameterized Complexity_. Berlin: Springer, 2013.
  * [28] N. Donselaar, “Probabilistic parameterized polynomial time,” in _International Conference on Current Trends in Theory and Practice of Informatics_. Springer, 2019, pp. 179–191.
  * [29] J. A. Montoya and M. Müller, “Parameterized random complexity,” _Theory of Computing Systems_ , vol. 52, no. 2, pp. 221–270, 2013.
  * [30] E. S. Ristad, _The Language Complexity Game_. MIT Press, 1993.
  * [31] I. van Rooij, M. Blokpoel, J. Kwisthout, and T. Wareham, _Cognition and Intractability: A Guide to Classical and Parameterized Complexity Analysis_. Cambridge, UK: Cambridge University Press, 2019.

##  Appendix A Proofs of Results

All of our intractability results will be derived using polynomial-time
reductions from the following problems:

3-Satisfiability (3SAT) [18, Problem LO2]
Input: A set $U$ of variables, a set $C$ of disjunctive clauses over $U$ such
that each clause $c\in C$ has $|c|=3$.
Question: Is there a satisfying truth assignment for $C$?

Dominating set [18, Problem GT2]
Input: An undirected graph $G=(V,E)$ and a positive integer $k$.
Question: Does $G$ contain a dominating set of size $k$, i.e., is there a
subset $V^{\prime}\subseteq V$, $|V^{\prime}|=k$, such that for all $v\in V$,
either $v\in V^{\prime}$ or there is at least one $v^{\prime}\in V^{\prime}$
such that $(v,v^{\prime})\in E$?

For each vertex $v\in V$, let the complete neighbourhood $N_{C}(v)$ of $v$ be
the set composed of $v$ and the set of all vertices in $G$ that are adjacent
to $v$ by a single edge, i.e., $v\cup\\{u~{}|~{}u~{}\in
V~{}\rm{and}~{}(u,v)\in E\\}$. We assume below for each instance of Dominating
set an arbitrary ordering on the vertices of $V$ such that
$V=\\{v_{1},v_{2},\ldots,v_{|V|}\\}$; we also assume analogous orderings for
the variables and clauses of each instance of 3-Satisfiability such that
$U=\\{u_{1},u_{2},\ldots,u_{|U|}\\}$ and $C=\\{c_{1},c_{2},\ldots,c_{|C|}\\}$.

For technical reasons [17, 18], our intractability results are initially
derived relative to decision versions of our problems, i.e., problems whose
answers are either “Yes” or “No”. Problem TeamEnvVer is already phrased as a
decision problem. The decision versions of ContDes:Sres and TeamDesLSres
(denoted by ContDesLSres,D and TeamDesLSres,D, respectively) ask if the
structures requested in each problem ($(C,T)$ and $T$, respectively) exist or
not. The following lemma (based on the observation that any algorithm for non-
decision problem X can be used to solve XD) will be useful below in
transferring results from decision problems to their associated non-decision
problems.

###### Lemma 1

If XD is not solvable in polynomial time relative to conjecture C then X is
not solvable in polynomial time relative to conjecture C.

All proofs of results given here are modifications of proofs given previously
in [9, 10, 11]. It is thus appropriate to describe the nature of these
modifications, starting with an overview. All previous proofs had $|Q|=1$, so
modifications were made to ensure

  1. 1.

the initial values of $h$ and $\\#_{ec}$ stated in the first part of each
result, and

  2. 2.

the higher values of $|Q|$, $h$, and $\\#_{ec}$ stated in the second part of
each result.

As previous work was done with respect to distributed constructions tasks,
$\\#_{ec}=0$ was easy to ensure in the modifications once we introduced our
task model (which was new to the current paper) and focused on tasks that
involve achieving specified environment / agent / state configurations — we
just needed robots to get to specified environmental positions where we
previously had them also building structures in those positions. Other result-
specific modifications are summarized below in text prior to each result’s
proof.

We now present the proofs of results in our paper, starting with Result A.
Part (1) of this result requires relatively straightforward modification of
environments and FSR transitions to convert the previous $h=5$ construction in
[11] to the needed $h=1$ construction while preserving $|Q|=1$. Part (2)
requires relatively straightforward modification to allow $|Q|>1$ and $h>1$.

Result A (Modified from Lemma 4 in the Supplementary Materials of [11]):
TeamEnvVer is not exact polynomial-time tractable when $|Q|=1$, $h=1$, and
$\\#_{ec}=0$ unless $P=NP$. Moreover, this intractability holds for any
version of TeamEnvVer when $h>1$, $|Q|>1$, and $\\#_{ec}>0$ for fixed values
of $h$, $|Q|$, and $\\#_{ec}$.

Proof: Lemma 4 in [11] gives a reduction from 3SAT to problem ContEnvVersyn,
which is our problem TeamEnvVer without the restriction on $\\#_{ec}$. In the
instance of ContEnvVersyn created by this reduction, all possible truth
assignments to the variables in $U$ in the given instance of 3SAT are
generated one at a time by the movements of a team of $|U|+3$ single-state FSR
of 4 types (Variable, Carry, CarryDetect, and CarrySignal) in an environment
$E$ such that the values in the truth assignment are encoded in the positions
of the Variable FSRs. Each of these assignments is in turn checked by a
single-state FSR of a fifth type (Evaluate) against the collection of clauses
$C$ in the given instance of 3SAT, and if the truth assignment satisfies the
conjunction of the clauses in $C$, the Evaluate FSR moves one square to the
east and modifies the type of the square it is now placed on.

We modify this reduction to create a reduction from 3SAT to TeamEnvVer when
$h=1$ and $\\#_{ec}=0$ as follows:

  1. 1.

As each of the five types of FSR can only occupy very specific non-overlapping
areas in $E$, change the types of the squares in the areas occupied by each
type of FSR to new FSR-type-specific square-types $e_{Var}$, $e_{Car}$,
$e_{CD}$, $e_{CS}$, and $e_{Evl}$.

  2. 2.

Create a single-state universal FSR that can simulate all five types of FSR by
combining modified versions of all transitions in the five types of FSR, where
a transition of the form $\langle,q,f,m,q\rangle$ in FSR-type $t$ is changed
to $\langle q,enval(e_{t},(0,0))\wedge f,m,q\rangle$, i.e., FSR-type-specific
transitions can only trigger if the FSR is in the environmental area
associated with FSR of that type.

  3. 3.

Replace the environmental modification made by the Evaluate FSR with $*$.

Note that this reduction creates an FSR team in which $h=1$, $|Q|=1$, and
$\\#_{ec}=0$. If we then make the target configuration of the task associated
with this instance of TeamEnvVer be that an FSR is positioned immediately to
the west of the initial position of the Evaluate FSR, the proof of reduction
correctness in Lemma 4 also establishes that this reduction is correct. As
3SAT is $NP$-hard [18, Problem L02], this reduction establishes that
TeamEnvVer is also $NP$-hard when $h=1$, $|Q|=1$, and $\\#_{ec}=0$ and and
hence not polynomial-time solvable under these restrictions unless $P=NP$.

We now need to establish the $NP$-hardness of TeamEnvVer for fixed values of
$h>1$, $|Q|>1$, and $\\#_{ec}>0$. Observe that the instance of TeamEnvVer
constructed above makes no environmental modifications and hence trivially
makes at most $\\#_{ec}$ environmental modifications for any fixed value of
$\\#_{ec}$. In the case of $h$ and $|Q|$, construct a modified instance of
TeamEnvVer above in which there is an additional “holding area” in $E$
consisting of $h-1$ squares enclosed by obstacle-squares such that no FSR
inside this area can leave it. Populate this holding area with arbitrary FSR
based on $|Q|$ states such that none of these FSR makes an environmental
modification. As $h$ and $|Q|$ are of fixed value, this modified instance of
TeamEnvVer can still be constructed in time polynomial in the size of the
given instance of 3SAT. Moreover, the reduction from 3SAT to this modified
instance of TeamEnvVer shows the $NP$-hardness of TeamEnvVer for the specified
values of $h$ and $|Q|$.

Part (1) of Result B requires relatively straightforward modification to use
transition-template library $L$, which actually ends up simplifying the
original proof in [11]. As the team consisted of a single FSR, this trivially
gives $h=1$. Part (2) requires complex and decidedly non-trivial modification
to allow $|Q|>1$ and relatively straightforward modification of the
construction used in the proof of part (2) of Result A to allow $h>1$.

Result B (Modified from Lemma 5 in [11]): ContDesLSres is not exact
polynomial-time tractable unless $P=NP$when $|Q|=1$, $h=1$, and $\\#_{ec}=0$.
Moreover, this intractability holds for any version of ContDesLSres when
$h>1$, $|Q|>1$, and $\\#_{ec}>0$ for any fixed values of $h$, $|Q|$, and
$\\#_{ec}$..

Proof: Lemma 5 in [11] gives a reduction (based on a reduction in [12]) from
Dominating set to a problem ContDes${}^{fast}_{D,syn}$ which is essentially
ContDesLSres,D in which selection from a library $L$ of transition-templates
is simulated by specifying bounds in the problem input on $|f|$, the maximum
length of any transition trigger-formula. This reduction creates a somewhat
complex environment for a team composed of a single-state FSR [12, Figure
2(b)]. In order to force the transitions in such a robot to encode a candidate
dominating set of size $k$ in the graph $G$ in the given instance of
Dominating set, the robot has to navigate from the southwestmost corner of the
environment to the top of the $(k+1)$st column in subgrid SG1 [12, Figure
2(c)]. From there, the robot navigates the $|V|$ columns of subgrid SG2 [12,
Figure 2(d)], where each column represents the vertex neighbourhood of a
particular vertex in $G$ and the robot progresses eastward from one column to
the next if and only if that robot has a transition corresponding to a vertex
in the neighbourhood encoded in the first column. Subgrid SG2 thus checks if
the robot encodes an actual dominating set of size $k$ in $G$, such that the
robot enters the northeastmost square of the environment and builds the
requested structure there if and only if the $k$ east-moving transitions in
the robot encode a dominating set of size $k$ in $G$.

Given the above, consider the following reduction from Dominating set to
ContDesLSres,D, Given an instance $\langle G=(V,E),k\rangle$ of Dominating
Set, construct an instance $\langle
E^{\prime},E^{\prime}_{T},|T|,p_{I},|L|,r,|Q|,d,h,\\#_{ec}\rangle$ of
ContDesLSres,D as follows: Let $E^{\prime}$ be the environment constructed in
Lemma 5 of [11] with the northwest $e_{N}$-based and SG1 subgrids removed,
$E^{\prime}_{T}$ be the version of $E_{T}$ in that same lemma,
$p^{\prime}_{I}=E^{\prime}_{1,1}$, $L=\\{\langle
q,enval(y,(0,0)),*,goEast,q^{\prime}\rangle~{}|~{}y\in\\{e_{1},\ldots,e_{|V|}\\}\\}\cup\\{\langle
q,*,*,goNorth,q^{\prime}\rangle\\}$, $|T|=|Q|=h=1$, $r=0$, $d=k+1$, and
$\\#_{ec}=0$. Finally, let the target configuration of the task associated
with this instance be an FSR positioned in the northeastmost square in
$E^{\prime}$; let us call this position $p_{F}$. This instance of
ContDesLSres,D can be constructed in time polynomial in the size of the given
instance of Dominating Set.

Observe that the use of $L$ means that we no longer need subgrid SG1 and the
restrictions on $|f|$ posited in Lemma 5 mentioned above to force the created
FSR to have $k$ east-moving transitions corresponding to a candidate
dominating set of $k$ distinct vertices in $G$. Hence, by slight
simplifications and modifications of the proof of correctness of the reduction
in Lemma 5 mentioned above, it can be shown that there is a dominating set of
size $k$ in graph $G$ in the given instance of Dominating set if and only if
there is an FSR with the structure specified in the constructed instance of
ContDesLSres,D such (1) the lone FSR in $T$ can progress to $p_{F}$ if this
FSR starts at $p_{I}$ and (2) the $k+1$ transitions in this FSR are $k$ east-
moving transitions from $L$ whose activation-formula predicates correspond to
the vertices in a dominating set of size $k$ in $G$, and the final transition
in $L$. As each transition in this FSR has an activation-formula consisting of
either $*$ or a single predicate evaluating if that square has a particular
square-type, there can be at most one transition enabled at a time and the
operation of this FSR in $E^{\prime}$ is deterministic. As the single robot in
$T$ can only move north or east and does one of either in each move, the
number of transitions executed in this construction task is the Manhattan
distance from $p_{I}$ to $p_{F}$ in $E$. This distance is
$(|V|+1)+(|V|^{2}+1)<|E^{\prime}|=c_{1}|E^{\prime}|^{c_{2}}<c_{1}(|E^{\prime}|+|Q|)^{c_{2}}$
when $c_{1}=c_{2}=1$, which means that this navigation task is
$(1,1)$-completable. As Dominating set is $NP$-complete [18, Problem GT2], the
reduction above establishes that ContDesLSres,D is $NP$-hard; our main result
then follows from Lemma 1.

We now need to establish the $NP$-hardness of ContDesLSres,D for fixed values
of $h>1$, $|Q|>1$, and $\\#_{ec}>0$. The case of $\\#_{ec}$ can be handled by
mechanisms and logic analogous to those used for ContEnvVer in the proof of
Result A. The modifications required in the case of $h$ and $|Q|$ are more
complex. Let us first consider the modifications required when $h=1$ and
$|Q|>1$:

  1. 1.

Add square-types $e_{q1},e_{q2},\ldots,e_{q(|Q|-1)},e_{E}$ to $E_{T}$.

  2. 2.

Add transition-templates $\\{\langle
q,enval(e_{qi},(0,0)),goSouth,q^{\prime}\rangle,\langle
q,enval(e_{qi},(0,0)),$ $goEast,q^{\prime}\rangle~{}|~{}1\leq
i\leq|Q|-1\\}\cup\\{\langle q,enval(e_{E},(0,0)),goEast,q^{\prime}\rangle\\}$
to $L$.

  3. 3.

Create a subgrid area $SGA$ of $|Q|$ columns, where column $j$ includes a
subcolumn consisting of the square-types $e_{q1},e_{q2},\ldots,e_{qj}$ if
$1\leq j\leq\min(k+2,|Q|-1)$ and a subcolumn consisting of the square-types
$e_{q1},e_{q2},\ldots,e_{q(k+2)},q_{ej}$ if $\min(k+2,|Q|-1)+1\leq
j\leq|Q|-1$. Let the bottom square of the subcolumn in column $j$ be
immediately to the west of the top square of the subcolumn in column $j+1$ for
$1\leq j\leq|Q|-2$, and there be a single square of type $e_{E}$ immediately
to the east of the bottom square in the subcolumn in column $|Q|-1$,

  4. 4.

Create a new environment $E^{\prime\prime}$ by positioning the subgrid $SGA$
specified in (3) in $E^{\prime}$ to the immediate east of SG2 such that the
square of type $e_{E}$ is to the immediate east of the southwestmost square of
SG2. Fill all previously unspecified squares in $E^{\prime\prime}$ with the
square-type $e_{N}$.

  5. 5.

Reset $k^{\prime}$ to $k+2$.

  6. 6.

Reset $p_{I}$ to the northwestmost square of $SGA$ in $E^{\prime\prime}$.

In order to progress eastward from column $j$ to $j+1$, $1\leq j\leq|Q|-1$,
the current state $q$ of the FSR must have the transition $\langle
q,enval(e_{qj},(0,0)),goEast,q^{\prime}\rangle$. However, to progress
southwards from the beginning of the subcolumn to the bottom of the subcolumn
in column $j$, the current state $q$ of the FSR must also have the transitions
in the set $\\{\langle q,enval(e_{qi},(0,0)),goSouth,q\rangle~{}|~{}1\leq
i<\min(j-1,k+1)\\}$. Recall that the rules of deterministic FSR operation
require that there cannot be more than one transition enabled at a time.
Hence, an FSR that can successfully navigate the first $|Q|-1$ columns of
subgrid $SGA$ must have a different state for each column of $SGA$. To then
enter and traverse SG2, the state $q_{|Q|-1}$ of the FSR must have the
transitions $\langle q_{|Q|-1},enval(e_{E},(0,0)),goEast,q_{|Q|-1}\rangle$,
$\langle q_{|Q|-1},*,goNorth,q_{|Q|-1}\rangle$, and $k^{\prime\prime}\leq
k^{\prime}$ eastward-moving transitions that correspond to a dominating set of
size $k^{\prime\prime}\leq k$ in graph $G$ in the given in stance of
Dominating set. In the case when $h>1$, further modify $E^{\prime\prime}$ to
include an additional “holding area” consisting of $h-1$ squares enclosed by
obstacle-squares such that no FSR inside this area can leave it, and populate
this holding area with arbitrary FSR based on $|Q|$ states such that none of
these FSR makes an environmental modification. As $h$ and $|Q|$ are of fixed
value, the modified instance of ContDesLSres,D can still be constructed in
time polynomial in the size of the given instance of Dominating set. Moreover,
the reduction from Dominating set to this modified instance of ContDes:Sres,D
shows the $NP$-hardness of ContDesLSres,D for the specified values of $h$ and
$|Q|$.

Result C requires a very straightforward modification of the algorithm
presented previously in [9] to incorporate our new task model.

Result C (Modified from Result A in [9]): TeamDesLSres is exact polynomial-
time tractable when $|Q|\geq 1$, $h=1$, and $\\#_{ec}\geq 0$.

Proof: The algorithm in the proof of Result A in [9], which tests for each
controller $c$ in $L$ whether a team based entirely on $c$ can construct $X$
at $p_{X}$ in at most $c_{1}(|E|+|Q|)^{c_{2}}$ timesteps, operates in
polynomial time. We need only modify that algorithm such that after each
timestep of robot team operation we check if the target configuration
associated with the task is in $E$, which can also be done in polynomial time.

Part (1) of Result D requires relatively straightforward modification to
environment and FSR transitions to force $h$ to be exactly 2 (as the proof
presented previously in [10] only required that $h\leq 2$). Part (2) requires
relatively straightforward modification of the construction used in the proof
of part (2) of Result A to allow $h>2$ and $|Q|>1$.

Result D (Modified from Lemma A.7 in [10]): TeamDesLSres is not exact
polynomial-time tractable when $|Q|=1$, $h=2$, and $\\#_{ec}=0$. Moreover,
this intractability holds for any version of TeamDesLSres when $h>2$, $|Q|\geq
1$, and $\\#_{ec}\geq 0$.

Proof: Lemma A.7 in [10] gives a reduction from 3SAT to a version of problem
SelAlgD which is for all intents and purposes our problem TeamDesLSres,D
without the restriction on $\\#_{ec}$. In the instance of SelAlgD created by
this reduction, a team $T$ of $|U|$ single-state FSR chosen from an FSR
library of size 2 encodes a truth assignment to the variables in $U$ in the
given instance of 3SAT and each of the southmost $|C|$ rows in environment $E$
encodes a clause in $C$ in this same instance. All robots in the team can
reach the northmost row of $E$ and deposit a strip of $|T|$ squares of type
$e_{X}$ if and only the truth-assignment to the variables encoded in $T$
satisfies all clauses in $C$; moreover, as at least one robot moves north in
each timestep and strictly less than $|E|$ such moves can be made if $T$
encodes such a valid truth-assignment, the task is $(1,1)$-completable.

We modify this reduction to create a reduction from 3SAT to TeamDesLSres,D
when $h=2$ and $\\#_{ec}=0$ as follows:

  1. 1.

Add new clauses $c_{T}$ and $c_{F}$ to $C$ and new variables $u_{T}$ and
$u_{F}$ to $U$ such that clause $c_{T}$ ($c_{F}$) is satisfied if and only if
variable $u_{T}$ ($u_{F}$) is assigned value $True$ ($False$).

  2. 2.

Add two extra columns to $E$ corresponding to new clauses $c_{T}$ and $c_{F}$
and add two to the size of the robot team.

  3. 3.

In each of two robots in $L$, replace the transition that modifies the type
$e_{B}$ of a square in the northmost row of $E$ to $e_{X}$ with a transition
that stays on a square of type $e_{B}$.

Note that this reduction creates an FSR team in which $h=2$, $|Q|=1$, and
$\\#_{ec}=0$. If we then make the target configuration of the task associated
with this instance of TeamDesLSres,D be $|T|$ FSRs positioned in the central
squares of the northmost row of $E$, the proof of reduction correctness in
Lemma A.7 also establishes that this reduction is correct. As 3SAT is
$NP$-hard [18, Problem L02], this reduction establishes that TeamDesLSres,D is
also $NP$-hard when $h=1$, $|Q|=1$, and $\\#_{ec}=0$; our main result then
follows from Lemma 1.

To complete the proof, observe that to establish the $NP$-hardness of
TeamDesLSres,D and polynomial-time unsolvability unless $P=NP$ of TeamDesLSres
for fixed values of $h>2$, $|Q|>1$, and $\\#_{ec}>0$, we can use appropriately
modified versions of the constructions and logic that we used to show similar
results for TeamEnvVer in the proof of Result A above.

Finally, in Result E, inapproximability in senses (1) and (2) follows from the
previously presented proofs in [10] and the $NP$-hardness of decision versions
of our verification and design problems shown in Results A, B, and D above.
Inapproximability in sense (3) is new to the current paper but follows in a
very straightforward manner from a known class inclusion result in
computational complexity theory and the inapproximability proof for sense (2).

Result E (Modified from Results A.4 and A.5 in [10]): None of the intractable
versions of TeamEnvVer, ContDesLSres, or TeamDesLSres described in Results A,
B, and D are polynomial-time approximable in senses (1–3) unless $P\neq BPP$
and $P=NP$.

Proof: That approximate polynomial-time solvability in sense (1) for any of
the listed problems implies $P=NP$ follows from the $NP$-hardness of the
decision versions of each of these problems (which is established in the
proofs of Result A, B, and D) and Corollary 2.2 in [21].

With respect to approximate polynomial-time solvability in sense (2), it is
widely believed that $P=BPP$ [20, Section 5.2] where $BPP$ is considered the
most inclusive class of decision problems that can be efficiently solved using
probabilistic methods (in particular, methods whose probability of correctness
is $\geq 2/3$ and can thus be efficiently boosted to be arbitrarily close to
one). Hence, if any of the listed problems has a probabilistic polynomial-time
algorithm which operates correctly with probability $\geq 2/3$ then the
decision version of that problem is by definition in $BPP$. However, if
$BPP=P$ and we know that each of the decision versions of the listed problems
is $NP$-hard by the proofs of Results A, B, and D, this would then imply by
the definition of $NP$-hardness that $P=NP$.

With respect to approximate polynomial-time solvability in sense (3), it is
known that $ZPP\subseteq BPP$, where $ZPP$ is the class of decision problems
that can be always be solved correctly by algorithms with expected polynomial
runtime [23]. Hence, if any of the listed problems is approximately solvable
in sense (3) then the decision version of that problem is by definition in
$ZPP$ as well as $BPP$. However, as each of the decision versions of the
listed problems is $NP$-hard by the proofs of Results A, B, and D, this would
then imply by the definition of $NP$-hardness and the widely-believed
conjecture $P=BPP$ that $P=NP$.