Stationary excitation waves and multimerization in arrays of quantum emitters

We explore the features of an equally-spaced array of two-level quantum emitters, that can be either natural atoms (or molecules) or artificial atoms, coupled to a field with a single continuous degree of freedom (such as an electromagnetic mode propagating in a waveguide). We investigate the existence and characteristics of bound states, in which a single excitation is shared among the emitters and the field. We focus on bound states in the continuum, occurring in correspondence of excitation energies in which a single excited emitter would decay. We characterize such bound states for an arbitrary number of emitters, and obtain two main results, both ascribable to the presence of evanescent fields. First, the excitation profile of the emitter states is a sinusoidal wave. Second, we discuss the emergence of multimers, consisting in subsets of emitters separated by two lattice spacings in which the electromagnetic field is approximately vanishing.


Introduction
Cooperative effects play an important role in the behavior and time evolution of quantum systems. The emission of light by well-separated atoms and molecules is uncorrelated, so that the states of a single quantum emitter determines the main characteristics of radiation. On the other hand, when quantum interference effects become important, typically in closely packed ensembles of emitters, cooperative effects take over and drastically modify the features of the emitted radiation. Typically, such effects compete with the dephasing induced by the dipoledipole interactions.
In this scenario, among the most interesting phenomena there are certainly super-and subradiance, by which the spontaneous emission of radiation in a transition between two atomic levels leads to coherent emission by an atomic ensemble, enhancing [1][2][3] or reducing [4] the decay rate. Clearly, these effects bear a profound influence on the formulation of the state of the radiating system. Since super-and sub-radiance are general phenomena, that can be observed also in large quantum systems, such as artificial atoms, quantum dots [5] and superconducting circuits [6], we will keep our discussion as general as possible, and refer generically to quantum "emitters" of radiation.
On the other hand, light scattering by cold atomic clouds can induce photon-mediated effective long-range interactions between atoms, leading to cooperative effects even at low atomic densities. While superradiance has been extensively studied since Dicke's seminal proposal [1], subradiance in large cold atom clouds has been observed only rather recently [7,8]. Although these phenomena take place at light wavelengths that are typically much larger than interatomic distances [9,10], a number of interesting quantum resonance effects appear at shorter wavelengths, comparable to the distance among atoms.
One of the most interesting features of such 1D systems is the presence of single-excitation bound states in the continuum (BICs) [29,30]. Such states represent an extreme case of subradiance, where the excitation is shared in a stable way between the emitters and the field, so that radiation is completely trapped by the emitters even though the energy would be sufficient to yield photon propagation. The physical characteristics of quantum emitters in waveguides has been studied in a variety of situations, for single [13,20,31,32] as well as double and multiple emitters [4,26,. In all cases, the quantum correlations between the emitters are due to the 1D photonic field.
The case of n = 3 and n = 4 emitters, analyzed via a non-perturbative treatment based on the Friedrichs-Lee model [63][64][65][66], was studied in [67]. The role of non-perturbative photonmediated interactions between emitters was shown to be crucial to find the correct bound states.
In this article, we will apply this formalism to find an explicit expression for the BICs for an arbitrary number n of emitters. We will focus on the global state of the emitters and will show that their excitation amplitude profile follows a sinusoidal law: single-excitation BICs are therefore excitation waves, conceptually analogous to the stationary spin waves emerging in 1D magnetic models [68][69][70]. We will also unearth the presence of multimerization: under suitable conditions, the excitation amplitude displays a modular repetition of quasi-decoupled components. One might expect that multimers arise independently of each other: by contrast, we shall see that the evanescent fields that are present between adjacent multimers play a crucial role. Our procedure will be based on a nearest-neighbor approximation for the inter-emitter interactions and will hold under very general assumptions.
The paper is structured as follows. In Section 2 we present the model, introduce its propagator matrix, and discuss its role in the evaluation of the BICs. In Section 3 we determine the structure of the BICs for an arbitrary number of emitters. Finally, in Section 4, we discuss the emergence of multimerized states and find a condition for constructing new BICs by arranging together smaller building blocks.

The model
In this section, after first introducing in Subsection 2.1 the basic features of the model to be investigated, we show in Subsection 2.2 how the BICs in the single-excitation sector can be conveniently described via the propagator matrix. The method outlined in this section will be then applied in Section 3.

Generalities
We consider a system of n identical two-level emitters, with excitation energy ε, equally spaced at a distance d (see Fig. 1), with ground and excited states |g j and |e j , coupled to a bosonic field with energy profile ω(k). The Hamiltonian reads where Figure 1: Schematic representation of the physical system: an array of two-level atoms, labelled by the index j, are characterized by a ground state |g j and an excited state |e j , separated by the energy ε. The atoms are placed at equal distances d and coupled to a linear waveguide mode.
is the free Hamiltonian and the interaction. In the above formulas, σ − j = |g j e j | and σ + j = |e j g j | are the lowering and raising operators of the j-th emitter, and b † (k), b(k) are the photon creation and annihilation operators, satisfying the canonical commutation relations The interaction has a rotating-wave form, with F (k) being the form factor that describes the coupling of the emitters with a boson of momentum k. 1 In order to proceed, some assumptions about the functions ω(k) and F (k) will be needed. A concrete example is represented by an array of two-level atoms coupled with a single transverse mode of the electromagnetic field in a linear waveguide. We assume a dispersion relation of the form expressed in natural units. In the case of a waveguide with rectangular cross section, whose sides are related by L z < L y , the "photon mass" m is proportional to the inverse smaller size L −1 z [71]. The form factor in this case reads where γ > 0 is a coupling constant. We shall focus henceforth on this case. However, as shown in Appendix B, the results that we shall discuss in the following turn out to be largely independent of the particular form of the coupling F (k) and the dispersion relation ω(k). We remark that, in order to simplify the computation, it is a common choice to replace F (k) with a constant form factor, and assign linear dispersion relations to left-and right-propagating photons (see, e.g., Ref. [41,72]). Although such an approximation is able to describe bound states in the continuum, it fails to capture non-Markovian effects [73], such as those which constitute the most relevant phenomenology presented in this work.
Observe that the rotating-wave form of the interaction Hamiltonian entails that the total number of excitations is conserved, [H, N ] = 0. In the single-excitation sector N = 1, the state vectors read where |vac is the field vacuum state, and a ∈ C n , ξ ∈ L 2 (R) are constrained by the normalization condition n j=1 Here, the vector a = (a 1 , . . . , a n ) is the excitation amplitude profile of the emitters and ξ(x) the photon wavefunction.

Propagator and bound states
The single-excitation sector was investigated for n = 2, 3, 4 atoms in Refs. [36,37] and shown to contain, for generic dispersion relations and form factors (in particular, for those given in Eqs. (4)- (5)) and for selected interatomic distances, nontrivial atom-photon bound states. Our main objective is to extend these results to arbitrary n and unearth genuine collective effects. Bound states are given by the solutions of the equation for some real energy E, with |Ψ as in Eq. (7) and normalized as in Eq. (9), for some a and ξ(x).
In particular, bound states with energy E > ω(0) are bound states in the continuum (BICs), since they are embedded in the continuous set of the field frequencies.
As shown in Appendix A, the content of Eq. (10) can be rephrased as follows. We consider the propagator matrix of the model, which we compute in Appendix B: where k(E) is the positive solution of ω(k) = E, i.e. the momentum of a photon with energy E, and β j (E) a real-valued function satisfying Interestingly, all terms β j (E) would vanish if we would chose a linear dispersion relation and a constant form factor, namely ω(k) = ω 0 + v(k − k 0 ) and F (k) = F 0 . However, any more accurate physical model with a correction to the linear-dispersion-relation approximation and/or the flat-coupling approximation will give nonvanishing values of β j (E). This is a crucial point: despite these terms being (for j = 0) exponentially small in md, they will play a fundamental role in determining the very structure of the BICs, as we will see in the next sections.
Then, a real energy E is an eigenvalue in Eq. (10) if this energy will corresponds to a BICs if E > ω(0). As shown in Appendix A and C, a corresponding eigenstate |Ψ is such that • the vector of atomic excitation amplitudes a = (a 1 , a 2 , . . . , a n ) is a solution of the matrix equation with the additional constraint • the wavefunction ξ(x) is where with the correction η(x) being real and satisfying Eq. (16) shows that the photon wavefunction ξ(x) corresponding to a BIC with energy E is a linear combination of n copies of ξ 1 (x), each centered at one of the emitters, x j = (j − 1)d: the contribution of each emitter is weighted by the corresponding component of a. Again, the term η(x), which physically corresponds to a small evanescent field, would again vanish if we would take a linear approximation of the dispersion relation. All relevant quantities are summarized in Table 1.

Quantity Description Definition ε
Excitation energy of the emitters d Distance between consecutive emitters ω(k) Dispersion relation of the boson continuum (4) m Lower bound of the continuum (mass cutoff) Propagator of the model, and its (j, )th element (11)  Summing up, in order to compute the BICs emerging in our model, we need to solve Eqs. (13) and (14), with the propagator matrix being given by Eq. (11); once we solve this problem, thus finding the admissible energies E of the BICs and the corresponding excitation profiles a, the associated photon wavefunction is given by Eq. (16). The next section will be devoted to this problem.

Structure of the bound states in the continuum
The evaluation of BICs for an array of quantum emitters was performed, under less general hypotheses, in Ref. [36] for n = 2 emitters and in Ref. [67] for n = 3, 4 emitters. When n > 4, the eigensystem becomes more and more involved because of the presence of the terms β j− (E) in the expression of the propagator, whose number grows with n.
At the physical level, these terms may interpreted as (field-mediated) inter-emitter interaction terms. Even though these terms are exponentially suppressed in md, and hence "small" in the physically interesting regime md 1, their role turns out to be fundamental: without these contributions the eigensystem would be degenerate and no preferred set of stable emitter configurations can be found, while, when taking into account the additional terms, such degeneracy is lifted.
This phenomenon can be better understood when looking at the photon wavefunction associated with a given emitter configuration: as shown in Eq. (17), the latter will be the sum of a dominant sinusoidal contribution plus an evanescent contribution η(x), given by Eq. (18). Albeit small, such tails are not confined and induce field-mediated coupling among all emitters, which causes instability unless particular configurations are chosen. An example is displayed in Fig. 2: emitter configurations in which only two neighboring emitters resonate, with the other ones "switched off", are not allowed even though the state would be stable in the absence of the other emitters.
Therefore, a generic BIC in a multi-atom array is expected to be a collective state: the excitation on the BICs must be distributed in a very specific way along the whole array, and must be chosen in such a way that the effects of all such fields even out, yielding a stable state.
In the absence of evanescent field one would have a manifold of degenerate local solutions with the same energy E, consisting of radiation trapped between pairs of emitters, that could be pasted together in an arbitrary way. However, in the presence of evanescent fields, most of these states are unstable, and only a few selected ones emerge as stable configurations. The role of the evanescent fields is to select collective stable states. In particular, even those emitters (if any) which have zero amplitude and thus are decoupled from the field must be placed at specific positions along the chain (namely, the nodes of the field), and, in this sense, they also cooperate to yield the overall stability of the state. We will see this selection mechanism at work in stabilizing multimer configurations in Sec. 4. Figure 2: An example: configuration with n = 3 emitters in positions x = 0, d, 2d, in which a 1 = a 2 , a 3 = 0 and the oscillating part of the field is confined between the atoms that share the excitation. If there were no evanescent contributions, the field would be exactly confined between the first two emitters and the resulting two-emitter state would be stable, regardless of the presence of the third emitter at x = 2d. However, the presence of an evanescent field, which couples with the third emitter, causes instability, as it tends to transfer the excitation to the third emitter, and eventually to decay. The heuristic discussion above will be supported by calculations. Collective configurations will emerge as a natural consequence of the structure of the propagator: we will prove that the atomic excitation amplitudes in a stable configuration are interpolated by (either exact or "deformed") sinusoidal functions, providing an interesting analogy with the stationary spin waves that emerge in the Heisenberg model [68][69][70].
The present section is organized as follows: • in Subsection 3.1 we present the idea on which the computation of BICs is based; • in Subsection 3.2 we will sum up the evaluation of the energies of all possible BICs and their corresponding excitation profiles a, obtained as solutions of Eq. (14), with all details being reported in Appendix D; • in Subsection 3.3, after computing the photon wavefunctions ξ(x) corresponding to such profiles via Eq. (16), we discuss our findings.
A particular feature of some BICs that comes out as a direct consequence of our discussion, namely the presence of multimerized BICs, will be investigated in Section 4.

A nearest-neighbor approximation
We shall search for bound states in the continuum with values of energy (close to) E = E ν , where we define the νth resonant energy by This choice is motivated by the following observation: if we were to discard all terms β j− (E) with j = in the propagator (11), Eq. (13) would be satisfied if and only if E = E ν for some ν, and thus all BICs would be found at such values of energy. Correspondingly, the system would exhibit an (n − 1)-dimensional degenerate space of BICs with energy E = E ν . This phenomenon is extensively discussed in Ref. [67] and, for completeness, revised in Appendix D.3. As we have seen above, this is the case for a linear dispersion relation and a flat form factor, when all β j (E) for j = 0 vanish. The degenerate situation outlined above is drastically modified when one takes into account the full structure of the propagator, which is the case when one goes beyond the linear approximation. This problem was already analyzed for the cases of n = 3 and n = 4 emitters [67], where it was shown that the presence of nonvanishing off-diagonal terms β j− , no matter how small, lifts the degeneracy: the available BICs, with energy either equal or close to E ν , emerge for distinct, albeit close, values of ε. This means that the system exhibits at most one truly stable BIC. Besides, as discussed in Appendix D.2, BICs turn out to have a well-defined parity.
When n > 4, the analytic study of the full structure of the propagator becomes unfeasible, and therefore we need a way to extract information about the resonant degeneracy breaking patterns, while keeping calculations viable. The basic idea is the following: since the off-diagonal terms with |j − | > 1 are of higher order in e −md than those with |j − | = 1, and therefore, in order to study the degeneracy lifting at large md, it is enough to consider the first-order approximation Indeed, we will now show that the presence of the largest of such contributions, namely β 1 , is sufficient to fully remove the n−1 degeneracy and thus determine the nondegenerate BICs in the large md regime. Physically, since β j− is related to a field-mediated interaction involving the jth and th emitter, the replacement (21) may be regarded as a nearest-neighbor approximation.
Higher-order corrections in md (that is, terms β j− (E) with |j − | > 1) would only contribute by O(e −2md ) corrections, without changing the qualitative picture that we will outline.

Exact and deformed excitation waves
In the approximation discussed above, the propagator (11) of the model, evaluated at a resonant energy E ν (19), can be written in a convenient form. By directly substituting E = E ν in (11), where is the zeroth order approximation of the propagator, that is, the one that would be obtained by discarding all terms β j− with j = . By Eq. (11), it reads and corresponds to the adjacency matrix (discrete Laplacian) of a one-dimensional regular graph [75,76] with Dirichlet boundary conditions. Its main properties are recalled in Appendix E.
Consequently, when taking into account the role of the terms β 1 (E) in the propagator (physically, nearest-neighbor photon-mediated interactions between the emitters), we are perturbing the original propagator G −1 0 (E) with a multiple of the matrix ∆ n . Despite the latter term being exponentially small in md (recall Eq. (12)), it does crucially affect the solutions of Eq. (13), since it lifts the degeneracy of the original problem and allows us to individuate an unambiguous set of solutions each corresponding to a BIC emerging at an energy (close to) E ν when the excitation energy ε has a certain value.
This problem is studied in Appendix D.4, to which we refer for details. We will present here the results. First of all, for all j = 1, . . . , n, define the quantities χ (j) and a (j) via where N j is a normalization constant. These are the eigenvalues and eigenvectors of the matrix ∆ n , as shown in Appendix E; in particular, the eigenstates of ∆ n can be obtained by a sampling at the equally-spaced points x = ( − 1)d of sine functions with wavelength We will refer to the eigenvectors a (j) of ∆ n as excitation waves. In analogy with the two branches of phonon modes in lattices, we will label excitation waves with low and high frequency as acoustic (consecutive amplitudes in phase) and optical modes (consecutive amplitudes in phase opposition), respectively. The BICs of the model are then characterized as follows: • n 2 BICs can emerge at an energy E = E ν : each of them corresponds exactly to one of the sinudoidal waves a (j) , with j ranging on either even or odd integers as described in Table 2. Each of these BICs is "switched on" when the excitation energy ε of the emitters has exactly the value x d • n 2 BICs can emerge at an energy E ≈ E ν : the corresponding excitation profile is a "deformed" excitation waveã (j) = a (j) + δ (j) , with j ranging on remaining integers, the deformation to be evaluated numerically (see Figs. 8-9 in Appendix D.4). Each of these BICs is "switched on" when the excitation energy ε of the emitters has a value to be evaluated numerically as well.
The first result follows easily from Eqs. The physical situation can be described as follows. In order to steadily sustain a BIC, the n emitters must collectively share a part of the excitation following any of the sinusoidal amplitude waves with wavelength (26): depending on the values of ν and n, waves with a given parity (either even or odd) will be exact, and waves with the converse parity will be distorted in order to support a stationary excitation of the photon field. n even n odd j even j odd j even j odd ν even Table 2: Given a system of n emitters, n − 1 BICs at energy E = E ν (exact excitation waves) or E ≈ E ν (deformed excitation waves) can emerge, each for a specific value of the excitation energy ε of the emitters. The jth BIC, ordered by increasing frequency, will be exact or deformed depending on whether n, ν and j are even or odd. Notice that either the first (j = 1) or the last (j = n) sinusoidal excitation wave, for ν even or odd, respectively, will not have any deformed counterpart as eigenvector of the propagator: the corresponding eigenvector will be an unstable state whose eigenvalue has a large imaginary part, i.e. a superradiant state.

Discussion of the results
Once we have computed the amplitude waves that can sustain a BIC at energies (close to) E ν , Eqs. (16)- (17) give directly the associated photon wavefunction. Figures 3 and 4 display, for n = 30, and E = E 1 and E = E 2 , respectively, the simplest acoustic and optical atomic excitation waves that emerge in the emitter array, together with the associated field wavefunction ξ(x), in properly normalized units. Under our assumptions, the photon wavefunction (18) is largely dominated by the contribution of the single pole at k(E). By neglecting the small η term, it reads, in correspondence of the E = E ν resonance, where It is immediate to see that the contribution (29) to the photon amplitude vanishes identically outside the emitter chain. The value of each a (j) is related to the discontinuity of the first derivative of ξ(x) at the position of the th emitter; in particular, the wavefunction is smooth in a neighborhood of the th emitter if and only if a (j) = 0. In practice, each excited emitter induces a jump, proportional to the excitation amplitude, in the derivative of the photon wavefunction. All these results are compatible with the ones for n = 3, 4 reported in [67].

Multimerization
This final section will be devoted to a better understanding of an interesting phenomenon, which emerges in the structure of the BICs computed in Section 3. A detailed scrutiny of the results brings to light a multimerization effect: some BICs are characterized by a modular structure, in which the same excitation amplitude configuration repeats for a certain number of times along the chain, with each "monomer" separated from the adjacent one by an emitter in its ground state. In the upper panels of Fig. 5, the configuration is a dimer, with the central emitter working as a dynamical mirror [36], by forcing the field excitation amplitude to vanish at its position. This phenomenon is more general: BICs can, in fact, split into more than two parts. A trimer and a tetramer are shown in the lower panels of Fig. 5. All multimerized states appear to be composed of a number r of identical monomers, made up of h emitters, separated by a single emitter in its ground state. The number of components and the size of each component are related with n by n = rh + r − 1.
Moreover, each module is in itself a BIC for an array of h emitters. We will interpret this interesting phenomenon as a direct consequence of a simple mathematical property of the adjacency matrix ∆ h in (24), and thus of the particular structure (22) of the propagator of the model in the nearest-neighbor approximation.
In such a case, the action of the n-dimensional matrix ∆ n on the full array splits into the action of the h-dimensional matrix ∆ h on each block. As a direct consequence: • if a ∈ C h is an antisymmetric eigenvector of ∆ h , then a h = −a 1 and the vector A = (a, 0, a, 0, a, 0 is an eigenvector of ∆ n with the same eigenvalue; • if a ∈ C h is a symmetric eigenvector of ∆ h , then a h = a 1 and the vector is an eigenvector of ∆ n with the same eigenvalue. It is easy to show that the mathematical property outlined above, together with the expression (22) for the propagator of the system in the nearest-neighbor approximation, does explain the multimerization phenomenon outlined at the start of this section. As discussed in Subsection 3.2 (also see Table 2), an h-emitter excitation wave a (i.e. an eigenvector of ∆ h ) with the "right" parity, i.e. such that u ν · a = 0, satisfies provided Eq. (108) holds, and is thus associated with a BIC for a chain of h identical emitters. Now, the n-component vector A defined either via Eq. (39) or (40) (depending on the value of ν) is an eigenvector of ∆ n with the same eigenvalue and, by construction, satisfies the condition U ν · A = 0, with U ν being the n-component analogue of u ν : consequently, it is an eigenstate of the propagator and thus corresponds to a multimerized BIC for the n-emitter chain.

Discussion of the results
Our findings demonstrate that it is always possible to construct a multimerized BIC of n = rh + (r − 1) emitters from smaller (building) blocks of h emitters, separated by single emitters that act as dynamical mirrors. The rules for constructing such states are, nevertheless, very specific: the blocks must be equal and the way of connecting them to each other must follow one of the prescriptions (39)-(40), depending on their symmetry. These conditions are ultimately due to the small corrections to the propagator and the wavefunctions due to the evanescent fields.
In the absence of evanescent fields, multimerized configurations could be constructed arbitrarily among the infinite possibilities in a degenerate (n − 1)-dimensional eigenspace, by superposing "local" excitations involving few neighboring emitters: this is due to the fact that adjacent monomers would not interact with each other, since no photon field would be present between them, and our model only involves photon-mediated interactions. On the other hand, the presence of evanescent fields modifies the picture by inducing interactions between adjacent monomers, thus making most of multimerized one-excitation states unstable, and selecting only a few ones as stable configurations. Hence, the role of the evanescent fields is to induce the emergence of collective stable states. In particular, our findings show that one emitter in its ground state must be present between neighboring monomers in order to ensure stability: in this sense, even the emitter that does not share the excitation cooperates with the rest of the system to form the bound state.

Conclusions and outlook
We analyzed the emergence of BICs in a regular array of quantum emitters in a waveguide. BICs are present for an arbitrary number of emitters, and the excitation profile of the emitter states is a sinusoidal wave. We also discussed the presence of multimers, separated by two regions in which there is practically no electromagnetic field. A crucial role is played by the evanescent fields generated by quantum emitters in determining the physical structure of bound states in the continuum.
The techniques adopted in this article hinge upon an analysis of the propagator and are intrinsically non-perturbative. The main physical factor that limits the robustness of BICs is photon loss from the waveguide: we did not discuss here these losses, but their effect can be estimated by applying standard techniques based on the master equation.
The collective behaviour unearthed in our analysis does not depend on the specific geometry of the waveguide and the details of the dispersion relations, and is therefore valid for a number of physical implementations. The identification of (manifolds of) states in which a single excitation is coherently shared among distant artificial atoms enables one to analyze the features of longrange coherence, mediated by the photon field, and could possibly provide a tool to control qubits at an arbitrary distance.
Moreover, the manifolds of long-lived states, whose lifetime is determined by the waveguide losses and not by the atom-field coupling, can be used as quantum registers/memories [77] due to their dynamical stability, in particular against spatially separated decoherence sources [78].
These features can make them useful in hydrid situations [79,80].
Finally, we mention that additional properties emerge when the emitters are confined in a finite waveguide [81] or in 2-dimensional geometries [82], or when topological effects contribute to the robustness of the dressed states [83]. These situations and related aspects will be investigated in the future.

A Eigenvalue equation
In this appendix we derive the eigenvalue equations for BICs discussed in Sec. 2. Consider the eigenvalue equation where H is Hamiltonian (1)-(3) and |Ψ state (7), and project it onto the basis vectors |e j ⊗|vac and |g ⊗ b † (k) |vac . One gets The second equation gives and these yield, respectively, Eq. (16) by a Fourier transform, and the constraint (15), which is needed in order to make the integral in Eq. (16) well-defined. By plugging (44) into the first equation in (43) we get which, by defining G −1 (E) as in Eq. (11), is Eq. (14). The latter equation admits nontrivial solutions if and only if Eq. (13) holds.

B Calculation of the propagator
In this appendix we will compute the propagator for our model, thus showing the second part of Eq. (11). Let us introduce the self-energy matrix: so that the propagator can now be expressed as We will find an expression for the self-energy matrix (47), and thus for the propagator, under general assumptions on the dispersion relation ω(k) and the coupling function F (k). We make the following assumptions: (i) ω(k) and f (k) = |F (k)| 2 are even real-valued functions satisfying the integrability condition: (ii) ω (k) > 0 for all k > 0, so that ω(k) is an increasing function for k > 0, and ω(0) > 0. (iv) the following conditions hold: uniformly in the strip K, and (v) there is an open connected set A ⊂ C containing the real half line (ω(0), +∞) such that, for z ∈ A, the equation ω(κ) = z admits exactly two solutions in K; We remark that an expression for the self-energy can be found without condition (iii), and that the discussion that follows can be extended to a nonmonotonic ω(k), provided that one looks for the eigenvalues in an energy range (E 1 , E 2 ) such that the equation ω(k) = E has a single positive real solution for all E ∈ (E 1 , E 2 ). First of all, since both ω(k) and f (k) are even real-valued functions, i.e. satisfy, for all k ∈ R, their analytic continuations ω(κ), f (κ) to the strip K of the complex plane, with K ± = K ∩ C ± , satisfy the following symmetry properties: The self-energy matrix is defined for j, = 1, . . . , n and z ∈ C by and, because of the symmetry properties of the integrand, satisfies for all z ∈ C and j, = 1, . . . , n, so that we only need to compute it for j ≥ or, equivalently, we can substitute j − with |j − |.
We will evaluate the integral through a contour integration in the complex plane (see Fig. 6), by closing the integration over [−R, R] with a properly chosen curve and exploiting the residue theorem. The integrand in (58) has poles at the solutions of the equation ω(κ) = z. Let us consider z ∈ A: since ω(−κ) = ω(κ), by assumption (v) such solutions come in pairs ±κ(z). By the analyticity of ω(κ), they are locally continuous and, moreover, analytic away from the critical points of ω(κ), i.e. the points at which ω (κ) = 0. Aside from such (countably many) points, each solution κ(z) is a simple pole for the integrand. Choose R > 0 large enough so that the pair of poles is included in the rectangle and K ± R = K R ∩ C ± . We will close the integration contour on the boundary of K + R , consisting of two horizontal segments, [−R, R] and [−R, R] + im, and two vertical segments, ±R + i[0, m]. For j = , by hypothesis (iv) the integral on ∂K + converges and there is no contribution from the vertical segments as R → ∞, hence For j = , we have and thus, i.e. all off-diagonal contributions are finite as well and, moreover, exponentially suppressed in md.
The integral over the contour will equal the sum of the residua of all singularities enclosed by the contour, that is, all solutions of the equation ω(κ) = z with positive imaginary part. By our assumptions, for z ∈ A there is precisely one couple of singularities ±κ(z), with Im κ(z) > 0 for all z ∈ A ∩ C + , each being continuously linked with the corresponding real solution ±k(E) of the real equation ω(k) = E for E > ω(0). Indeed, by taking into account the equation and expanding around the solution k(E) of the equation thus implying that the pole contained in the integration contour is +κ(z), as shown in Fig. 6. By the residue theorem we thus have where and Res g j , κ(z) denotes its residue at the pole κ(z). For all noncritical values of z, we have hence lim κ→κ(z) and thus the residue is Thus we finally obtain the desired result: for all E > ω(0), with a residue term and a contour contribution Notice that β j (E) = β j (E) * is real, since by using the symmetry properties (55) it easy to see that the integrand function has the symmetry g j (k + im) * = g j (−k + im). Therefore, the contour term gives a contribution to the self-energy matrix Σ(E) which is real Hermitian, whereas the residue yields a contribution which is non-Hermitian.

C Calculation of the photon wavefunction
In this appendix we apply the results of Appendix B in order to compute the integral in (16), which gives the boson wavefunction ξ(x) in the position representation. From Eq. (16) we immediately get where By assuming that F (k) is real and has the same properties of f (k) in (i)-(iv) of Appendix B, we can write where can be computed by the residue theorem. The computation is a carbon copy of the computation of the self-energy in Appendix B, and by replacing f (k) with F (k)/ √ 2π in (70)-(72) one has with a residue term and a (real valued) contour contribution Analogously one gets since for z = E − iδ the pole with positive imaginary part is −κ(z) and ω (−k(E)) = −ω (k(E). Finally, by plugging (77) and (80) into (75) we obtain with η(x) real and of order O(e −m|x| ), that is the expression (17) in the main text.

D Calculation of the BICs
In this appendix we will evaluate explicitly the bound states in the continuum in our model, thus proving the results reported in Section 3. The appendix is organized as follows: • in Appendix D.1 we rephrase the problem in a convenient form; • in Appendix D.2 we show that, because of the presence of corrections β j− (E) to the propagator, BICs are predicted to have a definite parity; • in Appendix D.3 we briefly revise the problem in the approximation in which all contributions β j− (E), j = , to the propagator in Eq. (11) are neglected: in this approximation, BICs emerge at resonant energies E = E ν , ν = 0, 1, . . . , and for each of these values there is an (n − 1)-dimensional space of degenerate BICs that are "switched on" for the same value of the excitation energy ε of the emitters; • in Appendix D.4 we refine our analysis by taking into account the largest correction β 1 (E) to the propagator, showing that this correction alone is sufficient to lift the degeneracy: while BICs still emerge (nearly) at resonant energies E = E ν , we will now be able to distinguish n − 1 nondegenerate BICs, corresponding to distinct (albeit close) values of ε and to well-defined collective structures.

D.1 General strategy
Before getting started, let us rephrase Eqs. (13)- (14) in a more compact way which will turn largely useful for our purposes. We start by defining the two quantities with Z(E) as in Eq. (71). The propagator matrix G −1 (E) in (11) can be written as where we introduce the matrix A(θ, b) by for j, = 1, . . . , n, b = (b 1 , . . . , b n−1 ). It is immediate to show that, by Eq. (83), the two equations (13)- (14) correspond to an E-dependent eigenvalue and eigenvector problem for the matrix A(θ, b), namely In this way, the evaluation of the BICs can be performed as follows: 1. we search for the real eigenvaluesχ(θ, b) of the matrix A(θ, b), evaluate the corresponding eigenvectors a and check that the constraints (15), with k(E) = θ/d, are satisfied; 2. we obtain the energies E of the corresponding BIC by solving the equation with χ(E) as in Eq. (82); 3. finally, we obtain the photon amplitude of the eigenfunction by substituting the values of E and a in (17) .

D.2 Symmetry properties of the BICs
Even without explicitly solving the aforementioned problem, a fundamental property of BICs in our system can be inferred by the structure of the matrix A(θ, b): indeed, for all j, = 1, . . . , n, i.e., its entries are a function of the distance from the main diagonal. This is an immediate consequence of the fact that the Hamiltonian (3) is parity-invariant and that the emitters are equally spaced. Consequently, A(θ, b) is centrosymmetric, i.e. satisfies the property where J n is the exchange matrix, i.e. the matrix having ones on the counterdiagonal and which admits ±1 as eigenvalues, with corresponding eigenvectors (a 1 , a 2 , . . . , a n ) satisfying a n−j+1 = a j for j = 1, . . . , n, or a n−j+1 = −a j , for j = 1, . . . , n, respectively. As a consequence, J n and A always share a common basis of eigenvectors, which means that A always admits a basis of eigenvectors that are either centrally symmetric or antisymmetric. Moreover, if the eigenvalues of A are nondegenerate, then necessarily its corresponding eigenvectors will be centrally symmetric or antisymmetric. Physically, this is caused by the presence of tails of the photon wavefunctions, due to the fact that the dispersion relation is generally bounded from below and nonlinear. Though exponentially suppressed in the interatomic distance, these tails couple with all the emitters of the configuration, generally causing instability unless proper symmetry conditions are fulfilled [67]. Finally, we remark that the property discussed above is exact. While the explicit calculation of the BICs in Appendix D.4 is performed under a nearest-neighbor approximation (and does, indeed, yield eigenstates with definite symmetry), taking into account higher-order terms, or even the full structure of the propagator, would not spoil this property, which is thus more fundamental.

D.3 The degenerate case: resonance condition
As discussed in the main text, the terms b j− in the matrix (84) are exponentially suppressed in |j − |md, and are therefore small in the physically meaningful regime. As a first approximation, one may simply discard them and proceed with the study of the eigensystem of the matrix This matrix admits real eigenvalues only if θ = νπ for some ν ∈ N, when it becomes the rank-one matrix [67] A(νπ, 0) = i u ν u ν , with The spectrum of A(νπ, 0) is thus composed of • the simple eigenvalueχ = in, associated with the one-dimensional eigenspace spanned by u ν ; • the (n − 1)-degenerate eigenvalueχ = 0, whose eigenspace contains all the atomic amplitudes a orthogonal to u ν : u ν · a = 0.
While the eigenvector u ν corresponds to an unstable (superradiant) state of the system, the eigenvectors associated with the zero eigenvalue correspond to stable configurations of the system, which emerge if and only if The first equation is a resonance condition: the emitter spacing d must be an integer multiple of the half-wavelength π/k(E) corresponding to the energy E. Hence, if we define the resonant energies a BIC with energy E exists if and only if Observe that by Eq. (82) the second condition reads As a consequence, a BIC emerges in the spectrum whenever its energy E ν equals the excitation energy ε of the emitters (see Fig. 1) plus a coupling-dependent correction, that is small in the perturbative regime. We remark that condition (95) reads, explicitly, n j=1 (−1) jν a j = 0.
Due to degeneracy, no symmetry condition is imposed a priori onto the eigenvectors. Therefore, states in which the excitation is shared only by two consecutive emitters with opposite amplitudes, such as those depicted in Fig. 2, represent eigenstates as valid as the ones with central symmetry, provided the evanescent fields are neglected.

D.4 Nearest-neighbor approximation: degeneracy lifting and excitation waves
As discussed in Subsection 3.1, the degenerate situation outlined in Appendix D.3 is drastically modified when one takes into account the full structure of the self-energy and the propagator, and it suffices to only take into account the correction b 1 (E) in order to understand the "direction" in which the degeneracy is broken. In this case, an interesting structure emerges: adding the nearest-neighbor approximation to Eq. (93), we obtain where ∆ n is the n × n matrix in Eq. (24), whose properties are studied in Appendix E. Eq. (101) means that, when taking into account the role of b 1 , the eigenproblem studied in Appendix D.3 is perturbed by a term which is proportional to the adjacency matrix ∆ n . Since the unperturbed matrix A(νπ, 0) is rank-one, and thus highly degenerate, its eigensystem will be crucially affected by the structure of the perturbation, even for small b 1 .
Here, we are interested in the spectrum and eigenvectors of the matrix (101) in the limit b 1 → 0, since the nearest-neighbor approximation is physically valid when b 1 is small. It is worth noticing that the obtained result have an interesting universal character, being independent of the specific form of b 1 .
As shown in Appendix E, the eigenvalues χ (j) of the matrix ∆ n , with 1 ≤ j ≤ n, can be evaluated exactly, see Eq. (25) in the main text: these are the sinusoidal excitation waves.
By applying the propagator at resonance, θ = νπ, we get, The excitation waves a (j) and the vector u ν are related as follows: • if n is even, then u ν · a (j) = 0 for even j + ν (103) u ν · a (j) = 0 for odd j + ν (104) (a) j = 2, exact excitation wave (open orange circles) of the matrix A(νπ, b), for n = 100 emitters and even ν, are compared with the excitation waves a (j) (full blue circles), eigenstates of the nearest-neighbor adjacency matrix ∆ n . The upper and lower panels represent acoustic and optical waves, respectively. Eigenstates with even j coincide with the excitation waves, while those with odd j are deformed. The effect of deformation is more relevant for acoustic waves than for optical waves. Notice that the excitation wave with j = 1 corresponds to an eigenvalue with large imaginary part. The amplitude vectors a (j) andã (j) ν are normalized to their respective square norms.
• if n is odd, then u ν · a (j) = 0 for even j (105) u ν · a (j) = 0 for odd j Moreover, the products u ν · a (j) , even when nonvanishing, tend to become negligible for large (small) j, when ν is even (odd), as reported in Fig. 7. This is due to the fact neighboring excitation amplitudes in a (j) tend to have the same sign for small j and opposite sign for large j. Therefore, the scalar product with u ν for even ν, which is a constant vector, tends to vanish for states with large j, while the product with u ν for odd ν, which is a staggered vector of ±1, is minimal for small j.
As a consequence, • if n is even, then A(νπ, b) admits as eigenvectors the excitation waves with even j (respectively odd j) iff ν is even (respectively odd); • if n is odd, then A(νπ, b) admits as eigenvectors the excitation waves with even j, for both even and odd ν;  (open orange circles) of the matrix A(νπ, b), for n = 100 emitters and odd ν, are compared with the excitation waves a (j) (full blue circles), eigenstates of the nearest-neighbor adjacency matrix ∆ n . The upper and lower panels represent acoustic and optical waves, respectively. Eigenstates with odd j coincide with the excitation waves, while those with even j are deformed. The effect of deformation is more relevant for optical waves than for acoustic waves. Notice that the excitation wave with j = 100 corresponds to an eigenvalue with large imaginary part. The amplitude vectors a (j) andã (j) ν are normalized to their respective square norms.
in the cases listed above, the excitation waves are also eigenvalues of the at the resonance energy E ν given by (97).
Gathering all cases together, we observe, for every value of n and ν, that n 2 out of the n excitation waves are also eigenvectors of the propagator matrix, and thus correspond, for some excitation energy ε of the emitters, to a bound state in the continuum with resonant energy E = E ν . Specifically, an admissible BIC according to the aforementioned criterion is present in the spectrum of the Hamiltonian provided and thus if and only if the excitation energy ε of the emitters satisfies Notice that, since the χ (j) 's are all distinct, in principle only one among the possible BICs will actually be a stable state for the array, even though, in practice, all such states are expected to be long-lived, with an O(e −md ) lifetime. The remaining n 2 excitation waves are not eigenstates of the propagator, since they do not satisfy the condition u ν ·a (j) = 0. A numerical analysis of the spectrum of the matrix A(νπ, b) for small b 1 shows, besides the already discussed n 2 exact excitation waves, the following features: • One of the eigenvectors of A(νπ, b) correspond to a complex eigenvalue with a large imaginary part. This state, characterized by j = 1 for even ν, and j = n for odd ν, represents a perturbation of the unstable eigenvector u ν of the matrix A(νπ, 0), with eigenvalueχ = in.
• The remaining n 2 − 1 eigenvectors are deformed versions of the excitation waves that do not exactly satisfy the resonance condition. Their associated eigenvaluesχ ±,(j) have real parts independent of b 1 and close to the eigenvalues χ (j) in Eq. (25), while they also acquire an imaginary part (namely, a finite decay rate) that vanishes with b 1 .
Summarizing, as reported in Table 2 in the main text and Figs. 8-9, the (n − 1)-dimensional degenerate eigenspace associated with the zero eigenvalue of the approximate propagator A(νπ, 0) breaks into n − 1 nondegenerate eigenvectors, n 2 of which are the exact sinusoidal waves that are orthogonal to u ν (i.e., they satisfy the resonance condition at θ = νπ), while the remaining ones are "deformed" versions of the sinusoidal waves that have a small [O(b 1 )] but finite projection onto u ν . The jth deformed excitation wave will be written as with δ (j) ν being the deformation with respect to the jth exact sinusoidal wave. We also remark that, due to the behavior of the product u ν · a (j) , the size of the deformation of the wave crucially depends on j in an opposite way for the two cases (see Fig. 8): • for even ν the deformation δ (j) ν is larger for small j (i.e. for acoustic waves), while it becomes negligible for large j (i.e. for optical waves); • conversely, for odd ν, the deformation δ (j) ν is larger for large j (i.e. for optical waves), while it becomes negligible for small j (i.e. for acoustic waves). This is a natural consequence of the fact that, as shown in Fig. 7, for off-resonant excitation waves, the term u ν · a (j) monotonically decays (increases) with j for even (odd) ν. In other words, the deformation will be relevant only for acoustic waves if ν is even and for optical waves if ν is odd, for every value of n.
As already mentioned, the deformed waves described above correspond to eigenvalues with a finite imaginary part, and therefore do not represent stable states. However, it is possible to find an actual BIC with energy E E ν such that close to each deformed wave. Moreover, since the imaginary parts of the deformed wave amplitudes are negligible, as n increases they tend to coincide with very good approximation with the actual bound states, and the condition (108) at which they appear in the spectrum is still valid.

E Eigenvalues and eigenvectors of the adjacency matrix
Given n ∈ N, let ∆ n be the adjacency matrix (24). In this appendix we will discuss the interpretation of the eigenvalue problem for the adjacency matrix ∆ n as a discrete boundary value problem for a one-dimensional regular array, and evaluate its eigenvalues and eigenvectors. We intend to solve the eigenvalue equation where χ ∈ R and a = (a 1 , . . . , a n ) ∈ C n . For this purpose, let us consider an array of n + 2 elements by adding two fictitious nodes a 0 , a n+1 at the extrema of the array a. By construction, the eigenvalue problem for the matrix ∆ n is equivalent to a recurrence relation for the array a 0 , a 1 , . . . , a h , a h+1 with two boundary conditions at the extrema of the array: a +1 = χ a − a −1 , = 1, . . . , n, a 0 = a n+1 = 0, which can be interpreted as a discrete version of the Helmholtz equation with vanishing (Dirichlet) boundary conditions. The straightforward trigonometric identity sin( + 1)θ = 2 cos θ sin θ − sin( − 1)θ, implies that if χ = 2 cos θ, then a = sin θ is a solution of the recurrence relations in (112) for = 1, . . . , n. The boundary condition at = 0, a 0 = 0, is automatically satisfied, while the boundary condition at = n + 1, a n+1 = 0 fixes the admissible values of θ, and thus the eigenvalues of the original problem, namely sin(n + 1)θ = 0, that is θ (j) = jπ/(n + 1) with j = 1, . . . n.
Since the adjacency ∆ n is a symmetric matrix, with n real eigenvalues, we have obtained in this way its full spectrum: the eigenvalues χ (j) and the eigenvectors a (j) = (a where = 1, . . . , n, and N j is a suitable normalization coefficient. Plots in panels (a)-(b) are acoustic modes, interpolated by sine functions with frequency j/(n+1), while plots in panels (c)-(d) are optical modes, oscillating between a sine function with frequency (n + 1 − j)/(n + 1) and its opposite. The amplitude vectors a (j) are normalized to their square norms. Figure 10 shows the excitation waves a (j) for n = 100 emitters and selected values of j: while for small j the excitation profiles can be best understood as the sampling of a sinusoidal function with wavelength λ j , with in-phase consecutive amplitudes (acoustic waves), for high values of j the components of a (j) have alternating signs and values that lie alternatively on two opposite sinusoidal functions with period 2d(n + 1)/(n + 1 − j) (optical waves)