## 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 *r*_{0} separated by a distance *d*. Their relative permittivity is described by the Drude model, _{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* > 3*r*_{0} and *a*-*d* > 3*r*_{0} (*32*), each plasmonic nanoparticle can be treated as an electric dipole with an inverse polarizability *33*) and *k*_{0} = ω/*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

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 *p*_{A}) and B (*p*_{B}) along the *y* axis can be obtained from the multiple-scattering theory by solving the following equation*E*_{A(B)} = *E*_{0}*e*^{∓iqd/2}. Equation 1 can be expressed as **MP** = **E** for convenience. The expressions of the lattice sums *S _{ij}* are given in Eq. 7 (see Materials and Methods). Because of the inversion symmetry of the system,

*S*have the following symmetry properties:

_{ij}*S*

_{AA}(

*q*,ω) =

*S*

_{BB}(

*q*,ω),

*S*

_{AA}(−

*q*,ω) =

*S*

_{AA}(

*q*,ω), and

*S*

_{AB(BA)}(−

*q*,ω) =

*S*

_{BA(AB)}(

*q*,ω) as shown in Eq. 8. Below, we denote

*S*

_{AA(BB)}(

*q*,ω) by

*S*

_{00}(

*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

*S*.

_{ij}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 *q*. In the case of no gain or loss, i.e., **Μ**^{′}**P**^{′} = 0. The eigenfrequencies of the resonant modes are determined by the complex frequency solutions of *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 *r*_{0} = 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 *S*_{00}(−*q*,ω) = *S*_{00}(*q*,ω) and *S _{ij}*(−

*q*,ω) =

*S*(

_{ji}*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 *9*, *34*). Thus, the condition for BICs becomes the real-frequency solutions of

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 *S*_{AB}*e ^{iqd}* =

*S*

_{BA}

*e*

^{-iqd}. This makes the matrix

**Μ**

^{′}symmetric, which in turn makes the condition

*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, *f*) is plotted with red solid line. It can also be seen that, for a given value of *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 *S*_{AB}*e ^{iqd}*), and Im(−

*S*

_{BA}

*e*

^{−iqd}) 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.

Because _{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 *f*) shifts upward and moves the two propagating BICs toward the Γ point as *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**Μ**′ + *i*γσ* _{z}*)

**P**′ = 0, where σ

*is the*

_{z}*z*component of the Pauli matrix. The conditions for the existence of real-frequency solutions of Eq. 4 in conjunction with the condition

*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

*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

*S*

_{AB}

*e*) and Im(−

^{iqd}*S*

_{BA}

*e*

^{−iqd}) shown in Fig. 2A for the case of no gain or loss. The three intersection points of the functions

*S*

_{AB}

*e*), i.e., the nodal points of Im(

^{iqd}*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

*S*

_{BA}

*e*

^{−iqd}), or the nodal points of Im(

*h*) − γ.

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 (*q*_{0},ω_{0}) with *q*_{0},ω_{0}) with *pt*-BIC, we have *P*_{gain} > 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 *S _{ij}*(−

*q*,ω) =

*S*(

_{ji}*q*,ω) for

*i*≠

*j*.

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 *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*γσ

*in the determinant does not affect the symmetry of ω(*

_{z}*q*) = ω(−

*q*). In Fig. 3B, we also plot the analytical dispersion of a higher-frequency branch

*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., *q*_{0} *=* ±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 *q*_{0} *=* ±0.090 have completely different eigenstate profiles. The **E** field is spatially localized in the transverse plane for the state at *q*_{0} *=* 0.090 without far-field radiation, indicating a *pt*-BIC (black pentagon). However, an extended state is observed at *q*_{0} = −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* − *q*_{0} for the former case and *Q*^{−1} ∝ ω″(*q*) ∝ (*q* − *q*_{0})^{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 *q*_{0}, 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 *q*_{0} = −0.090, an unstable region exists between this lasing threshold and the accompanying *pt*-BIC at *q*_{0} = 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 *P*_{rad} and absorption power *P*_{abs} per unit cell from the eigenfields generated by COMSOL under the scattering boundary condition (see Materials and Methods). Here, a negative *P*_{abs} 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, *U*_{eff}. The fields **E** and **H** used to calculate *U*_{eff} 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} = *P*_{rad}/2*U*_{eff} and γ_{abs} = *P*_{abs}/2*U*_{eff}. 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.

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 *C*_{ext} due to the plane wave excitation shown in fig. S4B and the extinction power *P*_{ext} 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* _{n}* = λ

*and λ*

_{−}_{+}are two eigenvalues of the matrix

**M**in Eq. 1 with parameters

*q*and ω, and ∣

**P**

*〉 and*

_{n}**M**, i.e.,

**M**∣

**P**

*〉 = λ*

_{n}*∣*

_{n}**P**

*〉 and*

_{n}*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 (

*q*

_{0},ω

_{0}), i.e.,

*q*≅

*q*

_{0}, ω ≅ ω

_{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

*16*).

For the case of no gain and loss, because a BIC cannot be coupled to the external plane wave, the coupling strength **E**〉 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.,

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,

*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., *S*_{00}(−*q*,ω) = *S*_{00}(*q*,ω) and *S _{ij}*(−

*q*,ω) =

*S*(

_{ji}*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,

*q*= −0.090, we have

*q*= 0.090, implying a nonzero coupling strength

*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 (*p*_{A}, *p*_{B}) 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.5*a*, 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., *r*_{0} = *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 *E _{y}* 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*).

## 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: **E**_{ext} = **E**_{0}*e*^{i(kxx + qz)} and **E**_{0} = (0, *E*_{0},0). Only the transverse mode is considered, i.e., **P**_{A, B} = (0, *p*_{A, B},0). Using Bloch theorem, the coupled dipole equations for the dimerized chain can be written as*G* is the Green’s function for electric dipoles. The above coupled dipole equations can be reexpressed as**MP** = **E**, where α is the polarizability of the nanoparticles, and the lattice sums *S _{ij}* of the Green’s functions have the following forms

The lattice sums in Eq. 7 can be expressed in terms of certain special functions

*Li _{n}*(

*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 *P*_{rad}, we integrate the time-averaged Poynting vector *P*_{rad} = ∬_{s}**S**_{rad} ⋅ **n***dS*. For *P*_{abs}, we integrate the change of energy density due to gain or loss inside the two nanoparticles in a unit cell, namely, * _{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 *P*_{sca} and *P*_{abs} 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 *C*_{sca} = *P*_{sca}/*I*_{i} and *C*_{abs} = *P*_{abs}/*I*_{i}, where *I*_{i} is the incident irradiance. Their sum gives the extinction cross section, i.e., *C*_{ext} = *C*_{sca} + *C*_{abs}.

## 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

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

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