Near-Field Excitation of Bound States in the Continuum in All-Dielectric Metasurfaces through a Coupled Electric/Magnetic Dipole Model

Resonant optical modes arising in all-dielectric metasurfaces have attracted much attention in recent years, especially when so-called bound states in the continuum (BICs) with diverging lifetimes are supported. With the aim of studying theoretically the emergence of BICs, we extend a coupled electric and magnetic dipole analytical formulation to deal with the proper metasurface Green function for the infinite lattice. Thereby, we show how to excite metasurface BICs, being able to address their near-field pattern through point-source excitation and their local density of states. We apply this formulation to fully characterize symmetry-protected BICs arising in all-dielectric metasurfaces made of Si nanospheres, revealing their near-field pattern and local density of states, and, thus, the mechanisms precluding their radiation into the continuum. This formulation provides, in turn, an insightful and fast tool to characterize BICs (and any other leaky/guided mode) near fields in all-dielectric (and also plasmonic) metasurfaces, which might be especially useful for the design of planar nanophotonic devices based on such resonant modes.

Interestingly, BICs arising in metasurfaces are typically explored through far-field reflection/transmission; since they are not accessible under (propagating) plane wave illumination by definition, BICs are observed as the vanishing of the quasi-BIC in certain parameter space when the BIC condition is approached (quasi-BICs). Thus, it would be desirable to have means to directly address BICs through resonant mode characterization and near-field excitation, including local density of states (LDOS), rigorously dealing with complex interactions between particle and lattice resonances in planar infinite arrays. This would also shed light onto the various mechanisms precluding BIC radiation into the far field, symmetry protection and accidental degeneracy being particularly relevant to metasurfaces; see, e.g., References [42,45]. Full wave numerical calculations very often lack accuracy and speed when dealing with modes with diverging Q-factor, not to mention the fact that they cannot shed much light onto the underlying physics.
In this work, we extend a previously developed coupled electric and magnetic dipole (CEMD) formulation for plane wave excitation [19], to deal with the proper Green function of the scattering problems. This allows in practice to consider point-dipole near-field excitation, mode dispersion relation, and local density of electromagnetic states. Thereby, we analyze the emergence of BICs in a planar array of silicon nanospheres, characterizing their vectorial LDOS at a plane close to the unit cell. Such LDOS information is used to explore the proper conditions (position and polarization) for point dipole excitation of BICs, in turn revealed through the near-field pattern calculated across the metasurface plane, with symmetry protection being evidenced. Finally, the conclusions extracted from such rich near-field BIC properties are drawn.

Methods: Coupled Electric/Magnetic Dipole Model for Point Dipole Excitation
Let us first consider a monochromatic point dipole source located at at an arbitrary position r d = (x d , y d , z d ) with respect to a planar array of scatterers located in the plane z = 0 (see Figure 1a). Using the Weyl expansion, the field emitted by the dipole at the position r = (x, y, z), ψ d (ω, r − r d ) can be expressed as: with k = ω/c, where c and ω are the speed and frequency of light, P characterizes the polarization of the dipole (both electric and magnetic), and ← → M ω, Q x , Q y is the matrix Green function of a point (electric and magnetic) dipole.
Thus Equation (1) can be formally written as a sum of (integral over) plane waves: where the intensity and the vector polarization distribution, I d and v d , respectively, are The sign (±) of the exponentials that depend on z and z d is chosen depending on the sign of (z − z d ).
Within our CEMD model [19], for monochromatic plane wave illumination with wavevector k 0 = k xx + k yŷ + k zẑ and linear polarization v, the field scattered by the metasurface, ψ sca PW (ω, r, k 0 , a 1 , a 2 ), can be written as: where a 1 and a 2 are the primitive vectors that characterize the metasurface. In addition, ← → G ω, r, k x , k y , a 1 , a 2 is the lattice Green function, and ← → α ω, k x , k y , a 1 , a 2 is the renor-malized polarizability of the metasurface defined by the polarizability of the particles, ← → α (ω), and the lattice depolarization Green function, ← → G ω, k x , k y , a 1 , a 2 , as follows: Therefore, the field scattered by the metasurface when illuminated by a point source, ψ sca d (ω, r, r d ), can be obtained by replacing the Weyl expression of the incident dipole source, Equation (3), into the corresponding scattered field given by Equation (6), leading to: where the dependencies on primitive vectors, a 1 and a 2 , have been omitted. Recall that such electromagnetic field ψ sca d (ω, r, r d ) is formally equivalent to the Green function of the scattering geometry with argument r − r d within the dipole approximation. Incidentally, we do not refer to it as ← → G in order to avoid confusion with lattice Green functions within our CEMD model. Interestingly, note also that the LDOS of the metasurface (within the dipolar approximation) can be extracted from ψ sca d (ω, r, r d ) evaluated at the position of the dipole r = r d , as follows: The Green function given by Equation (8), including the LDOS expression (9) are our main analytical results herein, allowing us, in turn, to determine the electromagnetic field close to the metasurface (near field and LDOS) upon point source illumination, indeed crucial to unveil the excitation in the near field of resonant modes and BICs.
Since the Weyl expansion allows us to differentiate between propagating and evanescent components of the dipolar source, it is convenient to split the integral in Q-space over Q x , Q y in the form: where the subscripts pro and eva stand for propagating and evanescent waves, respectively. Therefore, it is possible to separately study the contributions of propagating and evanescent waves to the near field, and by extension, to the LDOS. This will be indeed crucial to unveil the evanescent character of BICs. As a final remark, note that the dipole intensity diverges at Q 2 x + Q 2 y = k 2 , so the integral in the Weyl expansion is interpreted as an improper integral. Although well defined, the numeric integration can present problems. To avoid them, a change of variable from Cartesian coordinates (Q x and Q y ) to pseudo-polar coordinates ( k 2 ± Q 2 and θ) is used. For the propagating components, Q 2 x + Q 2 y < k 2 : with 0 ≤ Q < k and −π ≤ θ < π. On the other hand, for the evanescent components, with 0 < Q < ∞ and −π ≤ θ < π. Since the Jacobian of the transformation is proportional to Q, the divergence disappears.

Si Nanosphere Metasurface
Let us apply our formulation to explore the LDOS and near field patterns of canonical BICs arising in a metasurface consisting of a square array of Si nanospheres of radius R = 100 nm with lattice constant a = b = 4R = 400 nm. A dispersionless dielectric permittivity of = 12.25 is used to simplify the theoretical analysis, bearing in mind both that it is actually similar to that of Si and other semiconductors (n = 3.5), and that losses might be very small for wavelengths above λ ∼ 500 nm (as considered here); see, e.g., Ref. [48] for polycrystalline Si, GaP, etc. The reflectance spectra for varying angles of incidence, calculated as in Ref. [19], are shown in Figure 1b-e. We focus on two symmetry-protected BICs arising at the Γ point from vertical (out-of-plane) dipole resonances (magnetic, MD, for TE polarization at λ = 709 nm in Figure 1b, and electric, ED, for TM polarization at λ = 552 nm in Figure 1e), which exhibit the expected quasi-BIC bands at off-normal angles with diverging Q-factors that become inaccessible from the far-field exactly at the BIC condition at θ = 0 • .

Vertical Electric-Dipole BIC
We first explore in Figure 2 the LDOS for the symmetry-protected ED-BIC shown in Figure 1e at λ = 552 nm, on a plane parallel to the metasurface in a region equivalent to one unit cell located at a subwavelength distance z d = a/3 ≈ 133 nm (z d ∼ λ/4). All the different electromagnetic contributions are displayed of either electric or magnetic dipolar character along the three Cartesian components (p x , p y , p z , m x , m y , m z ). Indeed, even though our formulation only retains the dipolar response of the nanospheres, such contributions correlate quite reasonably with the ED-resonant character of this BIC (basically dipolar indeed): the larger contribution to the LDOS by far is the vertical electric one, p z , in turn concentrated at the sphere center, since the ED resonance actually emerges from the electric field pointing perpendicular to the metasurface plane. The other electric in-plane electric contributions, though clearly observable, are substantially weaker, and so are all the magnetic terms. Interestingly, the magnetic vertical dipole m z exhibits a suppression of the LDOS in the entire unit cell. We next show in Figure 3a the dependence of the LDOS on the distance to the metasurface z d , with (x d , y d ) = (0, 0) where the LDOS is maximum (see Figure 2c), including explicitly the evanescent contribution; the limit z → 0 is avoided for the sake of validity of our dipolar approach. As the distance becomes smaller than z d /a ∼ 0.6 (z d ∼ λ ∼ 240 nm), the contribution to p z from the evanescent components grows exponentially, indicating that a strongly confined, evanescent mode (ED-BIC) is responsible for this enhanced LDOS. It is also evident in Figure 3a that if the p z dipole is displaced from the unit cell center to a position x d = 0.2a = 80 nm where the LDOS decreases, the enhanced evanescent contribution to the LDOS diminishes, in turn becoming nearly negligible at a position x d = 0.5a = 200 nm in the border of the unit cell where the p z -LDOS is minimum in Figure 2c. All this evidence supports the fact that such a strongly confined, evanescent mode corresponds to the expected ED-BIC. In addition, we show in Figure 3b the LDOS associated to two in-plane terms, p x and m y , at a position where their LDOS is close to the corresponding maxima (see Figure 2a,e). By contrast, these contributions show a small increase of the evanescent components upon approaching the metasurface plane, much weaker than that shown above for the ED-BIC, very likely stemming from the broad in-plane resonances overlapping with the ED-BIC as observed in Figure 1e. On the basis of the information provided by the LDOS, we now explore the near-field (NF) excitation of the BIC. According to Figures 2c and 3a, we place a vertical ED source at r d = (0, 0, a/3) (thus, z d ≈ 133 nm), where the p z -LDOS is very large. The resulting near-field map on a parallel plane at a fixed distance z = z d ≈ 133 nm from the metasurface is shown in Figure 4, separated into components. The total electric field amplitude in Figure 4a reveals a slightly higher concentration close to the projection of the ED source (center of the Si nanosphere metasurface), but, most importantly, the electric NF spreads over the metasurface plane without apparent decay. Indeed, if we analyze separately the contributions from the evanescent and propagating components (see Figure 4b,c, respectively), we observe that the propagating components radiate away from the metasurface center, wherein the dipole source is located, emerging as outgoing spherical waves with decaying amplitude. Conversely, we clearly identify the evanescent components in Figure 4b as responsible for the non-leaky mode across the metasurface plane expected for a BIC. Moreover, if we focus on the E z component shown in Figure 4d,e, we confirm the expected NF pattern of the corresponding vertical-ED BIC, spreading with identical phase in every unit cell, thus precluding radiation by symmetry protection. This, in turn, supports the correlation between this planar constant NF pattern in Figure 4a,b and the enhanced (evanescent) p z -LDOS in Figure 2c with the symmetry-protected ED-BIC.  To further confirm such correlation, we next show in Figure 5 the NF maps arising with other point-dipole sources located also at z d = z = a/3 ≈ 133 nm: in-plane electric (p x ) and magnetic (m x ), and vertical magnetic (m z ). Only in the unit cells corresponding to the location of the dipole source and nearest neighbors the NF pattern is significant, indicating that the electromagnetic field is basically scattered away from the metasurface, thus showing negligible excitation of the ED-BIC. In addition, it can be seen that the fields present much lower intensity than those arising when excited by the vertical electric dipole (p z , Figure 2). It is evident that none of them are able to efficiently couple to the ED-BIC, as predicted by the low LDOS stemming from such components; see Figure 2a,d,f. Finally, we also consider in Figure 6 the LDOS at the wavelength of the ED-BIC (λ = 552 nm) as a function of the nanosphere radius for all vectorial components. Clearly, the evanescent contribution to the p z term in Figure 6a shows a large LDOS maximum at the ED-BIC condition R/a = 1/4, which abruptly decreases with a slight mismatch in nanosphere size. Recall that this variation is nearly equivalent to a frequency mismatch for fixed nanosphere radius, so that this peak qualitatively correlates with the diverging Q-factor expected upon approaching the BIC, only limited here by the fact that the LDOS at a given distance from the metasurface remains finite. As expected, the other evanescent contribution to all other LDOS terms in Figure 6a remain structureless (and nearly negligible, except for a small contribution from the in-plane magnetic terms), and so do all the propagating contributions in Figure 6b.

Vertical Magnetic-Dipole BIC
For the sake of completeness, let us analyze the LDOS of the other symmetry-protected BIC supported by the Si nanosphere metasurface as shown in Figure 1b. The resulting LDOS maps are presented in Figure 7. First, the MD-BIC is neatly evident in Figure 7f through a large contribution from the vertical MD, m z , concentrated as expected at the vertical of the sphere center. Other moderately strong in-plane electric contributions are obtained (see Figure 7a,b), since the MD resonance actually emerges from the electric field circulating inside the nanospheres (recall that the Si nanospheres are non-magnetic). Interestingly, note also that the in-plane m x , m y contributions (and also as a significant p z term) yield relatively strong LDOS (see Figure 7d,e,c). This stems from the fact that the MD resonance is degenerate for a sphere due to symmetry and the proper MD-BIC coexists with a relatively strong in-plane magnetic contribution revealed as a broad background in Figure 1b. The reason that this degeneracy-induced coexistence is observed for the MDrelated BIC, but not for the ED-BIC shown above, most likely stems from the impact of the lattice: weaker (respectively, stronger) in-plane coupling of the vertical MD (respectively, ED) imposes a smaller (respectively, larger) frequency shift of the MD-BIC (respectively, ED) with respect to the bare MD (respectively, ED) nanosphere resonance. Let us next explore how the LDOS relates to the NF excitation shown in Figure 8. As expected, it is evident in Figure 8a that the vertical MD source (m z ) excites the MD-BIC spreading without losses across the metasurface plane. Indeed, this MD-BIC behavior stems from the evanescent contribution (see Figure 8b), with the propagating contribution revealing radiation away from the unit cell close to the MD source (see Figure 8c). Furthermore, the vertical magnetic field component is shown to be the predominant one, and its phase-matched condition across unit cells corroborates the symmetry-protection mechanism that precludes leakage of this MD-BIC (see Figure 8d,e).

Conclusions
In summary, we have developed a coupled electric and magnetic dipole analytical formulation to determine the Green function for an infinite planar array of electric/magnetic dipoles. This allows us to obtain any physical magnitude connected to the excitation by a point dipole source: particularly, we have focused on two relevant magnitudes: LDOS and near fields. Those magnitudes are of the utmost importance when dealing with so-called BICs, since these modes cannot be probed by definition from the far field (except indirectly). Thereby, we have investigated the LDOS associated to symmetry-protected BICs emerging in a metasurface consisting of an square array of Si nanospheres. The LDOS at a subwavelength distance from the metasurface clearly describe the nature of such BICs (vertical ED or MD, both with predominant evanescent contributions), in turn revealing the most appropriate location over the unit cell to effectively couple to the corresponding BIC mode. Such information is used to calculate the near field pattern when the metasurface is excited by a suitable point dipole source, leading to BIC excitation across the metasurface without leakage; the resulting pattern indeed supports their symmetry-protection mechanisms. Our formulation can be exploited to investigate BICs and any other guided/leaky modes supported by metasurfaces and metagratings of interest throughout the electromagnetic spectrum, as long as the predominant optical response of the meta-atoms is of electric and/or magnetic dipolar nature, which is the case of many plasmonic, all-dielectric, and/or hybrid sub-wavelength scatterers. In this regard, recall that our coupled electric/magnetic dipole model for far-field illumination has proven very accurate, not only qualitatively but even quantitatively, in describing the properties of a variety of metasurfaces in the microwave, THz, and visible domains [46][47][48]59,60]. Moreover, cutting-edge experiments for near-field excitation based on the model presented here are indeed underway.