Bound states in the continuum for an array of quantum emitters

TL;DR

This paper investigates bound states in the continuum in a system of equally spaced quantum emitters coupled to a one-dimensional boson field, analyzing both approximate and exact solutions for small numbers of emitters.

Contribution

It provides a detailed analysis of bound states in the continuum for quantum emitter arrays, including explicit solutions for three and four emitters, considering both distant and finite spacing effects.

Findings

Degenerate eigenspaces of bound states identified for distant emitters.

Explicit solutions obtained for systems with three and four emitters.

Finite spacing and dispersion effects influence the bound state energies.

Abstract

We study the existence of bound states in the continuum for a system of n two-level quantum emitters, coupled with a one-dimensional boson field, in which a single excitation is shared among different components of the system. The emitters are fixed and equally spaced. We first consider the approximation of distant emitters, in which one can find degenerate eigenspaces of bound states corresponding to resonant values of energy, parametrized by a positive integer. We then consider the full form of the eigenvalue equation, in which the effects of the finite spacing and the field dispersion relation become relevant, yielding also nonperturbative effects. We explicitly solve the cases n=3 and n=4.

Figures

Tables

Table 1. Table 1: Critical values of the distance d c subscript 𝑑 𝑐 d_{c} at which the nonperturbative eigenvalue pair between the resonant energies E 1 subscript 𝐸 1 E_{1} and E 2 subscript 𝐸 2 E_{2} appears, and corresponding energy E c subscript 𝐸 𝑐 E_{c} , for arrays with different number of equally spaced emitters.

$n$	4	6	8	10	12
$d_c$	0.05	0.18	0.26	0.30	0.33
$E_c$	101	28	20	16	15

Equations

$$H_0=\epsilon \sum _{j=1}^n|e_j⟩⟨e_j|+\int \mathrm{d}k\omega (k)b^{†}(k)b(k).$$

$$\mathcal{N}=\sum _{j=1}^n|e_j⟩⟨e_j|+\int \mathrm{d}k{b}^{†}(k)b(k),$$

$$H_{\mathrm{int}}=\sum_{j=1}^{n}\int\mathrm{d}k\Bigl{[}F_{j}(k)|e_{j}\rangle\langle g_{j}|b(k)+\mathrm{H.c.}\Bigr{]},$$

$$|G^{(n)}⟩=\underset{j=1}{\overset{n}{⨂}}|g_j⟩$$

$$|{\Psi }_1⟩=\sum _{j=1}^n{a}_j|E_j^{(n)}⟩\otimes |\mathrm{vac}⟩+|G^{(n)}⟩\otimes \int \mathrm{d}k\xi (k)b^{†}(k)|\mathrm{vac}⟩,$$

$$|E_j^{(n)}⟩=|e_j⟩\underset{\ell \ne j}{⨂}|g_{\ell }⟩.$$

$$F_j(k)=\sqrt{\frac{\gamma }{2\pi \omega (k)}}{\mathrm{e}}^{\mathrm{i}(j-1)kd},$$

$$U_{m}b(k)U_m^{†}=\frac{1}{\sqrt{m}}b(k)=:\tilde{b}\left(\frac{k}{m}\right),$$

$$H(m,\epsilon ,d,\gamma )=m{U}_{m}H\left(1,\tilde{\epsilon },\tilde{d},\tilde{\gamma }\right)U_m^{†},$$

$$\tilde{\epsilon }=\frac{\epsilon }{m};\tilde{d}=md;\tilde{\gamma }=\frac{\gamma }{m^2}.$$

$$\left\{\begin{array}[]{l}\displaystyle(\varepsilon-E)a_{j}=-\sqrt{\frac{\gamma}{2\pi}}\int\mathrm{d}k\,\frac{\mathrm{e}^{-\mathrm{i}(j-1)kd}}{\sqrt[4]{k^{2}+1}}\xi(k),\\ \\ \displaystyle\left(\sqrt{k^{2}+1}-E\right)\xi(k)=-\sqrt{\frac{\gamma}{2\pi}}\sum_{l=1}^{n}a_{\ell}\frac{\mathrm{e}^{\mathrm{i}(\ell-1)kd}}{\sqrt[4]{k^{2}+1}}.\end{array}\right.$$

$$\xi (k)=-\sqrt{\frac{\gamma }{2\pi }}\sum _{\ell =1}^n{a}_{\ell }\frac{{\mathrm{e}}^{\mathrm{i}(\ell -1)kd}}{\sqrt[4]{k^2+1}\left(\sqrt{k^2+1}-E\right)},$$

$$\sum _{\ell =1}^n{a}_{\ell }{\mathrm{e}}^{\pm \mathrm{i}(\ell -1)d\sqrt{E^2-1}}=0.$$

$$(\epsilon -E)a_j=\frac{\gamma }{2\pi }\int \mathrm{d}k\frac{{\sum }_{l=1}^n{a}_{\ell }{\mathrm{e}}^{\mathrm{i}(\ell -j)kd}}{\sqrt{k^2+1}\left(\sqrt{k^2+1}-E\right)},$$

$${\mathrm{G}}^{-1}(E)\mathbf{a}=0,$$

$${\mathrm{G}}^{-1}(z)=(\epsilon -z)\text{openone}-\mathrm{\Sigma }(z),$$

$${\Sigma }_{j\ell }(z)=\frac{\gamma }{2\pi }\int \mathrm{d}k\frac{{\mathrm{e}}^{-\mathrm{i}(j-\ell )kd}}{\sqrt{k^2+1}\left(\sqrt{k^2+1}-z\right)}.$$

$$\lim_{\delta\downarrow 0}\bigl{[}\mathrm{\Sigma}(E+\mathrm{i}\delta)-\mathrm{\Sigma}(E-\mathrm{i}\delta)\bigr{]}\neq 0.$$

$$\det {\mathrm{G}}^{-1}(E)=0$$

$${\Sigma }_{jl}(E\pm \mathrm{i}0)=\frac{\pm \mathrm{i}\gamma }{\sqrt{E^2-1}}\left({\mathrm{e}}^{\pm \mathrm{i}|j-l|d\sqrt{E^2-1}}\pm \mathrm{i}{b}_{|j-l|}(E)\right),$$

$$b_j(E)=\frac{\sqrt{E^2-1}}{\pi }{\int }_1^{\infty }\frac{{\mathrm{e}}^{-j\lambda d}}{\sqrt{{\lambda }^2-1}}\frac{E}{E^2+{\lambda }^2-1}\mathrm{d}\lambda ;$$

$$b_0(E)=-\frac{1}{\pi }\log \left(E-\sqrt{E^2-1}\right).$$

$$\frac{|b_j(E)|}{|b_0(E)|}\le \exp (-jd)\text{for }E>1,$$

$$\xi (x)=$$

$$=$$

$${\xi }_1(x)=$$

$$=$$

$$\eta (x)=\frac{1}{2\pi }\sqrt{\frac{E^2-1}{2E}}{\int }_1^{\infty }\mathrm{d}\lambda \frac{{\mathrm{e}}^{-|x|\lambda }}{\sqrt[4]{{\lambda }^2-1}}\frac{\sqrt{{\lambda }^2-1}-E}{E^2+{\lambda }^2-1}$$

$${\left[A_n(\theta ,{\beta }_0,\mathbf{\beta })\right]}_{j\ell }={\mathrm{e}}^{\mathrm{i}|j-\ell |\theta }+\mathrm{i}{\beta }_{|j-\ell |},j,\ell =1,\dots ,n,$$

$$\mathrm{G}^{-1}(E)=-\frac{\mathrm{i}\gamma}{\sqrt{E^{2}-1}}A_{n}\big{(}\theta(E),\chi(E),\bm{b}(E)\big{)},$$

$$\theta (E)=$$

Full text

I Introduction

The physics of effectively 1-dimensional systems is recently attracting increasing attention, thanks to the unprecedented possibilities offered by modern quantum technologies. A number of interesting and versatile experimental platforms are available nowadays, to implement an effective dimensional reduction and enable photon propagation in 1D. These schemes differ in scope and make use of diverse physical systems, such as optical fibers onedim3 ; onedim4 , cold atoms focused1 ; focused2 ; focused3 , superconducting qubits onedim5 ; onedim6 ; mirror1 ; mirror2 ; atomrefl1 ; leo5 , photonic crystals kimble1 ; kimble2 ; onedim1 ; onedim2 ; ck , and quantum dots in photonic nanowires semiinfinite1 ; semiinfinite2 , the list being far from exhaustive. Light propagation in these systems is characterized by different energy dispersion relations and interaction form factors, yielding novel, drastically dimension-dependent features, that heavily affect dynamics, decay and propagation cirac1 ; cirac2 .

Although the physics of single quantum emitters in waveguides is well understood focused1 ; mirror2 ; boundstates1 ; lalumiere ; threelevel , novel phenomena arise when two refereeA1 ; refereeA2 ; PRA2016 ; oscillators ; baranger ; baranger2013 ; NJP ; yudson2014 ; laakso ; pichler ; Fedorov1 or more pichler2 ; bello ; bernien ; dong ; fang ; fang14 ; goban ; ck ; gu ; guimond ; lalumiere ; lodahl ; paulisch ; ramos14 ; ramos ; cirac1 ; boundstates1 ; tsoi ; yudsonPLA ; yudson2008 ; calajo15 emitters are present, since the dynamics is influenced by photon-mediated quantum correlations. In this and similar contexts, sub- and super-radiant states often emerge. However, while standard (Dicke) superradiance effects occur at light wavelength much larger than typical interatomic distances Dicke ; SRreview ; Kaiser1 ; Kaiser2 , considering wavelengths comparable to the interatomic distance brings to light a number of interesting quantum resonance effects.

In this article, we will apply the resolvent formalism cohentannoudji to study the existence of single-excitation bound states in the continuum in a system of $n$ quantum emitters. In these states, the excitation is shared in a stable way between the emitters and the field, even though the energy would be sufficient to yield photon propagation. The case $n=2$ has already been considered, both in the one- and two-excitation sectors PRA2016 ; PRA2018 . Here, we extend the results to general $n$, under the assumption of large interatomic spacing compared to the inverse infrared cutoff of the waveguide mode. We will then consider how the corrections to such approximation crucially affect the physical picture of the system, by explicitly analyzing the cases $n=3$ and $n=4$ and briefly reviewing $n=2$.

The paper is structured as follows. In Section II we introduce the physical system, the interaction Hamiltonian and the relevant parameters. In Section III we outline the general properties of bound states in the continuum. In Section IV we analyze and discuss the eigenvalues in the continuum and the corresponding eigenspaces. In Section V we comment on the existence of nonperturbative eigenstates, that emerge when the interatomic spacing is smaller than a critical value, depending on the number $n$. In Section VI we summarize the result and outline future research.

II Physical system and Hamiltonian

We shall consider a system of $n$ two-level emitters, equally spaced at a distance $d$ and characterized by the same excitation energy $\varepsilon$. Henceforth, we shall occasionally refer to the emitters as “atoms”. The ground and excited state of each emitter will be denoted by $|g_{j}\rangle$ and $|e_{j}\rangle$, respectively, with $j=1,\dots,N$. The emitter array is coupled to a structured one-dimensional bosonic continuum (e.g., a waveguide mode), characterized by a dispersion relation $\omega(k)\geq 0$, with $k\in\mathbb{R}$, and represented by the canonical field operators $b(k)$ and $b^{\dagger}(k)$, satisfying $[b(k),b^{\dagger}(k^{\prime})]=\delta(k-k^{\prime})$. In absence of interactions, the Hamiltonian of the system reads

[TABLE]

When the total Hamiltonian $H=H_{0}+H_{\mathrm{int}}$ is considered, the interacting dynamics generally does not preserve the total number of excitations

[TABLE]

unless a rotating-wave approximation is applied. In this case, the interaction Hamiltonian reads

[TABLE]

where $F_{j}(k)$ is the form factor describing the strength of the coupling of the $j$th emitter with a boson of momentum $k$, and $H$ can be diagonalized in orthogonal sectors characterized by a fixed eigenvalue of $\mathcal{N}$. The system is sketched in Fig. 1.

The zero-excitation sector is spanned by the single state $|G^{(n)}\rangle\otimes|\mathrm{vac}\rangle$, coinciding with the ground state of $H_{0}$, with

[TABLE]

and $|\mathrm{vac}\rangle$ satisfying $b(k)|\mathrm{vac}\rangle=0$ for all $k$’s. In this Article, we will focus on the possibility to find bound states in the one-excitation sector, in which the state vectors can be expanded as

[TABLE]

with

[TABLE]

In particular, we will consider a continuum with a massive boson dispersion relation $\omega(k)=\sqrt{k^{2}+m^{2}}$, characterized by the form factors

[TABLE]

determined by the $\bm{p}\cdot\bm{A}$ interaction of QED cohentannoudji , with $\gamma$ a coupling constant with the dimensions of squared energy.

The Hamiltonian $H=H_{0}+H_{\mathrm{int}}$, defined by the massive dispersion relation and by the form factors in Eq. (7), depends on the four parameters $\varepsilon$, $m$, $d$ and $\gamma$, all with physical dimension. However, it is easy to show that $H$ can be recast in a form in which only dimensionless combinations of such parameters appear. Define $U_{m}$ as the unitary operation that acts on the field operators as

[TABLE]

while acting trivially on the atomic sector. Then the following identity holds:

[TABLE]

with the (dimensionless) parameters in the right-hand side defined as

[TABLE]

In the following, since the spectra of the two Hamiltonians appearing in (9) are identical up to a factor $m$, we will focus on the properties of $H(1,\tilde{\varepsilon},\tilde{d},\tilde{\gamma})$, dropping the tilde from the dimensionless parameters and measuring momentum $k$ and energy $E$ in units of $m$.

III Bound states in the continuum

By considering the expressions (1), (3), and (7), that define the Hamiltonian, and the expansion (5) of the state vector, the eigenvalue equation in the one-excitation sector reads

[TABLE]

From the second equation

[TABLE]

one infers that, since $\xi(k)$ must be normalizable for a bound state, the vanishing of the denominator, occurring at $k=\pm\sqrt{E^{2}-1}$ for $E>1$, must be compensated by the vanishing of the numerator at the same points. Therefore, the atomic excitation amplitudes and the energy eigenvalue of bound states in the continuum necessarily satisfy the following constraint:

[TABLE]

By using the expression (12), one obtains the relation

[TABLE]

involving only the atomic excitation amplitudes and the eigenvalue $E$. The equation above can be expressed in the compact form

[TABLE]

with $\bm{a}=\{a_{j}\}_{1\leq j\leq n}$ and $\mathrm{G}^{-1}$ the inverse propagator matrix in the single-atomic-excitation subspace, generally defined for a complex energy $z$ by

[TABLE]

where the self-energy matrix $\mathrm{\Sigma}$ has elements

[TABLE]

The self-energy and the inverse propagator are well defined only for non-real arguments and on the real half-line $(-\infty,1)$, and are characterized by a discontinuity for $z=E\in[1,\infty)$, where generally

[TABLE]

Therefore, the coincidence of the two limits is a necessary condition for (15) to be well defined and, a fortiori, for $E$ to be an eigenvalue. Finally, notice that Eq. (15) always admits a trivial solution, which correspond, due to (12), to the null vector. If $\mathrm{G}^{-1}(E)$ is well defined, the equation

[TABLE]

provides a necessary and sufficient condition for $E$ to be an eigenvalue with a nontrivial solution $\bm{a}\neq\bm{0}$, providing the atomic excitation amplitudes of the corresponding eigenstate.

The integrals that define the elements of the self-energy in (17) can be evaluated by analytic continuation in the complex plane for $z=E\pm\mathrm{i}0$ and $E>1$, yielding

[TABLE]

with the first term derives from integration around one of the poles at $k=\pm\sqrt{z^{2}-1}$ and the second one

[TABLE]

from integration around one of the branch cuts of the analytic continuation. Notice that the $b_{j}$ functions are real for $E>1$. In the case $j=0$, the integral can be evaluated analytically and yields

[TABLE]

In the general case, the cut contribution must be evaluated numerically. However, a relevant property follows from the definition (21),

[TABLE]

implying that, for a sufficiently large spacing $d$, the terms $b_{j>0}$ can be neglected as a first approximation. In the following, we will show that, interestingly, the inclusion of such terms in the analysis on one hand entails selection rules that remove the degeneracy of bound states in the continuum, on the other hand displaces by orders $O(\mathrm{e}^{-d})$ the energies, resonance distances and amplitudes that satisfy the constraint in Eq. (13).

The boson (photon) eigenfunction (12) in the position representation reads

[TABLE]

with

[TABLE]

where

[TABLE]

is the $O(\mathrm{e}^{-x})$ cut contribution. Notice that the principal value prescription is required to define the integral appearing in $\xi_{1}$ for $E>1$, while the integral in $\xi$ is regularized by the constraint (13).

IV Eigenvalues and eigenstates

IV.1 Block-diagonal representation of the propagator

Given the form (15) of the eigenvalue equation for the atomic amplitude vector $\bm{a}$ and the dependence of the propagator on the inter-atomic distance $d$ and the transition energy $\varepsilon$, it is convenient to introduce the matrix $A_{n}(\theta,\beta_{0},\bm{\beta})$, with $\bm{\beta}=\{\beta_{p}\}_{1\leq p\leq n-1}$, depending on $n+1$ real parameters and defined as

[TABLE]

in terms of which the propagator reads

[TABLE]

with

[TABLE]

and $b_{j>0}(E)$ as defined in Eq. (21).

The matrix $A_{n}$ can be recast in a block-diagonal form by exploiting the invariance of the Hamiltonian with respect to spatial reflections around the midpoint between the first and $n$-th emitter, transforming the local basis $|E_{j}^{(n)}\rangle$ with the unitary transformation

[TABLE]

The action of such transformation, that is also real and symmetric, on the components in the local basis can be expressed for even $n=2h$ and odd $n=2h+1$ in terms of the $h\times h$ identity matrix $\openone_{h}$ and “exchange” matrix $J_{h}$ (i.e. the matrix with ones on the counterdiagonal as the only nonvanishing elements) as

[TABLE]

and

[TABLE]

respectively. The transformation $U_{n}$ generalizes the change from the local basis to the Bell basis for $n=2$ PRA2016 . In the new representation, the self-energy and the propagator are block diagonal:

[TABLE]

where $A^{-}_{n}(\theta,\chi,\bm{b})$ is the $\left\lfloor n/2\right\rfloor\times\left\lfloor n/2\right\rfloor$ matrix acting on the antisymmetric space, and $A^{+}_{n}(\theta,\chi,\bm{b})$ is the $\left\lceil n/2\right\rceil\times\left\lceil n/2\right\rceil$ matrix acting on the symmetric space of the qubits. Therefore, the eigenvalue equation (15) can be reduced to the quest for nontrivial solutions of the two decoupled linear systems

[TABLE]

Eigenvectors with indefinite reflection symmetry are allowed only if the same energy $E$ is an eigenvalue for both systems (35) for the same set of parameters $\varepsilon$, $d$ and $\gamma$. Examples of eigenstates with definite symmetry, whose relevance will be discussed in the following, are shown in Fig. 2.

Throughout this section, we will first analyze bound states by neglecting $O(\mathrm{e}^{-d})$ terms in the self-energy, and then discuss the consequences of including all the $b_{j>0}$ terms in the cases $n=2,3,4$.

IV.2 Large spacing approximation

When $dm$ is large, the terms $b_{j}$, with $j>0$, in the self-energy are exponentially suppressed and will be neglected as a first approximation, namely $\bm{b}=\bm{0}$. Both matrices $A_{n}^{\pm}(\theta,\chi,\bm{0})$ are singular if and only if $\theta=\nu\pi$, with $\nu\in\mathbb{N}$, and $\chi=0$. The former condition selects the possibile eigenvalues in terms of the spacing $d$

[TABLE]

which will be called resonant energies in the following the, while the latter condition

[TABLE]

provides a constraint involving the excitation energy, the spacing and the order $\nu$ of the resonance. Equation (37) defines a discrete family of curves in the $(\varepsilon,d)$ plane, identifying the values $\varepsilon$ for which a bound state in the continuum exists. The emitter configurations associated to the eigenvalues (36) satisfy different conditions, derived from the constraint (13), according to the parity of the resonance. For even $\nu$, for all the eigenvectors, the atomic excitation amplitudes must sum to zero

[TABLE]

while for odd $\nu$ one obtains

[TABLE]

Hence, each eigenvalue $E_{\nu}(d)$ is characterized by an $(n-1)$-fold degeneracy. It is worth observing that, since both matrices $A^{\pm}_{n}$ are characterized by the same singularity conditions at this level of approximation, the same eigenvalue can occur in both the symmetric and antisymmetric sector. In such cases, the eigenstates are not characterized by a well-defined symmetry.

The boson wavefunction associated to the eigenstates can be derived according to Eq. (III), considering $E=E_{\nu}(d)$. Neglecting the $\eta$ contribution in (III), the single-emitter contribution to the field is given by the oscillating function

[TABLE]

whose half-wavelength coincides with $d/\nu$. The boson wavefunction in the same approximation thus reads

[TABLE]

for even $\nu$, and

[TABLE]

for odd $\nu$. In both cases, due to the conditions (38) and (39), respectively, the field vanishes identically for $x<0$ and $x>nd$, and is therefore confined inside the emitter array.

Finally, it is worth observing that all possible $n$-emitter eigenstates can be obtained as linear combinations of two-emitter eigenstates at different positions. However, we will show in the following that $O(\mathrm{e}^{-d})$ effects, however small, remove this degeneracy, and imply selection rules related to the reflection symmetry of the atomic eigenstates.

IV.3 Full form of the self-energy

The degeneracy observed by approximating the self-energy as discussed in the previous subsection is lifted by considering the terms $b_{j}$, with $j>0$. We now discuss in detail this phenomenon. The effect of these terms can be summarized in the following points:

• i)

At given $d$ and $E_{\nu}(d)$, only one of the two matrices $A_{n}^{\pm}(\nu\pi,\chi(E_{\nu}(d)),\bm{b}(E_{\nu}(d)))$, namely the one for which

[TABLE]

continues to be singular for some values of $\varepsilon$ and $\gamma$. The matrix satisfying the property (43) is the antisymmetric one for odd $n$ and the one with symmetry $(-1)^{\nu+1}$ for even $n$. Details on this general result are given in the Appendix.

• ii)

The values of $\chi(E)$ (and hence of $\varepsilon$, through Eq. (37)) corresponding to the eigenstates with energy $E_{\nu}(d)$ will depend on the eigenstate. For any fixed $\varepsilon$, only one stable state with energy $E_{\nu}(d)$ can generally be found, with the orthogonal states becoming unstable (although they can be long-lived).

• iii)

If $A_{n}^{\pm}(\nu\pi,0,\bm{0})$ does not satisfy condition (43), then $A_{n}^{\pm}(\nu\pi,\chi(E_{\nu}(d)),\bm{b}(E_{\nu}(d)))$ is in general no longer singular. However, the corresponding stable states do not entirely disappear, but undergo a slight change in their amplitude and energy, which is now displaced with respect to $E_{\nu}(d)$. Such states must be studied numerically.

Here, we will explicitly examine these effects in the three cases $n=2,3,4$. Moreover we shall focus on eigenstates connected by continuity to the resonant bound states discussed in the previous subsection, postponing comments on the emergence of strong-coupling eigenstates, characterized by energies $E\gtrsim 10^{2}$ distant from the resonant values, to the remaining part of this Article.

IV.3.1 $n=2$

With respect to the inclusion of the cut terms in the self-energy, $n=2$ represents an oversimplified case, since the linear systems $A^{\pm}_{n}(\theta,\chi,\bm{b})$ reduce to single equations, and the singularity conditions read

[TABLE]

corresponding to eigenstates in which the emitter excitation amplitudes exactly satisfy

[TABLE]

The peculiarity of $n=2$ lies in the fact that the condition $\theta=\nu\pi$, with odd $\nu$ in the symmetric sector and even $\nu$ in the antisymmetric sector, still holds for both symmetries. Therefore, eigenvalues will be fixed by the condition $\chi=b_{1}$, that generalizes Eq. (36), and the constraint on the emitter excitation energy thus reads

[TABLE]

In this case, the inclusion of $b_{1}=O(\mathrm{e}^{-d})$ in the self-energy does not shift energies away from the resonant values and does not remove any degeneracy, since the symmetric and antisymmetric eigenstates already occurred for different $\nu$’s PRA2016 .

IV.3.2 $n=3$

For a system of three emitters, the eigenvalue equation breaks down into a single equation for the antisymmetric sector and a system of two equations in the symmetric case. In the former case, the eigenvalues are determined by the solution of

[TABLE]

As in the $n=2$ case, the real part of the above equation is sufficient to ensure that the resonance condition $\theta=\nu\pi$, here with any $\nu\in\mathbb{N}$ is still valid, and the corresponding energy must be in the form (36). The constraint on $\varepsilon$ and $d$ for the existence of an antisymmetric eigenstate, with the atomic excitation proportional to $(|E_{3}^{(1)}\rangle-|E_{3}^{(3)}\rangle)/\sqrt{2}$, is now determined by the equation $\chi(E)=b_{2}(E)$.

Instead, in the symmetric sector, where the eigenenergies are determined by the equation

[TABLE]

it is possible to directly check that, after imposing $\theta=\nu\pi$ with an integer $\nu$, one can find no solution, as their existence would imply at least one of the conditions $b_{2}(E)=\pm 3\sqrt{b_{1}(E)^{2}\pm 2b_{1}(E)}$. Actually, the energy of the symmetric bound state in the continuum

[TABLE]

is shifted by an amount of $O(\mathrm{e}^{-md})$ with respect to the resonant value $E_{\nu}(d)$, corresponding to a shift $\delta\theta\simeq(-1)^{\nu}b_{1}(E_{\nu}(d))$ in the phase. The values of $(\varepsilon,d)$ at which the symmetric bound states occur can now be derived from the condition

[TABLE]

with $E$ given by (49). For the lowest-order resonances $\nu=1$, one can observe that the energy of the symmetric state is shifted downwards with respect to the value $E_{1}(d)$, that is exact for the antisymmetric state. This effect is evident in Fig. 3, in which the behavior of the eigenvalues corresponding to bound states in the continuum for both parity sector is represented in terms of $d$. The trajectories of the bound states are displayed in Fig. 4.

While the excitation amplitudes of antisymmetric bound states are constrained to the values

[TABLE]

the amplitudes of the symmetric states depend on the parameters and on the magnitude of the cut contributions. If the terms $b_{j>0}$ are neglected, the symmetric bound state is characterized by

[TABLE]

with the second value sensitive to $O(\mathrm{e}^{-d})$ corrections when the $b_{j}$’s are included. These states were pictorially represented in the top panels (a)-(b) of Fig. 2, for relevant values of the parameters $d$ and $\gamma$. In the following section, we will find that bound states with different amplitudes, not connected by continuity to the ones described above, can emerge in the case $\varepsilon\gg 1$, a regime in which, however, the validity of the quasi-one-dimensional QED on which our model is based becomes questionable.

A relevant parameter that characterizes the features of bound states in the continuum is the total probability of atomic excitations

[TABLE]

that “measures” how the single excitation is shared between the emitters and the field. In this case, the probabilities $p^{(3)\pm}_{\nu}$ for the symmetric $(+)$ and antisymmetric $(-)$ eigenstates read

[TABLE]

up to order $O(\mathrm{e}^{-d})$. As we found in the case $n=2$ PRA2016 , the emitter excitation decreases with coupling and distance and increases with energy. In Fig. 5 we show the probabilities for the symmetric and antisymmetric states with $\nu=1$, computed from the approximate expressions (54)-(55) with varying $d$ and $\gamma$. In the whole parameter range, the approximate expressions provide, even for small $d$, a very good estimate of the exact values, which differ by less than $10^{-3}$ in the symmetric case and less than $2.5\times 10^{-2}$ in the antisymmetric case.

IV.3.3 $n=4$

For a system made of $n=4$ emitters, the eigenvalues in both symmetry sectors are determined by the singularity conditions of the $2\times 2$ matrices

[TABLE]

If the cut contributions are neglected, the singularity conditions yield $\theta=\nu\pi$ and $\chi=0$, and two complementary pictures emerge according to the resonance parity. For even $\nu$, the three-dimensional subspace corresponding to the eigenvalue $E_{\nu}(d)$ is spanned by the whole antisymmetric sector and by the symmetric state with

[TABLE]

For odd $\nu$, the eigenspace of $E_{\nu}(d)$ is still three-dimensional, spanned by the whole symmetric sector and by the antisymmetric state with

[TABLE]

When the $b_{j>0}$ terms are included, it is still possible to find eigenstates with resonant energy $E_{\nu}(d)$ in the antisymmetric sector for even $\nu$ and in the symmetric sector for odd $\nu$. In the former case, such states occur when the parameters $(\varepsilon,d,\gamma)$ satisfy

[TABLE]

The conditions derived from the two branches of the above equation, quadratic in $\chi$, yield the two eigenstates characterized, at the lowest order in $b_{j}$, by the amplitudes

[TABLE]

and the atomic excitation probabilities

[TABLE]

In the case of odd $\nu$, if the model parameters satisfy

[TABLE]

one finds symmetric eigenstates with $E=E_{\nu}(d)$, amplitudes

[TABLE]

and atomic excitation probabilities

[TABLE]

These are the states that were pictorially represented in the lower panels (c)-(f) of Fig. 2, for relevant values of the parameters $d$ and $\gamma$. The atomic probabilities of the four classes of eigenstates defined by Eqs. (60)-(63) are shown in Fig. 6.

The states defined by the amplitudes (57)-(58) persist as eigenstates even after the introduction of the cut integration terms. However, their energies and the ratios between local amplitudes are shifted by a quantity $O(\mathrm{e}^{-md})$ with respect to $E_{\nu}(d)$ and to the values in Eqs. (57)-(58), respectively. Specifically, at a fixed distance $d$, the antisymmetric state with amplitudes connected by continuity to (58) is characterized by an eigenvalue $E<E_{1}(d)$, slightly smaller than the resonant value. The total atomic probabilities corresponding to states in this class reads

[TABLE]

with even $\nu$, for the symmetric state, and

[TABLE]

with odd $\nu$, for the antisymmetric one.

The numerical analysis of the determinant of the matrices (IV.3.3) reveals the existence of a new class of nondegerate eigenstates, characterized, in the distance range $2\lesssim d\lesssim 6$, by the amplitudes

[TABLE]

with energy close to $E_{\nu}(d)$ for odd $\nu$, and

[TABLE]

with energy close to $E_{\nu}(d)$ for even $\nu$. The energy of such states is shifted with respect to the resonant values. In particular, one of the symmetric states (67) is characterized by an eigenvalue slighlty smaller than $E_{1}(d)$, which makes it the lowest-energy bound state in the continuum for a system of $n=4$ emitters at a fixed spacing $d$, as can be observed in Fig. 7. The states (67) and (68) are characterized by the values

[TABLE]

and

[TABLE]

of the emitter excitation probability, respectively, with $E_{\nu}(d)$ the closest resonant energy to the actual eigenvalue. The behavior of the lowest-energy bound states in the continuum is shown in detail in Fig. 8.

V Pair formation of high-energy eigenstates

Condition (19), which determines the eigenvalues of the system, is a complicated equation in $E$, featuring the functions $\theta(E)$, $\chi(E)$ and $\phi(E)$. In the previous section, we have analyzed the solutions that can be connected by continuity to the resonant energies (36) in the limit $\mathrm{e}^{-d}\to 0$. However, the non-polynomial character of Eq. (19) can generally gives rise to new solutions at finite $d$, which are unrelated to the resonant eigenvalues and eigenspaces. In particular, this phenomenon is facilitated for very small $d$ (in units $m^{-1}$), when the magnitude of all the $b_{j>0}$ is relevant and comparable to that of $b_{0}$, and expanding the equations for small $\mathrm{e}^{-d}$ becomes immaterial.

Figures 9 and 10 display general features of such nonperturbative states, for $n=3$ and $n=4$, respectively. These features are confirmed for higher $n$. At a sufficiently high value of the distance, all the eigenvalues are connected by continuity to $E_{\nu}(d)$, with $\nu\in\mathbb{Z}_{+}$. When distance decreases, additional eigenvalues start appearing in the $(E,d)$ plane, between $E_{\nu}(d)$ and $E_{\nu+1}(d)$, immediately branching in two distinct eigenvalues, whose energy increases when distance is further decreased. The observed processes of pair formation in the cases $n=3,4$ occur roughly at the same value of $d$. To quantify the range in which the phenomenon occurs we define the critical distance $d_{c}^{(n)}$ as the value which marks the appearence of the first eigenstate of this class between $E_{1}(d)$ and $E_{2}(d)$. We obtain the values $d_{c}=0.063$ for the $n=3$ system and $d_{c}=0.052$ for $n=4$. Notice that no state of this kind is observed with energy below $E_{1}(d)$. The value of energy $E_{c}$ corresponding to the critical distance is $E_{c}\simeq 79$ for $n=3$ and $E_{c}\simeq 101$ for $n=4$. Thus, independently of the values of the parameters $\varepsilon$ and $\gamma$, the energy of such states exceeds the mass $m$ by at least two orders of magnitude, an energy range in which the validity of our model, at least in a waveguide QED context, is far from being ensured. However, as one can observe from Tab. 1, the critical energy decreases to an order $10$ for larger systems.

The nonperturbative eigenvalues always correspond to symmetric eigenstates, in which the field is characterized by a central half-wavelength that is far from multiple integers of the interatomic spacing, as can be observed in both Figs. 9–10. From the expression (12) one infers that, in such high-energy states, the field wavefunction is suppressed and the single excitation is almost entirely shared by the emitters. Finally, we observe that, for $n>4$, we have found the existence of more than one pair of nonperturbative eigenstates between $E_{\nu}$ and $E_{\nu+1}$.

VI Conclusions

We have studied the existence and main features of bound states in the continuum for a multi-emitter system in a one-dimensional configuration. We have found that, remarkably, finite-spacing non-Markovian effects can break the degeneracies typical of the Markovian approximation, affecting eigenstates, eigenvalues and the physical model that features specific bound states. Future research will be devoted to the study of degeneracy breaking and the subsequent collective effects in systems with a large number of emitters.

Acknowledgments

PF, DL, SP, and DP are partially supported by Istituto Nazionale di Fisica Nucleare (INFN) through the project “QUANTUM”. FVP is supported by INFN through the project “PICS”. PF is partially supported by the Italian National Group of Mathematical Physics (GNFM-INdAM).

Appendix

Appendix A General properties of the eigenvalue equation

The method used to characterize resonant bound states for a system of $n$ emitters in the case of general $n$ is based on the decomposition (35) in decoupled parity sectors. In Section IV.2, we proved that, neglecting the $b_{j>0}$ terms, the eigenvalue equation reduces to $\chi(E)=0$, yielding $(n-1)$-times degenerate eigenvalues $E_{\nu}(d)$, with $\nu\in\mathbb{Z}_{+}$, corresponding to eigenvectors whose atomic excitation amplitudes are constrained by (38) or (39) according to the sign $(-1)^{\nu}$. Here, we prove that the resonant energies $E_{\nu}(d)$ persist as exact eigenvalues even after the introduction of cut integration terms, for some value of the excitation energy $\varepsilon$.

The reduction to a block-diagonal form provided by the transformations (32) and (33) enables one to recast the eigenvalue equation into the decoupled problems

[TABLE]

For definiteness, let us first consider the case of even $n=2h$. Let us introduce for convenience the quantities

[TABLE]

and the real and symmetric matrices

[TABLE]

[TABLE]

and $\mathcal{C}^{\pm}_{q}$ as the $q\times q$ matrix characterized by the elements

[TABLE]

If $\nu$ is even, then

[TABLE]

and

[TABLE]

while, for odd $\nu$,

[TABLE]

and

[TABLE]

Fixing $E=E_{\nu}(d)$ and considering the expression of $\chi(E)$, Eq. (71) can be generally recast in the form

[TABLE]

implying that $E_{\nu}(d)$ is an eigenvalue of the system if and only if $\varepsilon$ is the real eigenvalue of some matrix $\mathcal{M}$. From the expressions (76)-(78), one can notice that, in the antisymmetric sector for even $\nu$ and in the symmetric sector for odd $\nu$, the matrix $\mathcal{M}$ is Hermitian, entailing the existence of $n$ values of $\varepsilon$, real and generally distinct, corresponding to physical systems in which a bound state with energy $E_{\nu}(d)$ is present. Those values of $\varepsilon$ collapse to a single degenerate value in the $\mathrm{e}^{-d}\to 0$ limit. In the cases (77)-(79), insteads, $\mathcal{M}$ is not Hermitian, its the eigenvalues are generally no longer real, and the bound state energies displace from the resonant values.

The case of odd $n=2h+1$ is slighlty different. There, for all resonance orders $\nu$, in the antisymmetric sector

[TABLE]

leading to a condition (80) with a Hermitian $\mathcal{M}$, which implies that all the $E_{\nu}(d)$ are eigenvalues corresponding to antisymmetric bound states for generally different physical systems. On the other hand, the matrix $\mathcal{M}$ corresponding to all resonances in the symmetric sector is never Hermitian, since it features an imaginary and symmetric contribution proportional to $\mathcal{C}_{h+1}^{\pm}$.

Appendix B Unstable states

The resolvent formalism, employed in the main text to evaluate the existence and properties of bound states, also provides information on the lifetime of unstable states. The step required to perform this kind of analysis in the analytic continuation of the self-energy to the second Riemann sheet

[TABLE]

where $z$ is a complex energy. The lifetimes of unstable states are determined by the solutions $z_{p}=E_{p}-\mathrm{i}\gamma_{p}/2$ of the equation

[TABLE]

with

[TABLE]

We are now going to consider the properties of the complex poles of the propagator.

n=3

The block-diagonalization procedure applied to a system of three emitters implies the singularity conditions:

[TABLE]

for symmetric states and

[TABLE]

for antisymmetric states, with

[TABLE]

Introducing the functions $R_{3}(\theta,\bm{b})=\mathrm{Re}(f_{3}(\theta,\bm{b}))$, $S_{3}(\theta,\bm{b})=\mathrm{Im}(f_{3}(\theta,\bm{b}))$, the real and imaginary part of roots of the complex poles for the two blocks read

[TABLE]

The behavior of the complex poles of the propagator for $n=3$ is reported in panel (a) of Fig. 11.

n=4

The singularity condition for the symmetric and antisymmetric blocks in the $n=4$ system read

[TABLE]

respectively, with

[TABLE]

where we have defined $R_{4}^{\pm}(\theta,\bm{b})=\mathrm{Re}(f_{4}^{\pm}(\theta,\bm{b}))$, $S_{4}^{\pm}(\theta,\bm{b})=\mathrm{Im}(f_{4}^{\pm}(\theta,\bm{b}))$. In this way approximate decoupled solutions are

[TABLE]

The behavior of the complex poles of the propagator for $n=4$ in the symmetric and antisymmetric sectors is reported in panels (b)-(c) of Fig. 11.

We finally comment on the phenomenon of emergence of nonperturbative eigenstates in the low-spacing regime. Such poles appear when one of the complex poles with negative imaginary part in the second Riemann sheet approaches the real axis (see Fig. 12). Due to the analytic properties of the resolvent, this pole actually merges on the real axis with a pole of the analytic continuation

[TABLE]

in the upper half-plane. Further decreasing the spacing, the two poles split on the real axis and increase their energy difference.