# Fileset

[bic_final.pdf](https://mdr.nims.go.jp/filesets/ec8f9fe2-30a0-42fe-8482-69d915a6160d/download)

## Creator

[Tetsuyuki Ochiai](https://orcid.org/0000-0003-2933-0014)

## Rights

[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Generation of multiple bound states in the continuum through doubly degenerate quasiguided modes](https://mdr.nims.go.jp/datasets/52695ba8-6468-4b8c-941e-5e516113751c)

## Fulltext

Generation of multiple bound states in the continuum through doubly degeneratequasiguided modesTetsuyuki OchiaiResearch Center for Electronic and Optical Materials,National Institute for Materials Science (NIMS), Tsukuba 305-0044, Japan(Dated: August 5, 2024)We present a detailed theoretical analysis of a peculiar generation of multiple bound states in thecontinuum (BICs) in two-dimensional periodic arrays of dielectric spheres. They emerge in high-symmetry lattices with the C6v and C4v point groups and involve doubly degenerate quasi-guidedmodes at the Γ point that can couple to external radiation. By tuning a system parameter, thedoubly degenerate modes can exhibit accidental BICs at a critical parameter. In the vicinity of thecritical parameter, the two bands originating from the degenerate mode exhibit multiple off-Γ BICs.They move and annihilate across the two bands by changing the parameter around the critical one.A ring-like high Q channel pinned with multiple BICs emerges particularly for the C6v case. Acrossthe critical coupling, the total vorticity of the multiple BICs is conserved. The k · p perturbationtheory explains some features of the phenomena reasonably well.I. INTRODUCTIONOptical bound states in the continuum (BIC) are lo-calized eigenstates embedded in the radiation continuum[1]. They have infinite quality factors or, in other words,vanishing decay rates, although they are inside the lightcone. Strong light confinement via the infinite qualityfactor enables us to investigate various applications suchas lasing [2], sensing [3], nonlinear optics [4], and so on,via BICs.The BICs are found typically in monolayers of spheresand photonic crystal slabs at the Γ point [5–8]. Off-ΓBICs are also found at generic k points [9, 10], via the for-mation of polarization vortices [11]. The former BICs aresymmetry-protected and irrelevant to physical parame-ters unless the relevant symmetry is unchanged. Thelatter are topologically protected, moving in momentumspace. Therefore, a parameter scan is necessary to findthe latter BICs at a given k point other than Γ. In thissense, the off-Γ BICs are sometimes called accidental.Recently, merging of an at-Γ BIC and off-Γ BICs in anon-degenerate isolated band attracts much interest asthe so-called super BIC [12, 13]. This BIC is obtained bytuning system parameters such that the off-Γ BICs movetoward the Γ point, where the symmetry-protected BICexists. The super BIC has superb properties as it involvesan extreme suppression of the decay rate in a broad re-gion of the momentum space around Γ. The decay ratebehaves as k6 (k8 for higher-order charges [14]), in strik-ing contrast to the k2 (k4 for higher-order charges) be-havior of the symmetry-protected BICs. Merging amongBICs is available also off the Γ point [15], resulting in a(k − kBIC)4 behavior.Another type of super BICs was found recently inmonolayers of spheres [16]. From now on, we call it thenext-to-super BIC. It is obtained at the Γ point by tuningsystem parameters and exhibits a k4 scaling of the decayrate at a critical parameter. In contrast to the ordinarysuper BIC, it involves doubly-degenerate eigenmodes atthe Γ point. Moreover, off-critical parameters result in aring-like high Q channel in the triangular lattice system.This high Q channel is either P or S-polarization-likeand can be switched by changing the parameter acrossthe critical value.Once we have a high Q channel, we have a continu-ous distribution of optical modes with strongly enhancedlight-matter interactions. The channel can be selectivelyexcited by the polarization of the incident light. As aresult, we can steer the coherent radiation continuouslyalong the ring. The output light has nearly the samefrequency, but the angle changes continuously. Such abeam steering paves the way for various applications.In this paper, we further investigate the next-to-superBICs and show that multiple off-Γ BICs are generated atoff-critical parameters nearby. We show more complexbehaviors than in Ref. [16] can be observed. In addi-tion to the ring-like channel, discrete multiple BICs canbe generated through the next-to-super BICs. In thiscase, the ring is formed asymmetrically in the parame-ter space. If we change the parameter across the criticalvalue, the ring disappears, and multiple discrete BICsare again generated. Across the critical parameter, thetotal vorticity of the BICs is conserved. In addition, theappearance of the ring depends on the lattice structure.One of the points behind the above phenomena is spa-tial symmetry. We show that the k · p perturbation the-ory based on spatial symmetry explains some propertiesof the phenomena reasonably well.The ring-like high Q channel can be explained in termsof the multipolar lattice [17, 18]. However, the lattice-structure dependence needs complementary approachesto understand the phenomena. The k · p perturbationprovides one reasonable scenario as shown in this paper.The multiple off-Γ BICs are also produced by break-ing spatial symmetries of photonic crystal (PhC) slabs[19–21]. There, degenerate or nondegenerate symmetry-protected BICs at Γ turn into multiple off-Γ BICs. Sincethe original symmetry-protected BICs exist irrespectiveof system parameters, the resulting multiple BICs aredeterministic, and parameter tuning is unnecessary.2This paper is organized as follows. In Secs. II andIII, we present numerical examples of the accidental at-ΓBICs and resulting multiple off-Γ BICs in the triangularand square lattices of dielectric spheres, respectively. InSec. IV, we analyze these phenomena via the k ·p pertur-bation theory of the Maxwell equation. Finally, in Sec.V, we summarize the results.II. TRIANGULAR LATTICELet us first consider the next-to-super BICs in a mono-layer of dielectric spheres arranged in a triangular lat-tice. The system has the C6v point group symmetry[22]. Therefore, the eigenmodes are classified accord-ing to the irreducible representations of C6v. At the Γpoint, the eigenmodes with the E1 representation of C6vcan couple to external radiation. The other eigenmodesof A1, A2, B1, B2 and E2 representations are uncoupled,provided their eigenfrequencies are below the diffractionthreshold. Thus, the former eigenmodes generally havefinite lifetimes, whereas the latter eigenmodes have infi-nite lifetimes and are symmetry-protected BICs.However, tuning system parameters can cause the for-mer eigenmodes to have infinite lifetimes. In this way,accidental at-Γ BICs can occur for the E1 modes.Figure 1 shows the resonance angular frequency ωk andwidth γk of the E1 modes at the Γ point (k = 0). Theyare evaluated with the photonic Korringa-Kohn-Rostoker(KKR) method [23–25] together with the curve fittingto the Breit-Wigner formula [26] of the scattering phaseshift δ as a function of angular frequency ω:e2iδ = e2iδk(ω − ωk − iγkω − ωk + iγk)n, (1)where δk is the background (frequency-independent)phase shift, and n is the degree of degeneracy of theresonant mode concerned. The scattering phase shift isderived from the S-matrix of the monolayer [27]. At sev-eral radii, the eigenmodes exhibit infinite quality factoraccidentally.Such an accidental t-Γ IC can be easily found theo-retically and experimentally by monitoring the changeof the resonance signal in the transmission spectrum ofthe normal incidence. For instance, Fig. 2 shows thechange of the transmission spectrum with the sphere ra-dius around the critical one of BIC1 in Fig. 1. Since themodes other than E1 are symmetry-protected, they donot affect the transmission spectrum. Solely the modesof the E1 representation emerge as an asymmetric res-onance signal of the Fano shape [28] in the spectrum.The spectrum changes with the radius, and we can finda narrowing of the resonance width toward the criticalradius. The vanishing resonance width corresponds tothe accidental at-Γ BIC.The above property of the next-to-super BIC presentsa marked contrast to the super BIC. In the latter case, 0.4 0.5 0.6 0.7 0.8 0.3  0.35  0.4  0.45  0.5BIC1BIC4BIC2BIC3BIC5(a)ω0a/2πcrs/a0.00.51.01.52.0 0.3  0.35  0.4  0.45  0.5(b)(☓10-3)γ 0a/2πcrs/az-even  z-odd FIG. 1. The resonance angular frequencies ω0 (a) and widthsγ0 (b) of the E1 modes at the Γ point in the triangular-latticemonolayer of dielectric spheres with the C6v point-group sym-metry, as a function of sphere radius rs. The modes are clas-sified according to the parity in the z direction, assuming thecenters of the spheres are on the z = 0 plane. The dielec-tric constant of the spheres is taken to be 12, and that of thebackground material is 1. The lattice constant is denoted bya. Accidental at-Γ BICs (next-to-super BICs) occur aroundthe modes indicated by solid circles. The inset in (a) standsfor the geometry of the monolayer.the BIC consists of a symmetry-protected BIC and topo-logically protected BICs, so it is impossible to observe thesuper BIC via the transmission spectrum of the normalincidence. It is available only through a detailed analy-sis of the momentum and parameter dependence of thespectrum.Figure 3 shows the real and imaginary photonic bandstructure around the critical radius of BIC1. Here, thereal photonic band structure is referred to as the reso-nance angular frequency ωk, and the imaginary one is tothe resonance width γk, as a function of Bloch momen-tum k. Although the real band structure does not changeso much in its shape, the imaginary band structure ex-hibits a clear contrast between the two bands of oppositeparities and among the three values of on- and off-criticalradii. Below the critical radius [Figs. 3 (a,b)], the odd-parity band exhibits the off-Γ BIC at about ka/2π = 0.1.At the critical radius [Figs. 3(c,d)], both the bands ex-hibit a flat region of nearly zero values in their imaginaryparts. This suppression of the decay rate resembles thatin the super BIC. However, now γk behaves as k4 insteadof the k6 behavior of the super BIC. Above the criticalradius [Figs. 3 (e,f)], the even-parity band exhibits theoff-Γ BIC.For comparison, the photonic bands originating from3 0 0.2 0.4 0.6 0.8 10.45 0.50 0.55 0.600.36ars=0.41aTransmittanceωa/2πcFIG. 2. Specular transmission spectra of the normal incidencein the triangular lattice monolayer of dielectric spheres. Thesphere radius varies from 0.36a to 0.41a with a 0.005a step.The other parameters are the same as in Fig. 1. The vanishingresonance signal of the red curve indicates the accidental at-ΓBIC of BIC1 in Fig. 1.the symmetry-protected BICs at Γ exhibit a quick blow-up of the imaginary parts with the k2 scaling as shownin Fig. 4. Near the Γ point, we do not observe the trendof decreasing the imaginary part toward a minimum. Weshould point out that the imaginary part is much largerthan that from the E1 mode, whose imaginary part issignificantly suppressed.By scanning the Bloch momentum in all the directionsaround the Γ point, the real band structure consists oftwo surfaces touched quadratically at the Γ point, asshown in Fig. 5. The upper and lower bands are S-and P-polarization-like, respectively.Figure 6 shows the Q-value map of the photonic bandmodes at the off-critical radii. It is remarkable that thering-like high Q channels found in Ref. [16] are formed.The channels lies in the S-polarization-like upper bandat rs < rc and in the P-polarization-like lower band forrs > rc. The other bands (the lower band at rs < rc andthe upper band at rs > rc) do not have such channels.These channels indicate that the minimum decay rateγk at the off-critical radii found in Fig. 3 emerges atnearly the same distance from the Γ point regardless ofmomentum orientation. As we move the radius rs towardrc, we can show the ring shrinks to the Γ point and acrossrs = rc, the ring moves from the upper band to the lowerband or vice versa.Looking closely at the rings, we find twelve off-Γ BICs:six on the equivalent ΓM axes and six on the equiva-lent ΓK axes. These axes are the mirror axes of thetriangular-lattice Brillouin zone. The polarization vor-0.5100.5120.5140.5160.5180.520  0.15 (M←) Γ (→K) 0.15  (☓10-5)(e)ωka/2πcka/2π 01234  0.15 (M←) Γ (→K) 0.15  (f)γ ka/2πcka/2π 0.5140.5160.5180.5200.522 (c)ωka/2πc01234(d)γ ka/2πc0.5200.5220.5240.5260.5280.530(a)ωka/2πc evenodd01234Γ KM(b)γ ka/2πcFIG. 3. Real and imaginary photonic band structures of thez-even parity originating from the E1 mode around the criticalradius rc of BIC1 in Fig. 1. The momentum is taken alongthe ΓM and ΓK directions. The band structure is classifiedaccording to the in-plane parity concerning the ΓM and ΓKaxes. The sphere radius is fixed as rs = 0.38a in (a) and (b),rs = 0.385a (≃ rc) in (c) and (d), and rs = 0.388a in (e)and (f). The other parameters are the same as in Fig. 1. Theinset in (b) represents the first Brillouin zone of the triangularlattice.tices of the BICs are very elongated along the rings. Thevortex charge diagram is shown in Fig. 7. The vortexcharges are opposite between ΓM and ΓK. Off the mirroraxes, the Q values in the ring are still very high, of order10−8. Thus, the ring-like high Q channel is of quasi-BICspinned with the true BICs on the mirror axes.In the above argument, the polarization ellipse in two-dimensional (2D) momentum space is for the complex 2D4 0.52 0.53 0.54 0.55 0  0.05  0.1A2E2A1(a)ωka/2πcka/2π (ΓK)even odd  0 1 2 3 4 5 0  0.05  0.1(☓10-3)(b)γ ka/2πcka/2π (ΓK)FIG. 4. Real and imaginary photonic band structure origi-nating from symmetry-protected BICs at the Γ point. Thez-parity is odd. The momentum is taken along the ΓK di-rection. The sphere radius is fixed as rs = 0.38a. The otherparameters are the same as in Fig. 1.-0.1 0 0.1kxa/2π-0.1 0 0.1kya/2π 0.52 0.53ωka/2πcFIG. 5. Real band structure originating from the E1 mode atthe Γ point near the accidental BIC of BIC1 in Fig. 1. Thesphere radius is 0.38a. The band structure with |ka/2π| ≤0.15 is plotted. The other parameters are the same as in Fig.1.polarization vector e+ defined bye+ = t+P0k̂ + t+s0k̂⊥, (2)E±(x) =∑g(t±PgP±g + t±SgSg)eiK±g ·x, (3)P±g = ±Γgq0k̂g − |kg|q0ẑ, Sg = (k̂g)⊥, (4)K±g = kg ± Γg ẑ, Γg =√q20 − (kg)2, (5)kg = k + g, q0 =ωc. (6)where E±(x) is the electric field of the photonic bandmode above (superscript +) and below (-) the monolayer(see Appendix B), k is a 2D Bloch momentum, and g isa 2D reciprocal lattice. In the far field, solely the g =0 component survives, so that the far-field polarizationis determined by t+P0 and t±S0. In addition, the vortexcharge q is defined byq =12π∮dθk, (7)where θk is the argument of the long axis of the polariza-tion ellipse, and the integration contour orbits the vortexcore.Around the BIC2 in Fig. 1, we have similar behav-ior in the ring-like quasi-BICs pinned with multiple trueBICs on the mirror axis. However, a slightly differentbehavior than in Fig. 6 is observed. Figure 8 showsthe evolution of the multiple BICs through the criticalradius of BIC2. Across the critical radius, the ring-likeBIC moves from the upper band to the lower band. Atrs = 0.452a(> rc), there are extra BICs on the ΓM axesother than the ring-like BICs of the accidental at-Γ BICorigin. The latter ring is deformed from a circular shapeto a hexagonal shape by a substantial perturbation. Ifwe further increase the radius, the extra BICs merge withthe true BICs on ΓM in the ring, and the ring is furtherdeformed.The BIC3 in Fig. 1 is found at rs = 0.495a, near theclose-packing condition of the triangular lattice. Withinthe photonic KKR method employed in this paper, thereis not enough parameter space to enlarge the sphere ra-dius beyond the close-packing condition. We do not con-sider this BIC here.The BIC4 has a relatively high resonance frequency,giving rise to a more complex multiple off-Γ BIC distri-bution around the critical radius. Figure 9 shows theevolution of the multiple off-Γ BICs through the criticalradius of BIC4. Now the ring is formed for the lowerband at rs = 0.45a(> rc). The upper band exhibits ad-ditional BICs on ΓM. This ring pinned with true BICson ΓM and ΓK, and the additional BICs shrink to theΓ point as rs → rc. Then, at rs = 0.44a(< rc), the ringdisappears in the lower band, and discrete BICs are gen-erated from the Γ point of the upper band. The totalvorticity is conserved across the critical radius.5(a)-0.10.00.1k ya/2π104 106 108 1010 1012Qk(c)-0.1 0.0 0.1kxa/2π-0.10.00.1k ya/2π(b)(d)-0.1 0.0 0.1kxa/2π0.0972 0.0982-0.0020.002-0.002 0.0020.09860.09960.095 0.096-0.0020.0020.09780.0988-0.002 0.002q=+1q=-1q=+1q=-1FIG. 6. Q-value map of the two bands originating from the E1 mode at the Γ point near the accidental BIC of BIC1 in Fig.1. (a,c) The upper (a) and lower (c) bands (in the real band structure) at rs = 0.38a(< rc). The real band structure at thisparameter is shown in Fig. 5. (b,d) The upper (b) and lower (d) bands at rs = 0.388a(> rc). Insets show the close-up viewsoverlaid by the polarization ellipse map of the photonic band with off-Γ BICs. The vortex charge is denoted by q.Γ KMq=-1q=+1FIG. 7. Schematic illustration of the vortex charge diagramin the ring. The charge distribution is common between thetwo rings of Fig. 6.The BIC5 also have a high resonance frequency, result-ing in a complex behavior of multiple off-Γ BICs affectedby a nearby band of a higher-order vortex charge at theΓ point. We do not consider this BIC here.III. SQUARE LATTICEA similar design of the accidental at-Γ BIC and subse-quent multiple off-Γ BIC generation is available for thesquare-lattice systems with the C4v point group. A fine-tuning of a system parameter results in the accidentalat-Γ BIC of a doubly degenerate E mode at the Γ point.Above and below the critical parameter, multiple off-ΓBICs emerge in the two bands originating from the Emode at Γ. The eigenmodes of A1, A2, B1, and B2 repre-sentations of C4v at the Γ point are symmetry-protectedBICs provided that there are no open diffraction channelsother than the specular one.Figure 10 shows the design of the accidental at-Γ BICin the square lattice of identical spheres. By changingthe sphere radius, we can find three critical radii of theaccidental at-Γ BICs at rs ≃ 0.415a, 0.485a, and 0.49afor the z-even, z-odd, and z-even parities, respectively.Considering the parameter regions around the criticalradius rc, we have multiple off-Γ BICs. Figure 11 shows6upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2104 106 1081010Qkupper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.445a rc 0.452a 0.455a-0.2 0 0.2 -0.25 0 0.25 -0.2 0 0.2FIG. 8. Evolution of multiple off-Γ BICs around the criticalradius rc ≃ 0.45a of BIC2 in Fig. 1. The upper two rows showthe Q value map of the upper and lower bands originatingfrom the E1 mode at the Γ point. The lower row shows thevortex charge diagram of the off-Γ BICs. The inset shows thereal band structure around the critical radius.upper bandlower bandBIC vorticityq=+1q=-1rs=0.44a rc 0.45a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.44a rc 0.45a-0.2 0 0.2 -0.2 0 0.2104 106 1081010Qkupper bandlower bandBIC vorticityq=+1q=-1rs=0.44a rc 0.45a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.44a rc 0.45a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.44a rc 0.45a-0.2 0 0.2 -0.2 0 0.2FIG. 9. Evolution of multiple off-Γ BICs around the criticalradius rc ≃ 0.445a of BIC4 in Fig. 1. The upper two rowsshow the Q value map of the upper and lower bands originat-ing from the E1 mode at the Γ point. The lower row showsthe vortex charge diagram of the off-Γ BICs. The inset showsthe real band structure around the critical radius. 0.4 0.5 0.6 0.7 0.35  0.4  0.45  0.5(☓10-3)(a)BIC1BIC2BIC3ω0a/2πcrs/a0.00.20.40.60.81.0 0.35  0.4  0.45  0.5(b)γ 0a/2πcrs/az-even    z-odd FIG. 10. Resonance angular frequencies ω0 (a) and widths γ0(b) of the E modes at the Γ point in the square-lattice mono-layer of dielectric spheres with the C4v point-group symmetry,as a function of sphere radius rs. The modes are classified ac-cording to the parity in the z direction, assuming the centersof the spheres are on the z = 0 plane. The dielectric constantof the spheres is 12, and that of the background material is 1.The lattice constant is denoted by a. Accidental at-Γ BICsoccur around the modes indicated by the solid circles. Insetin (a) stands for the geometry of the monolayer.the real and imaginary photonic band structures of thez-even parity around the critical radius of BIC1 in Fig.10. As in Fig. 3, the real band structure does not changeso much, while the imaginary band structure changes re-markably by changing the sphere radius. Concerning theΓM direction, a similar trend as in Fig. 3 is observed.Namely, above and below the critical radius, the bandthat exhibits the off-Γ BIC at a nonzero momentum is in-terchanged. Below the critical radius [Fig. 11(a,b)], theodd-parity band exhibits the BIC while above the crit-ical radius [Fig. 11(e,f)], the even-parity band exhibitsthe BIC in ΓM. The suppression of the decay rate atabout the critical radius is weaker in ΓM but is strongerin ΓX. Such an anisotropy is manifest also in the realband structure. As for the ΓX direction, both the even-and odd-parity bands exhibit the off-Γ BICs below thecritical radius, while the BICs are absent above the crit-ical radius.By scanning all the directions around the Γ point, thereal band structure is shown in Fig. 12. The anisotropy isevident in the band structure. The upper band is S-like,and the lower band is P-like.Figure 13 shows the Q-value maps of the two photonicbands originating from the E mode at the off-criticalradii of Fig. 11. Now, the ring-like BICs found in thetriangular-lattice system are absent. A similar pattern as70.4750.4800.4850.4900.495(☓10-5)(a)ωka/2πcevenodd 0 1 2 3 4 5Γ XM(b)γ ka/2πc0.4700.4750.4800.4850.490(c)ωka/2πc 0 1 2 3 4 5(d)γ ka/2πc0.4650.4700.4750.4800.4850.2 (X←) Γ (→M) 0.2(e)ωka/2πcka/2π 0 1 2 3 4 50.2 (X←) Γ (→M) 0.2(f)γ ka/2πcka/2πFIG. 11. Real and imaginary photonic band structure of thez-even parity, originating from the E mode around the criticalradius rc of BIC1 in Fig. 10. The momentum is scanned inthe ΓX and ΓM directions. The photonic bands are furtherclassified according to the parity in these directions. Thesphere radius is 0.41a (a,b), 0.415a(≃ rc) (c,d), and 0.42a(e,f). The other parameters are the same as in Fig. 10. Insetin (b) represents the first Brillouin zone of the square lattice.in Fig. 13 (a) was observed in Ref. [16]. At rs = 0.41a,the total vorticity vanishes in the upper band. However,the lower band sustains the additional BICs on the ΓXaxes. The net vorticity of the upper and lower bands isthus nonzero. All the BICs are discrete and move to-ward the Γ point as rs → rc. At rs = rc, the BICs arecollapsed there. Across the critical coupling, multiple off-Γ BICs are again generated now on the ΓM axes. Thetotal vorticity is conserved across the critical coupling.As for the BIC2 of Fig. 10, the evolution of the multi-ple off-Γ BICs across the critical radius is shown in Fig.14. Since the two bands intersect in ΓM, the Q-valuebecomes singular there. At rs = 0.48a, the upper bandexhibits the eight off-Γ BICs. They are elongated, andtheir polarizations are perpendicular and parallel to theΓX and ΓM axes, respectively. The lower band holds fouroff-Γ BICs on the ΓM axes. All the BICs moves towardthe Γ point as rs → rc. After the collapse, the four off-ΓBICs are generated as in BIC1.The BIC3 exists in a close vicinity of another E1 mode,-0.2-0.1 0 0.1 0.2kxa/2π -0.2-0.1 0 0.1 0.2kya/2π 0.48 0.49ωka/2πcFIG. 12. Real band structure originating from the E modeat the Γ point near the accidental BIC of BIC1 in Fig. 10.The sphere radius is 0.41a. The z parity is even. The bandstructure with |ka/2π| ≤ 0.2 is plotted. The other parametersare the same as in Fig. 10.resulting in mixed four bands. We do not consider thisBIC here.IV. k · p PERTURBATION THEORYLet us consider the phenomena obtained in Secs. IIand III in the k · p perturbation. The eigenmodes in themonolayer are determined by∇×(1ε(x)∇×H(x))=ω2c2H(x). (8)Here, H is the time-harmonic magnetic field with angu-lar frequency ω, and ε(x) is the dielectric function thatis periodic in the in-plane coordinate x∥. The inversedielectric function can be expanded as1ε(x)=∑geig·x∥ηg(z). (9)In addition, the radiation field is also expanded by theplane waves via Bloch’s theorem asH(x) =∑gei(k+g)·x∥hg(z). (10)Accordingly, the equation to be solved becomes∑g′(i(k + g) + ẑ∂z)× [ηg−g′(z)(i(k + g′) + ẑ∂z)× hg′(z)] =ω2c2hg(z),(11)8q=-1q=+1q=-1q=-1(a)-0.2-0.10.00.10.2k ya/2π104 106 108 1010Qk(c)-0.2 -0.1 0.0 0.1 0.2kxa/2π-0.2-0.10.00.10.2k ya/2π(b)(d)-0.2 -0.1 0.0 0.1 0.2kxa/2π0.153 0.163-0.0050.0050.04 0.050.040.050.160 0.170-0.0050.0050.04 0.050.040.05FIG. 13. Q-value map of the two bands originating from the E mode at the Γ point near accidental BIC of BIC1 in Fig. 10.(a,c): the upper (a) and lower (c) bands at rs = 0.41a(< rc). (b,d): the upper (b) and lower (d) bands at rs = 0.42a(> rc).The insets show the close-up view of the Q-value map overlaid with the polarization ellipse map.which is symbolically expressed asH|ψ⟩ = E|ψ⟩, E =ω2c2. (12)Suppose we have a nearly accidental BIC of a doublydegenerate eigenmode at the Γ point. We perform thek ·p perturbation starting from these states as the zerothorder ones. The effective Hamiltonian to be diagonalizedis given by [29]Heffab = ⟨ϕ0a|H(2)|ϕ0b⟩+∑n ̸=01ϵ0 − ϵn⟨ϕ0a|H(1)|ϕn⟩⟨ϕn|H(1)|ϕ0b⟩, (13)H(0)|ϕ0a⟩ = ϵ0|ϕ0a⟩, H(0)|ϕn⟩ = ϵn|ϕn⟩, (14)where H(i) is the operator in the left-hand side of Eq.(11) expanded by Bloch momentum k (i = 0, 1 and 2 im-ply k = 0, k-linear, and k-quadratic, respectively). Thestate |ϕ0a⟩ (a = 1, 2) represents the doubly degeneratenearly accidental BIC mode at Γ. Here, we assume thedegeneracy is not lifted in the first-order perturbation.That is,⟨ϕ0a|H(1)|ϕ0b⟩ = 0. (15)This property is verified explicitly by the C6v or C4v sym-metry. More generally, the matrix element ⟨ϕm|H(1)|ϕn⟩is nonzero only if (ϕm, ϕn) is attributed to (E1, E2),(E1, A1), (E1, A2), (E2, B1), (E2, B2), and vice versa forC6v, and to (E,A1), (E,A2), (E,B1), (E,B2), and viceversa for C4v (see Appendix A).The sum in the second term of Eq. (13) generally in-cludes the continuous radiation modes and discrete quasi-guided modes other than the mode concerned, as the in-termediate states. The continuous radiation modes givethe decay rate to the original eigenmode via Fermi’sgolden rule [30, 31]. However, if we start from themode of the nearly accidental at-Γ BIC, the decay rateis strongly suppressed.The continuous radiation modes belong to the sameirreducible representation as the original doubly degen-erate mode below the diffraction threshold. Therefore,9upper bandlower bandBIC vorticityq=+1q=-1rs=0.48a rc 0.495a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.48a rc 0.495a-0.2 0 0.2 -0.2 0 0.2104 106 1081010Qkupper bandlower bandBIC vorticityq=+1q=-1rs=0.48a rc 0.495a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.48a rc 0.495a-0.2 0 0.2 -0.2 0 0.2upper bandlower bandBIC vorticityq=+1q=-1rs=0.48a rc 0.495a-0.2 0 0.2 -0.2 0 0.2FIG. 14. Evolution of multiple off-Γ BICs around the criticalradius rc ≃ 0.485a of BIC2 in Fig. 10. The upper two rowsshow the Q-value map of the upper and lower bands originat-ing from the E mode at the Γ point. The lower row showsthe vortex charge diagram of the off-Γ BICs. The inset showsthe real band structure (at rs = 0.495a).the matrix elements vanish between the nearly acciden-tal at-Γ BIC and radiation modes. That is why the decayrate is strongly suppressed.In the C6v case, the effective hamiltonian is expressedasHeff = αk21 + β((k2x − k2y)σ3 + 2kxkyσ1), (16)α = α(2)E1E1+∑n∈E2|α(1)E1E2|2ϵ0 − ϵn+12∑n∈A1|α(1)E1A1|2ϵ0 − ϵn+12∑n∈A2|α(1)E1A2|2ϵ0 − ϵn, (17)β = β(2)E1E1+12∑n∈A1|α(1)E1A1|2ϵ0 − ϵn− 12∑n∈A2|α(1)E1A2|2ϵ0 − ϵn(18)where σi(i = 1, 2, 3) is the Pauli matrix and parame-ters α(2)E1E1and so on are given in Appendix A. The re-sulting effective parameters α and β are complex withIm[α] ≥ 0. Their imaginary parts come solely fromthe one in the eigenfrequency of the E1 mode, namely,ϵ0 = (ω0−iγ0)2/c2, provided that the intermediate modesare lossless (Im[ϵn] = 0). Therefore, if we start with theE1 mode of the perfect accidental at-Γ BIC (γ0 = 0),Im[α] = Im[β] = 0 so that we still have a vanishing de-cay rate even at finite k, within the second-order k · pperturbation. This is the case of the extreme suppressionof the decay rate around the Γ point of the accidental(next-to-super) BIC. Since the third-order terms vanishby the inversion symmetry, the imaginary frequencies ofthe two bands behave k4 as pointed out in Ref. [16].By diagonalizing the effective Hamiltonian, the eigen-frequencies of the two bands that stem from the E1 modeat Γ are given byω2c2=(ω0 − iγ0)2c2+ (α± β)k2, (19)H+k (x) ∝ kxH(01)E1(x) + kyH(02)E1(x), (20)H−k (x) ∝ −kyH(01)E1(x) + kxH(02)E1(x). (21)where (H(01)E1,H(02)E1) forms the E1 representation at Γ.If γ0 ≪ ω0, the decay rate γk is evaluated asγk ≃ γ0 −c22ω0Im[α± β]k2. (22)The radiation field of H+(−)k is P-like (S-like). The prop-erty of Im[α] > 0 promises a monotonic decreasing of thedecay rate γk with increasing |k| for at least one band,because either Im[α + β] or Im[α − β] is promised to bepositive. The behavior of the other band depends on therelative magnitude of Im[α] and Im[β]. Across Im[β] = 0,the band with fast decreasing decay rate with |k| is inter-changed. This condition is satisfied at the critical cou-pling. There, γ0 = 0 and thus Im[β] = 0. Moreover, sincethe dispersion relation is isotropic, we can have ring-likeBICs at the minimum of γk. These results explain Fig.6 reasonably well.More direct and quantitative comparison between theKKR calculation and the k · p perturbation needs spe-cial care. The latter requires a regularization of thenon-normalizable (due to the radiation loss) quasi-guidedmodes of the E1 representation. Therefore, the param-eters α(1)E1A1and so on become regularization-dependent.To avoid such an artifact, we consider the fitting of thereal and imaginary band structures around the Γ pointby complex parameters of α and β. The result is shownin Fig. 15. Since the k · p perturbation is valid onlyfor small |k|, it can describe the trend there reasonablywell. In particular, the fitting of the real band structureis perfect. However, the fitting of the imaginary bandstructure is not good enough. In the second-order k · pperturbation, the lower bound of γk (γk ≥ 0) is ignored,so that the fitted γk can go down beyond the lower bound.This is a drawback of the present approach.As for Figs. 8 and 9, the spatial anisotropy is verystrong, so that the k · p perturbation that predicts theisotropic dispersion is not so efficient.10 0.522 0.523 0.524 0.525 0.526 0  0.05  0.1(a)ωka/2πcka/2π (ΓK)KKR, evenKKR, oddk·p, evenk·p, odd0.00.51.01.52.0 0  0.05  0.1(b)(☓10-5)γ ka/2πcka/2π (ΓK)FIG. 15. Comparison in the complex band structures of Fig.3 (a,b) between the KKR result and the k · p perturbationof Eq. (19) with a numerical fitting. Fitted values are α =0.111 + i9.56× 10−5 and β = −0.132 + i1.20× 10−3.In the C4v case, the effective Hamiltonian becomesHeff = αk21 + β1(k2x − k2y)σ3 + 2β2kxkyσ1, (23)α = α(2)EE+12∑n∈A1|α(1)EA1|2ϵ0 − ϵn+12∑n∈B1|α(1)EB1|2ϵ0 − ϵn+12∑n∈A2|α(1)EA2|2ϵ0 − ϵn+12∑n∈B2|α(1)EB2|2ϵ0 − ϵn, (24)β1 = β(2)1:EE+12∑n∈A1|α(1)EA1|2ϵ0 − ϵn+12∑n∈B1|α(1)EB1|2ϵ0 − ϵn− 12∑n∈A2|α(1)EA2|2ϵ0 − ϵn− 12∑n∈B2|α(1)EB2|2ϵ0 − ϵn, (25)β2 = β(2)2:EE+12∑n∈A1|α(1)EA1|2ϵ0 − ϵn− 12∑n∈B1|α(1)EB1|2ϵ0 − ϵn− 12∑n∈A2|α(1)EA2|2ϵ0 − ϵn+12∑n∈B2|α(1)EB2|2ϵ0 − ϵn. (26)By diagonalizing the effective Hamiltonian, the eigenfre-quency becomesω2c2=(ω0 − iγ0)2c2+ αk2 ±√β21(k2x − k2y)2 + 4β22k2xk2y.(27)Again, Im[α] > 0, promising a monotonic decreasing ofthe imaginary part in Eq. (27) with increasing |k| for atleast one band branch. Now, Eq. (27) is not isotropic.Thus, the ring-like high Q channel is not formed. Thistrend is entirely consistent with the numerical results ob-tained in Sec. III.On the ΓX axis (ky = 0), the two eigenstates becomeω2c2=(ω0 − iγ0)2c2+ (α± β1)k2x, (28)H+k (x) = H(01)E (x), H−k (x) = H(02)E (x), (29)where (H(01)E ,H(02)E ) forms the E representation at Γ.Since the E representation behaves as (x, y), H+k is P-polarized and H−k is S-polarized. One of the eigenmodesis promised to have a decreasing decay rate with |kx|, andthe two eigenmodes are interchanged at Im[β1] = 0 viathe critical coupling condition Im[γ0] = 0.On the ΓM axis (kx = ky), the eigenstates becomeω2c2=(ω0 − iγ0)2c2+ 2(α± β2)k2x, (30)H±k (x) ∝ H(01)E (x)±H(02)E (x). (31)The eigenmode of superscript ”+” is P-polarized, and ”-” is S-polarized. One of the two eigenmodes is promisedto have a decreasing decay rate with |kx|, and the twoeigenmodes are interchanged at Im[β2] = 0 via the criti-cal coupling condition Im[γ0] = 0. Again, these featuresare consistent with Figs. 13 and 14.The symmetry-protected BICs are also available forC3v point-group systems. There, the eigenmodes of theA1 and A2 representations at the Γ point are symmetry-protected, provided there are no open diffraction chan-nels other than the specular one. In this case, it is pos-sible to have accidental at-Γ BICs for the eigenmodes ofthe E representation by tuning system parameters.The effective Hamiltonian around the nearly accidentalat-Γ BIC mode is given byHeff = αk21 + β((k2x − k2y)σ3 + 2kxkyσ1), (32)α = α(2)EE +∑n∈E′ |α(1)EE |2ϵ0 − ϵn+12∑n∈A1|α(1)EA1|2ϵ0 − ϵn+12∑n∈A2|α(1)EA2|2ϵ0 − ϵn, (33)β = β(2)EE+12∑n∈A1|α(1)EA1|2ϵ0 − ϵn− 12∑n∈A2|α(1)EA2|2ϵ0 − ϵn, (34)where the prime in the second term of Eq. (33) means theunperturbed E mode is excluded in the sum. In this case,a significant contribution to the decay rate emerges fromthe continuous radiation modes of the E representationvia Fermi’s golden rule. Therefore, suppressing of thedecay rate at nonzero |k| near the accidental BIC (at Γ)does not occur. Consequently, the multiple off-Γ BICsare absent in the k · p perturbation theory.11In the above arguments, we do not rely on the mono-layer of spheres presented in Secs. II and III, but onthe spatial symmetry of C6v, C4v, and C3v. Accordingly,the multiple BIC generation with the present scenario isavailable in other photonic membranes with C6v or C4v.We can show PhC slabs with the triangular or squarelattice of circular air holes exhibit multiple BIC genera-tions. The ring-like high Q channel is also formed in theC6v case.V. SUMMARYIn summary, we have presented a detailed theoreticalanalysis of the multiple BIC generation in the monolay-ers of spheres with C6v or C4v point group. A tuningof system parameters results in an accidental at Γ BICof doubly degenerate eigenmodes at the Γ point. ThisBIC, obtained at a critical parameter, is the next-to-super BIC with extremely suppressed decay rates aroundthe Γ point. Off-critical parameters yield multiple off-ΓBICs that are generated from the next-to-super BIC atthe Γ point.In the C6v point-group system, a ring-like high Q chan-nel pinned with multiple off-Γ BICs on the mirror axescan be formed in the two bands originating from thedoubly degenerate E1 mode at the Γ point. The spa-tial anisotropy destroys the ring in the C4v point-groupsystem, but multiple BICs are certainly formed. Acrossthe critical parameters, these BICs move from the upperband to the lower band or vice versa, conserving the totalvorticity of the BICs.We also show that these phenomena are owing to thespatial symmetry of C6v or C4v and the k ·p perturbationexplains the phenomena reasonably well. In addition, theC3v system does not support the multiple BIC generationthrough a similar design.ACKNOWLEDGMENTSThis work was supported by JSPS KAKENHI GrantNo. 22K03488.Appendix A: Matrix elements in the k · pperturbationWe summarize various matrix elements relevant tothe effective Hamiltonian. The symmetry relation underpoint-group operations is crucial. It is given byH(i)R1R2(k) = D†R1(A)H(i)R1R2(Ak)DR2(A), (A1)where the matrix element is the one between the modesof irreducible representations R1 and R2, and DR(A) isthe representation matrix of group element A in the irre-ducible representation R. The equation (A1) constrainsthe possible form of the matrix elements as follows.In the C6v case, The nonzero matrix elements in thefirst-order k · p perturbation are given byH(1)E1E2(k) = α(1)E1E2(kxσ1 + kyσ3), (A2)H(1)E1A1(k) = α(1)E1A1(kxky), (A3)H(1)E1A2(k) = α(1)E1A2(−kykx), (A4)H(1)E2B1(k) = α(1)E2B1(−kykx), (A5)H(1)E2B2(k) = α(1)E2B2(kxky), (A6)where α(1)E1E2and so on are complex. The vanishing ma-trix elements between E1 and B1(B2) are crucial for thek4 scaling of the decay rate in the symmetry-protectedBICs of the B1 and B2 representations, which have ahigher order vortex charge of q = −2.In the second order, the relevant matrix element be-comesH(2)E1E1(k) = α(2)E1E1k2 + β(2)E1E1((k2x − k2y)σ3 + 2kxkyσ1),(A7)where α(2)E1E1and β(2)E1E1are real.In the C4v case, the nonzero matrix elements in thefirst-order k · p perturbation are given byH(1)EA1(k) = α(1)EA1(kxky), (A8)H(1)EA2(k) = α(1)EA2(−kykx), (A9)H(1)EB1(k) = α(1)EB1(kx−ky), (A10)H(1)EB2(k) = α(1)EB2(kykx), (A11)where α(1)EA1and so on are complex. In the second order,the relevant matrix element becomesH(2)EE(k) = α(2)EEk21 + β(2)1:EE(k2x − k2y)σ3 + β(2)2:EE2kxkyσ1,(A12)where α(2)E1E1, β(2)1:E1E1, β(2)2:E1E1are real.In the C3v case, the nonzero matrix elements in thefirst-order k · p perturbation are given byH(1)EE(k) = α(1)EE(kxσ1 + kyσ3), (A13)H(1)EA1(k) = α(1)EA1(kxky), (A14)H(1)EA2(k) = α(1)EA2(−kykx), (A15)where α(1)EE and so on are complex. The matrix element ofEq. (A13) vanishes by the time-reversal symmetry if the12two E modes that form the matrix element are identical.In the second order, the relevant matrix element becomesH(2)EE(k) = α(2)EEk2 + β(2)EE((k2x − k2y)σ3 + 2kxkyσ1),(A16)where α(2)EE and β(2)EE are real.The parameters α(1,2)R1R2and β(2)R1R2are not determinedsolely from the symmetry relation [Eq. (A1)].Appendix B: Electric fields of quasi-guided modesWe employ the following method for evaluating the po-larization ellipse map. At a true BIC point, the kernelof a secular matrix gives the eigenmode, and its determi-nant vanishes on the real frequency axis.In the monolayer of spheres, the secular matrix S isgiven by [23, 32]S(Lβ)(L′β′) = δLL′δββ′ − tβl Gββ′LL′ , (B1)Gββ′LL′ =1l(l + 1)∑L1L2[(P β)†]LL1GL1L2· [P β′]L2L′ , (B2)GLL′ = 4π∑L′′il−l′−l′′⟨L|L′′|L′⟩SL′′ , (B3)SL =∑X ̸=0h(1)l (q0|X|)Y ∗L (X̂)eik·X . (B4)where L = (l,m) (|m| ≤ l) is the angular momentum in-dex, β(= M,N) is the index to classify two transversevector spherical waves, tβl is the so-called t-matrix ofthe sphere, P β is the transformation matrix from scalarspherical waves to vector spherical waves, ⟨L|L′′|L′⟩ isthe Clebsh-Gordan coefficient, h(1)l is the spherical hankelfunction of the first kind, YL is the spherical harmonics,and X is the 2D real-lattice vector. The structure con-stant SL is efficiently calculated via the Ewald method[33]. Through ψβL = Ker[S], the electric field outside themonolayer is expressed asE(x) =∑LXh(1)l (q0|x−X|)YL( ̂x−X)V Leik·X , (B5)V L =∑L′βP βLL′ψβL′ . (B6)Its plane-wave-expansion coefficient t±σg (σ = P, S) be-comest±σg =2πq0ΓgAUC∑L(−i)lYL(K̂±g )σ±g · V L, (B7)where AUC is the area of the unit cell.Off the BIC point, the determinant is nonzero on thereal frequency axis but can be zero in the complex fre-quency plane. This zero should be the resonance pole ofthe S-matrix, namely, Ωk = ωk−iγk. We assume that theright eigenstate of S with the least absolute eigenvaluecan be approximated as the kernel on the real frequencyaxis near the zero. Then, we can obtain the plane-wavecoefficients t±σg(ω) through Eq. (B7). We further assumethat it can be Taylor expanded around ω = Ωk, namely,t±σg(ω) ≃ t±σg(Ωk) + (ω − Ωk)C±σg (σ = P, S). (B8)Then, from the two points ωa and ωb on the real fre-quency axis, we can estimate t±σg(Ωk) ast±σg(Ωk) =(ωa − Ωk)t±σg(ωb)− (ωb − Ωk)t±σg(ωa)ωa − ωb. (B9)We confirm that the resulting t±σ0(Ωk) vanishes at theBIC points.[1] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos,M. Soljacic, and M. Soljačić, Nat. Rev. Mater. 1, 16048(2016).[2] A. Kodigala, T. Lepetit, Q. Gu, B. Bahari, Y. Fainman,and B. Kanté, Nature 541, 196 (2017).[3] S. Romano, G. Zito, S. Torino, G. Calafiore, E. Penzo,G. Coppola, S. Cabrini, I. Rendina, and V. Mocella,Photonics Res. 6, 726 (2018).[4] K. Koshelev, S. Kruk, E. Melik-Gaykazyan, J. H. Choi,A. Bogdanov, H. G. Park, and Y. Kivshar, Science 367,288 (2020).[5] H. Miyazaki and K. Ohtaka, Phys. Rev. B 58, 6920(1998).[6] P. Paddon and J. F. Young, Phys. Rev. B 61, 2090(2000).[7] T. Ochiai and K. Sakoda, Phys. Rev. B 63, 125107(2001).[8] S. H. Fan and J. D. Joannopoulos, Phys. Rev. B 65,235112 (2002).[9] C. W. Hsu, B. Zhen, J. Lee, S. L. Chua, S. G. Johnson,J. D. Joannopoulos, and M. Soljačić, Nature 499, 188(2013).[10] Q. Jiang, P. Hu, J. Wang, D. Han, and J. Zi, Phys. Rev.Lett. 131, 013801 (2023).[11] B. Zhen, C. W. Hsu, L. Lu, A. D. Stone, and M. Soljačić,Phys. Rev. Lett. 113, 257401 (2014).[12] J. Jin, X. Yin, L. Ni, M. Soljačić, B. Zhen, and C. Peng,Nature 574, 501 (2019).[13] M.-S. S. Hwang, H.-C. C. Lee, K.-H. H. Kim, K.-Y. Y.Jeong, S.-H. H. Kwon, K. Koshelev, Y. Kivshar, andH.-G. G. Park, Nat. Commun. 12, 4135 (2021).[14] M. Kang, L. Mao, S. Zhang, M. Xiao, H. Xu, and C. T.Chan, Light Sci. Appl. 11 (2022).[15] M. Kang, S. Zhang, M. Xiao, and H. Xu, Phys. Rev.Lett. 126, 117402 (2021).[16] A. S. Kostyukov, V. S. Gerasimov, A. E. Ershov, and13E. N. Bulgakov, Phys. Rev. B 105, 075404 (2022).[17] Z. Sadrieva, K. Frizyuk, M. Petrov, Y. Kivshar, andA. Bogdanov, Phys. Rev. B 100, 115303 (2019).[18] S. Gladyshev, A. Shalev, K. Frizyuk, K. Ladutenko, andA. Bogdanov, Phys. Rev. B 105, L241301 (2022).[19] T. Yoda and M. Notomi, Phys. Rev. Lett. 125, 053902(2020).[20] C. F. Doiron, I. Brener, and A. Cerjan, Nat. Commun.13, 7534 (2022).[21] X. Wang, J. Wang, X. Zhao, L. Shi, and J. Zi, ACSPhotonics 10, 2316 (2023).[22] T. Inui, Y. Tanabe, and Y. Onodera, Group Theory andits Applications in Physics (Springer, 1996).[23] K. Ohtaka, J. Phys. C Solid State Phys. 13, 667 (1980).[24] A. Modinos, Physica A 141, 575 (1987).[25] N. Stefanou, V. Karathanos, and A. Modinos, J. Phys.Condens. Matter 4, 7389 (1992).[26] L. D. Landau and E. M. Lifshitz, Quantum Mechanics:Non-relativistic Theory (Pergamon, Oxford, 1977).[27] K. Ohtaka, J. Inoue, and S. Yamaguti, Phys. Rev. B 70,35109 (2004).[28] U. Fano, Phys. Rev. 124, 1866 (1961).[29] T. Ochiai, Phys. Rev. B 86, 075152 (2012).[30] T. Ochiai and K. Sakoda, Phys. Rev. B 64, 045108(2001).[31] L. C. Andreani and D. Gerace, Phys. Rev. B 73, 235114(2006).[32] K. Ohtaka, Phys. Rev. B 19, 5057 (1979).[33] K. Kambe, Zeitschrift fur Naturforsch. - Sect. A J. Phys.Sci. 22, 322 (1967).