# Strongly subradiant states in planar atomic arrays

Ilya Volkov  School of Physics and Engineering, ITMO University, St.
Petersburg 197101, Russia  Nikita Ustimenko  School of Physics and
Engineering, ITMO University, St. Petersburg 197101, Russia  Danil Kornovan
School of Physics and Engineering, ITMO University, St. Petersburg 197101,
Russia  Roman Savelev  School of Physics and Engineering, ITMO University, St.
Petersburg 197101, Russia  Alexandra Sheremet  School of Physics and
Engineering, ITMO University, St. Petersburg 197101, Russia  Mihail Petrov
<EMAIL_ADDRESS>School of Physics and Engineering, ITMO University,
St. Petersburg 197101, Russia

###### Abstract

The optically trapped ensembles of atoms provide a versatile platform for
storing and coherent manipulation of quantum information. However, efficient
realization of quantum information processing requires long-lived quantum
states protected from the decoherence e.g. via spontaneous emission. Here, we
theoretically study collective dipolar oscillations in finite planar arrays of
quantum emitters in free space and analyze mechanisms that govern the
emergence of strongly subradiant collective states. We demonstrate that the
external coupling between the collective states associated with the symmetry
of the array and with the quasi-flat dispersion of the corresponding infinite
lattice plays a crucial role in the boost of their radiative lifetime. We show
that among different regular arrangements of the atoms the square atomic
arrays support eigenstates with minimal radiative losses that scale with the
total number of atoms $N_{tot}$ as $\propto N_{tot}^{-5}$.

Keywords: atomic lattices; subradiance; collective quantum excitataions;
external coupling; flat bands.

##  I Introduction

In the recent years, a significant progress has been achieved in the
development of artificial quantum interfaces based on arrays of cold atoms
Barredo2016Nov ; barredo_synthetic_2018 ; Gross2017 ; rui_subradiant_2020 .
Such ensembles of qubits being ordered advanced optical manipulation methods
in the free space rui_subradiant_2020 ; barredo_synthetic_2018  or in the
vicinity of nanophotonic structures Corzo2016 ; Prasad2020 ;
liedl_collective_2023  have already demonstrated potential in storing,
processing and transmitting quantum information Corzo2019 . However,
implementation of these technological achievements in practical devices
requires solving a number of fundamental problems such as storing and
protecting quantum information from decoherence.

One of the possible ways to overcome these problems is to suppress spontaneous
emission of quantum emitters, which can be achieved in the subradiant
collective states in the ensembles of entangled qubits. Such states have been
largely explored recently in one-dimensional atomic chains coupled to a
waveguide Sheremet2023Mar ; Kornovan2019Dec ; pennetta_observation_2022 , or
atomic clouds in free space guerin_subradiance_2016 . However, the advantages
of two-dimensional arrays such as methods of precise optical positioning with
high filling factor Barredo2016Nov , convenience of transmittance and
reflection spectrum measurements Facchinetti2018Feb ; Rui2020Jul , possible
convenient application of Rydberg blockade Saffman2010Aug , and many others,
resulted in outstanding progress in optical manipulation over qubit states
srakaew_subwavelength_2023 . Moreover, there have been a number of novel
theoretical approaches to photons manipulation with atomic lattices
ballantine_parity-time_2021 ; ballantine_cooperative_2021 ;
nefedkin_nonreciprocal_2023 ; pedersen_quantum_2023 .

The studies of the long-living states in atomic lattices were in the center of
several works already showing that in the square two-dimensional (2D) lattices
Facchinetti2016Dec  the eigenstates with out-of-plane dipole moment
orientation and in-phase coherence (near the $\Gamma$-point of the Brillouin
zone (BZ), see Fig. 1) was shown to have a lifetime proportional to $N^{2}$,
where $N$ is the number of atoms along one of the directions. Such states fall
in the type of bound states in continuum extensively studied recently in
photonics Hsu2016 ; koshelev_bound_2021 , and have been further a subject for
investigation in the number of papers Ballantine2020Apr ; de_paz_bound_2023 ;
gladyshev_bound_2022  where their appearance in lattices of different
symmetries was also discussed Sadrieva2019 .

Figure 1: (a) A two-dimensional finite square regular lattice of $N\times N$
with $\sigma_{z}$ transition and separated by distance $a$. The distribution
of the highly non-radiative mode of $B_{2}$ symmetry is shown with arrows. (b)
The polaritonic band diagram of the infinite dipolar lattice. The main
mechanisms responsible for the formation of non-radiative states are shown.

In this work, we consider planar finite arrays of two-level atoms arranged in
different regular structures with the focus on the square atomic arrays
schematically shown in Fig. 1(a). We identify the main factors affecting the
radiative losses suppression in such structures, including interference of the
out-of-phase oscillating neighbor dipoles in the far zone, external coupling
of the states associated with the symmetry of the structure, and “accidental”
external coupling that emerges due to quasi-flat polaritonic band dispersion
of the corresponding infinite lattice, see Fig. 1(b). We establish the
conditions under which these mechanisms are present and demonstrate the
interplay between them on the example of the square array. We show that the
radiative decay rate of the most subradiant states in square arrays decrease
as fast as $\propto N_{tot}^{-3}$ or $\propto N_{tot}^{-5}$ ($N_{tot}=N^{2}$
is the total number of atoms in the array), which leads to very long lifetimes
of the eigenstates even in compact arrays. Moreover, there exist optimal
geometrical parameters that provide resonant increase of the lifetime by a few
orders. We also demonstrate the possibility of selective excitation of the
most subradiant states by vector Bessel beams with a suitably chosen orbital
angular momentum.

##  II General formalism

Let us consider an array of identical two-level atoms in free space with a
transition frequency $\omega_{0}$ (transition wavelength $\lambda_{0}$) and
spontaneous decay rate
$\Gamma_{0}={\omega_{0}^{3}|d|^{2}}/({3\pi\hbar\varepsilon_{0}c^{3}})$, where
$\varepsilon_{0}$ is vacuum permeability and $d$ is the amplitude of the
transition dipole moment. We assume that the atoms are separated by the
distance more than $\sim 0.1\lambda_{0}$, $\lambda_{0}=2\pi c/\omega_{0}$,
thus neglecting the Casimir interaction between atoms in the ground state due
to its very fast spatial decay $\sim 1/r^{6}$ Buhmann2015Mar . Consequently,
the coupling between the atoms is determined by the dipole-dipole interaction,
governed by the electromagnetic free space Green’s tensor
$\mathbf{G}(\mathbf{R},\omega)$ Novotny2012Sep .

To describe such interaction, we introduce dipole moment operator of the
$i$-th atom as
$\hat{\mathbf{d}}_{i}=\mathbf{d}\hat{\sigma}_{i}+\mathbf{d}^{*}\hat{\sigma}_{i}^{{\dagger}}$,
where $\hat{\sigma}_{i}=|g_{i}\rangle\langle e_{i}|$,
$\hat{\sigma}_{i}^{{\dagger}}=|e_{i}\rangle\langle g_{i}|$ are the atomic
lowering and rising operators describing the transition from the excited state
$|e_{i}\rangle$ to the ground $|g_{i}\rangle$ state of the $i$-th atom and
vice versa, respectively. The interaction Hamiltonian of the system
$\hat{\mathcal{H}}_{\text{int}}$, describing both re-scaterring processes
governed by dipole-dipole interaction between $N_{tot}$ atoms and spontaneous
emission into the free space, can be obtained as the total interaction energy
between dipole moments and electric field:
$\hat{\mathcal{H}}_{\text{eff}}=-\sum_{i=1}^{N_{tot}}\hat{\mathbf{d}}_{i}\hat{\mathbf{E}}(\mathbf{r}_{i})$,
where $\hat{\mathbf{E}}(\mathbf{r}_{i})$ is the electric field operator at the
position $\mathbf{r}_{i}$. Exploiting classical principle of electric field
superposition and the connection between the electric field and Green’s tensor
Belov2005Aug ; Markel2007Feb ; AsenjoPRA2017 , one can directly obtain the
effective Hamiltonian reitz_cooperative_2021 ; ruostekoski_cooperative_2023 :

$\hat{\mathcal{H}}_{\text{eff}}=\hbar\Bigl{(}\omega_{0}-\frac{i}{2}\Gamma_{0}\Bigr{)}\sum\limits_{i=1}^{N_{tot}}\hat{\sigma}_{i}^{{\dagger}}\hat{\sigma}_{i}-\\\
-\sum\limits_{i\neq
j}\frac{3\pi\hbar\Gamma_{0}}{k_{0}}\mathbf{e_{d}^{*}}\cdot\mathbf{G}(\mathbf{R}_{ij})\cdot\mathbf{e_{d}}\
\hat{\sigma}_{i}^{{\dagger}}\hat{\sigma}_{j},$ (1)

where $\mathbf{e_{d}}=\mathbf{d}/d$ is the unit vector of the atomic
transition dipole moment. In Eq.(1), the first term describes the coupling of
atoms to free space, while the second one describes the atom-atom coupling.
The above equation is obtained within the Markovian approximation by
neglecting the frequency dependence of the Green’s tensor,
$\mathbf{G}(\mathbf{R},\omega)\approx\mathbf{G}(\mathbf{R},\omega_{0})$, which
is justified by a very narrow resonance of a single atom
$\Gamma_{0}\ll\omega_{0}$.

The eigenvalues and eigenstates of the Hamiltonian (1) define the complex
eigenfrequencies $\omega_{j}-i\Gamma_{j}/2$ and eigenfunctions $\Psi^{(j)}$ of
the collective dipolar excitations, respectively. Real and imaginary parts of
the the eigenfrequency correspond to the frequency and radiative decay rate of
the $j$-th collective eigenstate. For convenience throughout the paper we
define the frequency with respect to that of a single atom,
$\Delta\omega_{j}=\omega_{j}-\omega_{0}$.

Normally, in the absence of the external magnetic field, atoms are spherically
symmetric, and the simplest model of such an atom is the one with triply
degenerate excited state ($S\leftrightarrow P$ transition). In the presence of
the $z$-oriented magnetic field [normal to the $x-y$ plane of the array, see
Fig. 1 (a)] the states of an atom can be spectrally separated into three
states depending on the polarization of the transition dipole moment,
$\sigma_{-}$, $\sigma_{+}$, or $\sigma_{z}$. Consequently, the collective
eigenstates of the planar arrays of such atoms can be considered separately
for different polarizations of the atomic transitions.

##  III Square atomic array

It is known that the radiative losses are mostly suppressed in symmetric
structures, which possess fewer radiation channels than non-symmetric ones. In
this section, we focus on one of the most simple and high-symmetry case of 2D
atomic array $N\times N$ arranged is a square lattice with a period $a$. Such
illustrative example allows us to demonstrate the competition between
different mechanisms of the radiative loss suppression.

###  III.1 Classification of the collective eigenstates

Figure 2: The symmetry classification of the basis functions
$\psi^{(m_{x},m_{y})}$ with respect to the irreducible representations and
parity numbers $m_{x}$ and $m_{y}$ [see Eq. (2)]. Color shows the real part of
the wavefunction amplitude at each lattice site.

The eigenstates of the finite lattices can be classified and characterized
based on their distribution in the quasi-reciprocal space Kornovan2016 ;
Asenjo-Garcia2017Aug ; Sheremet2023Mar  which can be done by expanding the
dipole moments distribution of the $j$-th eigenstate $\Psi^{(j)}$ in the
standing (Bloch) waves basis:

$\Psi^{(j)}(n_{x},n_{y})=\sum_{m_{x},m_{y}}c^{(j)}_{m_{x},m_{y}}\psi^{(m_{x},m_{y})}=\\\
=\dfrac{2}{N+1}\sum_{m_{x},m_{y}}c^{(j)}_{m_{x},m_{y}}\sin(q_{0}m_{x}n_{x})\sin(q_{0}m_{y}n_{y}),$
(2)

where $c^{(j)}_{m_{x},m_{y}}$ are the expansion coefficients,
$q_{0}=\pi/(N+1)$ is the discretization quasi-wavevector, $m_{x,y}=1..N$ are
the numbers of the basis state in the quasi-reciprocal space, and
$n_{x,y}=1..N$ are the atoms positions in the array in the real space. Such
basis is straight-forwardly constructed as a tensor product of the two sets of
standing waves in $x$ and $y$ directions. Two indices, $m_{x}$ and $m_{y}$,
that specify the basis state $\psi^{(m_{x},m_{y})}$, show that the standing
wave distribution of the dipole moments changes the sign $m_{x}-1$ and
$m_{y}-1$ times in $x$ and $y$ directions, respectively. For example, the
basis state that is closest to the corner of the BZ [$M$-point in Fig. 1(b)]
in the reciprocal space is characterized by the indices $m_{x}=m_{y}=N$.

Alternative natural classification of the eigenstates is based on their
symmetry properties. The square array belongs to the $C_{4v}$ point symmetry
group and thus the eigenstates $\Psi^{(j)}$ should transform according to one
of the irreducible representations (irreps) $A_{1}$, $A_{2}$, $B_{1}$,
$B_{2}$, $E$. However, the basis given by Eq. (2) does not account for the
specific symmetry of the square array. Therefore, not all of the basis states
$\psi^{(m_{x},m_{y})}$ fall into the irreps of the $C_{4v}$ point symmetry
group and, thus, should be symmetrised. Then, the following correspondence
between the basis states and irreps, summarized in Fig. 2, can be established.
The basis states with odd value of $(m_{x}+m_{y})$ transform according to $E$
irreducible representation and have “azimuthal numbers” $\pm 1$ (or their
linear combination), i.e. phases of the neighbor corner dipoles differ by
$\pi/2$. The basis states with $m_{x}=m_{y}$ transform according to $A_{1}$
($B_{2}$) irreducible representation if $m_{x,y}$ is odd (even), and have
“azimuthal numbers” $0$ ($2$), i.e. phases of the neighbor corner dipoles are
the same (differ by $\pi$). If $m_{x}\neq m_{y}$ and $(m_{x}+m_{y})$ is even
the mode $\psi^{(m_{1},m_{2})}$ hybridises with the mode
$\psi^{(m_{2},m_{1})}$ and they form two linear combinations:

$\psi^{(m_{1},m_{2})^{\pm}}=\dfrac{1}{\sqrt{2}}\left(\psi^{(m_{1},m_{2})}\pm\psi^{(m_{2},m_{1})}\right).$
(3)

Antisymmetric modes transform according to $A_{2}$ or $B_{1}$, while symmetric
modes transform according to $A_{1}$ or $B_{2}$ irrep.

###  III.2 Mechanism of the radiation suppression

If the eigenstate $\Psi^{(j)}$ of the atomic array is mainly characterized by
a single dominant contribution of a basis state $\psi^{(m_{x},m_{y})}$, its
radiative losses $\Gamma_{j}$ are qualitatively determined by its location in
the BZ, which defines the efficiency of the destructive interference between
the neighbouring dipoles in the array. The losses are typically weaker for the
states that are closer to the $M$-point in the reciprocal space, i.e. when
$m_{x},m_{y}\rightarrow N$, due to the larger mismatch between the wavevector
of such state and the free-space wavevector.

The second factor affecting the radiative losses is associated with the
symmetry of the structure and the eigenstates. For instance, Bloch waves
propagating in in an infinite square lattice in directions near the $M$-point
and symmetric with respect to the $\Gamma M$-direction are degenerate owing to
the symmetry of the lattice, see Fig. 1. However, in a finite array the
degeneracy is lifted due to the external coupling of these waves according to
Friedrich-Wintgen mechanism Friedrich1985a  leading to the formation of
symmetric and antisymmetric states. The latter ones transform according to
$A_{2}$ or $B_{1}$ irreps and are antisymmetric with respect to the diagonal
vertical planes of symmetry $n_{x}=\pm n_{y}$, i.e. the excitation of the
corner dipole moments of such states is necessarily zero. Due the suppressed
scattering from the sharp edges of the structure the radiative losses of these
states can substantially decrease Asenjo-Garcia2017Aug ; Poddubny2020
WiersigPRL2006 . On the opposite, symmetric modes that transform according to
$A_{1}$ or $B_{2}$ with non-zero corner dipole moments are characterized with
the increased losses.

Figure 3: (a) Normalized radiative decay rate of two most subradiant
eigenstates corresponding to $B_{2}$ (green curves) and $A_{2}$ (violet
curves) symmetry in $12\times 12$ square array as a function of the normalized
period $\tilde{a}$. (b,c) Real part of the wavefunction plotted in real space
for the states shown in (a) with green and violet solid curves, respectively.
(d-g) Expansion of the eigenstates, Eq. (2), shown in (a) with (d,e) black and
(f,g) blue points; size of the circles indicate the amplitude of the
corresponding harmonics; wavy lines illustrate the coupling between different
harmonics.

Additional mechanism of the losses suppression may appear due to the external
coupling of the eigenstates unrelated to the symmetry of the structure, but
rather to the accidental degeneracy of the Bloch waves in the corresponding
infinite lattice. As we show further, such regime can be realised for small
enough periods, when the lattice dispersion becomes quasi-flat. In this case
two (or typically several) basis states interact via radiation continuum
resulting in additional suppression or increase of their radiative losses. The
overall radiative losses of an eigenstate are then determined as a result of
the competition between these three mechanisms of loss suppression. In order
to determine which of the eigenstates exhibits minimal losses for a given size
and period of the structure, we perform straight-forward numerical
calculations and the results are presented in the next subsection.

###  III.3 Numerical calculations of the finite array eigenstates

In Fig. 3(a) we show numerically calculated emission rates
$\Gamma/\Gamma_{0}$, normalized by the emission rate of a single atom, for two
collective states of the $12\times 12$ square array as a function of the
normalized period $\tilde{a}=a/\lambda_{0}$. These states are characterized by
the minimal emission rates among the others. Solid lines correspond to the
$\sigma_{z}$ dipoles, while the dashed lines to the $\sigma_{+}$ or
$\sigma_{-}$ ones. One can observe that for periods $\tilde{a}\gtrsim 0.32$
the radiative losses monotonously increase with the decrease of the period.
Standing wave decomposition Eq. (2) performed in this range of periods reveals
that all states are characterized by the dominant contribution either of a
single basis state or two basis states in the symmetrized form, Eq. (3). The
numerical analysis revealed that the most subradiant states (in the array with
even $N$) transform according to either $B_{2}$ or to $A_{2}$ irrep. The
distribution of the wavefunction of the $\sigma_{z}$ polarized states in real
space is shown in Figs. 3(b,c), while its representation in the reciprocal
space, i.e. the decomposition Eq. (2), for the period $\tilde{a}=0.4$ is shown
in Figs. 3(d,e). The distribution for the $\sigma_{\pm}$ look similar. As
expected, the $B_{2}$ state is the closest one to the $M$-point and it is well
approximated by the $\psi^{(N,N)}$ basis function. The dipole moments of this
state in the neighbour atoms are almost out-of-phase, which leads to their
efficient destructive interference.

The radiative losses of such state, however, are minimal only for the periods
$\tilde{a}\gtrsim 0.36$. For smaller periods, the $A_{2}$ state, which is the
second closest to the $M$-point, becomes more subradiant. This state has the
dominant components $m_{x}=N,m_{y}=N-2$ and $m_{x}=N-2,m_{y}=N$ in the
antisymmetric combination $\psi^{(N,N-2)^{-}}$, Eq. (3) (see also Fig. 1). In
accordance to the symmetry properties of the irrep $A_{2}$, the dipole moments
of such state are equal to zero on diagonals. This leads to small probability
that the excitation resides on the corners of the square, see Fig. 3(c),
leading to suppressed scattering from sharp edges and consequently weak
radiative losses of such state in the most range of the periods.

Such antisymmetric subradiant states can also be viewed as the result of the
external coupling between two eigenstates of the rectangular array that
appears under stretching/squeezing the lattice to the square one YangIOP2016 ,
see also Supplementary Materials, Sec. A for the details. Consequently, the
mechanism of the symmetry-induced external coupling is very sensitive to the
deformation of the lattice: in a rectangular lattice the efficiency of the
coupling drops very fast, resulting in the substantially increased losses of
the $\psi^{(N,N-2)^{-}}$ state even for very small detunings from the square
geometry of the order of 1%. On the other hand, the $\psi^{(N,N)}$ state, not
being associated with the specific symmetry of the structure, survives under
strong deformation of the lattice up to tens of percents.

For even smaller periods $\tilde{a}\approx 0.3$ the decay rates of both
$A_{2}$ and $B_{2}$ states for $\sigma_{z}$ polarization exhibit a minimum for
a given number of atoms, see Fig. 3(a). At these minima the wavefunctions of
the states in the real space remain similar to those in Figs. 3(b,c). However,
the expansion in the reciprocal space, Eq. 2, shown in Figs. 3(f,g), reveals
that the amplitudes of other harmonics are substantially increased. The
interaction between these harmonics leads to additional suppression of
radiative losses near the period $\tilde{a}\approx 0.3$. At the same time, the
calculations for the $\sigma_{\pm}$ dipole transitions, dashed curves in Fig.
3(a), show that radiative losses are monotonous function of the period in this
case and no resonant suppression of the losses via external coupling mechanism
is observed.

###  III.4 Infinite atomic lattices

Appearance of the external coupling near $\tilde{a}\sim 0.3$ for $\sigma_{z}$
polarized dipoles and its absence for the $\sigma_{\pm}$ ones can be linked to
the behaviour of the infinite lattice dispersion of the polaritonic states
near the $M$-point. To illustrate this, we performed calculations of the band
diagrams of the square lattice for both polarizations of the dipole
transitions and for different periods, see details of the calculations in the
Supplementary Materials, Sec. B. First, we consider the most interesting case
of the $z$-oriented dipoles. In Fig. 4(a), we plot the dispersion along the
main direction in the first BZ. Fig. 4(b) complements the panel Fig. 4(a) with
the 2D colorplot in the whole BZ. Above the light cone, shown with the black
dashed lines in Fig. 4(a) and with white dashed curve in Fig. 4(b), the Bloch
states are leaky and the eigenstates of the corresponding finite arrays are
characterized with large radiative losses. On the other hand, the modes below
the light line are guided (the eigenfrequency is purely real) and the
eigenstates of the corresponding finite arrays are subradiant.

As it was shown in the previous subsection, the eigenstates of the finite
arrays that are closest to the $M$-point, are the most subradiant due to the
largest mismatch with the wavevectors of the free space. One can observe in
Fig. 4(c) that the dispersion near the $M$-point turns out to be qualitatively
different for different periods. For the large periods $\tilde{a}\gtrsim 0.3$
the interaction between the dipoles is rather weak, so that the dispersion is
parabolic. For the lattice period near $\tilde{a}=0.294$ the dispersion
exhibits an inflection point and becomes quartic, as shown with the red curve
in Fig. 4(c). For even smaller periods $\tilde{a}\lesssim 0.294$, the
dispersion becomes nonmonotonous near the $M$-point along the $\Gamma M$
direction [green curve in Fig. 4(c)], which explains the external coupling of
the several states as now their energies match with each other, analogously to
the one-dimensional case Kornovan2021Dec . In contrast, dispersion in the case
of $\sigma_{+,-}$ polarized dipole transition, shown in Fig. 4(d) with dashed
curves, remains monotonic near the $M$-point for any subwavelength period.
This explains the absence of the external coupling in finite arrays for this
polarization.

Figure 4: (a,b) Dispersion of an infinite square lattice of $\sigma_{z}$
polarized atoms for the period $\tilde{a}=0.35$ (a) along high symmetry path
and (b) in the whole BZ. (c,d) Dispersion near the M point for periods
$\tilde{a}_{1}=0.35$ (blue curves), $\tilde{a}_{2}=0.294$ (red curves) and
$\tilde{a}_{3}=0.28$ (green curves) for (c) $\sigma_{z}$ polarized atoms and
(d) $\sigma_{\pm}$ polarized atoms.

##  IV Discussion

###  IV.1 Variation of the radiative losses with the size of arrays with
different geometries

The mechanisms that determine the formation of the subradiant states in atomic
arrays shown on the example of the square geometry have rather general nature
and are present in various configuration of atomic ensembles. Calculations
performed for a few other main geometries of the regular atomic arrays such as
rotated square, triangle, hexagon (see Appendix, Sec C. for the details) have
revealed that the square geometry is preferable from the point of view of the
state lifetime. To illustrate this, in Fig. 5(a) we show decay rate of the
four characteristic eigenstates (corresponding to different geometries of the
array) as a function of the total number of the atoms in the array $N_{tot}$,
calculated for the non-optimal period $\tilde{a}=0.4$ when there is no
accidental external coupling between the eigenstates caused by the quasi-flat
dispersion. The polarization of the dipoles was chosen to be along the $z$
direction. The $N_{tot}^{-5}\equiv N^{-10}$ scaling law of the $A_{2}/B_{1}$
state in a square array, shown with purple pentagrams in Fig. 5(a), turns out
to be unique. For instance, diagonal square or triangle arrays exhibit only
$\propto N_{tot}^{-1.5}$ scaling laws providing much shorter lifetimes than
the square geometry even for rather small $N_{tot}$.

Another type of the eigenstate of the square array, $A_{1}/B_{2}$ states with
the dominant $\psi^{(N,N)}$ contribution, exhibit only $\propto
N_{tot}^{-3}\equiv N^{-6}$ dependence of the decay rate. Such difference in
the scaling laws for this particular geometry can be qualitatively explained
in a simple manner, as elaborated in the Supplementary Materials, Sec. D.
Despite the slower decrease of the decay rate for $A_{1}/B_{2}$ states, such
asymptotic dependence is valid only for large number of atoms, while it does
not necessarily reflect the difference between the radiative losses for a
particular $N_{tot}$. One can observe in Fig. 5(a), that the antisymmetric
state becomes the most subradiant only for $N_{tot}\gtrsim 300$ for a given
period $\tilde{a}=0.4$.

Figure 5: (a,b) Radiative losses of subradiant states as a function of the
total number of atoms in the array $N_{tot}$ for different geometries for (a)
fixed distance between the neighbor atoms $\tilde{a}=0.4$, (b) optimal period
for each $N_{tot}$, shown in panel (c). Dashed lines indicate corresponding
polynomial functions.

In Fig. 5(b) we plot the calculated decay rates for the optimal periods
$\tilde{a}_{opt}$, for which the losses are minimal (i.e. period is optimized
for each value of $N_{tot}$), as indicated in Fig. 5(c). Here, the presence of
external coupling leads to a change of the scaling laws in some geometries
while remaining the same $\propto N_{tot}^{-3}$ and $\propto N_{tot}^{-5}$ in
the square geometry. This contrasts with the one-dimensional case, where
interaction between two standing waves leads to the decrease of the asymptotic
behaviour of the subradiant state decay rate from $\propto N^{-3}$ to $\propto
N^{-7}$ Kornovan2019Dec ; Kornovan2021Dec . Although, e.g. in a triangle
geometry decrease of the decay rate also becomes close to $\propto
N_{tot}^{-5}$, the absolute values of the decay rate still remain minimal for
the $A_{2}/B_{1}$ states of the square geometry.

###  IV.2 Excitation of the subradiant states

In order to demonstrate the possibility of excitation of the strongly
subradiant states appearing for the $\sigma_{z}$ polarization, we consider the
scattering of vector Bessel beams Richards1959Dec ; Kotlyar2019Oct  by a
square atomic array. The use of such specific excitation field has two
underlying reasons. First, we are considering the states characterized by
large wavevectors. The use of Bessel beams, which can have arbitrary angular
momentum, allows to partially match the azimuthal component of the wavevector
of the incident field to that of the eigenstate. Second, Bessel beams have
non-zero longitudinal electric field, which allows for efficient excitation of
the $\sigma_{z}$ polarized states. The considered highly inhomogeneous
incident field may excite various modes of the lattice with different
efficiency, therefore we have optimized orbital number of the beam $m$ to
maximize the overlap between external electric field and the desired
collective state, see the Supplementary Materials, Sec. E for details of the
scattering calculations).

In Fig. 6, we show numerically calculated scattering cross-section spectra of
$6\times 6$ square atomic array irradiated by normally incident vector Bessel
beam with a spin $s=1$ due to circular polarization, and orbital numbers $l=7$
or $l=9$ Kotlyar2019Oct . Total cross-section consists of various peaks
matching eigen modes of the lattice characterized with different dominant
basis states. The parameters of the beams were chosen to achieve the maximum
coupling with the eigen states with dominant contribution of basis states
$\psi^{(N,N-2)^{-}}$ and $\psi^{(N,N)}$ respectively. Lattice periods
correspond to optimal periods in terms of radiative losses of these states,
$\tilde{a}=0.268$ and $\tilde{a}=0.281$ in the case of $6\times 6$ array.
Cross-section spectra are normalized with a maximum value of a single atom
cross-section $\sigma_{0}$.

As one can see in Fig. 6 insets, the optimal incident field profile barely
overlaps the atomic array from the outside, setting $\pi$ phase shift between
neighboring atoms at the edges with almost no influence at internal atoms.
Therefore, a total angular momentum of the beam $l+s$ matches to either a half
of total number of edge atoms $2(N-1)$, see Fig. 6(a), or the same number but
excluding corner atoms $2(N-2)$, see Fig. 6(b). In the first case, the beam is
capable to excite the $B_{2}$ state with a dominant $\psi^{(N,N)}$
contribution. In the second case, where the field at the corners is absent, it
is optimal to excite $A_{2}$ state with a dominant $\psi^{(N,N-2)^{-}}$
contribution. Thereby a fine tuning of the incident beam frequency and angular
momentum allows to selectively excite the desired long-living eigen mode of
the array.

Figure 6: Normalized scattering cross section spectra of 6$\times$6 square
atomic lattice irradiated by Bessel beam with (a) azimuthal number $m=9$ and
period $\tilde{a}=0.281$ and (b) $m=7$ and $\tilde{a}=0.268$. Solid black
curves correspond to the total scattering cross section spectra, while colored
curves correspond to different eigenstates with dominant contributions of
basis states indicated in the plots. Insets show distributions of $z$
component of the incident electric field.

###  IV.3 Relation to high quality modes in nanophotonic cavities

Importantly, the generality of the reported results is also underlined by the
tight connection to the research area of engineering of planar dielectric
nanophotonic cavities, such as homogeneous dielectric cavities GuoIEEE2003 ;
LiIEEE2006 ; YangIEEE2007 ; YangPRA2007 ; BittnerPRE2013 ; CaoRevModPhys2015 ;
YangIOP2016  and periodic or quasi-periodic nanostructured cavities TsiaOE2004
; NohPRA2010 .

The resonators made of homogeneous dielectric are expected to have the highest
$Q$-factors for spherical or cylindrical geometries, which exhibit the highest
symmetry. However, the utilized fabrication methods or proposed technological
applications sometimes demand the use of the cavities of the other shapes,
e.g. cubes, squares, hexagons etc., which were also shown to support rather
high $Q$-factors and have potential in development of micro and nanolasers
YangIOP2016 ; ZhengZhengACSNano2018 ; TiguntsevaACSNano2020 . The dispersion
of homogeneous dielectric is, however, quasi-linear and the external coupling
due to quasi-flat band in this case is not possible. Unlike the homogeneous
cavitites, the dispersion in the nanostructured ones, i.e. photonic-crystal or
nanoparticle cavities, is more tunable. Consequently, quasi-flat and
nonmonotonous bands in such structures result in the quality factor boost in
the finite arrays NohPRA2010 . By proper translating all mechanisms of the
subradiant states formation reported in this work for the dipolar arrays onto
the nanophotonic platform could lead to the development of the novel designs
of the compact optical micro and nanocavities.

##  V Summary

To summarize, we have studied main factors that define the radiative loss
suppression in finite-size planar atomic arrays. We have shown that the square
arrays support eigenstates with minimal losses among different regular
arrangements of the atoms, that scale with the number of atoms as $\propto
N_{tot}^{-5}$. The two distinct types of external coupling associated with the
symmetry of the array and with the quasi-flat dispersion of the corresponding
infinite array play the dominant role in the strong boost of the radiative
lifetime of such states. We believe, our findings might provide useful insight
in designing of the cold atoms based qubit arrays for storing and coherent
manipulation of quantum information.

##  VI Acknowledgement

The work was supported by the Federal Academic Leadership Program Priority
2030. The numerical computations were supported by the Russian Science
Foundations grant No. 21-72-10107. We thank Yuri Kivshar and Ivan Iorsh for
significant discussions.

## References

  * (1) D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, “An atom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays,” Science, vol. 354, no. 6315, pp. 1021–1023, 2016.
  * (2) D. Barredo, V. Lienhard, S. De Léséleuc, T. Lahaye, and A. Browaeys, “Synthetic three-dimensional atomic structures assembled atom by atom,” Nature, vol. 561, pp. 79–82, Sept. 2018.
  * (3) C. Gross and I. Bloch, “Quantum simulations with ultracold atoms in optical lattices,” Science, vol. 357, pp. 995–1001, Sept. 2017.
  * (4) J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, “A subradiant optical mirror formed by a single structured atomic layer,” Nature, vol. 583, no. 7816, pp. 369–374, 2020.  arXiv: 2001.00795 Publisher: Springer US.
  * (5) N. V. Corzo, B. Gouraud, A. Chandra, A. Goban, A. S. Sheremet, D. V. Kupriyanov, and J. Laurat, “Large Bragg Reflection from One-Dimensional Chains of Trapped Atoms Near a Nanoscale Waveguide,” Physical Review Letters, vol. 117, pp. 133603–133603, 2016.  arXiv: 1604.03129.
  * (6) A. S. Prasad, J. Hinney, S. Mahmoodian, K. Hammerer, S. Rind, P. Schneeweiss, A. S. Sørensen, J. Volz, and A. Rauschenbeutel, “Correlating photons using the collective nonlinear response of atoms weakly coupled to an optical mode,” Nature Photonics, vol. 1911.09701, no. 12, pp. 719–722, 2019.  arXiv: 1911.09701 Publisher: Springer US.
  * (7) C. Liedl, S. Pucher, F. Tebbenjohanns, P. Schneeweiss, and A. Rauschenbeutel, “Collective Radiation of a Cascaded Quantum System: From Timed Dicke States to Inverted Ensembles,” Physical Review Letters, vol. 130, p. 163602, Apr. 2023.
  * (8) N. V. Corzo, J. Raskop, A. Chandra, A. S. Sheremet, B. Gouraud, and J. Laurat, “Waveguide-coupled single collective excitation of atomic arrays,” Nature, vol. 566, no. 7744, pp. 359–362, 2019.
  * (9) A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, “Waveguide quantum electrodynamics: Collective radiance and photon-photon correlations,” Rev. Mod. Phys., vol. 95, no. 1, p. 015002, 2023.
  * (10) D. F. Kornovan, N. V. Corzo, J. Laurat, and A. S. Sheremet, “Extremely subradiant states in a periodic one-dimensional atomic array,” Phys. Rev. A, vol. 100, no. 6, p. 063832, 2019.
  * (11) R. Pennetta, D. Lechner, M. Blaha, A. Rauschenbeutel, P. Schneeweiss, and J. Volz, “Observation of Coherent Coupling between Super- and Subradiant States of an Ensemble of Cold Atoms Collectively Coupled to a Single Propagating Optical Mode,” Physical Review Letters, vol. 128, p. 203601, May 2022.
  * (12) W. Guerin, M. O. Araújo, and R. Kaiser, “Subradiance in a Large Cloud of Cold Atoms,” Physical Review Letters, vol. 116, p. 083601, Feb. 2016.
  * (13) G. Facchinetti and J. Ruostekoski, “Interaction of light with planar lattices of atoms: Reflection, transmission, and cooperative magnetometry,” Phys. Rev. A, vol. 97, no. 2, p. 023833, 2018.
  * (14) J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, “A subradiant optical mirror formed by a single structured atomic layer,” Nature, vol. 583, pp. 369–374, 2020.
  * (15) M. Saffman, T. G. Walker, and K. Mølmer, “Quantum information with Rydberg atoms,” Rev. Mod. Phys., vol. 82, no. 3, pp. 2313–2363, 2010.
  * (16) K. Srakaew, P. Weckesser, S. Hollerith, D. Wei, D. Adler, I. Bloch, and J. Zeiher, “A subwavelength atomic array switched by a single Rydberg atom,” Nature Physics, vol. 19, pp. 714–719, May 2023.
  * (17) K. E. Ballantine and J. Ruostekoski, “Parity-time symmetry and coherent perfect absorption in a cooperative atom response,” Nanophotonics, vol. 10, pp. 1357–1366, Feb. 2021.
  * (18) K. E. Ballantine and J. Ruostekoski, “Cooperative optical wavefront engineering with atomic arrays,” Nanophotonics, vol. 10, pp. 1901–1909, May 2021.
  * (19) N. Nefedkin, M. Cotrufo, and A. Alù, “Nonreciprocal total cross section of quantum metasurfaces,” Nanophotonics, vol. 12, pp. 589–606, Feb. 2023.
  * (20) S. P. Pedersen, L. Zhang, and T. Pohl, “Quantum nonlinear metasurfaces from dual arrays of ultracold atoms,” Physical Review Research, vol. 5, p. L012047, Mar. 2023.
  * (21) G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, “Storing Light with Subradiant Correlations in Arrays of Atoms,” Phys. Rev. Lett., vol. 117, no. 24, p. 243601, 2016.
  * (22) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, “Bound states in the continuum,” Nature Reviews Materials, vol. 1, no. 9, pp. 16048–16048, 2016.
  * (23) K. Koshelev, Z. Sadrieva, A. Shcherbakov, Y. Kivshar, and A. Bogdanov, “Bound states in the continuum in photonic structures,” Physics-Uspekhi, Dec. 2021.  arXiv:2207.01441 [physics, physics:quant-ph].
  * (24) K. E. Ballantine and J. Ruostekoski, “Subradiance-protected excitation spreading in the generation of collimated photon emission from an atomic array,” Phys. Rev. Res., vol. 2, no. 2, p. 023086, 2020.
  * (25) M. B. De Paz and P. A. Huidobro, “Bound states in the continuum in subwavelength emitter arrays,” Physical Review Research, vol. 5, p. 033108, Aug. 2023.
  * (26) S. Gladyshev, A. Shalev, K. Frizyuk, K. Ladutenko, and A. Bogdanov, “Bound states in the continuum in multipolar lattices,” Physical Review B, vol. 105, p. L241301, June 2022.
  * (27) Z. Sadrieva, K. Frizyuk, M. Petrov, Y. Kivshar, and A. Bogdanov, “Multipolar origin of bound states in the continuum,” Physical Review B, vol. 100, p. 115303, Sept. 2019.  arXiv: 1903.00309 Publisher: American Physical Society.
  * (28) S. Y. Buhmann, Dispersion Forces I: Macroscopic Quantum Electrodynamics and Ground-State Casimir, Casimir–Polder and van der Waals Forces (Springer Tracts in Modern Physics, 247).  Berlin, Germany: Springer, 2015.
  * (29) L. Novotny and B. Hecht, Principles of Nano-Optics.  Cambridge, England, UK: Cambridge University Press, 2012.
  * (30) P. A. Belov and C. R. Simovski, “Homogenization of electromagnetic crystals formed by uniaxial resonant scatterers,” Phys. Rev. E, vol. 72, no. 2, p. 026615, 2005.
  * (31) V. A. Markel and A. K. Sarychev, “Propagation of surface plasmons in ordered and disordered chains of metal nanospheres,” Phys. Rev. B, vol. 75, no. 8, p. 085426, 2007.
  * (32) A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, “Atom-light interactions in quasi-one-dimensional nanostructures: A green’s-function perspective,” Phys. Rev. A, vol. 95, p. 033818, Mar 2017.
  * (33) M. Reitz, C. Sommer, and C. Genes, “Cooperative quantum phenomena in light-matter platforms,” PRX Quantum, vol. 3, p. 010201, 2021.  arXiv: 2107.02674 Publisher: American Physical Society.
  * (34) J. Ruostekoski, “Cooperative quantum-optical planar arrays of atoms,” Physical Review A, vol. 108, p. 030101, Sept. 2023.  arXiv:2309.08487 [cond-mat, physics:physics, physics:quant-ph].
  * (35) D. F. Kornovan, A. S. Sheremet, and M. I. Petrov, “Collective polaritonic modes in an array of two-level quantum emitters coupled to an optical nanofiber,” Physical Review B, vol. 94, pp. 245416–245416, Dec. 2016.  arXiv: 1608.03202v1.
  * (36) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, “Exponential Improvement in Photon Storage Fidelities Using Subradiance and ``Selective Radiance” in Atomic Arrays,” Phys. Rev. X, vol. 7, no. 3, p. 031024, 2017.
  * (37) H. Friedrich and D. Wintgen, “Interfering resonances and bound states in the continuum,” Physical Review A, vol. 32, no. 6, pp. 3231–3242, 1985.
  * (38) A. N. Poddubny, “Quasiflat band enabling subradiant two-photon bound states,” Physical Review A, vol. 101, no. 4, pp. 1–8, 2020.  Publisher: American Physical Society.
  * (39) J. Wiersig, “Formation of long-lived, scarlike modes near avoided resonance crossings in optical microcavities,” Phys. Rev. Lett., vol. 97, p. 253901, Dec 2006.
  * (40) Y.-D. Yang and Y.-Z. Huang, “Mode characteristics and directional emission for square microcavity lasers,” Journal of Physics D: Applied Physics, vol. 49, p. 253001, may 2016.
  * (41) D. F. Kornovan, R. S. Savelev, Y. Kivshar, and M. I. Petrov, “High-Q Localized States in Finite Arrays of Subwavelength Resonators,” ACS Photonics, vol. 8, no. 12, pp. 3627–3632, 2021.
  * (42) B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system,” Proc. R. Soc. London A - Math. Phys. Sci., vol. 253, no. 1274, pp. 358–379, 1959.
  * (43) V. V. Kotlyar, A. G. Nalimov, and S. S. Stafeev, “Exploiting the circular polarization of light to obtain a spiral energy flow at the subwavelength focus,” J. Opt. Soc. Am. B, JOSAB, vol. 36, no. 10, pp. 2850–2855, 2019.
  * (44) W.-H. Guo, Y.-Z. Huang, Q.-Y. Lu, and L.-J. Yu, “Modes in square resonators,” IEEE Journal of Quantum Electronics, vol. 39, no. 12, pp. 1563–1566, 2003.
  * (45) C. Li, L. Zhou, S. Zheng, and A. W. Poon, “Silicon polygonal microdisk resonators,” IEEE Journal of Selected Topics in Quantum Electronics, vol. 12, no. 6, pp. 1438–1449, 2006.
  * (46) Y.-D. Yang and Y.-Z. Huang, “Mode analysis and $q$-factor enhancement due to mode coupling in rectangular resonators,” IEEE Journal of Quantum Electronics, vol. 43, no. 6, pp. 497–502, 2007.
  * (47) Y.-D. Yang and Y.-Z. Huang, “Symmetry analysis and numerical simulation of mode characteristics for equilateral-polygonal optical microresonators,” Phys. Rev. A, vol. 76, p. 023822, Aug 2007.
  * (48) S. Bittner, E. Bogomolny, B. Dietz, M. Miski-Oglu, and A. Richter, “Experimental observation of localized modes in a dielectric square resonator,” Phys. Rev. E, vol. 88, p. 062906, Dec 2013.
  * (49) H. Cao and J. Wiersig, “Dielectric microcavities: Model systems for wave chaos and non-hermitian physics,” Rev. Mod. Phys., vol. 87, pp. 61–111, Jan 2015.
  * (50) K. K. Tsia and A. W. Poon, “Dispersion-guided resonances in two-dimensional photonic-crystal-embedded microcavities,” Opt. Express, vol. 12, pp. 5711–5722, Nov 2004.
  * (51) H. Noh, J.-K. Yang, I. Vitebskiy, A. Figotin, and H. Cao, “Giant resonances near the split band edges of two-dimensional photonic crystals,” Phys. Rev. A, vol. 82, p. 013801, Jul 2010.
  * (52) Z. Liu, J. Yang, J. Du, Z. Hu, T. Shi, Z. Zhang, Y. Liu, X. Tang, Y. Leng, and R. Li, “Robust subwavelength single-mode perovskite nanocuboid laser,” ACS Nano, vol. 12, no. 6, pp. 5923–5931, 2018.  PMID: 29746781.
  * (53) E. Tiguntseva, K. Koshelev, A. Furasova, P. Tonkaev, V. Mikhailovskii, E. V. Ushakova, D. G. Baranov, T. Shegai, A. A. Zakhidov, Y. Kivshar, and S. V. Makarov, “Room-temperature lasing from mie-resonant nonplasmonic nanoparticles,” ACS Nano, vol. 14, no. 7, pp. 8149–8156, 2020.  PMID: 32484650.