Strongly subradiant states in planar atomic arrays

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}$.


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 [1][2][3][4].Such ensembles of qubits being ordered advanced optical manipulation methods in the free space [2,4] or in the vicinity of nanophotonic structures [5][6][7] have already demonstrated potential in storing, processing and transmitting quantum information [8].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 [9][10][11], or atomic clouds in free space [12].However, the advantages of twodimensional arrays such as methods of precise optical positioning with high filling factor [1], convenience of transmittance and reflection spectrum measurements [13,14], possible convenient application of Rydberg blockade [15], and many others, resulted in outstanding progress in optical manipulation over qubit states [16].Moreover, there have been a number of novel theoretical approaches to photons manipulation with atomic lattices [17][18][19][20].
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 [21] the eigenstates with out-of-plane dipole moment orientation and in-phase coherence (near the Γ-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 [22,23], and have been further a subject for investigation in the number of papers [24][25][26] where their appearance in lattices of different symmetries was also discussed [27].
In this work, we consider planar finite arrays of twolevel 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 quasiflat 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 ∝ N −3  tot or ∝ N −5 tot (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 ω 0 (transition wavelength λ 0 ) and spontaneous decay rate Γ 0 = ω 3 0 |d| 2 /(3πℏε 0 c 3 ), where ε 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 ∼ 0.1λ 0 , λ 0 = 2πc/ω 0 , thus neglecting the Casimir interaction between atoms in the ground state due to its very fast spatial decay ∼ 1/r 6 [28].Consequently, the coupling between the atoms is determined by the dipoledipole interaction, governed by the electromagnetic free space Green's tensor G(R, ω) [29].
To describe such interaction, we introduce dipole moment operator of the i-th atom as di = dσ i +d * σ † i , where σi = |g i ⟩⟨e i |, σ † i = |e i ⟩⟨g i | are the atomic lowering and rising operators describing the transition from the excited state |e i ⟩ to the ground |g i ⟩ state of the i-th atom and vice versa, respectively.The interaction Hamiltonian of the system Ĥ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: Ĥeff = − Ntot i=1 di Ê(r i ), where Ê(r i ) is the electric field operator at the position r i .Exploiting classical principle of electric field superposition and the connection between the electric field and Green's tensor [30][31][32], one can directly obtain the effec-tive Hamiltonian [33,34]: where e d = 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, G(R, ω) ≈ G(R, ω 0 ), which is justified by a very narrow resonance of a single atom Γ 0 ≪ ω 0 .
The eigenvalues and eigenstates of the Hamiltonian (1) define the complex eigenfrequencies ω j −iΓ j /2 and eigenfunctions Ψ (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, ∆ω j = ω j − ω 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 ↔ 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, σ − , σ + , or σ 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 highsymmetry case of 2D atomic array N × 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.

A. Classification of the collective eigenstates
The eigenstates of the finite lattices can be classified and characterized based on their distribution in the quasi-reciprocal space [9,35,36] which can be done by expanding the dipole moments distribution of the j-th eigenstate Ψ (j) in the standing (Bloch) waves basis: mx,my sin(q 0 m x n x ) sin(q 0 m y n y ), (2) where c (j) mx,my are the expansion coefficients, q 0 = π/(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 ψ (mx,my) , 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 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 Ψ (j) should transform according to one of the irreducible representations (irreps) 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 ψ (mx,my) 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" ±1 (or their linear combination), i.e. phases of the neighbor corner dipoles differ by π/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 π).If m x ̸ = m y and (m x + m y ) is even the mode ψ (m1,m2) hybridises with the mode ψ (m2,m1) and they form two linear combinations: Antisymmetric modes transform according to A 2 or B 1 , while symmetric modes transform according to A 1 or B 2 irrep.

B. Mechanism of the radiation suppression
If the eigenstate Ψ (j) of the atomic array is mainly characterized by a single dominant contribution of a basis state ψ (mx,my) , its radiative losses Γ 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 → 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 Γ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 [37] 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 = ±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 [36, 38] [39].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.
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.

C. Numerical calculations of the finite array eigenstates
In Fig. 3(a) we show numerically calculated emission rates Γ/Γ 0 , normalized by the emission rate of a single atom, for two collective states of the 12 × 12 square array as a function of the normalized period ã = a/λ 0 .These states are characterized by the minimal emission rates among the others.Solid lines correspond to the σ z dipoles, while the dashed lines to the σ + or σ − ones.One can observe that for periods ã ≳ 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 σ 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 ã = 0.4 is shown in Figs.3(d,e).The distribution for the σ ± look similar.As expected, the B 2 state is the closest one to the M -point and it is well approximated by the ψ (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 ã ≳ 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 ψ (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 [40], 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 ψ (N,N −2) − state even for very small detunings from the square geometry of the order of 1%.On the other hand, the ψ (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 ã ≈ 0.3 the decay rates of both A 2 and B 2 states for σ 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 ã ≈ 0.3.At the same time, the calculations for the σ ± 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.

D. Infinite atomic lattices
Appearance of the external coupling near ã ∼ 0.3 for σ z polarized dipoles and its absence for the σ ± 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 ã ≳ 0.3 the interaction between the dipoles is rather weak, so that the dispersion is parabolic.For the lattice period near ã = 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 ã ≲ 0.294, the dispersion becomes nonmonotonous near the M -point along the Γ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 [41].In contrast, dispersion in the case of σ +,− 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.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 ã = 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 −5 tot ≡ 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 ∝ N −1.5 tot 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 ψ (N,N ) contribution, exhibit only ∝ N −3 tot ≡ 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 ≳ 300 for a given period ã = 0.4.In Fig. 5(b) we plot the calculated decay rates for the optimal periods ã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 ∝ N −3 tot and ∝ N −5 tot 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 ∝ N −3 to ∝ N −7 [10,41].Although, e.g. in a triangle geometry decrease of the decay rate also becomes close to ∝ N −5 tot , the absolute values of the decay rate still remain minimal for the A 2 /B 1 states of the square geometry.

B. Excitation of the subradiant states
In order to demonstrate the possibility of excitation of the strongly subradiant states appearing for the σ z polarization, we consider the scattering of vector Bessel beams [42,43] 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 σ 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 × 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 [43].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 ψ (N,N −2) − and ψ (N,N ) respectively.Lattice periods correspond to optimal periods in terms of radiative losses of these states, ã = 0.268 and ã = 0.281 in the case of 6 × 6 array.Cross-section spectra are normalized with a maximum value of a single atom cross-section σ 0 .
As one can see in Fig. 6 insets, the optimal incident field profile barely overlaps the atomic array from the outside, setting π 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 ψ (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 ψ (N,N −2) − contribution.Thereby a fine tuning of the incident beam frequency and angular momentum allows to selectively excite the desired longliving eigen mode of the array.

C. 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 [40,[44][45][46][47][48][49] and periodic or quasi-periodic nanostructured cavities [50,51].
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 [40,52,53].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 [51].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 ∝ N −5  tot .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.main text) fall into the same irreducible representation.In Figs.S1(a,b) eigenstates with the dominant contribution of the ψ (N,N ) and ψ (N −2,N ) harmonics, respectively, are shown.
For large enough periods and not very large anisotropy of the lattice, the state with dominant contribution of ψ (N,N ) harmonic does not interact with any other state upon the change of the a y /a x ratio around a y /a x = 1, since this state is not associated with the specific symmetry of the square lattice.Consequently, the radiative decay rate in Fig. S1(a) shown with the point curve changes only slightly even for relatively large detunings of the a y /a x from one up to several percents.On the other hand, the state with dominant contribution of ψ (N −2,N ) harmonic hybridizes with the ψ (N,N −2) harmonic when the periods are close to be equal, i.e. the lattice is almost square.Such symmetry-induced external coupling leads to the strong suppression of losses of such state in the square lattice.Reciprocally, one can say, that the antisymmetric states in the square lattice are rather sensitive to the deformation, resulting in the substantially increased losses even for a very small detunings ∆(a y /a x ) ≲ 1% from the square geometry, see Fig. S1(b), point curve.
For even larger detuning (a y /a x ≈ 0.9 (≈ 0.94) for B 2 (A 2 ) state for the given period ãx = 0.31) there appears additional external coupling between with the ψ (N,N −2) harmonic caused by the quasi-flat band dispersion of the lattice.This also leads to substantial suppression of losses with the same level as in the square lattice with optimal period, shown with solid curves in Figs.S1(a,b).Fig. S1(c,d) summarizes the effects of stretching/squeezing of the square lattice on the decay rates of the subradiant states.

B. CALCULATION OF THE INFINITE LATTICE DISPERSION
Here, we provide the details of the calculations of the dispersion of the various twodimensional dipole lattices in a free space.According to the Bloch theorem, it is possible to present a eigenstate of an arbitrary Bravais lattice as |ψ k ⟩ = j e ik•R j σ j |0⟩, where R j is a position vector of atom j = 1 . . .∞.Such Bloch state is fully defined by the quasimomentum k, which can always be chosen to be within the first Brillouin zone.Plugging a Bloch state into the Hamiltonian, Eq. ( 1) from the main text, for N → ∞ and multiplying both sides of it by ⟨ψ † k |, we obtain the following general relations for the eigenfrequency and the decay rate [1]:

∆ω(k) Γ
where • e d e ik•R ij , e d is the unit vector parallel to the dipole moment transition, and G(R ij ) is the electromagnetic Green's tensor of a free space [2]:

FIG. 1 .
FIG. 1.(a) A two-dimensional finite square regular lattice of N × N with σz transition and separated by distance a.The distribution of the highly non-radiative mode of B2 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.

FIG. 2 .
FIG. 2. The symmetry classification of the basis functions ψ (mx,my ) with respect to the irreducible representations and parity numbers mx and my [see Eq. (2)].Color shows the real part of the wavefunction amplitude at each lattice site.
FIG. 3. (a) Normalized radiative decay rate of two most subradiant eigenstates corresponding to B2 (green curves) and A2 (violet curves) symmetry in 12 × 12 square array as a function of the normalized period ã. (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.

FIG. 5 .
FIG. 5. (a,b) Radiative losses of subradiant states as a function of the total number of atoms in the array Ntot for different geometries for (a) fixed distance between the neighbor atoms ã = 0.4, (b) optimal period for each Ntot, shown in panel (c).Dashed lines indicate corresponding polynomial functions.

FrequencyFIG. 6 .
FIG. 6. Normalized scattering cross section spectra of 6×6 square atomic lattice irradiated by Bessel beam with (a) azimuthal number m = 9 and period ã = 0.281 and (b) m = 7 and ã = 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.
radiative decay rate of most subradiant B 2 (a,c) and A 2 (b,d) states in rectangular 12 × 12 array.In (a) and (b) period ãx along x-axes is fixed whereas period ãy along y-axis is changing.In colormaps (c) and (d) both periods ãx and ãy are varied.In (a) and (c) solid lines correspond to ãx = 0.29, dotted lines -to ãx = 0.31, whereas in (b) and (d) solid lines correspond to ãx = 0.301, dotted lines -to ãx = 0.31.

Figure S2 .
Figure S2.Table containing the schematics of the considered arrays and their main characteristics.