Research ArticleOPTICS

Coexistence of a new type of bound state in the continuum and a lasing threshold mode induced by PT symmetry

See allHide authors and affiliations

Science Advances  19 Aug 2020:
Vol. 6, no. 34, eabc1160
DOI: 10.1126/sciadv.abc1160


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.


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 (217). Various mechanisms have been proposed to produce BICs (713). 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 (1315).

On the other hand, systems with parity-time (PT) symmetry have drawn considerable interest recently (1828). 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 (2931). 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.


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, ε(ω)=1ωp2/ω(ω+iγ0), where ωp is the plasma frequency and γ0 is the collision frequency. In the study of BICs, we are only interested in the existence of a bound state that is decoupled from the far field so that the radiative Q factor diverges (6). Like previous works of BICs, we ignore the intrinsic loss γ0 and focus on the radiation loss. In the case of d > 3r0 and a-d > 3r0 (32), each plasmonic nanoparticle can be treated as an electric dipole with an inverse polarizability α01(ω)=1r03ωp23ω2ωp22i3k03, where the imaginary part denotes the radiation loss (33) and k0 = ω/c, with c being the speed of light in vacuum. The PT-symmetric gain and loss can be added to particles A and B by coating them with thin layers of gain and loss media, respectively. The inverse polarizability can now be written as αA1(ω)=α01(ω)+iγ and αB1(ω)=α01(ω)iγ, where γ > 0 represents the amounts of gain and loss added to particles A and B, respectively.

Fig. 1 BICs in dimerized chains and their phase diagram.

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 ω¯p, and the solid red line denotes the nodal line of Im(f). Three intersection points of the nodal lines of Re(f) and Im(f) can give rise to three BICs. Two propagating BICs meet and are annihilated at q = 0 when a critical ω¯pc is reached, above which the propagating BICs disappear suddenly. (F) Phase diagram for the existence of propagating BICs.

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 Eext=E0êyei(kxx+qzωt), the transverse polarizations of particles A (pA) and B (pB) along the y axis can be obtained from the multiple-scattering theory by solving the following equation[αA1SAASABSBAαB1SBB][pApB]=[EAEB](1)as derived in Materials and Methods, where EA(B) = E0eiqd/2. Equation 1 can be expressed as MP = E for convenience. The expressions of the lattice sums Sij are given in Eq. 7 (see Materials and Methods). Because of the inversion symmetry of the system, Sij have the following symmetry properties: SAA(q,ω) = SBB(q,ω), SAA(−q,ω) = SAA(q,ω), and SAB(BA)(−q,ω) = SBA(AB)(q,ω) as shown in Eq. 8. Below, we denote SAA(BB)(q,ω) by S00(q,ω). The analytical study of BICs is made possible by expressing these lattice sums in terms of polylogarithm and Hurwitz-Lerch transcendent functions. Note that the material properties and geometrical parameters are contained separately in αA(B)1 and Sij.

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 P=(pAeiqd/2, pBeiqd/2)T to reflect the explicit phase difference between the two particles in the unit cell for the mode q. In the case of no gain or loss, i.e., αA1=αB1=α01, the periodic part of the eigenvector, P=(pA, pB)T, satisfies the following equations(α01S00)pASABeiqdpB=0SBAeiqdpA+(α01S00)pB=0(2)or simply ΜP = 0. The eigenfrequencies of the resonant modes are determined by the complex frequency solutions of det(Μ)=(α01S00)2SABSBA=0 denoted by ω(q) = ω(q) + iω(q). The function ω(q) gives the dispersion relation of the resonant modes, and the ratio ω(q)/2ω(q) gives the Q factor. The zero of ω(q) indicates the position of BICs at which the Q factor diverges. For any given q, there also exists other complex frequency solutions of det(Μ) = 0. Here, we are only interested in the leaky channel along which ω(q) is the smallest and vanishes at certain values of q when a BIC occurs.

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 pB=pA (9, 34). Thus, the condition for BICs becomes the real-frequency solutions off(q,ω)=def(α01S00)+SABeiqd=0h(q,ω)=defSBAeiqd+(α01S00)=0(3)

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 det(Μ)=[(α01S00)+SABeiqd][(α01S00)SABeiqd]=0 satisfied. Thus, although f (q,ω) = 0 and det(Μ) = 0 have different solutions in general, they do share the same real-frequency solutions at BICs.

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, ω¯p=ωp(a/2πc), and the nodal line of Im(f) is plotted with red solid line. It can also be seen that, for a given value of ω¯p, the nodal lines of Re(f) and Im(f) intersect each other at three points, giving rise to three BICs: one symmetry-protected BIC at q = 0 and two propagating BICs located symmetrically on the two sides of the Γ point. To see this more explicitly, in Fig. 2A, we plot the functions Im(α01S00), Im(−SABeiqd), and Im(−SBAeiqd) evaluated at the frequency of Re(f) = 0. This figure shows clearly that these three curves intersect at three BIC points at which Im(f) = 0, as indicated by vertical dashed lines.

Fig. 2 Analysis of the elements of linear response matrix.

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, γ=2×103/r03, and γ=6×103/r03. (A) Without PT-symmetric perturbation, the ordinary BICs correspond to the intersection points of three matrix elements: α01S00=SABeiqd=SBAeiqd or equivalently Im(f) = Im (h) = 0. (B) As γ < γc, the intersection points of Im(α01S00+iγ) [or Im(α01S00iγ)] and Im( − SABeiqd) [or Im( − SBAeiqd)] give rise to Im(f) = −γ [or Im(h) = γ] with a shift in q. (C) As γ > γc, only one intersection point exists for Im(α01S00+iγ) and Im( − SABeiqd) in the negative q space, giving rise to Im(f) = −γ.

Because Im(α01)=2k03/3 is independent of the material parameter ωp, the nodal line of Im(f) depends only on the geometrical parameters, while the nodal line of Re(f) depends on the value of ωp. In Fig. 1E, we plot the nodal lines of Re(f) at different values of ωp, where the dimensionless plasma frequency ω¯p=ωp(a/2πc) is indicated. The nodal line of Re(f) shifts upward and moves the two propagating BICs toward the Γ point as ω¯p increases. These propagating BICs merge with the symmetry-protected BIC when some critical ω¯pc is reached and then disappear when ω¯p>ω¯pc. To obtain a phase diagram for the existence of propagating BICs, we perform the similar study but vary the geometrical parameters as well. In Fig. 1F, we show a critical surface of ω¯pc as a function of parameters a and d in the region where dipole approximation is valid. Such a surface represents a boundary separating two phases. Below this boundary, propagating BICs are guaranteed and robust against small perturbations in both material and geometrical parameters. Above the boundary, there are no propagating BICs. Thus, our system provides an analytically solvable model for studying the phase diagram of propagating BICs without resorting to a tedious numerical search for divergent Q factors through parameter tuning.

Pairing of pt-BIC and lasing threshold mode

When gain and loss are added to the system, Eq. 2 becomes(α01S00+iγ)pASABeiqdpB=0SBAeiqdpA+(α01S00iγ)pB=0(4)or (Μ′ + iγσz)P′ = 0, where σz is the z component of the Pauli matrix. The conditions for the existence of real-frequency solutions of Eq. 4 in conjunction with the condition pB=pA for the cancellation of radiation fields become f (q,ω) = −iγ and h(q,ω) = iγ. Because h(q,ω) = f *(q,ω), we only need to consider the real-frequency solutions of f (q,ω) = −iγ. For the case of ω¯p=0.95, in Fig. 3A, we plot the nodal line of Re(f) with a dashed line and the nodal lines of Im(f) + γ with colored solid lines for different values of γ. For each value of γ, the intersection of the solid curves with the dashed curve produces three pt-BICs: one in the negative q space and two in the positive q space when 0 < γ < γc. They emerge from the three original BICs when the system has no gain or loss, but they are shifted in the q space. The shift of each pt-BIC can be explicitly seen from Fig. 2B, in which we choose γ=2×103/r03 and plot the functions Im(α01S00)+γ and Im(α01S00)γ with black dashed and yellow dashed lines, respectively, in addition to the two curves Im(−SABeiqd) and Im(−SBAeiqd) shown in Fig. 2A for the case of no gain or loss. The three intersection points of the functions Im(α01S00)+γ and Im( − SABeiqd), i.e., the nodal points of Im(f) + γ, show the positions of pt-BICs, which are indicated by the vertical dashed lines. The same vertical dashed lines are obtained from the intersection of Im(α01S00)γ and Im(−SBAeiqd), or the nodal points of Im(h) − γ.

Fig. 3 pt-BICs and lasing threshold modes in PT-symmetric dimerized chains.

(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 ω¯p=0.95 give three pt-BICs: one pt-BIC in the q < 0 region, which is shown more clearly in the inset, and two pt-BICs in the q > 0 region. At γc, the latter two pt-BICs are annihilated. They disappear suddenly when γ > γc, and the system is left with only one pt-BIC in the q < 0 region. (B) Dispersion of the leaky channel simulated by COMSOL Multiphysics (yellow solid line). The two red dotted lines denote the dispersions calculated analytically by solving det(Μ′ + iγσz) = 0 of Eq. 4. The lower-frequency branch, which agrees excellently with the COMSOL result, represents the leaky channel that has much higher Q values than that of the upper frequency branch. Three pairs of real-eigenfrequency modes obtained numerically are marked by polygons. (C) The imaginary part of the eigenfrequency ω(q) for the leaky channel obtained numerically (yellow solid line) and analytically (red dotted line). (D) Simulated field profiles of |E| in the transverse plane for the real-eigenfrequency states. The field patterns of the three BICs shown by black polygons are spatially localized, whereas the field patterns of the three lasing threshold modes marked by red polygons are spatially extended.

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 ω¯p is varied.

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 (q00) with pB=pA implies another solution of Eq. 4 at the point (−q00) with p¯B=SBA(q0,ω0)eiq0dSAB(q0,ω0)eiq0dp¯A. Different from the pt-BIC, we have p¯B<p¯A for this mode as observed from Fig. 2, which implies a net gain of power, i.e., Pgain > 0. Because this mode has a real frequency, the net gain must be exactly balanced by the radiation loss, implying a lasing threshold. The numerical verification of the balance between net gain and radiation loss at the lasing threshold will be given in the next subsection. It is interesting to note that at a lasing threshold, the vector P¯(p¯A,p¯B) actually corresponds to the left eigenvector of the equation P¯(Μ+iγσz)=0 as can be proved directly using the identity Sij(−q,ω) = Sji(q,ω) for ij.

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 ω¯p=0.95 and γ=1.7×103/r03, the results of ω′(q) and ω″(q) for the leaky channel calculated by COMSOL are shown by the yellow solid curve in Fig. 3 (B and C), respectively. The two red dotted curves illustrate the analytical results obtained from the condition det(Μ′ + iγσz) = 0 of Eq. 4. Excellent agreements are found between the theory and simulations. It is easy to show that the presence of the term iγσz in the determinant does not affect the symmetry of ω(q) = ω(−q). In Fig. 3B, we also plot the analytical dispersion of a higher-frequency branch ω>(q) (red dotted curve). The corresponding imaginary part ω>(q) is much larger in value and does not vanish. We are only interested in the dispersion of the leaky channel, ω′(q), along which pt-BICs and lasing threshold modes can occur at the zeros of ω″(q).

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) ∝ qq0 for the former case and Q−1 ∝ ω″(q) ∝ (qq0)2 for the latter. The quadratic behavior of Q−1 is a generic result of all known BICs (3537), 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.

Fig. 4 Decay rates calculated from COMSOL Multiphysics.

(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 γ=1.7×103/r03.

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 equationP=n1λnPnP¯nE=nWn(q,ω)λnPn(5)where λn = λ and λ+ are two eigenvalues of the matrix M in Eq. 1 with parameters q and ω, and ∣Pn〉 and P¯n are, respectively, the right and left eigenvectors of M, i.e., MPn〉 = λnPn〉 and P¯nM=λnP¯n, and Wn(q,ω)=P¯nE (33). The dispersion ω(q) obtained from the condition Re(λ) = 0 corresponds to the leaky channel we are interested in, i.e., the lower-frequency branch shown in Fig. 3B. The dispersion obtained from Re(λ+) = 0 corresponds to the higher-frequency branch in Fig. 3B. Close to a BIC at the point (q00), i.e., qq0, ω ≅ ω0, λ is vanishingly small and the response of the system ∣P〉 is mainly determined by the right eigenvector ∣P〉 and the quasi-mode coupling strength W(q,ω)=P¯E (16).

For the case of no gain and loss, because a BIC cannot be coupled to the external plane wave, the coupling strength P¯E is zero at BICs. If we rewrite the left eigenvector P¯ as (p¯Aeiqd/2,p¯Beiqd/2) and the external plane wave ∣E〉 as (EAeiqd/2,EBeiqd/2)T, we have P¯E=0 or p¯B=p¯A as ∣E′〉 ∝ (1,1)T. Because the matrix M′ in Eq. 2 is symmetric at the BICs as shown in Fig. 2A, its right eigenvector has the same entries as the left eigenvector, i.e., P=(P¯)T. This also implies a vanishing radiation loss at BICs, i.e., EP=0. The equivalence between zero coupling strength and zero radiation loss is expected for all ordinary BICs.

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, P, satisfying Eq. 4 has an antisymmetric form, i.e., pB=pA, the corresponding left eigenvector, P¯, will not be antisymmetric, i.e., p¯Bp¯A, implying P¯E0. This shows that a pt-BIC can be excited by external fields, although it does not radiate.

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., P0.090=12{1,1}T. From the symmetry relations S00(−q,ω) = S00(q,ω) and Sij(−q,ω) = Sji(q,ω) shown in Eq. 8, it is easy to show that the left eigenvector of the pt-BIC has the same entries as the right eigenvector of the accompanying lasing threshold at q = −0.090, namely, P¯q=PqT. Because it has been shown that the amplitude of the dipole moment of particle A is larger than that of particle B at q = −0.090, we have p¯B<p¯A for q = 0.090, implying a nonzero coupling strength P¯E at the pt-BIC. A nonzero coupling strength W makes the response vector ∣P〉 of Eq. 5 divergent when a pt-BIC is approached as λ→0. This is a typical result of the excitation precisely at the frequency of a resonant mode, which has zero spectrum width. Thus, the PT-symmetric perturbation changes the singular behavior of a pt-BIC, which not only reduces the Q factor divergence rate but also violates the equivalence between zero coupling strength and zero radiation loss.

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 ε1+iε1A for particle A and ε1+iε1B for particle B. The relations between the amount of gain or loss, i.e., ε1A<0 or ε1B>0, and the value of γ or −γ in Eq. 4 are given in section S4.

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 ε1iη, ε2 + iη, ε2iη, 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).

Fig. 5 Splitting of an ordinary BIC in a PT-symmetric photonic crystal fiber.

(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.


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.


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 aspA=αA(EA+n0G(na)pAeiqna+nG(na+d)pBeiqna),pB=αB(EB+n0G(na)pBeiqna+nG(nad)pAeiqna)(6)where G is the Green’s function for electric dipoles. The above coupled dipole equations can be reexpressed as[αA1SAASABSBAαB1SBB][pApB]=[EAEB]or in a more compact form MP = E, where α is the polarizability of the nanoparticles, and the lattice sums Sij of the Green’s functions have the following formsSAA=SBB=k03Σn0(1k0na+i1k02na21k03na3)eik0naeiqna,SAB=k03Σn(1k0na+d+i1k02na+d21k03na+d3)eik0na+deiqna,SBA=k03Σn(1k0nad+i1k02nad21k03nad3)eik0nadeiqna(7)

The lattice sums in Eq. 7 can be expressed in terms of certain special functionsSAA=SBB=k03(L1k0a+iL2k02a2L3k03a3),SAB=k03(Φ1,1k0a+iΦ1,2k02a2Φ1,3k03a3Φ0),SBA=k03(Φ1,1k0a+iΦ1,2k02a2Φ1,3k03a3Φ0)whereLn=Lin(ei(k0q)a)+Lin(ei(k0+q)a),Φm,n=Φ(ei(k0+q)a,n,mda)eimk0d+Φ(ei(k0q)a,n,mda)eimk0d,Φ0=(1k0d+i1k02d2+1k03d3)eik0d

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 relationsSAA(BB)(q,ω)=S(AA)BB(q,ω)SAB(q,ω)=SBA(q,ω)(8)or equivalently,M(q,ω)=MT(q,ω)(9)

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 givesM(q,ω)=σxM(q,ω)σxwhich gives rise to Eqs. 8 and 9.

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 Srad=c8πRe(Esca×Hsca) on a cylindrical surface enclosing a unit cell, i.e., Prad = ∬sSradndS. For Pabs, we integrate the change of energy density due to gain or loss inside the two nanoparticles in a unit cell, namely, Pabs=18πi=A,BViωIm(εi)Ei2dV. Im(εi) (i = Α,Β) are the imaginary parts of the relative permittivities of particles A and B, which include both the added gain/loss and the intrinsic loss of the plasmonic particles. The permittivities εA,B we used in simulation are related to the PT-symmetric parameter γ defined in Eq. 4, which is derived in section S4.

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 material for this article is available at

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.


Acknowledgments: D.H. thanks A. Sadreev, X. D. Zhang, and L. Shi for the helpful discussions. Funding: This work was supported by the National Natural Science Foundation of China (grant nos. 91750102, 11727811, and 11771248), the Hong Kong Research Grants Council (under grant nos. AoE/P-02/12, 16303119, and C6013-18G-A), the Projects of President Foundation of Chongqing University (grant no. 2019CDXZWL002), and the Natural Science Foundation of Xinjiang Autonomous Region (grant no. 2019D01C026). Author contributions: D.H. and C.T.C. conceived the original idea. Q.S., D.H., and Z.Q.Z. derived the theory and wrote the manuscript. Q.S., J.H., and S.D. performed the theoretical and numerical analysis. J.H. and C.Z. provided a rigorous mathematical proof for a theorem. All authors contributed to the analysis and the discussion of results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article