Bulk and surface bound states in the continuum

The eigenvalue structure of a general quantum system with a finite potential well consists usually of bound states and unbound continuum states, whose energies are positive if we set the asymptotic value of the potential at infinity as zero. While the bound states are localized in space and square integrable, the continuum modes are extended and non-normalizable. In 1929, however, Wigner and von Neumann suggested an exception to this picture by constructing explicitly a bound state in the continuum (BIC), that is, an eigenstate with energy above the continuum threshold but localized and square integrable [1]. Their basic idea consisted of imposing a modulation of an otherwise sinusoidal profile with a decaying envelope that is square integrable. From that, a suitable local potential can be designed such that the proposed wavefunction is an eigenstate of this potential. Both the potential and the wavefunction thus proposed oscillate in space and decay as a power law. The idea was taken up again by Albeverio who showed that the use of a suitable potential could generate a discrete embedded spectrum [2]. A few years later, Stillinger [3] and Herrick [4] suggested that BICs might be found in certain atomic and molecular systems. Later on, they suggested the use of superlattices to construct potentials that could support BICs [5, 6]. Subsequent experiments with semiconductor heterostructures provided the direct observation of electronic bound states above a potential well localized by Bragg reflections [7]. A different approach for the design of potentials that can support BICs comes from the concept of resonant states in quantum mechanics. Resonant states are localized in space but with energies in the continuum, and they eventually decay, i.e., they have a finite lifetime. Under certain conditions, the interference between resonances can lead to a resonant state of zero width. In other words, the localized state decouples from the continuum becoming a BIC. One example of this is the case of a hydrogen atom in a uniform magnetic field, modeled as a system of coupled Coulombic channels, where interference between resonances belonging to different channels leads to the ocurrence of BICs [8, 9]. More recently, BICs have been shown to occur in mesoscopic electron transport and quantum waveguides [10–17], and in quantum dot systems [18–22]. Here, the existence of BICs can be traced back to the destruction of the discrete-continuum decay channels by quantum interference effects. Another recent mechanism generates optical bound states in the continuum in an optical waveguide array by decoupling from the continuum by means of symmetry only [23]. A very recent type of BICs has been proposed consistent with BICs of interacting particles [24–27].

On the other hand, it has been admitted that the BIC phenomenon relies on interference, and thus is inherent to any wave-like system, besides quantum mechanics, such as optical systems described by the paraxial wave equation. In fact, the analogy between these two realms has gained much attention recently, and has given rise to experimental observations of many phenomena that are hard to observe in a condensed matter setting [28]. Examples of this include dynamic localization [29], Bloch oscillations [30, 31], the Zeno effect [32, 33] and Anderson localization [34–37]. The main appeal of using optical systems is that experiments can be designed to focus on a particular aspect without the need to deal with the presence of many other effects commonly present in quantum solids, such as many-body effects. In optics it is also possible to steer and manage the propagation of excitations and to tailor the optical medium. Thus, it is not surprising that there have been also a number of recent works on BICs in classical optical systems [38–43]. Finally, from the purely mathematical point of view the topic is being actively investigated [44].

In this work we extend Wigner and von Neumannʼs original concept to demonstrate that BICs can exist as bulk modes or surface modes of a discrete linear lattice, even in the presence of first and second nearest neighbor interactions. We show explicitly that the wavefunction thus generated decays faster than a power law at long distances from the mode center, while the local potential decays as a power law. The system seems structurally stable, and thus an experimental verification of these ideas in the optical domain (waveguide arrays) looks promising.

Let us consider a linear, one-dimensional lattice in the presence of a site energy distribution $\{{{\epsilon }_{n}}\}$. In optics, this could correspond to a set of weakly coupled optical waveguides, each of them characterized by a propagation constant ${{\epsilon }_{n}}$ and centered at ${{x}_{n}}=na$. In the coupled-mode approach, we expand the electric field $E(x,z)$ as a superposition of the fundamental modes centered at each waveguide, $E(x,z)={{\sum }_{n}}{{C}_{n}}(z)\psi (x-na)$, where $\psi (x)$ is the waveguide mode. After positing ${{C}_{n}}(z)={{C}_{n}}{\rm exp} ({\rm i}\lambda z)$ and after inserting this into the paraxial wave equation, one obtains the stationary equations for the mode amplitudes

$$\begin{eqnarray}&&(-\lambda +{{\epsilon }_{n}}){{C}_{n}}+\mathop{\sum }\limits_{m\ne n}{{V}_{nm}}{{C}_{m}}=0,\end{eqnarray} \tag{ 1 }$$

where ${{V}_{nm}}$ is the coupling between the nth and mth waveguides. From equation (1), it is possible to formally express

$$\begin{eqnarray}&&{{\epsilon }_{n}}=\lambda -\mathop{\sum }\limits_{m\ne n}{{V}_{nm}}({{C}_{m}}/{{C}_{n}}).\end{eqnarray} \tag{ 2 }$$

Now, for a homogeneous system, we have ${{\epsilon }_{n}}=0$ and C_n and λ are computed by solving the simpler eigenvalue equation. For instance, for an infinite (and also semi-infinite) lattice with interactions to first nearest neighbors only, ${{C}_{n}}\equiv {{\phi }_{\,n}}\sim \{{\rm sin} (kn),{\rm cos} (kn)\}$ and ${{\lambda }_{\phi }}=2V{\rm cos} (k)$, while for the case of an infinite lattice with first and second nearest neighbors, ${{\phi }_{\,n}}\sim \{{\rm sin} (kn),{\rm cos} (kn)\}$ and ${{\lambda }_{\phi }}=2{{V}_{1}}{\rm cos} (k)+2{{V}_{2}}{\rm cos} (2k)$. For the case of a simultaneous presence of long range interactions and boundaries, a simple closed form for ${{\phi }_{\,n}}$ is not readily available, so it is better to resort to a numerical ${{\phi }_{\,n}}$. Hereafter, we will keep the general notation ${{\phi }_{\,n}}$ for the unmodulated wavefunction.

Let us now consider the problem of building a localized eigenmode but with eigenenergy inside the lattice of allowed energy bands. Following the prescription of Wigner and von Neumann [1] we select one of the eigenstates ${{\phi }_{\,n}}$, with eigenvalue ${{\lambda }_{\phi }}$, and proceed to modulate its envelope in such a way that the wave thus modulated is an eigenstate of the inhomogeneous system with eigenvalue ${{\lambda }_{\phi }}$:

$$\begin{eqnarray}&&C_{n}^{0}={{\phi }_{\,n}}{{f}_{n}},\end{eqnarray} \tag{ 3 }$$

where f_n is a decaying and normalizable envelope: ${{f}_{n}}\to 0$ as $n\to \infty$ and ${{\sum }_{n}}|\;{{f}_{n}}{{|}^{2}}\lt \infty$. In order for this state to be an eigenstate of the system, we need to introduce a 'local potential' ${{\epsilon }_{n}},$ whose shape will be adjusted to render C_n⁰ as an eigenstate.

After inserting this ansatz into equation (2), we obtain

$$\begin{eqnarray}&&{{\epsilon }_{n}}={{\lambda }_{\phi }}-\mathop{\sum }\limits_{m\ne n}{{V}_{nm}}\left( \frac{{{f}_{m}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{m}}}{{{\phi }_{\,n}}} \right).\end{eqnarray} \tag{ 4 }$$

Let us now apply these ideas to four simple systems. We begin with the case of an infinite lattice with interactions to first neighbors only. Equation (4) reads

$$\begin{eqnarray}&&{{\epsilon }_{n}}={{\lambda }_{\phi }}-V\left( \frac{{{f}_{n+1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}} \right)-V\left( \frac{{{f}_{n-1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right).\end{eqnarray} \tag{ 5 }$$

(As an aside note, this expression is reminiscent of the discrete version of the ${{\phi }^{4}}$ field theory of Peierls coupling. See also [45].) We take a monotonically decreasing envelope f_n around some site n₀. To the right of n₀ we have

$$\begin{eqnarray}&&\left( \frac{{{f}_{n+1}}}{{{f}_{n}}} \right)=(1-{{\delta }_{n}}),\end{eqnarray} \tag{ 6 }$$

where ${{\delta }_{n}}\lt 1$. From this, we can solve formally for f_n:

$$\begin{eqnarray}&&{{f}_{n}}={{\mathop{\prod }\limits_{m={{n}_{0}}}^{n-1}}^{{}}}(1-{{\delta }_{m}})n\gt n_{0}^{{}}\end{eqnarray} \tag{ 7 }$$

with ${{f}_{{{n}_{0}}}}\equiv 1$.

To the left of n₀ we have

$$\begin{eqnarray}&&\left( \frac{{{f}_{n}}}{{{f}_{n+1}}} \right)=(1-{{\delta }_{n}})\end{eqnarray} \tag{ 8 }$$

from which we obtain

$$\begin{eqnarray}&&{{f}_{n}}=\mathop{\prod }\limits_{m=n}^{{{n}_{0}}-1}(1-{{\delta }_{m}})\quad n\lt n_{0}^{{}}.\end{eqnarray} \tag{ 9 }$$

Now, we must choose an appropriate form for ${{\delta }_{n}}$ such that the local and decreasing 'potential' $\{{{\epsilon }_{n}}\}$ remains bounded. From equation (5) we see that possible problems will emerge at the zeroes of ${{\phi }_{\,n}}$. To avert that, we choose ${{\delta }_{n}}$ in the form

$$\begin{eqnarray}&&{{\delta }_{n}}=\frac{a}{1+|n-{{n}_{0}}{{|}^{b}}}{{N}^{2}}\phi _{\,n}^{2}\phi _{n+1}^{2},\end{eqnarray} \tag{ 10 }$$

where $n\geqslant 1$ and N is the length of the chain, ${{\sum }_{n}}|{{\phi }_{\,n}}{{|}^{2}}=1$. The presence of the lattice size N is due to the fact that $\phi =O(1/\sqrt{N})$, and we want $N|{{\phi }_{\,n}}{{|}^{2}}=O(1)$. Note that other choices for ${{\delta }_{n}}$ are possible, with the only condition that the envelope be localized.

Let us now show that, with the choices (7), (9) and (10), the site energy distribution ${{\epsilon }_{n}}$ will not diverge at the zeroes of ${{\phi }_{\,n}}$.

For $n\gt {{n}_{0}}$, and using equation (7), ${{\epsilon }_{n}}$ can be written as

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{n}} & = & {{\lambda }_{\phi }}-V(1-{{\delta }_{n}})\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}} \right)-V(\frac{1}{1-{{\delta }_{n-1}}})\left( \frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right). \\ \end{array}\end{eqnarray} \tag{ 11 }$$

Now after using equation (10) and going to the limit ${{\phi }_{\,n}}\to 0$, one obtains

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{n}} & = & {\overbrace {{{\lambda }_{\phi }}-V\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}}+\frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right)}^{0}}+\frac{Va}{1+|n-{{n}_{0}}{{|}^{b}}}{{N}^{2}}{{\phi }_{\,n}}\phi _{n+1}^{3} \\ {} & {} & -\frac{Va}{1+|n-1-{{n}_{0}}{{|}^{b}}}{{N}^{2}}{{\phi }_{\,n}}\phi _{n-1}^{3}\cdots \\ \end{array}\end{eqnarray} \tag{ 12 }$$

and thus ${{\epsilon }_{n}}\to 0$ when ${{\phi }_{\,n}}\to 0$, for $n\gt {{n}_{0}}$.

For the case $n\lt {{n}_{0}}$, and using equation (9), ${{\epsilon }_{n}}$ can be written as

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{n}} & = & {{\lambda }_{\phi }}-V\left( \frac{1}{1-{{\delta }_{n}}} \right)\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}} \right)-V(1-{{\delta }_{n-1}})\left( \frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right). \\ \end{array}\end{eqnarray} \tag{ 13 }$$

Now after using equation (10) and going to the limit ${{\phi }_{\,n}}\to 0$, one obtains

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{n}} & = & {\overbrace {{{\lambda }_{\phi }}-V\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}}+\frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right)}^{0}}-\frac{Va}{1+|n-{{n}_{0}}{{|}^{b}}}{{N}^{2}}{{\phi }_{\,n}}\phi _{n+1}^{3} \\ {} & {} & +\frac{Va}{1+|n-1-{{n}_{0}}{{|}^{b}}}{{N}^{2}}{{\phi }_{\,n}}\phi _{n-1}^{3}+\cdots \\ \end{array}\end{eqnarray} \tag{ 14 }$$

the same limit as in equation (12).

For $n={{n}_{0}}$, we have ${{f}_{{{n}_{0}}}}=1,{{f}_{{{n}_{0}}+1}}=1-{{\delta }_{{{n}_{0}}}}$ and ${{f}_{{{n}_{0}}-1}}=1-{{\delta }_{{{n}_{0}}-1}}$. Therefore

$$\begin{eqnarray}&&{{\epsilon }_{{{n}_{0}}}}={{\lambda }_{\phi }}-V(1-{{\delta }_{{{n}_{0}}}})\left( \frac{{{\phi }_{{{n}_{0}}+1}}}{{{\phi }_{{{n}_{0}}}}} \right)-V(1-{{\delta }_{{{n}_{0}}-1}})\left( \frac{{{\phi }_{{{n}_{0}}}}-1}{{{\phi }_{{{n}_{0}}}}} \right).\end{eqnarray} \tag{ 15 }$$

Using now equation (10) and going to the limit ${{\phi }_{\,n}}\to 0$, one obtains

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{{{n}_{0}}}} & = & {\overbrace {{{\lambda }_{\phi }}-V\left( \frac{{{\phi }_{{{n}_{0}}+1}}}{{{\phi }_{{{n}_{0}}}}}+\frac{{{\phi }_{{{n}_{0}}-1}}}{{{\phi }_{{{n}_{0}}}}} \right)}^{0}}+Va{{N}^{2}}{{\phi }_{{{n}_{0}}}}(\phi _{{{n}_{0}}+1}^{3}-\phi _{{{n}_{0}}-1}^{3}) \\ \end{array}\end{eqnarray} \tag{ 16 }$$

therefore if ${{\phi }_{{{n}_{0}}}}\to 0$, ${{\epsilon }_{{{n}_{0}}}}$ will approach zero. Thus we conclude that the value of site energy distributions at those points where ${{\phi }_{\,n}}$ is zero, is also zero and divergences are averted.

The particular choice of parameters a and b will determine the decay rate of both the wavefunction envelope and the potential. On one hand, for practical purposes we need a potential that decays relatively fast, but we also need this potential not to be so abrupt that it can be indistinguishable from an impurity potential, since in that case our candidate state could be pushed out of the band. Using equation (10), the asymptotic decay of the envelope at large n values can be estimated in closed form, using the Euler–Maclaurin formula, to be

$$\begin{eqnarray}&&{{f}_{n}}\sim {\rm exp} (-{{\alpha }_{n}}|n-{{n}_{0}}{{|}^{1-b}}),\quad n\to \infty ,\end{eqnarray} \tag{ 17 }$$

where ${{\alpha }_{n}}={{B}_{n}}a/(1-b)$ with B_n an oscillating factor of order $O(1)$ on average. Thus, for $0\lt b\lt 1$ equation (17) gives a faster decay rate than the power law obtained by Wigner and von Neumann [1]. On the other hand, based on equation (10), we can obtain the asymptotic behavior of the potential equation (5) to be:

$$\begin{eqnarray}&&{{\epsilon }_{n}}\sim \frac{{{A}_{n}}Va}{|n-{{n}_{0}}{{|}^{b}}},\quad n\to \infty ,\end{eqnarray} \tag{ 18 }$$

where A_n is an oscillating factor of order $O(1)$ on average. Therefore, the envelope of the oscillating potential decays in space as a power law, qualitatively similar to the Wigner–von Neumann case. From equations (17) and (18) we see that if we make our BIC more localized by choosing a small value for b, then the site energy distribution will decay more slowly, and vice versa. A useful parameter to quantify the degree of localization of a state is its participation ratio R, defined by $R\equiv {{({{\sum }_{n}}|{{C}_{n}}{{|}^{2}})}^{2}}/{{\sum }_{n}}|{{C}_{n}}{{|}^{4}}$. For localized modes, $R\approx 1$ while for extended states $R\approx N$, where N is the number of sites in the lattice. We will use R as a measure of localization of all eigenstates in the presence of the $\{{{\epsilon }_{n}}\}$. In general we would like to have only one, or a few BICs generated by our procedure, keeping the rest as extended modes. This is yet another condition to be imposed on our selection of $a,b$. As the numerical results will show, these requirements can be met more easily when the system has only interactions to first nearest neighbors, rather than longer range couplings.

The numerical tools we will employ are quite basic and consist of diagonalization of medium size matrices, and the evaluation of long products and sums, usually carried out with the help of Mathematica [46].

Figure 1 shows results for a lattice of N = 1065 sites and nearest neighbour interactions, where we have chosen $V=1,\lambda =1.565\;24$ well inside the linear spectrum and used the trial function (10) with $a=1/2$, $b=3/4$ and ${{n}_{0}}=532$. While the BIC mode decays in space rather quickly ($\sim {\rm exp} (-\alpha {{n}^{1/4}})$) (17), the site energy distribution decays at a slower pace (${{n}^{-3/4}}$) (18). Now, once in the possession of the distribution $\{{{\epsilon }_{n}}\}$, we proceed to find all the eigenvalues and eigenvectors of the array, but this time including the site energy distribution ${{\epsilon }_{n}}$, and plot their participation ratio R as a function of the state eigenvalue. That way one can see which states fall inside the continuum of scattered states and which fall outside the band, becoming impurity-like states. In this example we have 11 states that fall outside the band (some of them very close to the band). This is shown in figure 1(d). Our BIC is the one with the smallest R, namely $R\approx 6$ with the next higher one being $R\sim 100$. This means that we have a single BIC while the rest of the states are extended. Examination of all eigenvalues reveals that a small number of them (1%) fall outside of the band $[-2V,2V]$, and become impurity-like states. This is a typical behavior observed in all cases examined in this work. The states outside the band were created as a consequence of the relative sharpness of the modulating potential chosen, as shown in figure 1(c), where the middle part of ${{\epsilon }_{n}}$ looks like an impurity potential. An examination of the density of states $(1/N){{\sum }_{n}}\delta (\lambda -{{\lambda }_{n}})$ of the modulated system shows no change from the unmodulated case, except for the presence of some few states outside the band. This means our BIC is not inside some minigap, but is indeed inserted in a single band. Another good case is obtained for $a=1/2,b=0.9$ which gives rise to a well-defined single BIC, plus nine modes outside the band. For some values of a, b it is possible to have more than one BIC (e.g., for n = 1065, $a=1/2$, b = 0.25). In these cases, the number of impurity-like states outside the band is also observed to increase.

Zoom In Zoom Out Reset image size

Figure 2 shows the states inside the band that are closest in energy to the embedded state. The embedded state is the only state inside the band whose amplitude decreases to zero at large distance from the mode center, while all the rest of the band states are extended. The presence of the modulation distorts the band states but does not change their extended character. The two states closest in energy to the embedded state look quite similar, at least at the scale of the graph, but a close examination reveals that they are not identical, and could not be since they are orthogonal. An examination of all modes revealed that the presence of the embedded mode seems to divide the other modes into two parts: to the left and to the right of the position of the center of the embedded mode. On each portion the mode oscillates differently.

Zoom In Zoom Out Reset image size

Let us now consider the case of a surface BIC, i.e., one where the BIC is centered at one of the ends of the lattice, say at n = 1. This case was recently considered by us [43], but here we include it for completeness and also give a more general treatment. The general approach is the same as in the bulk case, but with an envelope that decreases away from the surface

$$\begin{eqnarray}&&{{f}_{n}}=\mathop{\prod }\limits_{m=1}^{n-1}(1-{{\delta }_{m}})\end{eqnarray} \tag{ 19 }$$

for $n\gt 1$, while ${{f}_{1}}=1$, and where we choose

$$\begin{eqnarray}&&{{\delta }_{n}}=\frac{a}{{{n}^{b}}}{{N}^{2}}\phi _{\,n}^{2}\phi _{n+1}^{2}.\end{eqnarray} \tag{ 20 }$$

The site energy distribution is now given by

$$\begin{eqnarray}&&{{\epsilon }_{1}}={{\lambda }_{\phi }}-V\left( \frac{{{f}_{2}}}{{{f}_{1}}} \right)\left( \frac{{{\phi }_{2}}}{{{\phi }_{1}}} \right)\end{eqnarray} \tag{ 21 }$$

at the surface, and

$$\begin{eqnarray}&&{{\epsilon }_{n}}={{\lambda }_{\phi }}-V\left( \frac{{{f}_{n+1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}} \right)-V\left( \frac{{{f}_{n-1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right)\end{eqnarray} \tag{ 22 }$$

for $1\lt n\lt N$, and

$$\begin{eqnarray}&&{{\epsilon }_{N}}={{\lambda }_{\phi }}-V\left( \frac{{{f}_{N-1}}}{{{f}_{N}}} \right)\left( \frac{{{\phi }_{N-1}}}{{{\phi }_{N}}} \right)\end{eqnarray} \tag{ 23 }$$

at n = N. As in the previous case, different choices of a and b lead to different decay rates for the wavefunction ${{C}_{n}}={{\phi }_{\,n}}{{f}_{n}}$ and its associated local potential ${{\epsilon }_{n}}$. Figure 3 shows a good trial wavefunction obtained with $a=2/5,b=9/10$. As in the bulk case, the local potential distorts most of the band modes, without altering their extended character. A small percentage of the initial modes get pushed outside the band, becoming impurity-like states, whose spatial extent depends on their proximity to the band. In the case of figure 3, we have a single BIC, with $R\sim 3$ while the next higher R is about 213 and only three of the modes were pushed out of the band (two of them very close to the band). Another good case is obtained for $a=0.5,b=0.9$, where we also have a single BIC, and there are only four states outside the band. It is also possible to have more than one BIC. For instance, for $a=1/3,b=1/2$ we have two BICs of similar participation ratios R. In this case, there are 13 states outside the band. These features have also been observed in a recent work on surface BICs, using a modulation of the couplings [47]. Sometimes, the extra BICs thus produced are not quite localized, but consist of a well-localized central peak surrounded by a very small but delocalized tail. They constitute resonance-like states rather than BICs. This is the case observed for $a=1/3,b=1/4$. The number of states pushed out of the band could in principle be reduced by taking a smoother ${{\epsilon }_{n}}$; the price to pay would be a local potential with a more extended range.

Zoom In Zoom Out Reset image size

Let us now come back to the bulk BIC, but this time including interactions to first and second nearest neighbors. The site energy distribution is now given by

$$\begin{eqnarray}&&{{\epsilon }_{1}}={{\lambda }_{\phi }}-{{V}_{1}}\left( \frac{{{f}_{2}}}{{{f}_{1}}} \right)\left( \frac{{{\phi }_{2}}}{{{\phi }_{1}}} \right)-{{V}_{2}}\left( \frac{{{f}_{3}}}{{{f}_{1}}} \right)\left( \frac{{{\phi }_{3}}}{{{\phi }_{1}}} \right)\end{eqnarray} \tag{ 24 }$$

for n = 1,

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{2}} & = & {{\lambda }_{\phi }}-{{V}_{1}}\left( \frac{{{f}_{1}}}{{{f}_{2}}} \right)\left( \frac{{{\phi }_{1}}}{{{\phi }_{2}}} \right)-{{V}_{1}}\left( \frac{{{f}_{3}}}{{{f}_{2}}} \right)\left( \frac{{{\phi }_{3}}}{{{\phi }_{2}}} \right)-{{V}_{2}}\left( \frac{{{f}_{4}}}{{{f}_{2}}} \right)\left( \frac{{{\phi }_{4}}}{{{\phi }_{2}}} \right) \\ \end{array}\end{eqnarray} \tag{ 25 }$$

for n = 2,

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{n}} & = & {{\lambda }_{\phi }}-{{V}_{1}}\left( \frac{{{f}_{n+1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n+1}}}{{{\phi }_{\,n}}} \right)-{{V}_{1}}\left( \frac{{{f}_{n-1}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n-1}}}{{{\phi }_{\,n}}} \right) \\ {} & {} & -{{V}_{2}}\left( \frac{{{f}_{n+2}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n+2}}}{{{\phi }_{\,n}}} \right)-{{V}_{2}}\left( \frac{{{f}_{n-2}}}{{{f}_{n}}} \right)\left( \frac{{{\phi }_{n-2}}}{{{\phi }_{\,n}}} \right) \\ \end{array}\end{eqnarray} \tag{ 26 }$$

for $2\lt n\lt N-1$,

$$\begin{eqnarray}\begin{array}{rcl} {{\epsilon }_{N-1}} & = & {{\lambda }_{\phi }}-{{V}_{1}}\left( \frac{{{f}_{N}}}{{{f}_{N-1}}} \right)\left( \frac{{{\phi }_{N}}}{{{\phi }_{N-1}}} \right)-{{V}_{1}}\left( \frac{{{f}_{N-2}}}{{{f}_{N-1}}} \right)\left( \frac{{{\phi }_{N-2}}}{{{\phi }_{N-1}}} \right)-{{V}_{2}}\left( \frac{{{f}_{N-3}}}{{{f}_{N-1}}} \right)\left( \frac{{{\phi }_{N-3}}}{{{\phi }_{N-1}}} \right) \\ \end{array}\end{eqnarray} \tag{ 27 }$$

for $n=N-1$, and

$$\begin{eqnarray}&&{{\epsilon }_{N}}={{\lambda }_{\phi }}-{{V}_{1}}\left( \frac{{{f}_{N-1}}}{{{f}_{N}}} \right)\left( \frac{{{\phi }_{N-1}}}{{{\phi }_{N}}} \right)-{{V}_{2}}\left( \frac{{{f}_{N-2}}}{{{f}_{N}}} \right)\left( \frac{{{\phi }_{N-2}}}{{{\phi }_{N}}} \right)\end{eqnarray} \tag{ 28 }$$

for n = N, the end of the chain.

Parameters V₁ and V₂ denote the couplings to first and second nearest neighbors, respectively.

The envelope f_n is still given by equations (7) and (9), but to avoid possible divergences and to keep ${{\epsilon }_{n}}$ bounded, we now choose ${{\delta }_{n}}$ in the form

$$\begin{eqnarray}&&{{\delta }_{n}}=\frac{a}{1+|n-{{n}_{0}}{{|}^{b}}}{{N}^{4}}\phi _{n-1}^{2}\phi _{\,n}^{2}\phi _{n+1}^{2}\phi _{n+2}^{2}\end{eqnarray} \tag{ 29 }$$

for $1\lt n\leqslant N-2$. For n = 1 we choose

$$\begin{eqnarray}&&{{\delta }_{1}}=\frac{aN}{1+|1-{{n}_{0}}{{|}^{b}}}\phi _{1}^{2}\end{eqnarray} \tag{ 30 }$$

while at the other extreme, we choose

$$\begin{eqnarray}&&{{\delta }_{N-1}}={{\delta }_{N}}=\frac{a{{N}^{2}}}{1+|N-{{n}_{0}}{{|}^{b}}}\phi _{N-1}^{2}\phi _{N}^{2},\end{eqnarray} \tag{ 31 }$$

where as before, ${{\phi }_{\,n}}$ is the (normalized to unity) extended mode we want to modulate (with eigenvalue ${{\lambda }_{\phi }}$) and n₀ is the position of the center of the localized state. Now, in spite of the more complex form of ${{\delta }_{n}}$, the asymptotic form for C_n and ${{\epsilon }_{n}}$ is the same as in the case of coupling to nearest neighbors only. This is due to the extended character of ${{\phi }_{\,n}}$ which does not affect the form of the decaying envelope of ${{\delta }_{n}}$ in equations (29)–(31) but only produces oscillations. Results are shown in figure 4 for $N=1065,\lambda =1.693\;57,a=3.5,b=0.99$. In this example we take ${{V}_{1}}=1={{V}_{2}}$ which corresponds to a geometrical arrangement where the guides are forming a 'zigzag' structure. The presence of the extra coupling seems to increase the localization length of the BIC as compared to the case with only nearest neighbor coupling. This seems natural since an extra coupling would favor an increase of the spatial extent of the wavefunction. In this case we also seem to have several BICs, although only the one we built has a well-defined decaying rate (slow though); the other ones consist of well-localized centers with non-decaying tails and seem to constitute resonances rather than bona fide BICs. This time 20 states were pushed outside of the band. We attribute this increase with respect to the nearest-neighbor case to the steep shape of ${{\epsilon }_{n}}$, which can be visualized as a high impurity-like potential surrounded by a much smaller but bristlier potential (see figure 4(c)). The combination of these two factors gives rise to a larger number of impurity-like states outside the band.

Zoom In Zoom Out Reset image size

Finally, let us consider the case of a surface BIC in the presence of first and second nearest-neighbor interactions. In this case it is not possible to obtain a simple, closed-form expression for the unmodulated eigenfunction ${{\phi }_{\,n}}$. The presence of interactions longer than nearest neighbors complicates the use of the image method technique [48]. Thus, ${{\phi }_{\,n}}$ and λ have to be extracted directly from the initial diagonalization. The site energy distribution is the same as for the bulk case with first and second nearest interactions, equations (24)–(28). The envelope function is given by equation (19). The difference now lies in the decay function ${{\delta }_{n}}$:

$$\begin{eqnarray}&&{{\delta }_{n}}=\frac{a}{{{n}^{b}}}{{N}^{2}}\phi _{n-1}^{2}\phi _{\,n}^{2}\phi _{n+1}^{2}\phi _{n+2}^{2}\end{eqnarray} \tag{ 32 }$$

for $n\gt 2$, where we define ${{\delta }_{1}}=0$. With this choice, ${{\epsilon }_{n}}\to 0$ at those sites where ${{\phi }_{\,n}}\to 0$. The envelope function is then defined as in equation (19).

Figure 5 shows results for the case $N=533,{{V}_{1}}=1,{{V}_{2}}=0.5,\lambda =1.5655,a=5/2,b=3/4$. We notice that the site potential is more 'bristling' than in the case of nearest neighbors only. This lack of smoothness plus its relatively high value at the surface causes that more states are now pushed outside the band (30 in the case of figure 5). From figure 5(d) we notice about five states that could correspond to BICs. A close examination of them reveals that only three of them constitute true BICs while the other two are resonances. The rest of the modes with larger R are extended modes.

Zoom In Zoom Out Reset image size

Finally, we performed a numerical structural stability analysis of each modulated system by adding some random noise to the site energy distribution: ${{\epsilon }_{n}}\to (1+{{\xi }_{n}}){{\epsilon }_{n}}$, where ${{\xi }_{n}}\in [-0.05,0.05]$ is a random quantity. Some results are shown in figure 6 for the cases of BICs with nearest-neighbor coupling. In all cases examined it was observed that the perturbation did not affect the envelope of the BIC appreciably, although it had an effect on the rest of the modes, especially when the BIC was a bulk mode. In this case, the rest of the modes showed a general tendency towards localization, as evidenced by an overall decrease of their participation ratio (figure 6(b)). However, when the BIC was a surface mode, the effects were quite small (not discernible at the scale of figure 6(d)). Of course, in all cases, the position of the BICs in the band was slightly shifted, depending on the random realization. Also, the range of the interaction did not seem to play a significant role. We conclude that the BIC is structurally stable against small perturbations. This stability could be important when one tries to carry out an experimental realization of the systems described here; small errors in the building of the local potential ${{\epsilon }_{n}}$ are unavoidable, and thus the stability of the system gives us hope that an experimental realization (perhaps in optics) can be soon carried out.

Zoom In Zoom Out Reset image size

We have examined bound states whose eigenenergies are embedded in the continuous spectrum of a one-dimensional periodic lattice. We have demonstrated an explicit procedure to generate square-integrable, bulk and surface localized modes embedded in the continuum, considering interactions to first and second nearest neighbors. We have also given a prescription to generate the local bounded potential that gives rise to such modes. The envelope of these modes can be chosen to decrease in space faster than a power law, although the local potential decreases as a power law. The procedure of BIC generation can give rise to other BICs, as well as to resonances, and impurity modes located outside the continuous band. The BICs thus constructed are structurally stable against small perturbations.

This work was supported in part by Fondo Nacional de Ciencia y Tecnología (Grant 1120123), Programa Iniciativa Científica Milenio (Grant P10-030-F), and Programa de Financiamiento Basal (Grants FB0824 and FB0807).