Generation of multiple bound states in the continuum through doubly degenerate quasi-guided modes
Tetsuyuki Ochiai

 


Research Center for Electronic and Optical Materials, National Institute for Materials Science (NIMS), Tsukuba 305-0044, Japan



We present a detailed theoretical analysis of a peculiar generation of multiple bound states in the continuum (BICs) in two-dimensional periodic arrays of dielectric spheres. 
They emerge in high-symmetry lattices with the C_6v and C_4v point groups and involve doubly degenerate quasi-guided modes at the Γ point that can couple to external radiation. 
By tuning a system parameter,  the doubly degenerate modes can exhibit accidental BICs at a critical parameter. 
In the vicinity of the critical 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 C_6v case. Across the critical coupling, the total vorticity of the multiple BICs is conserved.   The *k·*p perturbation theory explains some features of the phenomena reasonably well. 






Generation of multiple bound states in the continuum through 
doubly degenerate quasi-guided modes
    Tetsuyuki Ochiai
    April 29, 2024
===================================================================================================





§ INTRODUCTION


Optical bound states in the continuum (BIC) are localized eigenstates embedded in the radiation continuum <cit.>.  They have  infinite quality factors or, in other words, vanishing decay rates,   although they are inside the light cone. 
Strong light confinement via the infinite quality factor enables us to investigate various applications such as lasing <cit.>, sensing <cit.>, nonlinear optics <cit.>, and so on, via BICs. 

The BICs are found typically in monolayers of spheres and photonic crystal slabs at the Γ point  <cit.>.  Off-Γ BICs are also  found at generic  *k  points  <cit.>, via the formation of polarization vortices <cit.>. The former BICs are symmetry-protected and irrelevant to physical parameters unless the relevant symmetry is unchanged. The latter are  topologically protected, moving in momentum space.  Therefore, a parameter scan is necessary to find the latter BICs at a given *k point other than Γ. In this sense, the off-Γ BICs are sometimes called accidental. 


Recently, merging of an at-Γ
 BIC and off-Γ BICs in a non-degenerate isolated band 
attracts much interest as the so-called super BIC <cit.>.  
This BIC is obtained by tuning system parameters such that the off-Γ BICs move toward the Γ point, where 
the symmetry-protected BIC exists.  
The super BIC has superb properties as it involves an extreme suppression of the decay rate in a broad region of the momentum space around Γ. The decay rate  behaves as |*k|^6, in  striking contrast to the |*k|^2 behavior of the ordinary symmetry-protected BICs. 


Another type of super BICs has been found recently in monolayers of spheres <cit.>. From now on, we call it the next-to-super BIC.  It is obtained at the Γ point by tuning system parameters and exhibits a |*k|^4 scaling of the decay rate at a critical parameter. In contrast to the ordinary super BIC, it 
involves doubly-degenerate eigenmodes at the Γ point. 
Moreover, off-critical parameters result in a ring-like high Q channel in the triangular lattice system. This high Q channel is either P or S-polarization-like and can be switched by changing the parameter across the critical value. 


Once we have a high Q channel, we have a continuous distribution of optical modes with strongly enhanced light-matter interactions. The channel can be selectively excited by the  polarization of the incident light.    
As a result, we can steer the coherent radiation continuously along the ring.  The output light has nearly the same frequency, but the angle changes continuously.   
Such a beam steering paves the way for various applications. 




In this paper, we further investigate the next-to-super BICs and  show that multiple off-Γ BICs are generated at off-critical parameters nearby. 
We show more complex behaviors than in Ref. <cit.> can be observed. 
In addition to the ring-like channel, discrete multiple BICs can be  generated through the next-to-super BICs. In this case, the ring is formed asymmetrically in the parameter space. If we change the parameter across the critical value, the ring disappears, and 
multiple discrete BICs are again generated. Across the critical parameter, the total vorticity of the BICs is conserved. Besides, the appearance of the ring depends on the lattice structure.    

One of the points behind the above phenomena is spatial symmetry. We show that the *k·*p perturbation theory based on spatial symmetry explains some properties of the phenomena reasonably well.   
 


The ring-like high Q channel can be explained in terms of the multipolar lattice <cit.>.   However, the lattice-structure dependence needs complementary approaches to understand the phenomena. 
The *k·*p perturbation provides one reasonable scenario as shown in this paper. 

    

The multiple off-Γ BICs are also produced by breaking spatial symmetries of  photonic crystal  (PhC) slabs  <cit.>. There, degenerate or nondegenerate symmetry-protected BICs at Γ turn into multiple off-Γ BICs. Since the original symmetry-protected BICs exist irrespective of system parameters, the resulting multiple BICs are deterministic and parameter tuning is  unnecessary.  






 


This paper is organized as follows. In Secs. II and III, we present numerical examples of  the accidental at-Γ BICs and resulting multiple off-Γ BICs in the triangular and square lattices of dielectric spheres, respectively. In Sec. IV, we analyze these phenomena via the *k·*p perturbation theory of the Maxwell equation. Finally, in Sec. V, we summarize the results.  
 
  



   
 
 



§ TRIANGULAR LATTICE
 

Let us first consider the next-to-super BICs in a monolayer of dielectric spheres arranged in a triangular lattice.  The system has the C_6v point group symmetry <cit.>.  
Therefore, the eigenmodes are classified according to the irreducible representations of C_6v.

At the Γ point, the eigenmodes with the E_1 representation of C_6v can couple to external radiation. The other eigenmodes of A_1,A_2,B_1,B_2 and E_2 representations are uncoupled, provided that their eigenfrequencies are below the diffraction threshold. Thus, 
the former eigenmodes generally have finite lifetimes, whereas  the latter eigenmodes have infinite lifetimes and are symmetry-protected BICs. 



However, tuning system parameters can cause the former eigenmodes to have  infinite lifetimes. In this way, accidental BICs can occur for the E_1 modes. 


Figure <ref> shows the resonance angular frequency ω_k and width γ_k of the E_1 modes at the Γ point (k=0). 



They are evaluated with the photonic Korringa-Kohn-Rostoker (KKR) method <cit.> together with the curve fitting to the Breit-Wigner formula  <cit.> of the scattering phase shift δ as a function of angular frequency ω: 

    ^2δ=^2δ_k(ω-ω_k-γ_k/ω-ω_k+γ_k)^n,
 
where δ_k is the background (frequency-independent) phase shift, and n is the degree of degeneracy of the resonant mode concerned.  
The scattering phase shift is derived from the S-matrix of the monolayer  <cit.>. 
At several radii, the eigenmodes exhibit infinite quality factor accidentally. 





Such an accidental BIC can be easily found theoretically and experimentally by monitoring the change in the resonance signal in the  transmission spectrum of the normal incidence.  
For instance, Fig. <ref> shows the change of the transmission spectrum with the sphere radius around the critical one of BIC1 in  Fig. <ref>.  



Since the modes other than E_1 are symmetry-protected, they do not affect the transmission spectrum. Solely the modes of the E_1 representation emerge as an asymmetric resonance signal of the Fano shape <cit.> in the spectrum. 
The spectrum changes with the radius, and we can find a narrowing of the resonance width toward the critical radius. The vanishing resonance width corresponds to the accidental BIC.  


The above property of the next-to-super BIC presents a marked contrast to the super BIC. In the latter case, the BIC consists of a symmetry-protected BIC with topologically protected BICs, so it is impossible to observe the super BIC via the transmission spectrum of the normal incidence.  It is available only through a detailed analysis of the momentum and parameter dependence of the spectrum. 
 


  
Figure <ref> shows the real and imaginary photonic band structure around the critical radius of BIC1. 



Here, the real photonic band structure is referred to as the resonance angular  frequency ω_k, and the imaginary one is to the resonance width γ_k, as a function of Bloch momentum k.    
Although the real band structure does not change so much in its shape, the imaginary band structure  exhibits a clear contrast between the two bands of opposite parities
and among the three values of on- and off-critical radii.  
Below the critical radius [Figs. <ref> (a,b)], the odd-parity band exhibits the off-Γ BIC at about 
ka/2π=0.1. 
At the critical radius [Figs. <ref>(c,d)], both the bands exhibit 
a flat region of nearly zero values in their imaginary parts. This suppression of the decay rate resembles that in the super BIC. However, now γ_k behaves as k^4 instead of the k^6 behavior of the super BIC.    
Above the critical radius [Figs. <ref> (e,f)], the even-parity band  exhibits the off-Γ BIC. 


For comparison, the photonic bands originating from the symmetry-protected BICs at Γ exhibit a quick blow-up of the imaginary parts as shown in Fig. <ref>. 



Near the Γ point, we do not observe the trend of decreasing the imaginary part toward a minimum. 
We should point out that the imaginary part is much larger than  that  from the E_1 mode, whose imaginary part is significantly  suppressed. 





By scanning the Bloch momentum in all the directions around the Γ point, 
the real band structure consists of two surfaces touched quadratically at the Γ point, as shown in Fig. <ref>.  



The upper and lower bands are S- and P-polarization-like, respectively. 


Figure <ref> shows the Q-value map of the photonic band modes at the off-critical radii.  



It is remarkable that the ring-like high Q channels found in Ref. <cit.> are formed.  
The channels lies in the S-polarization-like upper band at r_s<r_c and in the P-polarization-like lower band for r_s>r_c. The other bands (the lower band at r_s<r_c and the upper band at r_s>r_c) do not have such channels. 
These channels indicate that the minimum decay rate γ_k at the off-critical radii found in Fig. <ref> emerges at nearly the same distance from the Γ point regardless of momentum orientation.  
As we move the radius r_s toward r_c, we can show the ring shrinks to the Γ point and across r_s=r_c, the ring moves from the upper band to the lower band or vice versa. 



Looking closely at the rings, we find twelve off-Γ BICs: six on the equivalent ΓM axes and six on the equivalent ΓK axes. 
These axes are the mirror axes of the triangular-lattice Brillouin zone. 
The polarization vortices of the BICs are very elongated along the rings.  
The vortex charge diagram is shown in Fig.  <ref>. 



The vortex charges are opposite between  ΓM and ΓK.  
Off the mirror axes, the Q values in the ring are still very high, of order 10^-8. Thus,  the ring-like high Q channel is of 
quasi-BICs pinned 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 2D polarization vector *e^+ defined by 

    *e^+=t_P*0^+ *̂k̂ + t_s*0^+ *̂k̂_⊥, 
       *E^±(*x)=∑_*g (t_P*g^±*P_*g^± + t_S*g^±*S_*g)^*K_*g^±·*x,
       *P_*g^±=±Γ_*g/q_0*̂k̂_*g-|*k_*g|/q_0ẑ,  *S_*g=(*̂k̂_*g)_⊥,
       *K_*g^±=*k_*g±Γ_*gẑ,  Γ_*g=√(q_0^2-(*k_*g)^2), 
       *k_*g=*k+*g,   q_0=ω/c.
	
where *E^±(*x) is the electric field of the photonic band mode above (superscript +) and below (-) the monolayer (see Appendix B), *k is a 2D Bloch momentum, and *g is a 2D reciprocal lattice. 
In the far field, solely the *g=0 component survives, so that the far-field polarization is determined by t_P*0^+ and t_S*0^±. 
Besides, the vortex charge q is defined by 

    q=1/2π∮θ_k,

where θ_k is the argument of the long axis of the polarization ellipse, 
and the integration contour orbits the vortex core.  












Around the BIC2 in Fig. <ref>, we have similar behavior in the ring-like quasi-BICs pinned with multiple true BICs on the mirror axis.  
However, a slightly different behavior than in Fig. <ref> is observed.   
Figure <ref> shows the evolution of the multiple BICs through the critical radius of BIC2.  



Across the critical radius, the ring-like BIC moves from the upper band to the lower band. 
At r_s=0.452a>r_c,  there are extra BICs on the ΓM axes other than the ring-like BICs of the accidental BIC origin.  The latter ring is deformed from a circular shape to a hexagonal shape by a substantial *k·*p perturbation.  
If we further increase the radius, the extra BICs merge with the true BICs on ΓM in the ring, and the ring is further deformed. 



The BIC3 in Fig. <ref> is found at r_s=0.495a,  near the close-packing condition of the triangular lattice. Within the photonic KKR method employed in this paper, there is not enough parameter space to enlarge the sphere radius beyond the close-packing condition. We do not consider this BIC here.  


The BIC4 has a relatively high resonance frequency,  giving rise to  a more complex multiple BIC distribution  around the critical radius.   
Figure <ref> shows the evolution of the multiple BICs through the critical radius of BIC4. 



Now the ring is formed for the lower band at r_s=0.45a (>r_c). The upper band exhibits additional BICs on ΓM. This ring pinned with true BICs on ΓM and ΓK, and the additional BICs shrink to the Γ point as r_s→ r_c. Then, at r_s=0.44a (<r_c), the ring disappears in the lower band, and discrete BICs are generated from the Γ point of the upper band. 
The total vorticity is conserved across the critical radius. 


 







§ SQUARE LATTICE



A similar design of the accidental BIC and subsequent multiple BIC generation is available for the square-lattice systems with the C_4v point group. A fine-tuning of a system parameter  results in the accidental BIC   
of a doubly degenerate E mode at the Γ point.  Above and below the critical parameter, multiple BICs emerge in the two bands originating from the E mode at Γ. The eigenmodes of A_1, A_2, B_1, and B_2 representations of C_4v at the Γ point are symmetry-protected BICs provided that there are no open diffraction channels other than the specular one. 



Figure <ref> shows the design of the accidental BIC in the square lattice of identical spheres. 



By changing the sphere radius, we can find two critical radii of the accidental BICs at  r_s≃ 0.415a and 0.485a for the z-even and odd parities, respectively. 


Considering the parameter regions around the critical radius r_c, we have multiple off-Γ BICs. 
Figure <ref> shows the real and imaginary photonic band structures of the z-even parity around the critical radius of BIC1 
in  Fig. <ref>.  



As in Fig. <ref>, the real band structure does not change so much, while the imaginary band structure changes remarkably by changing the sphere radius. Concerning the ΓM direction, a similar trend as in Fig. <ref> is observed. Namely, above and below the critical radius, the band that exhibits the off-Γ BIC at a nonzero momentum, is interchanged. 
 Below the critical radius [Fig. <ref>(a,b)], the odd-parity band exhibits the BIC while above the critical radius [Fig. <ref>(e,f)], the even-parity band exhibits the BIC in ΓM. 
 The suppression of the decay rate at about the critical radius is weaker in ΓM but is stronger in ΓX. Such an anisotropy is  manifest also in the real band structure. As for the ΓX direction, both the even- and odd-parity bands exhibit the off-Γ BICs below the critical radius, while the BICs are absent above the critical radius.     
 





By scanning all the directions around the Γ point, the real band structure 
is shown in Fig. <ref>. 



The anisotropy is evident in the band structure.  
The upper band is S-like, and  the lower band is P-like. 


Figure <ref> shows the Q-value maps of the two photonic bands originating from the E mode at the off-critical radii of Fig. <ref>. 



Now, the ring-like BICs found in the triangular-lattice system are absent. A similar pattern as in Fig.  <ref> (a) was observed in Ref. <cit.>. 
At r_s=0.41a, the total vorticity vanishes in the upper band. However, the lower band sustains the additional BICs on the ΓX axes. The net vorticity of the upper and lower bands is thus nonzero. 
All the BICs are discrete and move toward the Γ point as r_s to r_c. At r_s=r_c, the BICs are collapsed there. Across the critical coupling, multiple BICs are again generated now on the ΓM axes.  The total vorticity  is conserved across the critical coupling.  
  

As for the BIC2 of Fig. <ref>, the evolution of the multiple BICs across the critical radius is shown in Fig. <ref>. 



Since the two bands intersect in ΓM, the Q-value becomes singular there. At r_s=0.48a, the upper band exhibits the eight off-Γ BICs. 
They are elongated, and their polarizations are perpendicular and parallel to the ΓX and ΓM axes, respectively.  
The lower band holds four off-Γ BICs on the ΓM axes. 
All the BICs moves toward the Γ point as r_s→ r_c. 
After the collapse, the four off-Γ BICs are generated as in BIC1. 





 


 
 


§ *K·*P PERTURBATION THEORY
 

Let us consider the phenomena obtained in Secs. II and III in the *k·*p perturbation.  
The eigenmodes in the monolayer are determined by 

    ×(1/ε(*x)×*H(*x))=ω^2/c^2*H(*x).

Here, *H is the time-harmonic magnetic field with angular frequency ω, and ε(*x) is the dielectric function that is periodic in the in-plane coordinate *x_.
The inverse dielectric function  can be expanded  as 

    1/ε(*x)=∑_*g^*g·*x_η_*g(z).
 
Besides, the radiation field is also expanded by the plane waves via Bloch's theorem as 

    *H(*x)=∑_*g^(*k+*g)·*x_*h_*g(z).

Accordingly, the equation to be solved becomes  

    ∑_*g'((*k+*g)+ẑ∂_z) 
       5pt ×[  η_*g-*g'(z) ((*k+*g')+ẑ∂_z)
    		×*h_*g'(z)]=ω^2/c^2*h_*g(z),

which is symbolically expressed as 

    H|ψ⟩= E|ψ⟩,    E=ω^2/c^2.




Suppose we have a nearly accidental BIC of a doubly degenerate eigenmode at the Γ point.    
We perform the  *k·*p perturbation starting from these states as the zeroth order ones. 
The effective Hamiltonian to be diagonalized is given by <cit.>

    H_ab^eff=⟨ϕ_0a| H^(2) |ϕ_0b⟩
       20pt +∑_n 01/ϵ_0-ϵ_n⟨ϕ_0a| H^(1) |ϕ_n⟩⟨ϕ_n| H^(1) |ϕ_0b⟩, 
        H^(0) |ϕ_0a⟩=ϵ_0  |ϕ_0a⟩,    H^(0) |ϕ_n⟩=ϵ_n  |ϕ_n⟩,

where H^(i) is the operator in the left-hand side of Eq. (<ref>) 
expanded by Bloch momentum *k (i=0,1 and 2 imply *k=0, *k-linear, and *k-quadratic, respectively).  
The state |ϕ_0a⟩ (a=1,2) represents the doubly degenerate nearly accidental BIC mode at Γ. 
Here, we assume the degeneracy is not lifted in the first-order perturbation. That is, 

    ⟨ϕ_0a| H^(1) |ϕ_0b⟩=0.
 
This property is verified explicitly by the C_6v or C_4v symmetry. 
More generally, the matrix element ⟨ϕ_m| H^(1) |ϕ_n⟩ is nonzero only if 
(ϕ_m,ϕ_n) is attributed to (E_1,E_2), (E_1,A_1), (E_1,A_2), (E_2,B_1), (E_2,B_2), and vice versa for C_6v, and to (E,A_1), (E,A_2),
 (E,B_1), (E,B_2), and vice versa for C_4v (see Appendix A).  


The sum in the second term of Eq. (<ref>) generally includes the continuous radiation modes and discrete quasi-guided modes other than the mode concerned, as the intermediate states. The continuous radiation  modes give the decay rate to the original eigenmode via  Fermi's golden rule <cit.>.  However, if we start from the mode of the nearly accidental BIC, the decay rate is strongly suppressed. 


  
The continuum radiation modes belong to the same irreducible representation as the original doubly degenerate mode below the diffraction threshold.  Therefore, the matrix elements vanish between the nearly accidental BIC and radiation modes. That is why the decay rate is strongly suppressed. 


In the C_6v case, the effective hamiltonian is expressed as 

    H^eff = α*k^2 1 +β( (k_x^2-k_y^2)σ_3+2k_xk_yσ_1 ),
       α = α_E_1E_1^(2) + ∑_n∈ E_2|α_E_1E_2^(1)|^2 /ϵ_0-ϵ_n
       10pt +1/2∑_n∈ A_1|α_E_1A_1^(1)|^2/ϵ_0-ϵ_n
    +1/2∑_n∈ A_2|α_E_1A_2^(1)|^2/ϵ_0-ϵ_n,
       β = β_E_1E_1^(2)
       10pt +1/2∑_n∈ A_1|α_E_1A_1^(1)|^2/ϵ_0-ϵ_n
    -1/2∑_n∈ A_2|α_E_1A_2^(1)|^2/ϵ_0-ϵ_n

where σ_i(i=1,2,3) is the Pauli matrix, and parameters α and β are complex with [α]≥ 0. Their imaginary parts come solely from the one in the eigenfrequency of the E_1 mode, namely,  ϵ_0=(ω_0-γ_0)^2/c^2, provided that the intermediate modes are lossless (Im[ϵ_n]=0).
Therefore, if we start with the E_1 mode of the perfect accidental 
BIC (γ_0=0), [α]=[β]=0 so that
we still have a vanishing decay rate even at finite *k, within the second-order *k·*p perturbation.  This is the case of the extreme suppression of the decay rate around the Γ point of the accidental (next-to-super) BIC.   
Since the third-order terms vanish by the inversion symmetry, the imaginary frequencies of the two bands behave |*k|^4 as pointed out in Ref. <cit.>.





By diagonalizing the effective Hamiltonian, the eigenfrequencies of the two bands that stem from the E_1 mode at Γ are given by 

    ω^2/c^2=(ω_0-γ_0)^2/c^2 + (α±β)*k^2,  
       *H_*k^+(*x)∝ k_x *H_E_1^(01)(*x)+k_y*H_E_1^(02)(*x),
       *H_*k^+(*x)∝ -k_y *H_E_1^(01)(*x)+k_x*H_E_1^(02)(*x).

where (*H_E_1^(01),*H_E_1^(02)) forms the E_1 representation at Γ. 
If γ_0≪ω_0, the decay rate γ_k is evaluated as 

    γ_k≃γ_0-c^2/2ω_0[α±β]*k^2.

The radiation field of *H_*k^+(-) is P-like (S-like).   
The property of  [α]>0 promises a monotonic decreasing of the decay rate γ_k with increasing |*k| for at least one band, because either [α+β] or [α-β] is promised to be positive. The behavior of the other band depends on the relative magnitude of [α] and [β].  
Across Im[β]=0, the band with fast decreasing decay rate with |*k| is interchanged.  This condition is satisfied 
at the critical coupling. There, γ_0=0 and thus [β]=0.  
Moreover, since the dispersion relation is isotropic, we can have ring-like BICs at the minimum of γ_k. 
These results explain Fig. <ref> reasonably well. 

As for Figs. <ref> and <ref>, the spatial anisotropy is very strong, so that the *k·*p perturbation that predicts the isotropic dispersion is not so efficient. 








In the C_4v case, the effective Hamiltonian becomes 

    H^eff = α*k^2 1 +β_1  (k_x^2-k_y^2)σ_3+2β_2 k_xk_yσ_1,
       α = α_EE^(2)
       10pt 
    +1/2∑_n∈ A_1|α_EA_1^(1)|^2/ϵ_0-ϵ_n +1/2∑_n∈ B_1|α_EB_1^(1)|^2/ϵ_0-ϵ_n
       10pt 
    +1/2∑_n∈ A_2|α_EA_2^(1)|^2/ϵ_0-ϵ_n
    +1/2∑_n∈ B_2|α_EB_2^(1)|^2/ϵ_0-ϵ_n
    ,
       β_1 = β_1:EE^(2)
       10pt 
    +1/2∑_n∈ A_1|α_EA_1^(1)|^2/ϵ_0-ϵ_n +1/2∑_n∈ B_1|α_EB_1^(1)|^2/ϵ_0-ϵ_n
       10pt 
    -1/2∑_n∈ A_2|α_EA_2^(1)|^2/ϵ_0-ϵ_n
    -1/2∑_n∈ B_2|α_EB_2^(1)|^2/ϵ_0-ϵ_n
    ,
       β_2 = β_2:EE^(2)
       10pt 
    +1/2∑_n∈ A_1|α_EA_1^(1)|^2/ϵ_0-ϵ_n -1/2∑_n∈ B_1|α_EB_1^(1)|^2/ϵ_0-ϵ_n
       10pt 
    -1/2∑_n∈ A_2|α_EA_2^(1)|^2/ϵ_0-ϵ_n
    +1/2∑_n∈ B_2|α_EB_2^(1)|^2/ϵ_0-ϵ_n.

By diagonalizing the effective Hamiltonian, the eigenfrequency becomes 

    ω^2/c^2=(ω_0-γ_0)^2/c^2 + 
    	α*k^2 ±√(β_1^2(k_x^2-k_y^2)^2+4β_2^2 k_x^2k_y^2).

Again, [α]>0, promising a monotonic decreasing of the imaginary part in Eq. (<ref>) with increasing |*k| for at least one band branch.
Now, Eq. (<ref>) is not isotropic. Thus, the ring-like high Q channel is not formed. This trend is entirely consistent with the numerical results obtained in Sec. III. 

On the ΓX axis (k_y=0), the two eigenstates become

    ω^2/c^2=(ω_0-γ_0)^2/c^2 + (α±β_1)k_x^2,  
       *H_*k^+(*x)=*H_E^(01)(*x),   *H_*k^-(*x)=*H_E^(02)(*x),
  
where (*H_E^(01),*H_E^(02)) 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 eigenmodes is promised to have a decreasing decay rate with 
|k_x|, and the two eigenmodes are interchanged at [β_1]=0 
via the critical coupling condition [γ_0]=0. 

On the ΓM axis (k_x=k_y), the eigenstates become

    ω^2/c^2=(ω_0-γ_0)^2/c^2 + 2(α±β_2)k_x^2,  
       *H_*k^±(*x)∝*H_E^(01)(*x)±*H_E^(02)(*x).
  
The eigenmodes of superscript "+" is P-polarized, and "-" is  S-polarized. 
One of the two eigenmodes is promised to have a decreasing decay rate with 
|k_x|, and the two eigenmodes are interchanged at [β_2]=0 
via the critical coupling condition [γ_0]=0. 
Again, these features are consistent with Figs. <ref> and   <ref>. 




The symmetry-protected BICs are also available for C_3v point-group systems. There, the eigenmodes of the A_1 and A_2 representations at the Γ point are symmetry-protected, provided there are no open diffraction channels other than the specular one. In this case, it is possible to have accidental at-Γ BICs for the eigenmodes of the E representation by tuning system parameters. 

The effective Hamiltonian for the nearly accidental BIC mode is given by 

    H^eff = α*k^2 1 +β( (k_x^2-k_y^2)σ_3+2k_xk_yσ_1 ),
       α = α_EE^(2) + ∑_n∈ E '|α_EE^(1)|^2 /ϵ_0-ϵ_n
       10pt +1/2∑_n∈ A_1|α_EA_1^(1)|^2/ϵ_0-ϵ_n
     	+1/2∑_n∈ A_2|α_EA_2^(1)|^2/ϵ_0-ϵ_n,
       β = β_EE^(2)
       10pt +1/2∑_n∈ A_1|α_EA_1^(1)|^2/ϵ_0-ϵ_n
     	-1/2∑_n∈ A_2|α_EA_2^(1)|^2/ϵ_0-ϵ_n,
 
where the prime in the second term of Eq. (<ref>) means the unperturbed E mode is excluded in the sum.   
In this case, a significant contribution to the decay rate emerges from the continuous radiation modes of the E representation via Fermi's golden rule. Therefore, suppressing of the decay rate at nonzero |*k| near the accidental BIC (at Γ) does not occur. Consequently, the multiple off-Γ BICs are absent in the *k·*p perturbation theory.  

In the above arguments, we do not rely on the monolayer of spheres presented in Secs. II and III, but on the spatial symmetry of C_6v, C_4v, and C_3v. Accordingly, the multiple BIC generation with the present scenario is available in other photonic membranes with C_6v or C_4v.  We can show PhC slabs with the triangular or square lattice of circular air holes exhibit 
the multiple BIC generation. The ring-like high Q channel is also  formed in the C_6v case. 



§ SUMMARY
 
In summary, we have presented a detailed theoretical analysis of the multiple BIC generation in the monolayers of spheres with C_6v or C_4v point group. A tuning of system parameters results in an accidental BIC of doubly degenerate eigenmodes at the Γ point.   This BIC, obtained at a critical parameter, is the next-to-super BIC with extremely suppressed decay rates around the Γ point. Off-critical parameters yield multiple off-Γ BICs that are generated from the next-to-super BIC at  the Γ point.  

In the C_6v point-group system, a ring-like high Q channel pinned with multiple  off-Γ BICs on the mirror axes can be formed in the two bands originating from the doubly degenerate E_1 mode at the Γ point.   
The spatial anisotropy destroys the ring in the C_4v point-group system, but multiple BICs are certainly formed. 
Across the critical parameters, these BICs move from the upper band to the lower band or vice versa, conserving the total vorticity of the BICs. 

We also show that these phenomena are owing to the spatial symmetry of C_6v or C_4v and the *k·*p perturbation explains the phenomena reasonably well. In addition, the C_3v system does not support the multiple BIC generation through a similar design. 

 
 
 
   

This work was supported by JSPS KAKENHI Grant No.  22K03488. 






§ MATRIX ELEMENTS IN THE *K·*P PERTURBATION

We summarize various matrix elements relevant to the effective Hamiltonian. The symmetry relation under point-group operations is crucial. 
It is given by 

    H_R_1R_2^(i)(*k)=D_R_1^†(A) H_R_1R_2^(i)(A*k)D_R_2(A),

where the matrix element is the one between the modes of irreducible representations  R_1 and R_2, and D_R(A) is the representation matrix of group element A in the irreducible representation R.  
The equation (<ref>) constrains the possible form of the matrix elements. 


In the C_6v case, 
The nonzero matrix elements in the first-order *k·*p perturbation are given by 

    H_E_1E_2^(1)(*k)=α_E_1E_2^(1)(k_xσ_1+k_yσ_3),
        H_E_1A_1^(1)(*k)=α_E_1A_1^(1)( k_x
     k_y),
        H_E_1A_2^(1)(*k)=α_E_1A_2^(1)( -k_y
     k_x),
        H_E_2B_1^(1)(*k)=α_E_2B_1^(1)( -k_y
     k_x),
        H_E_2B_2^(1)(*k)=α_E_2B_2^(1)(  k_x
     k_y).

In the second order, the relevant matrix element becomes 

    H_E_1E_1^(2)(*k)=α_E_1E_1^(2)*k^2 +β_E_1E_1^(2)((k_x^2-k_y^2)σ_3+2k_xk_yσ_1),

where α_E_1E_1^(2) and β_E_1E_1^(2) are real. 

In the C_4v case, 
the nonzero matrix elements in the first-order *k·*p perturbation  are given by 

    H_EA_1^(1)(*k)=α_EA_1^(1)( k_x
     k_y),
        H_EA_2^(1)(*k)=α_EA_2^(1)( -k_y
     k_x),
        H_EB_1^(1)(*k)=α_EB_1^(1)( k_x
     -k_y),
        H_EB_2^(1)(*k)=α_EB_2^(1)(  k_y
     k_x).

In the second order, the relevant matrix element becomes 

    H_EE^(2)(*k)=α_EE^(2)*k^2 1  +β_1:EE^(2)(k_x^2-k_y^2)σ_3+ β_2:EE^(2)2k_xk_yσ_1,

where α_E_1E_1^(2), β_1:E_1E_1^(2), β_2:E_1E_1^(2) are real. 



In the C_3v case, the nonzero matrix elements in the first-order *k·*p perturbation  are given by 

    H_EE^(1)(*k)=α_EE^(1)( k_xσ_1+k_yσ_3),
        H_EA_1^(1)(*k)=α_EA_1^(1)( k_x
     k_y),
        H_EA_2^(1)(*k)=α_EA_2^(1)( -k_y
     k_x).

The matrix element of Eq. (<ref>) vanishes by the time-reversal symmetry if the two E modes that form the matrix element are identical.   
In the second order, the relevant matrix element becomes 

    H_EE^(2)(*k)=α_EE^(2)*k^2 +β_EE^(2)((k_x^2-k_y^2)σ_3+ 2k_xk_yσ_1).





§ ELECTRIC FIELDS OF QUASI-GUIDED MODES
 
We employ the following method for evaluating the polarization ellipse map.  At a true BIC point, the kernel of a secular matrix gives the eigenmode, and its determinant vanishes on the real frequency axis.  




In the monolayer of spheres, the secular matrix S is given by <cit.>

    S_(Lβ)(L'β')=
    δ_LL'δ_ββ'-t_l^β G_LL'^ββ', 
       G_LL'^ββ'=1/l(l+1)∑_L_1L_2[(*P^β)^†]_LL_1G_L_1L_2·[*P^β']_L_2L',
       G_LL'=4π∑_L”^l-l'-l”⟨ L|L”|L'⟩ S_L”,
       S_L=∑_*X 0h_l^(1)(q_0|*X|)Y_L^*(*̂X̂)^*k·*X.

where L=(l,m) (|m|≤ l) is the angular momentum index, β(=M,N) is the index to classify two transverse vector spherical waves, t_l^β is  the so-called t-matrix of the sphere, *P^β is the transformation matrix from scalar spherical waves to vector spherical waves, ⟨ L|L”|L'⟩ is the Clebsh-Gordan coefficient, h_l^(1) is the spherical hankel function of the first kind, Y_L is the spherical harmonics, and *X is the 2D real-lattice vector.   
The structure constant S_L is efficiently calculated via the Ewald method <cit.>.  
Through ψ_L^β = Ker[ S], the electric field outside the monolayer is expressed as 

    *E(*x)=∑_L*Xh_l^(1)(q_0|*x-*X|)Y_L(*x-*X)*V_L ^*k·*X,
       *V_L=∑_L'β*P_LL'^βψ_L'^β.
	
Its plane-wave-expansion coefficient t_σ*g^± (σ=P,S) becomes 

    t_σ*g^±=	2π/q_0Γ_*gA_UC∑_L (-)^l Y_L(*̂K̂_*g^±)*σ_*g^±·*V_L,

where A_UC is the area of the unit cell. 


     
  
Off the BIC point, the determinant is nonzero on the real axis but can be zero in the complex frequency plane. This zero should be 
the resonance pole of the S-matrix, namely, Ω_k=ω_k-γ_k. 
We assume that  the right eigenstate of S with the least absolute eigenvalue can be approximated as the kernel on the real frequency axis near the zero. Then, we can obtain the plane-wave  coefficients t_σ*g^±(ω) through Eq.  (<ref>). 
We further assume that it can be Taylor expanded around ω=Ω_k, namely, 

    t_σ*g^±(ω)≃ t_σ*g^±(Ω_k)+(ω-Ω_k)C_σ*g^±   (σ=P,S).
 
Then, from the two points ω_a and ω_b on the real frequency axis, we can estimate t_σ*g^±(Ω_k) as 

    t_σ*g^±(Ω_k)=(ω_a-Ω_k)t_σ*g^±(ω_b)-(ω_b-Ω_k)t_σ*g^±(ω_a)/ω_a-ω_b.
 
We confirm that the resulting t_σ*0^±(Ω_k) vanishes at the BIC points. 
 
 










31
	fxundefined [1]
		ifx#1
	
	fnum [1]
		#1firstoftwo
		secondoftwo
		
	fx [1]
		#1firstoftwo
		secondoftwo
		
	
	
	
	
	
	noop [0]secondoftwo
	
	ref[1]@startlink#1@href
	href[1]#1@endlink
	anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
	startlink[1]
	endlink[0]
	
	rl [1]href #1
	
	
	
	
	
	
	
	
	
	
	
	
	@bib@innerbibempty
	
	[Hsu et al.(2016)Hsu,
		Zhen, Stone, Joannopoulos,
		Soljacic, and Soljačić]hsu2016bound
	author author C. W. Hsu, author B. Zhen, author A. D. Stone, author
		J. D. Joannopoulos, author
		M. Soljacic,  and author
		M. Soljačić, 10.1038/natrevmats.2016.48 journal journal Nat.
				Rev. Mater. volume 1, pages 16048
		(year 2016)NoStop
	[Kodigala et al.(2017)Kodigala, Lepetit, Gu, Bahari, Fainman, and Kanté]Kodigala2017
	author author A. Kodigala, author T. Lepetit,
		author Q. Gu, author
		B. Bahari, author Y. Fainman,  and author B. Kanté, 10.1038/nature20799
	journal journal Nature volume 541, pages 196 (year
		2017)NoStop
	[Romano et al.(2018)Romano,
		Zito, Torino, Calafiore,
		Penzo, Coppola, Cabrini,
		Rendina, and Mocella]Romano2018
	author author S. Romano, author G. Zito,
		author S. Torino, author G. Calafiore, author
		E. Penzo, author G. Coppola, author S. Cabrini, author I. Rendina,  and author V. Mocella, 10.1364/prj.6.000726
	journal journal Photonics Res. volume 6, pages 726 (year
		2018)NoStop
	[Koshelev et al.(2020)Koshelev, Kruk, Melik-Gaykazyan,
		Choi, Bogdanov, Park, and		Kivshar]Koshelev2020
	author author K. Koshelev, author S. Kruk,
		author E. Melik-Gaykazyan,
		author J. H. Choi, author A. Bogdanov, author
		H. G. Park,  and author
		Y. Kivshar, 10.1126/science.aaz3985 journal journal Science
				 volume 367, pages 288
		(year 2020)NoStop
	[Miyazaki and Ohtaka(1998)]Miyazaki1998a
	author author H. Miyazaki and author K. Ohtaka, @noop  journal journal Phys.
				Rev. B volume 58, pages 6920
		(year 1998)NoStop
	[Paddon and Young(2000)]Paddon:Y::61:p2090-2101:2000
	author author P. Paddon and author J. F. Young, @noop  journal journal Phys.
				Rev. B volume 61, pages 2090
		(year 2000)NoStop
	[Ochiai and Sakoda(2001a)]Ochiai:S::63:p125107:2001
	author author T. Ochiai and author K. Sakoda, 10.1103/PhysRevB.63.125107 journal journal Phys. Rev. B volume
			63, pages 125107 (year
		2001a)NoStop
	[Fan and Joannopoulos(2002)]Fan:J::65:p235112:2002
	author author S. H. Fan and author J. D. Joannopoulos, @noop  journal journal
			Phys. Rev. B volume 65, pages 235112
		(year 2002)NoStop
	[Hsu et al.(2013)Hsu,
		Zhen, Lee, Chua,
		Johnson, Joannopoulos, and Soljačić]Hsu2013b
	author author C. W. Hsu, author B. Zhen, author J. Lee, author
		S. L. Chua, author
		S. G. Johnson, author
		J. D. Joannopoulos,  and author M. Soljačić, 10.1038/nature12289 journal journal
			Nature volume 499, pages 188
		(year 2013)NoStop
	[Jiang et al.(2023)Jiang,
		Hu, Wang, Han, and		Zi]Jiang2023b
	author author Q. Jiang, author P. Hu, author J. Wang, author
		D. Han,  and author
		J. Zi, 10.1103/PhysRevLett.131.013801 journal journal
			Phys. Rev. Lett. volume 131, pages
		013801 (year 2023)NoStop
	[Zhen et al.(2014)Zhen,
		Hsu, Lu, Stone, and		Soljačić]Zhen2014
	author author B. Zhen, author C. W. Hsu,
		author L. Lu, author
		A. D. Stone,  and author
		M. Soljačić, 10.1103/PhysRevLett.113.257401 journal journal
			Phys. Rev. Lett. volume 113, pages
		257401 (year 2014)NoStop
	[Jin et al.(2019)Jin,
		Yin, Ni, Soljačić, Zhen, and Peng]Jin2019
	author author J. Jin, author X. Yin, author L. Ni, author
		M. Soljačić, author
		B. Zhen,  and author
		C. Peng, 10.1038/s41586-019-1664-7 journal journal
			Nature volume 574, pages 501
		(year 2019)NoStop
	[Hwang et al.(2021)Hwang,
		Lee, Kim, Jeong,
		Kwon, Koshelev, Kivshar,		and Park]Hwang2021a
	author author M.-S. S.			Hwang, author H.-C. C.			Lee, author K.-H. H.			Kim, author K.-Y. Y.			Jeong, author S.-H. H.			Kwon, author K. Koshelev, author Y. Kivshar, 		and author H.-G. G. Park,	10.1038/s41467-021-24502-0 journal journal Nat. Commun. volume 12, pages 4135 (year 2021)NoStop
	[Kostyukov et al.(2022)Kostyukov, Gerasimov, Ershov, and		Bulgakov]Kostyukov2022
	author author A. S. Kostyukov, author V. S. Gerasimov, author A. E. Ershov,  and author E. N. Bulgakov, 10.1103/PhysRevB.105.075404 journal journal Phys. Rev. B volume
			105, pages 075404 (year 2022)NoStop
	[Sadrieva et al.(2019)Sadrieva, Frizyuk, Petrov, Kivshar, and Bogdanov]Sadrieva2019
	author author Z. Sadrieva, author K. Frizyuk,
		author M. Petrov, author Y. Kivshar,  and author A. Bogdanov, 10.1103/PhysRevB.100.115303 journal journal
			Phys. Rev. B volume 100, pages
		115303 (year 2019) NoStop
	[Gladyshev et al.(2022)Gladyshev, Shalev, Frizyuk, Ladutenko, and Bogdanov]Gladyshev2022a
	author author S. Gladyshev, author A. Shalev,
		author K. Frizyuk, author K. Ladutenko,  and author A. Bogdanov, 10.1103/PhysRevB.105.L241301 journal journal
			Phys. Rev. B volume 105, pages
		L241301 (year 2022) NoStop
	[Yoda and Notomi(2020)]Yoda2020
	author author T. Yoda and author M. Notomi,	10.1103/PhysRevLett.125.053902 journal
		journal Phys. Rev. Lett. volume 125,		pages 053902 (year 2020) NoStop
	[Doiron et al.(2022)Doiron,
		Brener, and Cerjan]Doiron2022
	author author C. F. Doiron, author I. Brener, 		and author A. Cerjan, 10.1038/s41467-022-35246-w journal journal Nat. Commun. volume 13, pages 7534 (year 2022)NoStop
	[Wang et al.(2023)Wang,
		Wang, Zhao, Shi, and		Zi]Wang2023a
	author author X. Wang, author J. Wang, author X. Zhao, author
		L. Shi,  and author
		J. Zi, 10.1021/acsphotonics.2c01522 journal journal ACS
				Photonics volume 10, pages 2316
		(year 2023)NoStop
	[Inui et al.(1996)Inui,
		Tanabe, and Onodera]inui1996group
	author author T. Inui, author Y. Tanabe,  and		author Y. Onodera, @noop
	 title Group Theory and its Applications in Physics	(publisher Springer, year 1996)NoStop
	[Ohtaka(1980)]Ohtaka1980
	author author K. Ohtaka, 10.1088/0022-3719/13/4/022 journal journal J. Phys. C Solid State Phys. volume 13, pages 667 (year
		1980)NoStop
	[Modinos(1987)]MODINOS::141:p575-588:1987
	author author A. Modinos, @noop  journal journal
			Physica A volume 141, pages 575
		(year 1987)NoStop
	[Stefanou et al.(1992)Stefanou, Karathanos, and Modinos]Kristensson2016
	author author N. Stefanou, author V. Karathanos,  and author A. Modinos, 10.1049/SBEW524E journal
		journal J. Phys. Cond. Matt. volume
			4, pages 7389 (year 1992)NoStop
	[Landau and Lifshitz(1977)]landau2013quantum
	author author L. D. Landau and author E. M. Lifshitz, @noop  title Quantum Mechanics:
				Non-relativistic Theory (publisher Pergamon, address Oxford, year 1977)NoStop
	[Ohtaka et al.(2004)Ohtaka,
		Inoue, and Yamaguti]Ohtaka:I:Y::70:p035109:2004
	author author K. Ohtaka, author J. Inoue, 		and author S. Yamaguti,	@noop  journal journal Phys. Rev. B		volume 70, pages 35109 (year
		2004)NoStop
	[Fano(1961)]Fano1961
	author author U. Fano, 10.1103/PhysRev.124.1866 journal
		journal Phys. Rev. volume 124,		pages 1866 (year 1961)NoStop
	[Ochiai(2012)]PhysRevB.86.075152
	author author T. Ochiai, 10.1103/PhysRevB.86.075152 journal journal Phys. Rev. B volume
			86, pages 075152 (year 2012)NoStop
	[Ochiai and Sakoda(2001b)]Ochiai2001d
	author author T. Ochiai and author K. Sakoda, 10.1103/PhysRevB.64.045108 journal journal Phys. Rev. B volume
			64, pages 045108 (year
		2001b)NoStop
	[Andreani and Gerace(2006)]Andreani2006
	author author L. C. Andreani and author D. Gerace, 10.1103/PhysRevB.73.235114 journal journal Phys. Rev. B volume
			73, pages 235114 (year 2006) NoStop
	[Ohtaka(1979)]Ohtaka1979
	author author K. Ohtaka, 10.1103/PhysRevB.19.5057 journal
		journal Phys. Rev. B volume 19,		pages 5057 (year 1979)NoStop
	[Kambe(1967)]Kambe1967
	author author K. Kambe, 10.1515/zna-1967-0305 journal
		journal Zeitschrift fur Naturforsch. - Sect. A J. Phys. Sci.		volume 22, pages 322 (year
		1967)NoStop


