Modelling Purcell enhancement of metasurfaces supporting quasi-bound states in the continuum

3.1 Lossless q-BIC on TiO₂ metasurface

We can first get a basic understanding of the behavior of the PLE of q-BIC by considering an ideal, lossless metasurface that supports a symmetry-protected BIC. The detailed design of the metasurface is described in the supplementary materials [41]. The angle-dependent transmissivity T and reflectivity R of the metasurface were numerically simulated and plotted in Figure 2(a) and (b). We can clearly observe the symmetry-protected BIC at incident angle θ = 0° and wavelength λ = 612.6 nm, and as detuning θ was introduced, the q-BIC emerges with a ±θ symmetry. The nearfield in the dye medium was also recorded and the numerically predicted PLE was calculated directly by Eq. (1) and plotted in Figure 2(c). The T and R of each incident angle were fitted to obtain the spectral parameters ω ₀, Γ_tot and Γ_rad,in, with Eq. (S4) and (S5), respectively [41]. The extinction coefficient κ of the dye medium was modulated to calibrate the factor Γ_abs,dye/κ, which is the slope of Γ_abs against κ [41].

Figure 2:

The numerically simulated (a) transmissivity, (b) reflectivity and (c) PLE of the lossless TiO₂ metasurface. The (d) resonant wavelength and (e) decay rates Γ_tot, Γ_rad,in and Γ_abs are plotted as a function of the incident angle θ. (f) The numerical PLE and predicted PLE are plotted as a function of the emission angle θ. The PLE predicted with a small artificial Γ_abs are plotted as dashed lines.

Since we focus on analyzing the behavior of the metasurface within a small detuning from the BIC, we can model the change in spectral parameters by a polynomial expansion, specifically up to the quadratic term. The ω ₀ and the decay rates are plotted in Figure 2(d) and (e) with their polynomial fit. As illustrated in Figure 2(e), Γ_abs remains zero for all incident angles while Γ_tot and Γ_rad,in increases quadratically from the BIC point. On the other hand, we can see that the factor Γ_abs,dye/κ remain mostly constant throughout the fitted range with a slight increase at larger θ (see Fig. S1(f)), which indicates similar nearfield confinement at different θ. As a result, the total Q factor diverges at the BIC and decreases following the inverse-squared law of the detuning (θ) [42], [43], [44]. We then predict the PLE with the best fit of the parameters and compare with the numerical PLE in Figure 2(f). The PLE increases rapidly as we approach the BIC and appears to diverge at the BIC. Since non-radiative decay rates remain zero while Γ_tot and Γ_rad,in both varies quadratically against the incident angle, the out-coupling efficiency Γ_rad,in/Γ_tot remains constant over all incident angles and did not influence (or reduce) the PLE.

However, as illustrated in Figure 2(f), the behavior is significantly different once we artificially introduce a small Γ_abs, which models a non-zero absorption loss and/or random scattering loss. Since the Γ_abs increases Γ_tot slightly while not affecting Γ_rad,in, the out-coupling efficiency Γ_rad,in/Γ_tot is no longer a constant but drops sharply around the BIC, influencing the PLE to hit a maximum and drop sharply around the BIC. We can see this effect even with very small artificial Γ_abs equivalent to a linewidth broadening of 0.0001 nm, which is difficult to observe with conventional spectrometers. This prompts us to explore in more detail the effect of absorption and random scattering loss.

3.2 Lossy q-BIC on Si metasurface

In the previous section, we revealed that even a very small absorption loss from the metasurface would inevitably change the characteristic behavior of PLE from BIC. However, material loss is inevitable in realistic metasurfaces [45]; nanofabrication also introduces inherent roughness into the metasurface structure, which introduces random scattering loss in application [46]. This prompts us to take a closer look into BIC supported on realistic materials with loss instead of an ideal, lossless model. Therefore, we consider a symmetry-protected BIC supported on a bipartite Si metasurface, as described in the supplementary materials [41]. The symmetry-protected BIC at λ = 704.5 nm was detuned by displacing one of the Si nanoparticles by d, and the d-dependent transmissivity T, reflectivity R, absorptivity A = 1 − T − R and PLE were numerically simulated and plotted in Figure 3(a)–(d). The T, R and A were fitted to obtain the spectral parameters ω ₀, Γ_tot and Γ_rad,in, while Γ_abs,dye/κ was fitted by modulating κ of the dye medium [41].

Figure 3:

The numerically simulated (a) transmissivity, (b) reflectivity, (c) absorptivity and (d) PLE of the bipartite Si metasurface. (e) The decay rates Γ_tot, Γ_rad,in and Γ_abs are plotted as a function of the detuning d. (f) The numerical PLE and predicted PLE are plotted as a function of the detuning d.

As shown in Figure 3(e), Γ_abs remains almost constant against d with a small decrease at increasing d. This decrease is due to the gradual redshift of ω ₀ as d increases and Si being less absorptive at longer wavelengths. On the other hand, Γ_rad,in = Γ_rad is zero at the BIC and increases quadratically from there. This is like what we observed in the lossless case, and we expect a quadratic dependency due to the limited range and the symmetry about ±d. Therefore, Γ_tot, the sum of Γ_abs and Γ_rad, roughly follows the same shape of Γ_rad but is displaced upwards by Γ_abs. Since d only shifts the position of the Si nanoparticles and does not increase or reduce the amount of dye medium, the factor Γ_abs,dye/κ also remains mostly constant against d, as illustrated in Fig. S2(e) [41]. This allows us to assume Γ_abs,dye/κ is constant and use the value at d = 20 nm for the prediction of PLE. As illustrated in Figure 3(f), the numerically simulated PLE aligned with the PLE predicted by Eq. (5). First, the PLE increases when reducing d and approaching the BIC due to the increase in Q factor. However, as we get closer to the BIC configuration, the decrease in the out-coupling efficiency overwhelms the improvement in Q factor and the PLE drops sharply near the BIC.

Interestingly, the maximum in PLE occurs very close to the point at which Γ_abs and Γ_rad intersects. This coincides with the critical coupling condition, which states that light–matter interaction is maximized when Γ_abs = Γ_rad [32]. This also presents a new perspective on how to optimize Purcell enhancement with realistic metasurfaces. First, Γ_abs should be as small as possible, which can be achieved by using low-loss materials for the metasurface [45], or reducing plasmonic loss with nonlocal resonances [43], [47]. The Γ_abs sets a lower limit to Γ_tot and thereby also imposes an upper bound to the Q factor. Then, Γ_rad should be tuned to equal Γ_abs to achieve the critical coupling condition and maximize light–matter interaction. If Γ_rad is too small, the system is dominated by the absorptive loss and most of the enhanced emission would be absorbed by the metasurface instead of out-coupled. While for large Γ_rad, the increase in out-coupling efficiency is not worth the decrease in Q factor and therefore the PLE is reduced.

3.3 Asymmetric TiO₂ metasurface supporting accidental q-BIC

In the following sections, we explore the behavior of an off-Γ accidental q-BIC on experimental and simulated metasurfaces. We fabricated asymmetric TiO₂ metasurfaces that support accidental q-BICs by glancing angle deposition of Ti followed by rapid thermal annealing (RTA) on symmetric TiO₂ metasurfaces [41]. A numerical model is built to replicate the behavior of the experimental asymmetric TiO₂ metasurface, which the design is described in the supplementary materials [41]. Since the TiO₂ was formed through RTA of Ti, we expect some roughness in the metasurface and the random scattering loss is modelled by adding an extinction coefficient κ Ti O 2 = 0.005 to the TiO₂ material model. While the numerical simulation only simulates the coherent and periodic part of the EM field, the randomly scattered photons are decoherent from the resonant field. Therefore, we choose to model the random scattering loss as absorption as both absorbed photons and randomly scattered photons would not be part of the simulated field. We also assumed a uniform scattering density within the TiO₂ material, which is intended to account for both the surface roughness of the material as well as structural inhomogeneity within the material. As we can see in the angle-dependent T, R, A and PLE plotted in Figure 4(a)–(d), an accidental BIC is observed at an oblique incident angle θ = 0.8° and λ = 614 nm. The T, R and A were fitted to obtain the spectral parameters ω ₀, Γ_tot and Γ_rad,in, and Γ_abs,dye/κ was also fitted [41].

Figure 4:

The numerically simulated (a) transmissivity, (b) reflectivity, (c) absorptivity and (d) PLE of the asymmetric TiO₂ metasurface. (e) The decay rates Γ_tot, Γ_rad, Γ_abs, Γ_rad,↑ and Γ_rad,↓ are plotted as a function of the incident angle θ. (f) The numerical PLE and predicted PLE are plotted as a function of the incident angle θ.

The decay rates of the accidental q-BIC on the asymmetric TiO₂ metasurface behave markedly differently from that of the symmetric cases. As illustrated in Figure 4(e), Γ_abs is largest at normal incident and decreases following a linear style when increasing |θ|. The Γ_rad is also symmetrical over ±θ and remains non-zero at all θ. On the other hand, the asymmetry in the metasurface splits the point where Γ_rad,↑ and Γ_rad,↓ touch zero, corresponding to the ports at the superstrate (↑) and substrate (↓) sides, towards positive and negative incident angles respectively. Γ_rad,↑ and Γ_rad,↓ now touch zero at 0.92° and −0.53° respectively. To use the parameters to predict the PLE, the decay rates were fitted with fourth-order polynomials. The polynomial fits of the parameters were then substituted into Eq. (5) to derive a prediction of the PLE. As illustrated in Figure 4(f), Eq. (5) predicted the numerical PLE with good accuracy (see Fig. S3(e) for the simulated substrate side PLE) [41]. The PLE at the superstrate side shows a minimum at 0.9° while that at the substrate side shows a minimum at −0.5°. A local maximum in PLE is observed at −2.5° at the superstrate side and the model predicted it at −2.2°. On the substrate side, a local maximum is observed at 3.2° followed by a very slow decrease in PLE, while the model predicts a local maximum at 2.5°. Overall, we see larger discrepancies between the analytical and numerical PLE near the end of the fitted range. Due to the nature of polynomial regression, the edge of the data is usually less well represented compared to the data near the center (normal incident), as the higher-order terms are truncated and they are only significant for larger |θ|. Therefore, in the trade-off between overfitting and poor fitting accuracy, we proceeded with fourth-order polynomials as they fitted the decay rates better when compared to second-order polynomials, while further increasing the polynomial order did not improve significantly upon that. While this choice is unlikely to lead us to overfitting, this also limits our accuracy, especially for the data near the edges of the fitted range.

It is worth noting that R remained symmetric over ±θ, even on an asymmetric structure. We were also able to predict the PLE at both the superstrate side and the substrate side, only with the spectra incident from the superstrate side. Both features are because Lorentz reciprocity ensures the in-coupling and out-coupling process are symmetric, which ensures a symmetric R. Also, by simultaneously fitting T and R of ±θ, we were able to determine the in-coupling constants of all ports and thus model the PLE at both the superstrate side and the substrate side.

3.4 Asymmetric TiO₂ metasurface fabricated by glancing angle deposition

In this section, we experimentally investigate the behavior of accidental q-BIC on the asymmetric TiO₂ metasurfaces. The asymmetric metasurface was fabricated by glancing angle deposition on a symmetric TiO₂ metasurface. A Ti thin film (100 nm) was deposited on a SiO₂ glass substrate by electron beam deposition, and a resist (TU7, Obducat) was spin-coated on top. A resist array with square nanoparticles of sides 160 nm arranged in a square lattice of period 370 nm was then fabricated by nanoimprinting. The pattern was transferred to the Ti thin film by O₂ ashing (RIE-10NR) and etching with Cl₂ and CHF₃ (NLD-570). The Ti metasurface was then oxidized by RTA (900 °C, 10 min) under O₂ atmosphere to fabricate a symmetric TiO₂ metasurface. A 20 nm thin film of Ti was deposited onto the TiO₂ metasurface from a glancing angle by tilting the substrate in an electron beam deposition apparatus. The substrate was tilted by 30° along the xz-plane. The newly deposited Ti was also oxidized with the same RTA process to create the asymmetric TiO₂ metasurface. As highlighted in the SEM image of the asymmetric metasurface shown in the inset of Figure 5(d), voids formed only on the left-hand side of the TiO₂ nanoparticles, at the shadow of the nanoparticles during glancing angle deposition. Finally, a PMMA layer (450 nm) was spin-coated onto the TiO₂ metasurface. The PMMA is doped with an organic dye (Lumogen F 305 Red) at 1 wt% (weight percent), and the PLE on this Lumogen dye is evaluated. The T, R and PLE (I/I ₀) were measured on the rotation stages illustrated in the supplementary materials [41]. The measured T, R and PLE (I/I ₀) are plotted in Figure 5(a)–(c). We focus on the branch showing an accidental BIC at θ = 1.2° and λ = 710 nm. The spectral parameters ω ₀, Γ_tot and Γ_rad,in were fitted from the measured spectra directly. However, to calibrate the nearfield confinement with Γ_abs,dye/κ, we need to use a dye that absorbs light at around λ = 700–800 nm. In place of the PMMA layer doped with Lumogen, we spin-coated a similar PMMA layer with 0.25 wt% of IR780 iodide onto the asymmetric TiO₂ metasurface. The T and R of the metasurface with IR780 are then measured to fit for (mainly) Γ_abs and thus we derive Γ_abs,dye/κ [41].

Figure 5:

The measured (a) transmissivity, (b) reflectivity, and (c) PLE (I/I ₀) of the asymmetric TiO₂ metasurface fabricated by glancing angle deposition. The arrows indicate the mode of interest. The (d) resonant wavelength and (e) decay rates Γ_tot, Γ_rad, Γ_abs, Γ_rad,↑ and Γ_rad,↓ are plotted as a function of the incident angle θ. The inset in (d) highlights the asymmetry in the TiO₂ metasurface. The void formed under the shadow of TiO₂ nanoparticles during glancing angle deposition are outlined by the dashed lines. (f) The re-normalized PLE and predicted PLE are plotted as a function of the emission angle θ.

The decay rates behaved like the simulated asymmetric TiO₂ metasurfaces. The Γ_abs is largest at normal incident and decreases slightly at larger |θ|. As illustrated in Figure 5(e), Γ_rad remains symmetric over ±θ, while Γ_rad,↑ and Γ_rad,↓ touches zero at 1.3° and −1.2° respectively. The parameters also show a good fit when approximated with fourth-order polynomials. We then use the polynomial fits to derive the predicted PLE through Eq. (5). While our analytical model only predicts enhancement due to the Purcell effect, enhancement in the absorption by the Lumogen dye due to the presence of the metasurface is not described, yet is included in the measured I/I ₀. To minimize the influence of the absorption enhancement and the uneven excitation profile associated with higher absorption enhancement, the incident angle of the excitation laser was chosen at −31.5°, and the measured PLE was re-normalized against the baseline to isolate the Purcell effect. We can now compare the predicted PLE against the re-normalized PLE. As shown in Figure 5(f), the experimental Purcell enhancement is predicted by our analytical model. We also apply Lorentz reciprocity to predict the backside PLE, which our experimental measurements show good agreement with (see Fig. S4(c) for the measured backside PLE) [41].