Refractive index sensing with optical bound states in the continuum

We consider refractive index sensing with optical bounds states in the continuum (BICs) in dielectric gratings. Applying a perturbative approach we derived the differential sensitivity and the figure of merit of a sensor operating in the spectral vicinity of a BIC. Optimisation design approach for engineering an effective sensor is proposed. An analytic formula for the maximal sensitivity with an optical BIC is derived. The results are supplied with straightforward numerical simulations.


FIG. 1: (a)
The dielectric grating assembled of rectangular bars on glass substrate. The plane of incidence y0z is shaded grey. The magenta arrow shows the electric vector of the incident wave. The geometric parameters are w = 0.5h, b = 0.5h. (b, c) BIC eigenmode profiles for Si bars as E y in the y0z-plane with (n (0) c = 1.333) and without (n (0) c = 1) cladding, correspondingly.

I. INTRODUCTION
Bound states in the continuum (BICs) have revolutionized nanophotonics offering the opportunity to realize a new class of high throughput sensing devices and improving the control of the interaction between light and matter at the nanometer scale [1][2][3][4][5][6][7][8][9][10] . Strictly speaking, a BIC can be considered as a resonant mode with an infinite quality factor (Q-factor) in an open-cavity. This special mode has a frequency in the radiation continuum but does not lose energy because of symmetry mismatch with radiative waves. The BICs are hosted by a leaky band of high-quality resonances with Q-factor diverging in the Γ-point 11 . The leaky band with diverging Q-factor, in turn, induces a collapsing Fano feature in the transmittance spectrum [12][13][14][15] . Experimentally, this results in sharp spectral features in the transmission spectrum. Generally, the position of these extremely narrow Fano resonances is sensitive to the refractive index of the surrounding medium allowing to engineer optical sensors with a good sensitivity, S and an excellent figure of merit (FOM ). Among the different geometries of structures supporting BICs, all-dielectric resonant planar structures, including photonic crystal slabs, periodic arrays and metasurfaces, have been taken as an alternative to conventional plasmonic sensors due to the possibility to realize high-performance sensing systems in loss-free media with a remarkable sensitivity to a small volume of the analyte [16][17][18][19][20] . While S is affected by the spatial overlap between the nonradiating evanescent field and the surrounding cladding, FOM is proportional to the Q-factor and ultimately represents the sensor capability to follow tiny changes in the environment refractive index 21,22 . In this paper we propose the derivation of the analytical expressions for S and FOM for sensors based on dielectric grating (DGs), one of the major set-ups for studying optical BICs both theoretically [23][24][25][26][27][28][29] and experimentally 30,31 . By using a perturbative approach we find the differential sensitivity and the figure of merit of a DG-sensor operating in the spectral vicinity of a BIC. The reported analytic results could pave a way for an optimisation design approach for engineering an effective sensor based on BICs and provide a tool for a better understanding of the physics underlying this sensing mechanism.

II. BOUND STATES IN THE CONTINUUM
Here we consider the system shown in Fig. 1 (a). It is a dielectric grating (DG) assembled of rectangular bars. The bars are periodically placed on the glass substrate with period h. We assume that the incident field is TE polarized, i.e. the electric vector is aligned with the bars as shown in Fig. 1 (a). The propagation of the TE-modes is controlled by the Helmholtz equation for the x component of the electric field: where k is the vacuum wavenumber, and ǫ is the dielectric function. In what follows we assume that the bars can be made of different dielectric materials such as polycrystalline silicon, Si; and titanium dioxide, Ti0 2 . The refractive index (RI) of the substrate n 0 = 1.5. The RI sensor operating principle is detection of the shift in the wavelength of an optical resonance, λ res in response to the change of the RI of the cladding fluid, n c on top of the DG. As the reference value of the cladding RI we take that of water n (0) c = 1.333. Our finite-difference time-domain (FDTD) simulations demonstrate that for the set of geometric parameters specified in the caption to Fig. 1 a Si DG supports a BIC with kh = 2.5312. The BIC is observed as a specific point of the leaky band dispersion with infinite life-time of the resonant mode in the Γ-point (see e.g. 15 ). The field profile of the BIC in the Si DG with cladding fluid is shown in Fig. 1 (b). The BIC clearly falls into the symmetry protected type as it is symmetrically mismatched to the zeroth order diffraction channel 15 . It is remarkable that the BIC is robust with respect to variations of n c as it can be seen from Fig. 1 (c), where we demonstrate the field profile of the BIC at kh = 2.5699 with no cladding fluid. Notice that the mode profiles are indistinguishable to the naked eye. The simulations were also run for a DG made of Ti 2 0. The obtained BIC field profiles are almost identical to those in Fig. 1 (b) and Fig. 1 (c), and, hence, are not shown in the paper. The numerical data for both DG materials are collected in Table I.

III. REFRACTIVE INDEX SENSING
Although the BIC proper is totally decoupled from the outgoing waves, it is spectrally surrounded by a leaky band with with a diverging Q-factor. If the DG is illuminated from the far-zone, this leaky band induces a collapsing Fano feature in the transmittance spectrum. Below we investigate into application of the collapsing Fano resonance for RI sensing. Two quantities are of major importance for engineering an effective sensor: the differential sensitivity defined as and the figure of merit (FOM) given by where W is the width of the operating resonance in terms of wavelength.
Here we aim at deriving analytic expressions for both quantities of interest using a perturbative approach. The perturbative approach utilizes two small parameters: the first is the increment of the dielectric function ∆ǫ c due to the change of the RI of the cladding, and the second is the angle of incidence θ defined in Fig. 1 (a). First of all let us write the solution of Eq. (1) as where β = k sin(θ) is the propagation constant along the y-axis. After substituting the above into Eq. (1) and taking into account ∆ǫ c one finds where the perturbation operators are given by Notice that L c and L θ are independently parametrized with respect to n c and θ. It means that the contributions of the two operators into the perturbed eigenfrequency are additive in the first perturbation order. First let us see the effect of L c . In general, the application of perturbative approaches to resonant states in not trivial, since the eigenfields diverge in the far-zone 32-34 . One of the possible solutions is the use of the Wigner-Brillouin (WB) perturbation theory 35,36 . The application of the WB approach requires a very specific normalization condition for the resonant states involving both "volume" convolution with dielectric function, and "flux" term as a surface integral over the outer interface of the elementary cell. It can be argued, however, that the BIC proper is a specific case when the flux term can be dropped off due to the eigenfield vanishing in the far-zone 37 . Thus, for a single BIC mode the WB perturbation theory yields vacuum wavenumber where k (0) is the vacuum wavenumber of the BIC with the reference value of the cladding RI and with S c as the area of the elementary cell occupied by the cladding fluid, and E (BIC) x as the BIC eigenfield satisfying the following normalization condition where S is the total area of the elementary cell. Notice that although the integration domains are infinite both integrals converge due to the eigenfield vanishing in the far-zone. The normalization condition (10) naturally invites applying the standard Rayleigh-Schrödinger (RS) perturbative approach following 38 . Remarkably, unlike the WB method the first order RS approach does not involve any other states rather than the unperturbed state under consideration 39 . The first order RS solution reads We stress again that although the solution (8) does not require smallness of ∆ǫ c , it is still approximate since it neglects contribution of all unperturbed modes other than the BIC. In contrast the first order non-degenerate RS solution (11) only involves a single mode but requires smallness of ∆ǫ. Notice though, that Eq. (8) and Eq. (11) coincide up to O(∆ǫ 2 c ). Thus, the Wigner-Brillouin approach truncated to a single eigenstates produces an exact result up to O(∆ǫ 2 c ). Unfortunately, unlike L c the effect of L θ can not be described by the single state WB approach. This is because the perturbed solution is a radiating state and, thus, even in the first order approximation can not be written through the non-radiating unperturbed BIC. Physically, the application of L θ generates the dispersion of the leaky band hosting the BIC in the Γ-point. Since the band is symmetric with respect to β → −β the dispersion about the Γ-point reads The complex-valued parameter α can be found by fitting to numerical or experimental data. As it was already mentioned the corrections due to L c and L θ are additive in the leading order because both θ and ∆n c are regarded as small quantities. Thus, using λ = 2π/k, Eq. (8), and Eq. (12) we obtain the resonant wavelength, λ res as Notice that Eq. (13) explicitly predicts blue shift of the resonant wavelength with positive ∆ǫ c . Applying the definition of the sensitivity, Eq. (2) one finds The above expression for sensitivity exactly coincides with that previously derived by Mortensen, Xiao, and Pedersen 38 for guided modes in a bulk photonic crystal. Remarkably, Eq. (14) can also be derived with a straightforward application of the Hellmann-Feynman (HF) theorem 40 . Notice, though, that the HF theorem strictly requires the problem to be described by Hermitian operators. In general the radiating boundary conditions break the hermiticity. However, in the case of the BIC proper the eigenfield is localized about the DG and, thus, is not affected by the boundary condition in the far-zone.
To derive the FOM we recall that the width W of a high-Q resonance can be written as where Q is the quality factor Hence, for the FOM we have Finally, let us discuss the upper limit of the differential sensitivity. By comparing the integral (9) against Eq. (10) one writes where S DG is the area of the elementary cell occupied by the DG. The quantity I c reaches its maximal value, when the integral in Eq. (18) equals to zero. Thus, the maximal differential sensitivity in the spectral vicinity of an isolated BIC is given by

IV. NUMERICAL RESULTS
To verify the above findings we run FDTD simulations for computing the transmittance, T as a function of both k and θ. In Fig. 2 (a) we demonstrate the blue shift of the Fano feature with the increase of the cladding RI. In Figs. 2 (c-d) we plotted the shift of the resonant wavelength ∆λ res = λ (0) res − λ res in comparison against Eq. (13) at three different values of θ. On can see from Figs. 2 (c-d) that the data match to a good accuracy with deviation only becoming visible at a larger angle of incidence, θ = 4, (deg), when the higher perturbation orders in θ come into play.
In Table II we present the numerical values of differential sensitivity obtained from Figs. 2 (c-d) in comparison against the analytic result Eq. (14). To assess the accuracy of Eq. (14) we calculated the relative error where S theor stands for the theoretical value of the sensitivity found from Eq. (14), while S num is used for numerical data. Again one can see a good coincidence with δ > 1% only for θ = 4, (deg) with a TiO 2 DG. Otherwise the sensitivity remains constant about the normal incidence as predicted by Eq. (14).
In Fig. 3 (a) we show the Fano feature in transmittance collapsing on approach to normal incidence. The width of the resonances is extracted from the data in Fig. 3 (a) and plotted against θ 2 in Fig. 3 (b). One can see from Fig. 3 (b) that the numerical data comply with Eq. (15) and Eq. (16) which predict W ∝ θ 2 . Finally, in Fig. 3 (c) we plot the FOM against the angle of incidence together with fitting curves FOM ∝ θ −2 . One can see that the numerical results are in accordance to Eq. (17).

V. SUMMARY AND CONCLUSIONS
The equations derived provide a cue to designing an efficient sensor with an optical BIC. The expression for sensitivity obtained in this paper is identical to that derived by Mortensen, Xiao, and Pedersen 38 for bulk photonic crystal underlying the unique properties of BICs among the leaky modes supported by the grating. According to Eq (14), and Eq. (18) the obvious approach for maximising the sensitivity is either manufacturing the DG of a dielectric with small ǫ or choosing a BIC with the eigenfield predominantly occupying the free space above the DG. Both design paradigms can be implemented by solving the eignevalue problem with FDTD method with variation of the geometric parameters.
On the other hand, the FOM can be optimized by operating at near normal incidence. At this point some comments are due in regard to FOM, Eq. (17) diverging with the decrease of the angle of incidence. Experimentally, a zero linewidth can not obtained in a realistic set-up due to three major factors: (i) material losses due to absorption in the dielectric 41,42 ; (ii) fabrication inaccuracies 43 ; and (iii) finiteness of the structure 41,44 , as the BIC proper can not be supported by systems finite in all spatial dimensions 45 .
The effect of all three factors is saturation of FOM on approach to normal incidence. Which of the three will dominate in the staturation is determined by a specific design of the DG.
Another important aspect is the finite waist of the Gaussian beam illuminating the DG. Equations (14), and (17) can only be strictly applied to a plane wave, whereas a Gaussian beam is always a continuum of plane waves propagating to slightly different directions. Miniaturization of the sensor would imply a tightly focused beam with a broad momentum distribution in the Fourier plane. This would naturally compromise the FOM, as the scattering channels with larger angles of incidence become more populated.
Finally, let us discuss whether the limit, Eq. (19) can be overcome. First of all, Eq. (19) is derived under assumption of smallness of the perturbation operators. Therefore, the sensitivities larger than Eq. (19) are not strictly forbidden. For example, a larger sensitivity has been recently reported with surface plasmon-polaritons 46 . The drawback, however, is the drop of the quality factor due to strong coupling of the eigenmode to the outgoing waves. Thus, the increase of sensitivity with departure from the BIC proper can be only obtained at the cost of the sensor FOM. More promising is optimizing the sensitivity while keeping the high FOM intact. This can be done by breaking the second condition in derivation of Eq. (19), namely, the single mode approach. One possible route is looking for a BIC in the exceptional points [47][48][49] , which unavoidably emerge as a result of an intricate interplay between two eigenmodes. The other possible route is application of BIC existing just below the first Wood's anomaly in which case the contribution of the evanescent channels corresponding to the first diffraction order becomes extremely sensitive to the system's parameters 21 . The destruction of a BIC in crossing the Wood's anomaly inherently implies application of a multimodal perturbative approach since the contribution of states with finite live-time is required. We speculate that operating in the vicinity of the Wood's anomaly may result in a significant gain in sensitivity while keeping diverging the FOM and quality factor.