# Strongly interacting photons in 2D waveguide QED

M. Tečer  Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8,
I-35131 Padova, Italy  M. Di Liberto  Dipartimento di Fisica e Astronomia “G.
Galilei”, via Marzolo 8, I-35131 Padova, Italy  Padua Quantum Technologies
Research Center, Universitá degli Studi di Padova  Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy  P. Silvi
Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131
Padova, Italy  Padua Quantum Technologies Research Center, Universitá degli
Studi di Padova  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di
Padova, I-35131 Padova, Italy  S. Montangero  Dipartimento di Fisica e
Astronomia “G. Galilei”, via Marzolo 8, I-35131 Padova, Italy  Padua Quantum
Technologies Research Center, Universitá degli Studi di Padova  Istituto
Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy
F. Romanato  Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8,
I-35131 Padova, Italy  Padua Quantum Technologies Research Center, Universitá
degli Studi di Padova  CNR-IOM Istituto Officina dei Materiali, Trieste, Italy
G. Calajó  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova,
I-35131 Padova, Italy  Dipartimento di Fisica e Astronomia “G. Galilei”, via
Marzolo 8, I-35131 Padova, Italy  Dipartimento di Fisica e Astronomia “G.
Galilei”, via Marzolo 8, I-35131 Padova, Italy  Padua Quantum Technologies
Research Center, Universitá degli Studi di Padova  Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy  Dipartimento
di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131 Padova, Italy
Padua Quantum Technologies Research Center, Universitá degli Studi di Padova
Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131
Padova, Italy  Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo
8, I-35131 Padova, Italy  Padua Quantum Technologies Research Center,
Universitá degli Studi di Padova  Istituto Nazionale di Fisica Nucleare
(INFN), Sezione di Padova, I-35131 Padova, Italy  Dipartimento di Fisica e
Astronomia “G. Galilei”, via Marzolo 8, I-35131 Padova, Italy  Padua Quantum
Technologies Research Center, Universitá degli Studi di Padova  CNR-IOM
Istituto Officina dei Materiali, Trieste, Italy  Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy

(August 28, 2024)

###### Abstract

One dimensional confinement in waveguide Quantum Electrodynamics (QED) plays a
crucial role to enhance light-matter interactions and to induce a strong
quantum nonlinear optical response. In two or higher dimensional settings,
this response is reduced since photons can be emitted within a larger phase
space, opening the question whether strong photon-photon interaction can be
still achieved. In this study, we positively answer this question for the case
of a 2D square array of atoms coupled to the light confined into a two-
dimensional waveguide. More specifically, we demonstrate the occurrence of
long-lived two-photon repulsive and bound states with genuine 2D features.
Furthermore, we observe signatures of these effects also in free-space atomic
arrays in the form of weakly-subradiant in-band scattering resonances. Our
findings provide a paradigmatic signature of the presence of strong photon-
photon interactions in 2D waveguide QED.

The capability of engineering effective photon-photon interactions is a
compelling requirement in many quantum technology applications for the
generation and manipulation of complex states of light Chang  _et al._ (2014).
This task has been successfully achieved in light-matter interfaces, such as
cavity QED Reiserer and Rempe (2015), Rydberg atomic ensembles Firstenberg
_et al._ (2016, 2013); Liang  _et al._ (2018) and waveguide QED Sheremet  _et
al._ (2023); Roy  _et al._ (2017), where non-linear quantum optical effects
are strongly enhanced. In particular, waveguide QED has recently emerged as a
versatile and experimentally implementable platform where the light confined
in a 1D channel, either at optical Lodahl  _et al._ (2015); Hood  _et al._
(2016); Corzo  _et al._ (2019); Tiranov  _et al._ (2023) or microwave
frequencies Astafiev  _et al._ (2010); Brehm  _et al._ (2021); Mirhosseini
_et al._ (2019); Kannan  _et al._ (2023), couples to one or multiple quantum
emitters. In this scenario, light confinement leads to strong photon
correlations, resulting in anti-bunched Shen and Fan (2007a); Schrinski and
Sørensen (2022a) or bunched output photons Shen and Fan (2007a, b); Zheng  _et
al._ (2011); Mahmoodian  _et al._ (2018); Prasad  _et al._ (2020); Mahmoodian
_et al._ (2020); Le Jeannic  _et al._ (2022). Among the causes of the bunching
phenomena is the emergence of multi-photon bound states Shen and Fan (2007a);
Zheng  _et al._ (2011); Mahmoodian  _et al._ (2020); Tomm  _et al._ (2023)
where two or more photonic excitations propagate jointly in the system with
spatially correlated positions. Remarkably, the lifetime and the dispersion
properties of these states can be enhanced when ordered atomic arrays are
coupled to the waveguide Zhang and Mølmer (2020); Poddubny (2020); Bakkensen
_et al._ (2021); Calajó and Chang (2022) or to waveguide networks Marques  _et
al._ (2021); Schrinski  _et al._ (2023). In this setting, the realized bound
state can be interpreted as a propagating excitation, which acts as a moving
defect in the otherwise periodic medium, dragging other photons with it. The
recent progress in building scalable microwave resonator arrays coupled to
superconducting qubits Zhang  _et al._ (2022); Scigliuzzo  _et al._ (2022);
Gong  _et al._ (2021), in interfacing 2D photonic crystals with cold atoms Yu
_et al._ (2019) or quantum dots Lodahl  _et al._ (2015); Hauff  _et al._
(2022) and quantum emitters interacting with atomic matter waves Krinner  _et
al._ (2018); Kim  _et al._ (2023) have created new experimentally viable
opportunities for studying QED in 2D confined geometries. This scenario is
particularly fascinating when combined with the possibility to induce photon-
photon interactions and to build strongly correlated many-body quantum phases
of light Noh and Angelakis (2016).

Figure 1: (a) A square array of two level atoms is coupled to the light
confined into a 2D photonic waveguide. (b) Single photon dispersion relation
for $k_{0}d=\frac{\pi}{2}$. Upper right panel: first Brillouin zone and
$\Gamma,X,M$ symmetry points. Lower right panel: Bandgap energy as a function
of the interatomic distance $k_{0}d$. The gap closes at $k_{0}d>\pi$. The
dashed orange lines indicate the divergence of the dispersion relation at
$|\mathbf{k}|=k_{0}$.

In this letter, we show that photonic states exhibiting both spatial repulsion
and binding can occur by coupling a 2D atomic array to light confined into a
2D photonic structure (see Fig. 1). The emergence of such strong interactions
in this 2D environment represents a nontrivial outcome. Indeed, unlike the 1D
scenario, where photon-mediated interactions among emitters are infinite-
range, in 2D, they decay as the square root of the emitters separation. We
show that the interacting photonic states occurring in the system possess a
distinct 2D structure, exhibiting binding and repulsion along different
directions, and that their existence is made possible by an interplay between
the coherent and dissipative parts of the emitters collective interaction.
Furthermore, we identify meta-stable in-band localized scattering resonances.
Notably, these states also manifest in free-space atomic arrays, expanding the
range of platforms where our findings can be explored Rui  _et al._ (2020);
Srakaew  _et al._ (2023). Finally, we demonstrate how to excite the photon
bound states via a dynamical relaxation process. Our results provide a
realistic pathway for realizing strongly-correlated states of light in 2D
confined geometries in the sub-wavelength regime.

Model. We consider a two-dimensional square array of $N$ two-level atoms with
ground and excited states $|g\rangle$, $|e\rangle$ and lattice constant $d$
perfectly coupled (see Ref. Sup ) to the light confined into a 2D waveguide,
as shown in Fig. 1. To simplify the description of this complex light-matter
interacting system we integrate out the photonic degrees of freedom employing
a Born-Markov approximation. This procedure is valid as long as the atom-
photon dynamics, which is set by the atomic decay rate $\gamma$, occurs on a
time-scale slower than a photon freely propagating through the whole array
(see Ref. Sup ). We obtain an effective spin model for the atoms described by
the non-hermitian Hamiltonian ($\hbar=1$) Asenjo-Garcia  _et al._ (2017a, b);
Sheremet  _et al._ (2023)

$\hat{H}_{\rm
eff}=\sum_{ij}G\left(k_{0}|\mathbf{x}_{i}-\mathbf{x}_{j}|\right)\hat{\sigma}_{+}^{i}\hat{\sigma}_{-}^{j},$
(1)

where $k_{0}=2\pi/\lambda$ is the photon wavevector whose corresponding
frequency is resonant with the atomic transition frequency $\omega_{a}$ and
$\hat{\sigma}_{+}^{i}=|e\rangle\langle g|$ and
$\hat{\sigma}_{-}^{i}=(\hat{\sigma}_{+}^{i})^{\dagger}$ are the pseudospin
operators for the $i$-th atom located at the position $\mathbf{x}_{i}$. The
dissipative and coherent long-range photon-mediated interactions among the
atoms are encoded in the function
$G\left(k_{0}|\mathbf{x}_{i}-\mathbf{x}_{j}|\right)$, related to 2D
electromagnetic Green’s function of the nanophotonic structure Asenjo-Garcia
_et al._ (2017b). For confined photons with a quadratic isotropic dispersion
relation, it reads $G(z)=(\gamma/2)(\mathcal{Y}_{0}(z)-i\mathcal{J}_{0}(z))$
Sup ; González-Tudela  _et al._ (2015); González-Tudela and Cirac (2017);
Galve and Zambrini (2018); González-Tudela and Galve (2018), where
$\mathcal{J}_{0}(z)$ and $\mathcal{Y}_{0}(z)$ are respectively the zeroth
order Bessel functions of the first and second kind, which decay as
$z^{-1/2}$. The Hamiltonian in Eq. (1) provides a complete description of the
system within a fixed excitation sector in the absence of an external pumping
field. For a single excitation, the Hamiltonian is diagonalized by Bloch
waves, labelled by the wavevector $\mathbf{k}$, having eigenvalues
$E(\mathbf{k})=\epsilon^{(1)}(\mathbf{k})-i\gamma^{(1)}(\mathbf{k})/2$, where
$\epsilon^{(1)}(\mathbf{k})$ and $\gamma^{(1)}(\mathbf{k})$ provide the
single-excitation dispersion relation and collective decay rates,
respectively. The dispersion relation (Fig. 1) exhibits two distinct
polaritonic branches, along with a band gap near the atomic resonance
frequency, where propagation of excitations is forbidden. Differently from 1D
wQED Calajó and Chang (2022); Sheremet  _et al._ (2023); Poddubny (2020), the
bandgap closes for inter-atomic distances above the sub-wavelength regime (see
Fig. 1), $k_{0}d>\pi$, namely for $d>\lambda/2$ due to the occurrence of
higher diffraction orders. Note that, similar to the conventional 1D waveguide
scenario Albrecht  _et al._ (2019); Calajó and Chang (2022), the divergences
in the dispersion relation at the resonant wavevector, $|\mathbf{k}|=k_{0}$,
are associated with super-radiant modes Albrecht  _et al._ (2019); Zhang and
Mølmer (2019); Kumlin  _et al._ (2020).

Figure 2:  (a) Two-excitation spectrum as a function of the center of mass
momentum for $k_{0}d=\frac{\pi}{2}$ . Blue dots represent unbound states
energy continuum while the red line and the orange dot at the $\Gamma$ point
represent the bound states. The dashed black boxes indicate the regions where
the repulsive states lie. (b) Relative coordinate population distribution,
$p(r_{x},r_{y})$, of two kind of repulsive states at $\Gamma$ point for
$k_{0}d=0.5$. Similar distributions up to a phase modulation are observed at
the $M$ point. (c) Zoom of the bound states distribution at $\Gamma$ and M
point for $k_{0}d=0.5\pi,0.7\pi$. (d) Two-excitation spectrum at the $M$ point
(upper panel) and average relative distance between two excitations, $\langle
r\rangle$, (lower panel) as a function of the inter-atomic distance. In the
upper panel, the continuous line indicates the BS energy while the dashed one
indicates the scattering resonances branches discussed in the main text.
Panels (b),(c) ((a),(d)) were obtained considering an array made of
$N=30\times 30$ ($100\times 100$) atoms.

Interacting photons. To demonstrate that strong photon-photon interactions
occur in this 2D setting, we consider the two-body problem described by Eq.
(1) restricted to the two excitations subspace. By re-parameterizing the
arrangement of the atoms in terms of center-of-mass ($\mathbf{R}$) and
relative ($\mathbf{r}$) coordinates, we reduce the two-body problem to a
single particle one for the relative coordinate space and described by the
basis set $|\mathbf{K},\mathbf{r}\rangle$, where the parametric dependence on
the center of mass momentum $\mathbf{K}$ is explicitly indicated (see Ref. Sup
). The corresponding Hamiltonian matrix elements read:

$\langle\mathbf{K},\mathbf{r}|\hat{H}_{\textrm{eff}}|\mathbf{K},\mathbf{r}^{\prime}\rangle=\sum_{\epsilon=\pm}\cos\mathbf{K}\cdot\left(\mathbf{r}+\epsilon\mathbf{r}^{\prime}\right)G\left(k_{0}\left|\mathbf{r}+\epsilon\mathbf{r}^{\prime}\right|\right)\,.$
(2)

Numerical diagonalization returns the real part of the two-body spectrum,
$\epsilon^{(2)}(\mathbf{K})$, plotted in Fig. 2(a). We observe a continuum of
unbound two-particle states, whose energies correspond to the sum of two
single-particle excitations, separated by a bandgap hosting a dispersive
state. Most of the continuum states away from the high-symmetry points of the
Brillouin zone (BZ) are scattering states made of a pair of particles with
relatively high group velocity (see Fig. 1), thus experiencing a weak mutual
interaction Sheremet  _et al._ (2023). Instead, close to those high-symmetry
points, the relatively flat single-particle dispersion can provide an
enhancement of interactions. A convenient approach to study photon-photon
interactions can be formulated by mapping the matrix elements in Eq. (2) onto
a single-particle model for the relative coordinate $\mathbf{r}$ in the
presence of an impurity potential, described by the Hamiltonian:

$\hat{H}_{\mathbf{K}}=\sum_{\mathbf{r}\mathbf{r}^{\prime}}J^{\mathbf{K}}_{\mathbf{r}\mathbf{r}^{\prime}}|\mathbf{K},\mathbf{r}\rangle\langle\mathbf{K},\mathbf{r}^{\prime}|+U|\mathbf{K},0\rangle\langle\mathbf{K},0|\,,$
(3)

where the first term describes free particle propagation, while the second
term implements the short-range impurity potential with infinite strength,
$U\rightarrow\infty$, located at $\mathbf{r}=0$. Here the hopping coefficient
is
$J^{\mathbf{K}}_{\mathbf{r}\mathbf{r}^{\prime}}=(1/N)\sum_{\mathbf{q}}\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})e^{i(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{q}}$,
where $\mathbf{q}$ is the relative momentum and the dispersion is given by the
sum of individual photon energies,
$\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})=\epsilon^{(1)}(\mathbf{q})+\epsilon^{(1)}(\mathbf{K}-\mathbf{q})$.
The effective model of Eq. (3) can be formally obtained by reformulating the
Hamiltonian (1) in terms of hardcore bosons and then solving again the two-
body problem in the relative coordinate frame Sup . In this way, the hard-core
photon-photon interactions are mapped to the impurity term of Eq. (3), which
reflects the saturation of an atom after the absorption of one excitation.

The first paradigmatic consequence of strong photon interactions originating
from the hard-wall condition induced by the impurity potential is the presence
of two kinds of _repulsive_ scattering states indicated by the black dashed
squares in Fig. 2(a). Their population distribution in the relative coordinate
shows a smooth repulsion of the two excitations either with respect to the
$r_{x}=0$ and $r_{y}=0$ axis (type I) or with respect to the $r_{x}=r_{y}$
diagonal (type II), as shown in Fig. 2(b). The 1D analogue of these states
have been recently discussed both in free space Asenjo-Garcia  _et al._
(2017b) and wQED atomic arrays Zhang and Mølmer (2019); Albrecht  _et al._
(2019); Ostermann  _et al._ (2019); Needham  _et al._ (2019); Zhong  _et al._
(2020); Kornovan  _et al._ (2019); Schrinski and Sørensen (2022b) and are
commonly referred to as fermionic states due to the fact that their population
distribution resembles the one of noninteracting fermions.

Besides repulsive states we also identify the presence of photon-photon bound
states (BSs) Shen and Fan (2007a); Zheng  _et al._ (2011); Mahmoodian  _et
al._ (2020); Tomm  _et al._ (2023); Zhang and Mølmer (2020); Poddubny (2020);
Bakkensen  _et al._ (2021); Calajó and Chang (2022). Here the two excitations
are spatially correlated and their occurrence can be understood, within the
impurity model (3), in terms of a defect bound state localized around
$\mathbf{r}=\mathbf{0}$. The BSs energy, $E_{\rm BS}(\mathbf{K})$, indicated
by the red lines of Fig. 2(a), can be computed by using standard Green’s
function based methods (see Ref. Sup ) and it is given by the numerical
solution of $\mathcal{G}_{0}(0,E_{\rm BS}(\mathbf{K}))=1/U$ for
$U\rightarrow\infty$, with

$\mathcal{G}_{0}(\mathbf{r},E)=\frac{1}{(2\pi)^{2}}\int_{\rm
BZ}d^{2}\mathbb{q}\frac{e^{i\mathbb{q}\cdot\mathbb{r}}}{E-\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})}\,,$
(4)

being the free particle Green’s function. The BSs solutions need to respect
the constraint $E_{\rm
BS}(\mathbf{K})\neq\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})$, by
definition, meaning that their occurrence is conditioned to the existence of a
band gap in the relative coordinate space. This requirement is not matched for
$k_{0}d>\pi$ due to the closing of the single-excitation gap in 2D discussed
before (see Fig. 1), which prevents the gap formation also for the two-
particle spectrum. For $k_{0}d<\pi$, a gap in
$\epsilon^{(2)}_{\text{scat}}(\mathbf{K})$ exists for values of the center of
mass momentum in the range $|\mathbf{K}|>2k_{0}$ within the BZ. When $2k_{0}$
becomes equal to the $M$ point momentum, $|\mathbf{K}_{M}|=\sqrt{2}\pi/d$, the
condition $|\mathbf{K}|>2k_{0}$ cannot be satisfied within the BZ and the gap
is therefore closed (see Ref. Sup  and Fig. 2(d)). An additional gap occurs
exactly at the $\Gamma$ point due to a perfect vanishing of the density of
states Bakkensen  _et al._ (2021); Calajó and Chang (2022) and persists for
all inter-atomic distances satisfying $k_{0}d<\pi$.

Figure 3: (a) Decay rates $\gamma^{(2)}_{s}$ for the states discussed in the
main text as a function of inter-atomic distance $k_{0}d$. Here we fixed the
system size to a square lattice of $N=10\times 10$ atoms. (b) Finite size
scaling of the decay rates $\gamma_{s}$ with the system size $N$ in the 2D
waveguide case (filled symbols) compared to the free space scenario (empty
symbols) indicated by the grey scale curves. Here we fixed $k_{0}d=0.52\pi$
for the 2D waveguide while $k_{0}d=1.09\pi$ and $k_{0}d=0.73\pi$ for the
$\Gamma$ and $M$ point in the free space case. In both plots, we considered
the repulsive (II) state.

Once established the conditions of existence of the BSs, we can use the
impurity model to compute the bound states population distribution,
$p(\mathbf{r})=|\mathcal{G}_{0}(\mathbf{r},E_{\rm BS}(\mathbf{K}))|^{2}$, in
the relative coordinate space. An example is shown in Fig. 2(c) for the
$\Gamma$ and $M$ points at two different inter-atomic distances. These plots
clearly show a pronounced localization in the vicinity of $\mathbf{r}=0$,
meaning that the two excitations are bound and spatially correlated, as
expected for bound states. This binding, which we quantify by the average
distance between two excitations, $\langle
r\rangle=\sum_{\mathbf{r}}p(\mathbf{r})|\mathbf{r}|$, strongly depends on the
inter-atomic distance, $k_{0}d$, as shown in Fig. 2(d). The plot shows how the
bound states at $\Gamma$ and $M$ points are largely extended at short inter-
atomic distances and progressively become more localized. At the inter-atomic
distance of $k_{0}d\approx 0.52\pi$ the $M$ point bound state becomes most
localised while the $\Gamma$ point BS changes its relative coordinate
distribution from being peaked at next-nearest neighbour sites to being peaked
at the nearest neighbour ones (see first row of Fig. 2(c)). We observe that
this transition occurs when the nearest-neighbour coherent and dissipative
interactions become equal, suggesting that the 2D spatial structure of the BSs
depends on a nontrivial interplay between coherent and dissipative atom-atom
interactions. Finally, as anticipated, at $k_{0}d=\pi/\sqrt{2}$ the $M$ point
gap closes and the BS ceases to exist. Nevertheless, we can still identify the
presence of two localized states in the band, which can be interpreted as
scattering resonances Bakkensen  _et al._ (2021); Sup , whose energies and
localization are shown in Fig. 2(d)). Interestingly, close to $k_{0}d=\pi$, we
observe the formation of a region with low density of states (see insets in
Fig. 2(d)), which plays the role of an effective gap in a finite size system,
thus stabilizing them into quasi-bound states. This mechanism depends on the
shape of the dispersion relation and occurs also in free space atomic arrays,
as we later discuss.

Interacting states lifetime. The states described in the previous section for
an infinite atomic array have vanishing decay rates in the thermodynamic
limit. In realistic settings where the atomic lattice has a finite size, these
states acquire a finite life-time due to leakage of the excitations at the
array’s edges. The decay rate of these states is provided by the imaginary
part of the two-excitation sector eigenvalues, $\gamma^{(2)}_{s}$, obtained
via exact diagonalization of Eq. (1) Albrecht  _et al._ (2019); Yu  _et al._
(2019), where $s$ labels the different states. For both repulsive and bound
states, the decay rates are much smaller than the single atom emission rate,
$\gamma$, indicating a sub-radiant behavior. We find that, while the $\Gamma$
point BSs have a weak dependence with the inter-atomic distance, the $M$ point
BSs, the associated scattering resonance for $k_{0}d>\pi/\sqrt{2}$ and the
repulsive states strongly depend on that, with the shortest lifetime reached
at around $k_{0}d\sim 0.77\pi$. The largest lifetime for both $\Gamma$ and $M$
point bound states is instead achieved around $k_{0}d\approx 0.52$, where we
previously observed a pronounced localization of the states (see Fig. 2(d)).
At this atomic distance, the decay rates of the $\Gamma$ and $M$ point BSs
scale approximately as $\sim N^{-1.5}$, as shown in Fig. 3. The repulsive
states exhibit an even stronger sub-radiant behaviour with a scaling of $\sim
N^{-3}$. These findings demonstrate that long-lived states exhibiting photon-
photon interaction can occur even for finite size arrays and can potentially
be observed in the system dynamics, as we discuss later.

Free space array. To draw a comparison with the case of an atomic array
coupled to a 2D photonic structure, we consider a 2D square array of atoms
embedded in free space Bettles  _et al._ (2016); Shahmoon  _et al._ (2017);
Manzoni  _et al._ (2018). If we fix for simplicity the atomic polarisation to
be orthogonal to the photonic structure, the projected dyadic Green’s function
ruling the atomic interactions in Eq. (1) reads Asenjo-Garcia  _et al._
(2017a, b) (see Sup ):
$G(x)=(3\gamma_{0})(4(k_{0}x)^{3})e^{ik_{0}x}\Big{(}1-ik_{0}x-(k_{0}x)^{2}\Big{)}$,
where $\gamma_{0}=d_{eg}^{2}k_{0}^{3}/(3\pi\hbar\epsilon_{0})$ is the free
space spontaneous emission decay rate with $d_{eg}$ being the dipole moment
strength. In this scenario there is no bandgap in the single-excitation
subspace and therefore neither in the two-excitations subspace. Thus, no bound
states can form in this setting. However, as previously discussed for the 2D
waveguide, it is still possible to identify two-excitation scattering
resonances exhibiting a localised wavefunction in their relative coordinate
(see Sup ). For these states, we pinpoint regions characterized by a low
density of states, where they become meta-stable within the subradiant region
$k_{0}d<\sqrt{2}\pi$ (see Sup ). This is illustrated in Fig. 3 by the scaling
of sub-radiant decay rates at the $\Gamma$ and $M$ points, exhibiting
dependencies of $N^{-1.2}$ and $N^{-0.7}$, respectively. It is noteworthy that
these states can endure in the system for a timescale more than ten times
greater than that associated with the decay of a single atom. In addition to
these localized states, we also identify the presence of repulsive states
presenting a similar decay rate scaling as the ones for a 2D waveguide. All
these interacting excitations can potentially be observed in current
experimental implementations of subwavelength ordered atomic arrays Rui  _et
al._ (2020); Srakaew  _et al._ (2023).

Figure 4: (a) Two-excitation probability distribution in the relative (left
panel) and lab coordinate (right panel) space at times $\gamma t=0,20,300$ for
$\ell=1$ and $k_{0}d=0.52\pi$. (b) Time-evolution of the two-particle
correlator, $C_{\bm{\ell}}(t)$, for different inter-atomic distances $k_{0}d$
and different initial populations ($\ell=1,2$). A large overlap with a
localized bound state ($\ell=1$, as in panel $(a)$) appears as a persistent
correlation at long times only for $k_{0}d=0.52\pi$. For these figures, we
considered a $N=10\times 10$ atomic array.

Dynamical excitation. To observe the manifestation of bound states in the
system dynamics, we consider an experimentally feasible relaxation scenario
where two atoms are initially excited in a configuration that has a large
overlap with the target state. To be concrete, we consider the dynamical
excitation of the $M$ point BS by exciting two atoms in the bulk of the system
separated by a distance $\mathbf{x}_{i}-\mathbf{x}_{j}=\bm{\ell}=(\ell,\ell)$,
$|\psi(t=0)\rangle=\hat{\sigma}_{+}^{\mathbf{i}}\hat{\sigma}_{+}^{\mathbf{i+\bm{\ell}}}|0\rangle$.
For an inter-atomic distance of $k_{0}d\approx 0.52$ and an excitation
separation of $\ell=1$, this configuration has indeed a large initial overlap
with the population distribution of the $M$ point BS (see Ref. Sup ). In Fig.
4 we plot the two excitation probability distribution in the relative and real
space coordinate at different times for this initial configuration. Note that
we evolve the state via the time-dependent Schrödinger equation ruled by the
non-Hermitian Hamiltonian (1). With this approach, we avoid tracking the
excitations leaving the two-excitation subspace while still getting the
correct dynamics within this subspace (see Sup ). Since the dominant
contribution comes from the $M$ point BSs, the two excitations remain
localized in their relative position during the evolution , while the BSs
center of mass can diffuse across the lattice. As a measurable observable to
detect the dynamical excitation of the BS we define the two particle equal
time correlator
$C_{\bm{\ell}}(t)=\sum_{\mathbf{i}}\langle\psi(t)|\hat{n}_{\mathbf{i}}\hat{n}_{\mathbf{i}+\bm{\ell}}|\psi(t)\rangle/|\langle\psi(t)|\psi(t)\rangle|^{2}$,
with the vector index $\mathbf{i}=(i_{x},i_{y})$ spans the whole array,
$\hat{n}_{\mathbf{i}}=\hat{\sigma}_{+}^{\mathbf{i}}\hat{\sigma}_{-}^{\mathbf{i}}$
and the denominator normalizes this quantity with respect to the probability
of having two excitations in the system. This correlator indicates how likely
is to find two excitations at a relative distance $\bm{\ell}=(\ell,\ell)$ in
the entire array. Its time evolution is plotted in Fig. 4 for different
initial conditions. In particular, we consider the case $\ell=1$, the one
having a large overlap with the BS, and $\ell=2$. In the first case, the two
excitations stay bound together as signaled by $C_{\bm{\ell}}(t)$ relaxing to
a finite value, whereas in the second case, the correlator quickly drops to
zero. A large dropping is also observed, as expected, for an inter-atomic
distance of $k_{0}d\approx 1.2$, where BSs do not exist.

Conclusions and Outlook. We have discussed the emergence of two-photon
strongly-correlated and long-lived states, in particular repulsive states and
bound states, in a 2D array of atoms coupled to a two-dimensional waveguide.
Using a spin model formulation, we have characterized these states discussing
their dispersion, binding and decay properties and we have demonstrated their
manifestation in the system dynamics. The occurrence of these states crucially
relies on interfacing an array of two-level atoms with the light confined into
a two-dimensional photonic structure, but some of the discussed features could
persist also in free space atomic arrays. This work opens interesting
possibilities to study many-body physics with photons in this open and
dissipative system. Targets of ongoing investigation include exploring few-
body interacting topological states Ozawa  _et al._ (2019); Salerno  _et al._
(2018, 2020); Clark  _et al._ (2020); Léonard  _et al._ (2023) in atomic
arrays Perczel  _et al._ (2017a, b, 2020) and emergent many-body phases in
frustrated lattices Bergholtz and Liu (2013); Semeghini  _et al._ (2021).
Another interesting perspective is to investigate the occurrence of high
dimensional bound states in 3D atomic arrays Brechtelsbauer and Malz (2021);
Andreoli  _et al._ (2023) where the complete filling of the electromagnetic
environment could lead to strong non-linear optical effects.

This work was supported by the EU QuantERA2017 project QuantHEP, the EU
QuantERA2021 project T-NiSQ, the Quantum Technology Flagship project
PASQuanS2, the NextGenerationEU project CN00000013, the Italian Research
Center on HPC, Big Data and Quantum Computing, the Quantum Computing and
Simulation Center of Padova University, the INFN project QUANTUM and by the
Italian Ministry of University and Research via PRIN2022-PNRR project TANQU
and the Rita Levi-Montalcini program.

## References

  * Chang  _et al._ (2014) D. E. Chang, V. Vuletić, and M. D. Lukin, Nature Photonics 8, 685 (2014).
  * Reiserer and Rempe (2015) A. Reiserer and G. Rempe, Reviews of Modern Physics 87, 1379 (2015).
  * Firstenberg  _et al._ (2016) O. Firstenberg, C. S. Adams, and S. Hofferberth, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 152003 (2016).
  * Firstenberg  _et al._ (2013) O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gorshkov, M. D. Lukin, and V. Vuletić, Nature 502, 71 (2013).
  * Liang  _et al._ (2018) Q.-Y. Liang, A. V. Venkatramani, S. H. Cantu, T. L. Nicholson, M. J. Gullans, A. V. Gorshkov, J. D. Thompson, C. Chin, M. D. Lukin, and V. Vuletić, Science 359, 783 (2018).
  * Sheremet  _et al._ (2023) A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, Reviews of Modern Physics 95, 015002 (2023).
  * Roy  _et al._ (2017) D. Roy, C. M. Wilson, and O. Firstenberg, Reviews of Modern Physics 89, 021001 (2017).
  * Lodahl  _et al._ (2015) P. Lodahl, S. Mahmoodian, and S. Stobbe, Reviews of Modern Physics 87, 347 (2015).
  * Hood  _et al._ (2016) J. D. Hood, A. Goban, A. Asenjo-Garcia, M. Lu, S.-P. Yu, D. E. Chang, and H. Kimble, Proceedings of the National Academy of Sciences 113, 10507 (2016).
  * Corzo  _et al._ (2019) N. V. Corzo, J. Raskop, A. Chandra, A. S. Sheremet, B. Gouraud, and J. Laurat, Nature 566, 359 (2019).
  * Tiranov  _et al._ (2023) A. Tiranov, V. Angelopoulou, C. J. van Diepen, B. Schrinski, O. A. D. Sandberg, Y. Wang, L. Midolo, S. Scholz, A. D. Wieck, A. Ludwig, _et al._ , Science 379, 389 (2023).
  * Astafiev  _et al._ (2010) O. Astafiev, A. M. Zagoskin, A. Abdumalikov Jr, Y. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, Science 327, 840 (2010).
  * Brehm  _et al._ (2021) J. D. Brehm, A. N. Poddubny, A. Stehli, T. Wolz, H. Rotzinger, and A. V. Ustinov, npj Quantum Materials 6, 10 (2021).
  * Mirhosseini  _et al._ (2019) M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang, and O. Painter, Nature 569, 692 (2019).
  * Kannan  _et al._ (2023) B. Kannan, A. Almanakly, Y. Sung, A. Di Paolo, D. A. Rower, J. Braumüller, A. Melville, B. M. Niedzielski, A. Karamlou, K. Serniak, _et al._ , Nature Physics 19, 394 (2023).
  * Shen and Fan (2007a) J.-T. Shen and S. Fan, Physical Review A 76, 062709 (2007a).
  * Schrinski and Sørensen (2022a) B. Schrinski and A. S. Sørensen, New Journal of Physics 24, 123023 (2022a).
  * Shen and Fan (2007b) J.-T. Shen and S. Fan, Physical review letters 98, 153003 (2007b).
  * Zheng  _et al._ (2011) H. Zheng, D. J. Gauthier, and H. U. Baranger, Physical review letters 107, 223601 (2011).
  * Mahmoodian  _et al._ (2018) S. Mahmoodian, M. Čepulkovskis, S. Das, P. Lodahl, K. Hammerer, and A. S. Sørensen, Physical review letters 121, 143601 (2018).
  * Prasad  _et al._ (2020) A. S. Prasad, J. Hinney, S. Mahmoodian, K. Hammerer, S. Rind, P. Schneeweiss, A. S. Sørensen, J. Volz, and A. Rauschenbeutel, Nature Photonics 14, 719 (2020).
  * Mahmoodian  _et al._ (2020) S. Mahmoodian, G. Calajó, D. E. Chang, K. Hammerer, and A. S. Sørensen, Physical Review X 10, 031011 (2020).
  * Le Jeannic  _et al._ (2022) H. Le Jeannic, A. Tiranov, J. Carolan, T. Ramos, Y. Wang, M. H. Appel, S. Scholz, A. D. Wieck, A. Ludwig, N. Rotenberg, _et al._ , Nature Physics 18, 1191 (2022).
  * Tomm  _et al._ (2023) N. Tomm, S. Mahmoodian, N. O. Antoniadis, R. Schott, S. R. Valentin, A. D. Wieck, A. Ludwig, A. Javadi, and R. J. Warburton, Nature Physics , 1 (2023).
  * Zhang and Mølmer (2020) Y.-X. Zhang and K. Mølmer, Physical Review Letters 125, 253601 (2020).
  * Poddubny (2020) A. N. Poddubny, Physical Review A 101, 043845 (2020).
  * Bakkensen  _et al._ (2021) B. Bakkensen, Y.-X. Zhang, J. Bjerlin, and A. S. Sørensen, arXiv preprint arXiv:2110.06093 (2021).
  * Calajó and Chang (2022) G. Calajó and D. E. Chang, Physical Review Research 4, 023026 (2022).
  * Marques  _et al._ (2021) Y. Marques, I. Shelykh, and I. Iorsh, Physical Review Letters 127, 273602 (2021).
  * Schrinski  _et al._ (2023) B. Schrinski, J. A. Brimer, and A. S. Sørensen, arXiv preprint arXiv:2310.19115 (2023).
  * Zhang  _et al._ (2022) X. Zhang, E. Kim, D. K. Mark, S. Choi, and O. Painter, arXiv preprint arXiv:2206.12803 (2022).
  * Scigliuzzo  _et al._ (2022) M. Scigliuzzo, G. Calajò, F. Ciccarello, D. P. Lozano, A. Bengtsson, P. Scarlino, A. Wallraff, D. Chang, P. Delsing, and S. Gasparinetti, Physical Review X 12, 031036 (2022).
  * Gong  _et al._ (2021) M. Gong, S. Wang, C. Zha, M.-C. Chen, H.-L. Huang, Y. Wu, Q. Zhu, Y. Zhao, S. Li, S. Guo, _et al._ , Science 372, 948 (2021).
  * Yu  _et al._ (2019) S.-P. Yu, J. A. Muniz, C.-L. Hung, and H. Kimble, Proceedings of the National Academy of Sciences 116, 12743 (2019).
  * Hauff  _et al._ (2022) N. V. Hauff, H. Le Jeannic, P. Lodahl, S. Hughes, and N. Rotenberg, Physical Review Research 4, 023082 (2022).
  * Krinner  _et al._ (2018) L. Krinner, M. Stewart, A. Pazmino, J. Kwon, and D. Schneble, Nature 559, 589 (2018).
  * Kim  _et al._ (2023) Y. Kim, A. Lanuza, and D. Schneble, arXiv preprint arXiv:2311.09474 (2023).
  * Noh and Angelakis (2016) C. Noh and D. G. Angelakis, Reports on Progress in Physics 80, 016401 (2016).
  * Rui  _et al._ (2020) J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, Nature 583, 369 (2020).
  * Srakaew  _et al._ (2023) K. Srakaew, P. Weckesser, S. Hollerith, D. Wei, D. Adler, I. Bloch, and J. Zeiher, Nature Physics , 1 (2023).
  * (41) _See Supplementary Material._
  * Asenjo-Garcia  _et al._ (2017a) A. Asenjo-Garcia, J. Hood, D. Chang, and H. Kimble, Physical Review A 95, 033818 (2017a).
  * Asenjo-Garcia  _et al._ (2017b) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. Kimble, and D. E. Chang, Physical Review X 7, 031024 (2017b).
  * González-Tudela  _et al._ (2015) A. González-Tudela, C.-L. Hung, D. E. Chang, J. I. Cirac, and H. Kimble, Nature Photonics 9, 320 (2015).
  * González-Tudela and Cirac (2017) A. González-Tudela and J. I. Cirac, Physical Review A 96, 043811 (2017).
  * Galve and Zambrini (2018) F. Galve and R. Zambrini, Physical Review A 97, 033846 (2018).
  * González-Tudela and Galve (2018) A. González-Tudela and F. Galve, ACS Photonics 6, 221 (2018).
  * Albrecht  _et al._ (2019) A. Albrecht, L. Henriet, A. Asenjo-Garcia, P. B. Dieterle, O. Painter, and D. E. Chang, New Journal of Physics 21, 025003 (2019).
  * Zhang and Mølmer (2019) Y.-X. Zhang and K. Mølmer, Physical review letters 122, 203605 (2019).
  * Kumlin  _et al._ (2020) J. Kumlin, K. Kleinbeck, N. Stiesdal, H. Busche, S. Hofferberth, and H. P. Büchler, Physical Review A 102, 063703 (2020).
  * Ostermann  _et al._ (2019) L. Ostermann, C. Meignant, C. Genes, and H. Ritsch, New Journal of Physics 21, 025004 (2019).
  * Needham  _et al._ (2019) J. A. Needham, I. Lesanovsky, and B. Olmos, New Journal of Physics 21, 073061 (2019).
  * Zhong  _et al._ (2020) J. Zhong, N. A. Olekhno, Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Physical review letters 124, 093604 (2020).
  * Kornovan  _et al._ (2019) D. Kornovan, N. Corzo, J. Laurat, and A. Sheremet, Physical Review A 100, 063832 (2019).
  * Schrinski and Sørensen (2022b) B. Schrinski and A. S. Sørensen, New journal of physics 24, 123023 (2022b).
  * Bettles  _et al._ (2016) R. J. Bettles, S. A. Gardiner, and C. S. Adams, Physical review letters 116, 103602 (2016).
  * Shahmoon  _et al._ (2017) E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, Physical review letters 118, 113601 (2017).
  * Manzoni  _et al._ (2018) M. Manzoni, M. Moreno-Cardoner, A. Asenjo-Garcia, J. V. Porto, A. V. Gorshkov, and D. Chang, New journal of physics 20, 083048 (2018).
  * Ozawa  _et al._ (2019) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * Salerno  _et al._ (2018) G. Salerno, M. Di Liberto, C. Menotti, and I. Carusotto, Phys. Rev. A 97, 013637 (2018).
  * Salerno  _et al._ (2020) G. Salerno, G. Palumbo, N. Goldman, and M. Di Liberto, Phys. Rev. Res. 2, 013348 (2020).
  * Clark  _et al._ (2020) L. W. Clark, N. Schine, C. Baum, N. Jia, and J. Simon, Nature 582, 41 (2020).
  * Léonard  _et al._ (2023) J. Léonard, S. Kim, J. Kwan, P. Segura, F. Grusdt, C. Repellin, N. Goldman, and M. Greiner, Nature 619, 495 (2023).
  * Perczel  _et al._ (2017a) J. Perczel, J. Borregaard, D. E. Chang, H. Pichler, S. F. Yelin, P. Zoller, and M. D. Lukin, Physical review letters 119, 023603 (2017a).
  * Perczel  _et al._ (2017b) J. Perczel, J. Borregaard, D. E. Chang, H. Pichler, S. F. Yelin, P. Zoller, and M. D. Lukin, Physical Review A 96, 063801 (2017b).
  * Perczel  _et al._ (2020) J. Perczel, J. Borregaard, D. E. Chang, S. F. Yelin, and M. D. Lukin, Physical Review Letters 124, 083603 (2020).
  * Bergholtz and Liu (2013) E. J. Bergholtz and Z. Liu, International Journal of Modern Physics B 27, 1330017 (2013).
  * Semeghini  _et al._ (2021) G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Science 374, 1242 (2021).
  * Brechtelsbauer and Malz (2021) K. Brechtelsbauer and D. Malz, Physical Review A 104, 013701 (2021).
  * Andreoli  _et al._ (2023) F. Andreoli, B. Windt, S. Grava, G. M. Andolina, M. J. Gullans, A. A. High, and D. E. Chang, arXiv preprint arXiv:2303.10998 (2023).
  * Carmichael (1999) H. J. Carmichael,  _Statistical Methods in Quantum Optics 1_ (Springer Berlin Heidelberg, 1999).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University PressOxford, 2007).
  * Economou (2006) E. N. Economou,  _Green’s functions in quantum physics_ , Vol. 7 (Springer Science & Business Media, 2006).
  * Lehmberg (1970) R. Lehmberg, Physical Review A 2, 883 (1970).
  * Dung  _et al._ (1998) H. T. Dung, L. Knöll, and D.-G. Welsch, Physical Review A 57, 3931 (1998).

Strongly interacting photons in 2D Waveguide QED

– SUPPLEMENTARY INFORMATION –

Interacting photons in 2D waveguide QEDM. Tečer M. Di Liberto P. Silvi S.
Montangero F. Romanato G. Calajó

##  Appendix S1 2D Spin model derivation

Let us consider a 2D atomic ensemble interacting with the light confined into
a 2D photonic structure, as illustrated in Figure 1 of the main text. To
address the main physics originating from the 2D confinement we focus on the
ideal waveguide QED scenario assuming a perfect coupling of the atoms to the
2D structure. With this assumption we completely neglect spontaneous emission
of the atoms into the free space. While this approximation may seem strong, it
is frequently used in waveguide Quantum Electrodynamics (QED) literature
Sheremet  _et al._ (2023) and it accurately reflects the conditions
encountered in scenarios involving superconducting qubits coupled to arrays of
microwave resonators. In such cases, emission outside the resonators is
typically two to three orders of magnitude smaller compared to emissions into
other channels Scigliuzzo  _et al._ (2022). The Hamiltonian that describes
this coupled atom-photonic structure system, under the rotating-wave and
dipole approximations, is expressed as:

$\displaystyle\begin{split}&\hat{H}=\hat{H}_{0}+\hat{H}_{\rm int}\,,\\\
&\hat{H}_{0}=\hbar\sum_{i}\omega_{a}\hat{\sigma}_{+}^{i}\hat{\sigma}_{-}^{i}+\hbar\int
d^{2}{\mathbf{k}}\;\omega({\mathbf{k}})\hat{a}({\mathbf{k}})^{\dagger}\hat{a}({\mathbf{k}})\,,\\\
&\hat{H}_{\rm int}=\sum_{i}\int
d^{2}\mathbf{k}\left(g_{\mathbf{k}}(\mathbf{r}_{i})\hat{a}({\mathbf{k})}^{\dagger}\hat{\sigma}_{-}^{i}e^{i\mathbf{k}\cdot\mathbf{x}_{i}}+\text{h.c}.\right)\,,\end{split}$
(S1)

where $\hat{a}(\mathbf{k})^{\dagger}(\hat{a}(\mathbf{k})$) are the bosonic
creation (annihilation) operators of the photonic structure and

$g_{\mathbf{k}}(\mathbf{r}_{i})=\sqrt{\frac{\omega_{\mathbf{k}}}{2\hbar\epsilon_{0}}}\mathbf{d}^{i}_{eg}\cdot\mathbf{\Phi}_{\mathbf{k}}(\mathbf{x}_{i})\,,$
(S2)

is the atom-photon coupling constant with $\epsilon_{0}$ being the vacuum
permittivity, $\mathbf{d}^{i}_{eg}=\langle e|\hat{\mathbf{d}}^{i}|g\rangle$
the dipole moment of the i-th atom and
$\mathbf{\Phi}_{\mathbf{k}}(\mathbf{x}_{i})$ the spatial mode eigenfunction of
the field normalized according:

$\int
d^{2}\mathbf{x}\mathbf{\Phi}_{\mathbf{k}}(\mathbf{x})\cdot\mathbf{\Phi}_{\mathbf{k^{\prime}}}(\mathbf{x})=\delta_{\mathbf{k}\mathbf{k^{\prime}}}.$
(S3)

To simplify the description of this complex interacting system, we assume that
the standard conditions for employing the Born-Markov approximation are
satisfied. These conditions are guaranteed if the density of states of the
photonic bath is smooth and if time retardation effects can be neglected.
Consequently, we can eliminate the photonic degree of freedom, resulting in
the following master equation for the atoms: Carmichael (1999); Breuer and
Petruccione (2007):

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}\Big{[}\hat{H}_{\rm
eff}\hat{\rho}-\hat{\rho}\hat{H}_{\rm
eff}^{\dagger}\Big{]}+2\sum_{i,j}\Gamma_{ij}\hat{\sigma}_{-}^{i}\hat{\rho}\hat{\sigma}_{+}^{j}\,,$
(S4)

where we have defined the effective non-hermitian Hamiltonian:

$\hat{H}_{\rm
eff}=\sum_{i,j}G_{ij}\hat{\sigma}_{+}^{i}\hat{\sigma}_{-}^{j}\,.$ (S5)

Assuming a smooth position dependence of the atomic coupling trough the whole
atomic ensemble $g_{\mathbf{k}}(\mathbf{r}_{i})\simeq g_{\mathbf{k}}$ we
finally obtain the following expression:

$G_{ij}=\,\lim_{s\to
0}\int\frac{d^{2}\mathbf{k}}{2(2\pi)^{2}}\frac{|g_{\mathbf{k}}|^{2}}{is-(\omega_{a}-\omega_{\mathbf{k}})}e^{i\mathbf{k}\cdot\mathbf{x}_{ij}},$
(S6)

where $\mathbf{x}_{ij}=\mathbf{x}_{i}-\mathbf{x}_{j}$ .

Figure S1:  (a) The photonic structure’s dispersion relation in the vicinity
of a band edge is approximated to be quadratic and isotropic. The atomic
transition frequency is located in the band
($\Delta=\omega_{a}-\omega_{e}>0$). (b) Schematic representation of the ansatz
used to map the two-excitation sector into the relative coordinate frame. To
avoid double counting, we consider the relative positions of the green
highlighted atoms only with respect to the atoms in the green shaded area.

###  S1.1 Isotropic and quadratic dispersion

To simplify the discussion we assume that the photonic dispersion is isotropic
and quadratic. This choice is motivated by the common occurrence of such
dispersion in 2D confined photonic structures like 2D photonic crystals Yu
_et al._ (2019) or above the band edge of a square photonic lattice González-
Tudela and Cirac (2017) (as illustrated in Fig.S1(a)). Under this assumption,
the dispersion explicitly takes the form
$\omega(\mathbf{k})=\omega_{e}+A|\mathbf{k}|^{2}$, where $\omega_{e}$ is the
frequency at the band edge, and we consider $A>0$. If the atomic frequency
$\omega_{a}$ lies inside the band, the atoms are coupled to propagating
photons, and the collective interaction term is given by:

$G_{ij}=\lim_{s\to
0}\int_{0}^{k_{c}}\frac{k\,dk}{2(2\pi)^{2}}\frac{|g_{\mathbf{k}}|^{2}}{is-(\Delta-
Ak^{2})}\int_{0}^{2\pi}e^{ik|\mathbf{x}_{ij}|\cos\phi}d\phi=\lim_{s\to
0}\int_{0}^{k_{c}}\frac{k\,dk}{2(2\pi)}\frac{|g_{k}|^{2}}{is-(\Delta-
Ak^{2})}\mathcal{J}_{0}(k|\mathbf{x}_{ij}|)\,,$ (S7)

where $\Delta=\omega_{a}-\omega_{e}>0$ is the atom-band edge detuning, $k_{c}$
sets an higher frequency cutoff and $\mathcal{J}_{0}(x)$ is the zeroth order
Bessel function of the first kind. Performing the change of variable
$q=\sqrt{\frac{\Delta}{A}}k$ and assuming the coupling to be constant across
the integration domain, $g_{k}\approx g$, we finally obtain González-Tudela
_et al._ (2015):

$G_{ij}\approx\lim_{s\to
0}\frac{|g|^{2}}{A}\int_{0}^{\infty}\frac{qdq}{4\pi}\frac{\mathcal{J}_{0}\left(q|\mathbf{x}_{ij}|k_{0}\right)}{is-(1-q^{2})}=\frac{\gamma}{2}\left(\mathcal{Y}_{0}\left(k_{0}|\mathbf{r}_{ij}|\right)-i\mathcal{J}_{0}\left(k_{0}|\mathbf{r}_{ij}|\right)\right)\,.$
(S8)

where $\mathcal{Y}_{0}(x)$ is the zeroth order Bessel function of the second
kind, $k_{0}=\sqrt{\frac{\Delta}{A}}$ is the resonant wavevector and we
assumed $k_{c}\sqrt{\frac{A}{\Delta}}\gg 1$. In this form, $G_{ij}$ is
explicitly composed of two contributions, accounting for coherent (first term)
and dissipative (second term) collective interactions. We also define
$\gamma=\frac{|g|^{2}}{4A}$ as the spontaneous emission decay rate of an
individual atom coupled to a 2D photonic bath. This decay rate establishes the
scale of the system dynamics. Notably, the earlier assumption of neglecting
time retardation effects can be reformulated as $\gamma\ll c/L$, where
$c=2Ak_{0}$ is the photon group velocity at the atomic frequency and $L$ is
the size of the atomic ensemble.

###  S1.2 Single excitation dispersion

In the thermodynamic limit, the eigenstates of the effective Hamiltonian (S5)
with single excitations are spin waves characterized by well-defined
quasimomentum $\mathbf{k}$. To diagonalize the Hamiltonian, a Bloch wave
ansatz is employed, utilizing the Fourier transform of spin operators
$\hat{\sigma}_{-}^{n}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{x}_{n}}\hat{\sigma}_{-}^{\mathbf{k}}$
($\hbar=1$):

$\hat{H}_{\rm{eff}}=\sum_{n,m}G(k_{0}|\mathbf{x}_{m}-\mathbf{x}_{n}|)\hat{\sigma}_{+}^{m}\hat{\sigma}_{-}^{n}=\frac{1}{N}\sum_{n,m}G(k_{0}|\mathbf{x}_{m}-\mathbf{x}_{n}|)\sum_{\mathbf{k}^{\prime}\mathbf{k}}\hat{\sigma}_{+}^{\mathbf{k}^{\prime}}\hat{\sigma}_{-}^{\mathbf{k}}e^{i(\mathbf{k}\cdot\mathbf{x}_{n}-\mathbf{k}^{\prime}\cdot\mathbf{x}_{m})}\,.$
(S9)

The summation over atoms is then split into three cases: $m>n$, $m=n$, and
$m<n$, where the relations $><$ are defined on the 2D lattice, as illustrated
in Fig. S1(b) and explained in the next section. Introducing the relative
coordinate $\mathbf{r}=\mathbf{x}_{m}-\mathbf{x}_{n}$ for $m>n$ (and similarly
for $m<n$), the expression becomes:

$\displaystyle\hat{H}_{\rm{eff}}=\frac{1}{N}\sum_{n=1}^{N}\sum_{\mathbf{k},\mathbf{k^{\prime}}\in
BZ}\sum_{\mathbf{r}}G(k_{0}|\mathbf{r}|)\,e^{i\mathbf{x}_{n}\cdot(\mathbf{k}-\mathbf{k}^{\prime})}\Big{(}e^{-i\mathbf{k}^{\prime}\cdot\mathbf{r}}+e^{i\mathbf{k}\cdot\mathbf{r}}\Big{)}\hat{\sigma}_{+}^{\mathbf{k}^{\prime}}\hat{\sigma}_{-}^{\mathbf{k}}+G(0)\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
BZ}\hat{\sigma}_{+}^{\mathbf{k}^{\prime}}\hat{\sigma}_{-}^{\mathbf{k}}\,\cdot\frac{1}{N}\sum_{n=1}^{N}e^{i\mathbf{x}_{n}\cdot(\mathbf{k}-\mathbf{k}^{\prime})}\,.$
(S10)

Combining all contributions and using
$\lim_{N\to\infty}\frac{1}{N}\sum_{n}e^{i\mathbf{x}_{n}\cdot(\mathbf{k}-\mathbf{k}^{\prime})}=\delta(\mathbf{k}-\mathbf{k^{\prime}})$,
the final expression is:

$\hat{H}_{\rm{eff}}=\sum_{\mathbf{k}}\epsilon^{(1)}_{\mathbf{k}}\hat{\sigma}_{+}^{\mathbf{k}}\hat{\sigma}_{-}^{\mathbf{k}},$
(S11)

where the single-excitation dispersion relation is defined as:

$\epsilon^{(1)}_{\mathbf{k}}=-i\frac{\gamma}{2}+\frac{2}{N}\sum_{\mathbf{r}}G(k_{0}|\mathbf{r}|)\cos(\mathbf{k}\cdot\mathbf{r})\left(\sqrt{N}-|r_{x}|\right)\left(\sqrt{N}-r_{y}\right)\,,$
(S12)

with the sum over $\mathbf{r}$ taking values
$\mathbf{r}\in\Big{\\{}\big{(}r_{x}\in\big{[}1,\sqrt{N}-1\big{]},r_{y}=0\big{)}\Big{\\}}\,\cup\,\Big{\\{}\big{(}r_{x}\in\big{[}-\sqrt{N}+1,\sqrt{N}-1\big{]},r_{y}\in\big{[}1,\sqrt{N}-1\big{]}\big{)}\Big{\\}}$.

##  Appendix S2 Relative coordinate frame

In this section, we tackle the two-body problem outlined in Hamiltonian Eq.
(S5) within the relative coordinate framework. Our primary goal is to find
solutions for $H_{\rm eff}|\psi^{(2)}\rangle=E|\psi^{(2)}\rangle$, where the
two-excitation eigenstate generally adopts the form
$|\psi^{(2)}\rangle=\sum_{i<j}c^{(2)}{ij}|e_{i},e_{j}\rangle$, and
$|e_{i},e_{j}\rangle=\hat{\sigma}+^{i}\hat{\sigma}_{+}^{j}|0\rangle$. For
convenience, we re-parameterize the problem in terms of center-of-mass,
$\mathbf{R}=(R_{x},R_{y})$, and relative, $\mathbf{r}=(r_{x},r_{y})$,
coordinates.

In the case of an infinite system, we can take the thermodynamic limit and
apply Bloch’s theorem to express the center-of-mass wave function in relation
to the center-of-mass momentum $\mathbf{K}$. Compared to the analogous 1D wQED
case Bakkensen  _et al._ (2021); Zhang and Mølmer (2020); Calajó and Chang
(2022), the 2D relative coordinate ansatz is more involved, and it is
represented by the following expression, schematically illustrated in
Fig.S1(b),

$\displaystyle\begin{split}&\big{|}\mathbf{K},\mathbf{r}\big{\rangle}=\mathcal{N}_{\mathbf{r}}\theta(r_{x})\sum_{j_{y}=1}^{L-r_{y}}\sum_{j_{x}=1}^{L-r_{x}}\big{|}\phi\left(\mathbf{K},\mathbf{r},j_{x},j_{y}\right)\big{\rangle}+\mathcal{N}_{\mathbf{r}}\theta(-r_{x})(1-\delta(r_{y}))\sum_{j_{y}=1}^{L-r_{y}}\sum_{j_{x}=1+|r_{x}|}^{L}\big{|}\phi\left(\mathbf{K},\mathbf{r},j_{x},j_{y}\right)\big{\rangle},\\\
&\big{|}\phi\left(\mathbf{K},\mathbf{r},j_{x},j_{y}\right)\big{\rangle}=e^{i\frac{K_{x}d(2j_{x}+r_{x})}{2}+i\frac{K_{y}d(2j_{y}+r_{y})}{2}}\big{|}(j_{x},j_{y}),(j_{x}+r_{x},j_{y}+r_{y})\big{\rangle}\,.\end{split}$
(S13)

In the provided expression, the vector $|(i_{x},i_{y}),(j_{x},j_{y})\rangle$
denotes a state where the atoms at the sites $(i_{x},i_{y})$ and
$(j_{x},j_{y})$ are occupied. Here, $N=L\times L$ represents the system size,
and $\mathcal{N}_{\mathbf{r}}$ is the norm of the state
$\big{|}\mathbf{K},\mathbf{r}\big{\rangle}$. To ensure the comprehensive
counting of all two-excitation states without duplication, we assume that the
relative coordinate can assume values
$r_{y}\in[0,L-1],\,r_{x}\in[-(L-1),L-1]$. Applying Eq. (S5) to the state
(S13), we derive the following matrix elements:

$\langle\mathbf{K},\mathbf{r}|\hat{H}_{\textrm{eff}}|\mathbf{K},\mathbf{r}^{\prime}\rangle=\sum_{\epsilon=\pm}\cos\mathbf{K}\cdot\left(\mathbf{r}+\epsilon\mathbf{r}^{\prime}\right)G\left(k_{0}\left|\mathbf{r}+\epsilon\mathbf{r}^{\prime}\right|\right)\,.$
(S14)

In this manner, we have transformed the original two-body problem into a
single-body problem in the relative coordinate, which depended parametrically
from $\mathbf{K}$.

##  Appendix S3 Mapping to the impurity model

As explained in the main text, the spin model Hamiltonian presented in Eq.
(S5) can be transformed into a hardcore boson model. This transformation
involves promoting the spin operators to bosons, denoted as
$\hat{\sigma}^{i}_{-}\rightarrow\hat{b}_{i}$, and introducing an on-site
interaction with an infinite strength $U\rightarrow\infty$. This results in
the following expression ($\hbar=1$):

$\hat{H}_{\rm eff}=\hat{H}_{0}+\hat{H}_{\rm
int}=\sum_{ij}G_{ij}\hat{b}_{i}^{\dagger}\hat{b}_{j}+\frac{U}{2}\sum_{i}(\hat{b}_{i}^{\dagger})^{2}\hat{b}_{i}^{2},$
(S15)

where, the bosonic operators $\hat{b}_{i}$ and $\hat{b}_{i}^{\dagger}$
respectively annihilate and create excitations on the site $\mathbf{x}_{i}$.
In the two-excitation sector, we can again map the problem into the center-of-
mass, $\mathbf{R}$, and relative coordinates, $\mathbf{r}$, with the center-
of-mass momentum $\mathbf{K}$ being a well-defined quantum number in the
thermodynamic limit. The first term in Hamiltonian (S15) accounts for free
bosons, and its energy is straightforwardly the sum of individual photon
energies. This term can be expressed in a diagonal form with respect to the
relative momentum state $|\mathbf{K},\mathbf{q}\rangle$.

$\hat{H}_{0}=\sum_{\mathbf{q}}\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})|\mathbf{K},\mathbf{q}\rangle\langle\mathbf{K},\mathbf{q}|$
(S16)

where
$\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})=\epsilon^{(1)}(\mathbf{q})+\epsilon^{(1)}(\mathbf{K}-\mathbf{q})$.
We can Fourier transform this Hamiltonian with respect to the relative
coordinate using the relation
$|\mathbf{K},\mathbf{q}\rangle=(1/\sqrt{N})\sum_{\mathbf{r}}e^{i\mathbf{q}\mathbf{r}}|\mathbf{K},\mathbf{r}\rangle$.
By incorporating the interaction term mapped to the relative coordinate space,
we arrive at the impurity model Hamiltonian presented in the main text.

$\hat{H}_{\mathbf{K}}=\sum_{\mathbf{r}\mathbf{r}^{\prime}}J^{\mathbf{K}}_{\mathbf{r}\mathbf{r}^{\prime}}|\mathbf{K},\mathbf{r}\rangle\langle\mathbf{K},\mathbf{r}^{\prime}|+U|\mathbf{K},0\rangle\langle\mathbf{K},0|,$
(S17)

where
$J^{\mathbf{K}}_{\mathbf{r}\mathbf{r}^{\prime}}=(1/N)\sum_{\mathbf{q}}\omega_{\mathbf{K}}(\mathbf{q})e^{i(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{q}}$
and $|\mathbf{K},0\rangle=|\mathbf{K},\mathbb{r}=0\rangle$ is the site where
is placed the impurity.

As discussed in Economou (2006), it is possible to derive a closed expression
for the single-particle Green function in this particular problem:

$\mathcal{G}=\mathcal{G}_{0}+\mathcal{G}_{0}|\mathbf{K},0\rangle\frac{U}{1-U\mathcal{G}_{0}(0,E)}\langle\mathbf{K},0|\mathcal{G}_{0}\,,$
(S18)

where $\mathcal{G}_{0}$ is the Green function of the non interacting
Hamiltonian $\hat{H}_{0}$ and

$\mathcal{G}_{0}(\mathbf{r},E)=\langle\mathbf{K},\mathbf{r}|\mathcal{G}_{0}(E)|\mathbf{K},\mathbf{0}\rangle=\frac{1}{(2\pi)^{2}}\int_{\rm
BZ}d^{2}\mathbb{q}\frac{e^{i\mathbb{q}\cdot\mathbb{r}}}{E-\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})}\,.$
(S19)

The poles of the Green function, characterized by energies $E_{\rm BS}$,
correspond to the bound state solutions of the system. The equation governing
these bound states is expressed as:

$\mathcal{G}_{0}(0,0,E_{\rm BS})=\frac{1}{U}\rightarrow 0\,.$ (S20)

This equation can be numerically solved by seeking solutions in the gap for
each $\mathbb{K}$ value:

$\frac{1}{(N)^{2}}\sum_{\mathbb{q}\in_{\rm BZ}}\frac{1}{E_{\rm
BS}-\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})}=0\,,$ (S21)

where the sum over the relative momentum spans the Brillouin zone. Finally
with the Green function formalism, we can also express the solution for the
unnormalized bound states in the relative coordinate:

$|\psi_{\rm
BS},\mathbf{K}\rangle=\sum_{\mathbb{r}}\mathcal{G}_{0}(\mathbb{r},E_{\rm
BS})|\mathbf{K},\mathbb{r}\rangle.$ (S22)

Figure S2: (a)-(b)Single photon dispersion relation at various inter-atomic
distances both before and after the gap closing at $k_{0}d=\pi$. Panels
(c)-(e) illustrate the two-excitation dispersion of unbound photons, denoted
as $\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})$, as a function of the
relative momentum $\mathbf{q}$ with the center-of-mass momentum fixed at the
$M$-point,$\mathbf{K}=\left(\frac{\pi}{d},\frac{\pi}{d}\right)$. In the first
row (1), the dispersions are plotted within the first Brillouin zone, while in
the second row (2), they are unfolded with respect to the $\Gamma$, $X$, and
$M$ symmetry points. All plots are generated for a system size of $N=300\times
300$.

##  Appendix S4 Bandgap closing

As explained in the main text, the energy of bound states, denoted as $E_{\rm
BS}$, must differ from the sum of individual photon energies, given by:

$E_{\rm
BS}(\mathbf{K})\neq\epsilon^{(1)}(\mathbf{q})+\epsilon^{(1)}(\mathbf{K}-\mathbf{q})\,.$
(S23)

This condition cannot be satisfied in the absence of a gap in the single
excitation spectrum. In Fig.1 of the main text, we demonstrated the closure of
the single excitation bandgap above the sub-wavelength limit, $k_{0}d>\pi$.
This closure is visibly apparent by examining the single-particle dispersion
at different inter-atomic distances, as illustrated in Fig. S2(a)-(b). When
$k_{0}d>\pi$, the central polaritonic branch with a radius of
$|\mathbf{k}|=k_{0}$ reaches the edges of the 1BZ, and higher diffraction
orders prevent gap existence. Given that the condition given in Eq. (S23)
cannot be satisfied for $k_{0}d>\pi$, we focus on studying the unbound part of
the two-excitation spectrum
$\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})=\epsilon^{(1)}(\mathbf{q})+\epsilon^{(1)}(\mathbf{K}-\mathbf{q})$
for different values of $k_{0}d<\pi$ as a function of the relative momentum
$\mathbf{q}$. This allows us to determine the regions of validity for (S23).
At the $\Gamma$ point ($\mathbf{K}=\left(0,0\right)$), the gap closes
precisely at the same distance as for a single excitation because trivially
$\epsilon^{(2)}_{\text{scat}}(\mathbf{K},\mathbf{q})=2\epsilon^{(1)}({\mathbf{q}})$.
At the $M$ point ($\mathbf{K}=\left(\frac{\pi}{d},\frac{\pi}{d}\right)$), the
dispersion reveals multiple polaritonic branches with a radius of
$|\mathbf{q}|=k_{0}$, positioned in the middle and at the edges of the
Brillouin zone. The bandgap then closes when these polaritonic branches
intersect for $k_{0}d\geq\frac{\pi}{\sqrt{2}}$, as illustrated in Fig. S2.

Figure S3:  (a) The distribution of relative coordinates for $r_{x}=r_{y}=r$
is depicted for the $M$ point bound state (BS), marked by orange dots, along
with two scattering resonances at the $M$ point. SR1, represented by blue
triangles, denotes the scattering resonance found on the same energy branch as
the bound state when the bandgap closes. SR2, indicated by red squares,
represents the scattering resonance from the other energy branch located in
the low density of states region of the 2-excitation dispersion spectrum (see
Fig. 2 of the main text). (b) The inverse participation ratio is shown as a
function of system size for the same states as in panel (a). In all the plots
we considered a system size of $N=60\times 60$

##  Appendix S5 Localized scattering resonances

In the main text, we identified states within the band that exhibit
localization of the two excitations respect to their relative coordinates.
These states, lacking protection from the band gap, cannot be bound states;
instead, they are scattering resonances. Comprising a bound and unbound
component, these states are quantitatively distinct from bound states. This
distinction is evident in Fig. S3(a), where we plotted the spread of the
relative coordinate population distribution for $r_{x}=r_{y}$ for the two
scattering resonance branches discussed in the main text. The plot illustrates
that bound states exhibit an exponential decay with the relative distance,
whereas scattering resonances display a non-exponentially oscillating tail.
Note that a similar behavior has been discussed also for 1D waveguides
Bakkensen  _et al._ (2021). To further verify the nature of these states, in
Fig. S3(b) we present the inverse participation ratio defined as

$IPR=\frac{1}{\sum_{\mathbf{r}}\big{|}\psi(\mathbf{r})\big{|}^{4}}\,,$ (S24)

as a function of the system size. As observed, for bound states at
$k_{0}d=0.32\pi$, the IPR saturates to a constant value for $L\gtrapprox 40$.
Conversely, for scattering resonances at $k_{0}d=0.74\pi,0.9\pi$, the IPR
oscillates without reaching an asymptotic value, thereby confirming the bound
and scattering resonance nature of these states, respectively.

Figure S4: 2D free space array with atomic polarization orthogonal to the
array plane. $\mathrm{(a_{1}),(a_{2})}$ Single excitation dispersion for
$k_{0}d=1.09\pi,\,0.6\pi$. $\mathrm{(b_{1}),(b_{2})}$ Two excitation spectrum
at the $\Gamma$ and $M$ point of CM momenta as a function of inter-atomic
distance. Blue dots represent the sum of the free photon energies, while the
red (orange) line represents the energy of the scattering resonance. (c)
Probability distribution of the $\Gamma$ and $M$ point scattering resonance
for $k_{0}d=0.6\pi,\,1.09\pi$. (d) Decay rate of $\Gamma$ and $M$ point
scattering resonances as a function of $k_{0}d$ for system size $N=10\times
10$. In panels (a),(b), we considered a system size of $150\times 150$, while
in panel (c) we considered a system size $60\times 60$.

##  Appendix S6 Interacting states in 2D free space atomic arrays

The dyadic Green’s tensor, which characterizes the interaction between dipoles
in free space employed in the main text, is derived as a solution to the
dyadic equation:Lehmberg (1970); Dung  _et al._ (1998); Asenjo-Garcia  _et
al._ (2017a); Perczel  _et al._ (2017a):

$\frac{\omega_{A}^{2}}{c^{2}}\tilde{G}_{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})-\left(\partial_{\alpha}\partial_{\nu}-\delta_{\alpha\nu}\partial_{\eta}\partial_{\eta}\right)\tilde{G}_{\nu\beta}(\mathbf{r},\mathbf{r}^{\prime})=\delta_{\alpha\beta}\delta(\mathbf{r}-\mathbf{r}^{\prime})\,.$
(S25)

where the Green’s tensor, denoted as
$\tilde{G}_{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})$, characterizes the
component $\alpha={x,y,z}$ of the electromagnetic field at position
$\mathbf{r}$ emitted by a perfect dipole positioned at $\mathbf{r}^{\prime}$
and oriented in the $\beta={x,y,z}$ direction. The equation solution is
explicitly given by:

$\tilde{G}_{\alpha\beta}(\mathbf{r}-\mathbf{r}^{\prime})=\frac{e^{ik_{0}r}}{4\pi
k_{0}^{2}r^{3}}\Bigg{[}\left(1-i(k_{0}r)-(k_{0}r)^{2}\right)\delta_{\alpha\beta}+\left(1+\frac{3i}{k_{0}r}-\frac{3}{(k_{0}r)^{2}}\right)k_{0}^{2}r_{\alpha}r_{\beta}\Bigg{]},$
(S26)

where we have defined $r=|\mathbf{r}-\mathbf{r}^{\prime}|$, the resonant wave
vector $k_{0}=\frac{\omega_{a}}{c}$, and with $r_{\alpha}$
($\alpha=\\{x,y,z\\}$) we indicated the different components of the vector
$\mathbf{r}-\mathbf{r}^{\prime}$ in the Cartesian coordinates. In the main
text, we configure the atomic polarization to be oriented perpendicular to the
array’s plane. In this particular orientation, it is sufficient to consider
only the $\tilde{G}{zz}$ component of the Green’s tensor. This component
yields the collective interaction function used in the main text, denoted as
$G(x)=\frac{3\pi\gamma_{0}}{k_{0}}\tilde{G}{zz}$ with
$\gamma_{0}=d_{eg}^{2}k_{0}^{3}/(3\pi\hbar\epsilon_{0})$ being the free space
spontaneous emission decay rate and $d_{eg}$ the dipole moment strength. Once
identified the collective interaction term let’s now discuss the single and
two excitation spectrum. In contrast to the 2D confined scenario, the single
excitation spectrum of the free space array lacks a bandgap, as shown in Fig.
S4(a). This absence of a bandgap extends to the two-excitation spectrum,
implying the nonexistence of bound states in this system. Nonetheless,
subradiant repulsive and scattering resonances (SRs) persist in this setting
within the subradiant region enclosed in the first BZ, $k_{0}d<\sqrt{2}\pi$.
The former can be identified in a manner similar to the 2D confined case and
exhibit a comparable probability distribution of the two excitations. The
latter are more involved, as depicted in Fig. S4(b), where we illustrate the
two-excitation spectrum for the $\Gamma$ and $M$ points, highlighting the SRs
branches. Similar to the confined case, the existence of spectral regions with
a low density of states can result in highly bound SRs, as demonstrated in
Fig. S4(c) through the relative coordinate probability distribution of the two
SRs. Importantly, these SRs exhibit subradiant behavior even if they are
composed of two attractive excitations that can potentially scatter to each
other in free space. Their decay rate is strongly dependent on the inter-
atomic distances, as depicted in Fig. S4(d), showing pronounced minima,
respectively at $k_{0}d=1.09\pi$ and $k_{0}d=0.73\pi$ (highlighted in Fig.
S4(b) by the black arrow), where it becomes more than ten times smaller than
the single atom emission rate. Importantly, both critical distances are
moderately subwavelength making this regime potentially accessible in cold
atom arrays Rui  _et al._ (2020).

Note that the behaviour of the SR at the $\Gamma$ point is more subtle.
Indeed, this point is known to host subradiant states in single excitation
only for orthogonal atomic polarization. In this setting, all the atoms
oscillate in phase and orthogonally to the array, limiting the emission
strictly within that plane Asenjo-Garcia  _et al._ (2017a). Considering this,
the extended lifetime observed in the scattering resonance at the $\Gamma$
point appears to be connected to the remarkably flat dispersion at this point
(see Fig. S4 $\mathrm{(a_{1})}$). This flat dispersion induces localization in
the center of mass wavefunction, leading to a reduction in emissions at the
edges (see Poddubny (2020)).

Figure S5: (a) The projection of the state $|\psi(t)\rangle$ onto the
eigenstates in the two-excitation sector is displayed at times $\gamma
t=0,20,300$. The most prominent projection comes from the states with center-
of-mass momentum close to the $M$ point (labeled by $M$ BSs in the second
panel). The peak at $J_{s}\approx 2$ corresponds to the $\Gamma$ point bound
state, while the group of peaks at $J_{s}\approx 0$ corresponds to other long-
lived states. In all plots, the parameter is set as $k_{0}d=0.52\pi$, and an
array size of $N=10\times 10$ is considered.

##  Appendix S7 Dynamical excitation of the bound states

As detailed in the main text, we investigate the dynamical excitation of bound
states by solving the time-dependent Schrödinger equation governed by the
effective Hamiltonian (S5), projected in the double excitation sector:

$i\frac{\partial}{\partial t}|\psi^{(2)}(t)\rangle=H_{\rm
eff}|\psi^{(2)}(t)\rangle\,.$ (S27)

In this context, the initial state is set as
$|\psi(t=0)\rangle=\hat{\sigma}_{+}^{\bar{i}}\hat{\sigma}_{+}^{\bar{j}}|0\rangle$,
with the indices $\bar{i}$ and $\bar{j}$ indicating the positions of the
initially excited atoms separated by a distance $\bm{\ell}=(l,l)$. The
evolution of the state according to (S27) neglects the recycling term in the
master equation (S4). Consequently, the norm of the state is not conserved in
time, $\big{|}\langle\psi(t)|\psi(t)\rangle\big{|}^{2}<1$, as the population
of the single and zero excitation subspaces is not monitored during the state
relaxations. This approach yields accurate results as long as observables are
computed within the two-excitation subspace.

In the main text, we examine the excitation of the M point bound state
starting from an initial state with $\ell=1$, where two next-nearest neighbor
atoms (along the diagonal of the square plaquet) are excited. This initial
condition, when projected onto the (incomplete) two-excitation sector non-
Hermitian basis, significantly overlaps with the bound states having a center-
of-mass momentum $\mathbf{K}$ close to the $M$ point, as illustrated in the
first panel of Fig. S5(a). In this figure, an initial non-zero projection with
other sub-radiant (in blue) and super-radiant states (in orange) is also
observed. As the state evolves (second and third panel of Fig. S5(a)),
contributions from other states diminish, while those from the $M$ point bound
states persist over time. A weak excitation of other long-living states is
also noted. This includes the $\Gamma$ point bound state (less subradiant than
the M point) located at energy $J_{s}\approx 2$, where $s$ is an index
labeling the states and a contribution from other subradiant states located at
$J_{s}\approx 0$.