Abstract
Some photonic systems support bound states in the continuum (BICs) that have infinite lifetimes, although their frequencies and momenta are matched to vacuum modes. Using a prototypical system that can be treated analytically, we show that each of these BICs always splits into a pair of new type BIC and lasing threshold mode when a parity-time (PT)–symmetric perturbation is introduced. The radiation loss at the lasing threshold is exactly balanced by the net gain of the particles. These PT symmetry-induced BICs are different from ordinary BICs, as they can be excited by an external source but do not radiate, and they carry a different quality factor divergence rate from that of the ordinary BICs. While most of the attention of PT-symmetric systems is captured by the coalescence of modes at exceptional points, the splitting of ordinary BICs is a new phenomenon that illustrates the rich physics embedded in PT-symmetric systems.
INTRODUCTION
Bound states in the continuum (BICs) are a special kind of spatially localized states with the corresponding frequencies lying in the continuous spectrum. The existence of such states was first demonstrated in custom-constructed potentials in an electronic system (1) and later revealed to be a general phenomenon for both quantum and classical waves (2–17). Various mechanisms have been proposed to produce BICs (7–13). For periodic systems, they can generally be classified into two types: symmetry-protected BICs and propagating BICs. The symmetry-protected BICs are a result of their different symmetries with the other states in the environment so that no coupling can exist. A good example is the BIC appearing at the Brillouin zone (BZ) center of a periodic system (4). However, BICs can also appear at k ≠ 0 when the system has certain symmetries, which allow the far-field radiation to vanish at some special point on the leaky channel formed by the resonant modes inside the light cone. Such BICs are called propagating BICs and can be identified by the divergence of the quality (Q) factor through parameter tuning. The existence of propagating BICs is conventionally believed to be a result of field cancellation in the far field due to destructive interference of multiple resonant modes at the same frequency. The BICs in certain two-dimensional (2D) photonic crystal slabs carry some topological characteristics as they are the vortex centers of the polarization directions of far-field radiation (13–15).
On the other hand, systems with parity-time (PT) symmetry have drawn considerable interest recently (18–28). Such systems are non-Hermitian and can exhibit a PT phase transition at an exceptional point where the system’s eigenvalues change from real to complex pairs (22). The presence of exceptional points gives rise to many interesting phenomena, such as asymmetric power oscillation (23), unidirectional propagation (24, 25), Bloch oscillation (26), single-mode lasing (27), and rings of exceptional points (28).
The introduction of loss/gain into a system usually coalesces modes, characterized by the merging of two modes into one mode at the exceptional point. In this work, we show that PT-symmetric perturbation can actually induce a splitting of modes. We use a prototypical system that has analytic solutions to study the effects of PT-symmetric perturbation on BICs. More specifically, we first consider a periodic dimerized chain of identical plasmonic nanoparticles. By using the coupled dipole equations, we first study the phase diagram for the existence of propagating BICs in such a system without the PT-symmetric perturbation. Using analytical lattice sums, we are able to derive the condition for the existence of propagating BICs in such a system without resorting to a numerical search for a divergent Q factor. As a result, we find a well-defined surface in a 3D space spanned by material and geometrical parameters that separate one region in which propagating BICs are guaranteed to exist from another region in which they do not exist. Under a PT-symmetric perturbation, where equal amounts of gain and loss are added to two adjacent particles, we show rigorously that a BIC, which can be either symmetry protected or propagating, always splits into two modes along the leaky channel. One mode is a new type of BIC (which we call pt-BIC), and the other is a mode at the lasing threshold. Lasing can occur in the region between the pt-BIC and the lasing threshold. A remarkable property of the pt-BIC mode found here is that it can be excited by an external field but does not radiate itself. The appearance of lasing threshold is a result of the exact energy balance between the radiation loss and net gain obtained from the particles, which leads to a divergent Q factor. The divergence rate of the Q factor for the pt-BICs is half of that of the ordinary ones.
If the original periodic system has propagating BICs, the PT-symmetric perturbation can create two pt-BICs on one side of the BZ and two accompanying lasing thresholds on the other side of the BZ. Unlike the unstable lasing region found between a pt-BIC and a lasing threshold, the region between two pt-BICs or two lasing thresholds is a stable passive region, in which the net gain of a mode is less than its radiative loss. Thus, the presence of pt-BICs and lasing thresholds divides the system into regions of stable resonant modes and lasing modes. When the amount of gain and loss is increased to a critical value, the two pt-BICs on the same side of the BZ can be annihilated. Above the critical value, both pt-BICs disappear suddenly together with the stable region enclosed. The two lasing threshold modes on the other side of the BZ also simultaneously display the same behavior. These novel phenomena are numerically verified by using COMSOL Multiphysics.
It should be pointed out that the creation of pairs of pt-BICs and lasing threshold modes under a PT-symmetric perturbation, as well as the associated novel phenomena mentioned above, is guaranteed in the system of dimerized chain, independent of the material and geometrical parameters. It is also a generic phenomenon for other 1D systems with BICs. As an example, we have also applied the PT-symmetric perturbation to a system of photonic crystal fiber and find the splitting of a BIC into a pair of pt-BIC and lasing threshold mode with a lasing region in between. The experimental considerations are also discussed. We note that different types of BICs have been found in systems with PT symmetry in both broken-PT and unbroken-PT phases (29–31). These BICs are different from the pt-BIC found here, which originates from an ordinary BIC that already exists before applying a PT-symmetric perturbation. The novel phenomena found here manifest the new physics arising from the interplay between PT-symmetric perturbation and BICs, which is completely different from the known physics of PT symmetry associated with exceptional points.
RESULTS
Dimerized chain and phase diagram for ordinary BICs
The system of an infinite PT-symmetric dimerized periodic chain with period a is schematically shown in the upper inset of Fig. 1. In the absence of gain and loss, particles A and B in a unit cell are two identical spherical plasmonic nanoparticles of radius r0 separated by a distance d. Their relative permittivity is described by the Drude model,
Upper inset: Schematic of a dimerized chain. Gain and loss can be introduced into particles A and B, respectively. In the absence of gain and loss, the real and imaginary parts of the resonance frequency ω are plotted respectively in (A) and (B). The BICs are marked by circles. The real and imaginary parts of the function f(q,ω) defined in Eq. 3 are shown in (C) and (D), respectively. (E) The four dashed lines denote the nodal lines of Re(f) at different values of dimensionless plasma frequency
For a periodic chain along the z axis, a resonant state can be labeled by a Bloch momentum q inside the BZ [−π/a,π/a]. When the system is excited by an external plane wave of the form
To study the resonant modes of the system, we set E = 0 in Eq. 1 and look for those complex frequencies at which the matrix M has a zero eigenvalue, i.e., MP = 0. Thus, the resonant frequencies of the system can be determined from the condition det(M) = 0. It is convenient to express the polarization vector in the form
In this work, we choose the geometric parameters as r0 = 20 nm, d = 70 nm, and a = 190 nm. For a particular material parameter ωp = 0.95(2πc/a), the results of ω′(q) and ω″(q) are shown in Fig. 1 (A and B), respectively. Here, we only focus on the zeroth-order diffraction region, i.e., |q| < ω/c < 2π/a − |q|. From the expression of det(Μ′) and the symmetry relations S00(−q,ω) = S00(q,ω) and Sij(−q,ω) = Sji(q,ω), it is easy to observe the symmetry of the resonant frequency, i.e., ω(q) = ω( − q), as shown in Fig. 1 (A and B). The three zeros of ω″(q) indicate the existence of one symmetry-protected BIC at q = 0 and a pair of propagating BICs at some nonzero ±q.
The existence of BICs can also be explicitly determined via the following analytical approach. For a BIC with a real-frequency solution of Eq. 2, the nonradiative condition in the transverse plane requires a complete cancellation of far-field radiation, which in turn requires
However, the condition h(q,ω) = 0 is not independent of the condition f (q,ω) = 0 in the determination of BICs. The proof is given in section S1. Briefly, we show that the inversion symmetry of the system gives an identity f (−q,ω) = h(q,ω), while the time-reversal symmetry gives f *(q,ω) = f (−q,ω). Thus, we have h(q,ω) = f *(q,ω), which means that the real-frequency solutions of f (q,ω) = 0 are also the solutions of h(q,ω) = 0. Therefore, at BICs, we also have SABeiqd = SBAe-iqd. This makes the matrix Μ′ symmetric, which in turn makes the condition
In Fig. 1 (C and D), we plot the real and imaginary parts of the function f(q,ω) in the q-ω plane. It can be seen that Re(f) and Im(f) are, respectively, an even and an odd function of q, consistent with the identity f *(q,ω) = f (−q,ω) obtained previously at BICs. This identity holds in general in the q-ω space. A rigorous mathematical proof is given in section S6. Thus, the identity h(q,ω) = f *(q,ω) also holds in general. In Fig. 1E, the nodal lines of Re(f) are plotted with dashed lines at different values of dimensionless plasma frequency,
Imaginary parts of the matrix elements of Μ′ + iγσz as well as the functions Im(f) + γ and Im(h) − γ evaluated at the frequency of Re(f) = 0 are shown in (A) to (C) for the PT-symmetric parameters γ = 0,
Because
Pairing of pt-BIC and lasing threshold mode
When gain and loss are added to the system, Eq. 2 becomes
(A) The solid curves denote the nodal lines of Im(f) + γ at different values of the PT-symmetric parameter γ. When 0 < γ < γc, their intersection points with the nodal line of Re(f) (dashed curve) at
Figure 3A also shows the evolution of pt-BICs as a function of γ. As the value of γ is increased, the two pt-BICs in the positive q space move toward each other along the nodal line of Re(f), and finally, they are annihilated when γ = γc. The system is left with only one pt-BIC in the negative q space when γ > γc. Thus, it is interesting to see that the trajectory of pt-BICs follows the nodal line of Re(f) when γ is varied. On the contrary, the trajectory follows the nodal line of Im(f) + γ when the material parameter
Accompanying each pt-BIC, a mode at the lasing threshold emerges at the same frequency but in the opposite q space. The existence of such a lasing threshold mode can also be derived analytically from Eq. 4. It is easy to show that a real-frequency solution of Eq. 4 at a point (q0,ω0) with
In addition to the above analytical study, we have also performed numerical studies using COMSOL Multiphysics to obtain the quasi-eigenmodes and complex eigenfrequencies ω′(q) + iω″(q) under a scattering boundary condition. For the case of
The function ω″(q) has zeros at three pairs of discrete q points, i.e., q0 = ±0.090, ±0.400, and ± 0.442, indicated by colored polygons in Fig. 3B. By using COMSOL, we numerically simulate the cross-sectional field profiles of |E| for these three pairs of singularities and show the results in Fig. 3D. Taking the pair split from the original BIC at the Γ point as an example, the two states at q0 = ±0.090 have completely different eigenstate profiles. The E field is spatially localized in the transverse plane for the state at q0 = 0.090 without far-field radiation, indicating a pt-BIC (black pentagon). However, an extended state is observed at q0 = −0.090, indicating a lasing threshold mode (red pentagon) having far-field radiation, in agreement with the analytical prediction. Similar conclusions can be drawn for the other pairs of states marked by the triangles and squares in Fig. 3 (B and D). A quantitative analysis of these field patterns in terms of the powers of net gain/loss and radiation loss will be given in the next subsection.
It is important to note that the slope of ω″(q) at the zeros of ω″(q) is not zero as shown in Fig. 3C, unlike the slope at the zeros of ω″(q), which is also zero shown in Fig. 1B. The difference suggests two different behaviors of ω″(q) in the vicinity of a zero ω″, which in turn give two different Q factor divergence rates, i.e., Q−1 ∝ ω″(q) ∝ q − q0 for the former case and Q−1 ∝ ω″(q) ∝ (q − q0)2 for the latter. The quadratic behavior of Q−1 is a generic result of all known BICs (35–37), and as such, the linear behavior of Q−1 found here suggests that pt-BICs belong to a new class of BICs. An analytical discussion of the Q factor divergence rate as well its numerical verification is given in section S2. A direct consequence of the linear behavior of ω″(q) at pt-BICs is the creation of an unstable region on one side of q0, in which ω″(q) is positive and lasing can occur. In this region, the net gain obtained from the particles exceeds the radiation loss, and therefore, the energy of a mode grows exponentially in time. For the case of q0 = −0.090, an unstable region exists between this lasing threshold and the accompanying pt-BIC at q0 = 0.090 as shown in Fig. 3C. Together with the existence of two other pairs of pt-BICs and lasing threshold modes, the q space is divided into four stable passive regions and three unstable lasing regions. Besides the difference in the Q factor divergence rate, there exists another important difference between an ordinary BIC and a pt-BIC. It is well known that an ordinary BIC cannot be excited by an external field, nor can it radiate. However, this is no longer true for pt-BIC. As we will show in a later subsection that a pt-BIC can be excited by an external field even though it does not radiate itself.
Net gain/loss and radiation loss of the resonant modes
Here, via COMSOL simulation, we study numerically the decay or growth rate of the resonant modes and compare it with the analytical result of ω″(q) shown in Fig. 3C. In particular, we wish to numerically verify that ω″ vanishes at both the pt-BICs and the lasing threshold modes. To find this rate at any given q, we first evaluate separately the radiation power Prad and absorption power Pabs per unit cell from the eigenfields generated by COMSOL under the scattering boundary condition (see Materials and Methods). Here, a negative Pabs represents the power due to net gain in the system.
To find the decay/growth rate of a resonant mode, we further calculate the energy stored in the resonant mode, Ueff. The fields E and H used to calculate Ueff are the numerical results obtained from COMSOL simulations minus the analytical far-field radiation fields (see section S3). The decay rates due to both radiative and nonradiative processes can be obtained and denoted, respectively, by γrad = Prad/2Ueff and γabs = Pabs/2Ueff. A negative γabs denotes the growth rate. The results of −γabs and γrad are shown in Fig. 4B. In Fig. 4C, we plot their sum γtot = γrad + γabs (red dotted line) together with the analytically determined total decay rate (yellow solid line), i.e., −ω″(q), which is obtained directly from Fig. 3C. Excellent agreement between the theory and simulation is clearly seen. Although γtot is symmetric in the q space, its physics is very different. For example, in the positive q space, both γrad and γabs vanish at the two pt-BICs, i.e., q = 0.090 and 0.400. This confirms that, like ordinary BICs, pt-BICs have neither radiative loss nor nonradiative loss. In the negative q space, however, there exist two lasing threshold modes at q = −0.090 and −0.400, at which both γrad and γabs are nonzero although γrad + γabs = 0. This confirms that at lasing thresholds, the power obtained through the net gain from the particles is exactly balanced by the power loss due to radiation. Thus, the condition of ω″ = 0 for real eigenfrequency solutions cannot tell the difference between the two possible solutions: pt-BICs and lasing threshold modes. Their difference can only be seen by studying γrad and γabs separately.
(A) Dispersion of the leaky channel. (B) Absorption and radiation rate for the resonant modes. The curves of γrad and γabs intersect six times at which γrad = −γabs. Three of the intersection points with γrad = γabs = 0 correspond to pt-BICs, and the other three with finite γrad correspond to the lasing thresholds at which the radiation loss is exactly balanced by the net gain. (C) Total decay rate γtot = γrad + γabs. This rate is symmetric with respect to the Γ point and coincides with the −ω″(q) obtained analytically. The PT-symmetric perturbation used here is
In section S5, we also show two examples of excitation by external waves. The scattering and absorption for the excitation of a plane wave and point sources are given in figs. S4 and S5, respectively. It can be observed that both the extinction cross section Cext due to the plane wave excitation shown in fig. S4B and the extinction power Pext due to the excitation of point sources shown in fig. S5B exhibit the same behavior for q and −q modes, although the scattering and absorption behaviors are rather different for ±q. This is consistent with the result shown in Fig. 4 that γrad and γabs for the resonant modes are asymmetric in the q space, but the total decay rate γtot is symmetric.
Coupling of BICs and pt-BICs to external fields
When the system is excited by a plane wave ∣E〉 at a frequency ω with a wave vector q in the z direction, the response of the system can be obtained from the equation
For the case of no gain and loss, because a BIC cannot be coupled to the external plane wave, the coupling strength
However, this is no longer true for pt-BICs. We will see that a pt-BIC can be excited by an external field, but it does not radiate. This is because when gain and loss are introduced into the system, the matrix Μ′ + iγσz in Eq. 4 is no longer symmetric at the pt-BICs as shown in Fig. 2B. The left and right eigenvectors of a pt-BIC are hence different. Because the right eigenvector of a pt-BIC,
Let us take the state at q = 0.090 in Fig. 3B as an example. For the pt-BIC marked by the black pentagon, the right eigenvector has the form of zero radiation condition, i.e.,
To demonstrate the nonzero coupling strength at pt-BICs, we have simulated the dipole moments in a unit cell close to a pt-BIC excited by both point sources and plane waves as shown in fig. S6. For both types of sources, the magnitude of the simulated dipole moments (pA, pB) diverges at the same rate as 1/λ−, confirming the analytical prediction that the coupling strength is nonzero at pt-BICs.
A generic phenomenon
The splitting of an ordinary BIC into a pair of pt-BIC and lasing threshold mode under a PT-symmetric perturbation and the associated phenomena discussed above are consequences of the identity h(q,ω) = f *(q,ω), which makes the real-frequency solutions of Eq. 4 possible. Because this identity does not depend on the dielectric constant of nanoparticles, the phenomena also apply to dielectric nanoparticles as long as the dipole approximation is valid. Furthermore, because the identity holds independent of the choice of the system geometric parameters, as mathematically proved in the section S6, the phenomena are guaranteed and robust in our system. For the experimental realization, we can add gain and loss to particles A and B separately by coating them with a thin shell of dielectric material with relative permittivity
The novel phenomena found here are generic for other 1D PT-symmetric system having BICs. Here, we propose another system for the possible experimental realization of the found phenomena. The system considered is a photonic crystal fiber with periodic index modulations (38) shown in the inset of Fig. 5A, where the two components of the dielectric core are of equal thickness, i.e., d = 0.5a, and have relative permittivities ε1 = 4.9 and ε2 = 1. The radius of the core is chosen to be the same as the lattice constant, i.e., r0 = a. The real and imaginary parts of the eigenfrequencies of the transverse electric (TE) modes simulated by COMSOL Multiphysics are shown, respectively, in Fig. 5 (A and B). A symmetry-protected BIC can be found at the Γ point. The corresponding Ey field is spatially localized in the lateral direction as shown in the inset of Fig. 5B. To make the system PT symmetric, we divide each unit cell into four layers of equal thickness and introduce gain or loss into each sublayer in the order ε1 − iη, ε2 + iη, ε2 − iη, and ε1 + iη as shown in the inset of Fig. 5C so that the PT symmetric perturbation ε(z) = ε∗(−z) is satisfied. For the case of η = 0.2, the real and imaginary parts of the eigenfrequencies are shown, respectively, in Fig. 5 (C and D). It can be clearly seen that the PT-symmetric perturbation splits the original BIC into a pair of pt-BIC and lasing threshold mode located symmetrically on the two sides of the q space at which ω″ = 0. Although the pt-BIC and the lasing threshold mode have the same eigenfrequency, they have completely different eigenstates that are, respectively, localized and extended in the lateral direction as shown in the inset of Fig. 5D. To experimentally achieve the PT-symmetric modulation of the dielectric constant, one can place a periodic loss structure on the modulated dielectric core with gain, similar to a previous work (39).
(A) Dispersion of TE modes (real parts of eigenfrequencies). The inset shows the parameters of the photonic crystal fiber. (B) Imaginary parts of eigenfrequencies. The inset shows the simulated field profile of Ey in a unit cell for the eigenstate of the BIC. Under a PT-symmetric perturbation of strength η = 0.2, (C) and (D) show the real and imaginary parts of the eigenfrequencies, respectively. Inset of (C) shows the gain and loss profile in a unit cell. Inset of (D) shows the simulated field profiles of Ey for both the pt-BIC and the lasing threshold mode.
DISCUSSION
In summary, we have shown that a dimerized plasmonic chain provides an analytically solvable model for BICs. A phase diagram that separates one region with propagating BICs from another region without propagating BICs is obtained analytically without resorting to a tedious numerical search for divergent Q factors in a 3D parameter space. When a PT-symmetric perturbation is applied to the system, it is shown rigorously that each original BIC is split into a pair of pt-BIC and lasing threshold mode. The pt-BICs belong to a new class of BICs carrying a Q factor with a reduced divergence rate. Different from the ordinary BICs, this new class of BICs can be excited but themselves do not radiate. It should be pointed out that the reciprocity of the system shown in Eq. 9 guarantees the symmetry of eigenfrequency and decay rate between the states in the q and –q spaces (40). This is true even when a PT-symmetric perturbation is added to the system, consistent with the fact the lasing threshold mode emerges at –q accompanying each pt-BIC at q.
The existence of pt-BICs and lasing threshold modes divide the q space into regions of stable and unstable modes. A pair of pt-BICs or a pair of lasing threshold modes can also be annihilated under parameter changes. These novel phenomena manifest the new physics arising from the interplay between PT-symmetric perturbation and BICs. Such mode splitting phenomena are different from the physics of PT symmetry associated with the existence of exceptional points, which are characterized by mode merging. The fact that a pair of pt-BIC and lasing threshold mode corresponds to the right and left eigenvectors of the non–Hermitian matrix involving gain and loss suggests that the phenomena revealed here are generic and can apply to other photonic or acoustic 1D chains as well as 2D slabs supporting BICs. An example of photonic crystal fiber is demonstrated. Our work offers another perspective on the interesting connections between BICs, non–Hermitian physics, and lasing.
MATERIALS AND METHODS
Coupled dipole equation for dimerized chains
Let us consider a periodic dimerized lattice as specified in the upper inset of Fig. 1. For simplicity, a TE-polarized incidence with the electric field perpendicular to the chain is considered: Eext = E0ei(kxx + qz) and E0 = (0, E0,0). Only the transverse mode is considered, i.e., PA, B = (0, pA, B,0). Using Bloch theorem, the coupled dipole equations for the dimerized chain can be written as
The lattice sums in Eq. 7 can be expressed in terms of certain special functions
Lin(x) and Φ(z, s, a) are the polylogarithm function and the Hurwitz-Lerch transcendent function, respectively.
By changing the variable from n to –n in Eq. 7, it is straightforward to show that the lattice sums in Eq. 7 satisfy the following symmetry relations
Equation 9 is consistent with the inversion symmetry of the dimerized chain with no gain and loss. Under the inversion along the z direction (axis of the chain), z→−z. This yields particle A→particle B and leads to q→−q for the Bloch eigenstate. For the matrix M, the inversion gives
Equation 9 can also be understood from the reciprocity of the system (40). The reciprocity is always maintained in a linear and nonmagnetic optical system, which ensures that the resonant modes at q and –q have the same frequencies and decay rates, namely, ω′( − q) = ω′(q) and ω″( − q) = ω″(q). This relation holds in our system as manifested in Fig. 1 (A and B) for the case of no gain and loss, and in Fig. 3 (B and C) for the PT-symmetric case. However, Eq. 9 only gives an identity, f (−q,ω) = h(q,ω), between the states q and –q. One more identity, i.e., f *(q,ω) = f (−q,ω) proved in section S6, is needed to guarantee the existence of BICs and pt-BICs in this system.
Numerical simulation
Throughout this paper, all numerical studies are performed using COMSOL Multiphysics, in which the module EM waves, frequency domain is used. To calculate Prad, we integrate the time-averaged Poynting vector
It should be noted that the above definitions of Psca and Pabs are not only limited to the study of eigenstates of the system but also valid for the scattering problem when the system is excited by external sources. For the scattering of plane wave, the scattering and absorption cross sections can be respectively defined by Csca = Psca/Ii and Cabs = Pabs/Ii, where Ii is the incident irradiance. Their sum gives the extinction cross section, i.e., Cext = Csca + Cabs.
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/34/eabc1160/DC1
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
REFERENCES AND NOTES
- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).