[a,b]Deivid Ribeiro

# TeV Detection of the Extreme HSP Blazar RBS 1366 by VERITAS

###### Abstract

Extreme high-synchrotron-peak blazars (EHSPs) are postulated as the most
efficient and extreme particle accelerators in the universe but remain
enigmatic as a possible new class of TeV gamma-ray blazars. Blazars are active
galactic nuclei (AGNs) with jets of relativistic particles that generate non-
thermal emission pointed along the line-of-sight. Their spectral energy
distribution (SED) are characterized by synchrotron and inverse-Compton peaks,
indicating acceleration of leptonic and possibly hadronic particle populations
in the jet. EHSPs are characterized by a peak synchrotron frequency
$>${10}^{17}\text{\,}\mathrm{H}\mathrm{z}$$ with their Compton peak expected
to fall in the TeV range. Indeed, the handful of EHSPs detected by Imaging Air
Cherenkov Telescopes (IACTs) have presented challenges where some may be a
high-frequency extension of the blazar sequence while others peaking around 10
TeV may represent a different class of TeV emitters. Detections of the high-
energy and very-high-energy (HE; E > 100 MeV, VHE; E > 100 GeV) components of
the Compton peak will play an important role in constraining the acceleration
model derived from the SED. We present the discovery of TeV emission from RBS
1366, a candidate EHSP, by the VERITAS observatory. Using HE and VHE data from
the Fermi-LAT and VERITAS observatories, respectively, we characterize the
detection by providing an SED and model fit in the context of other EHSP
candidates. Our work confirms the status of RBS 1366 as an EHBL.

##  1 Introduction

Blazars are a subclass of active galactic nuclei with relativistic jets
pointed toward the observer, emitting gamma rays in the very-high-energy (VHE;
E>100 GeV) or "TeV" regime. The underlying mechanisms for the emission of
these gamma rays are evident in the observed spectral energy distribution
(SED), which is modeled to include the underlying particle populations, and
acceleration and cooling mechanisms in the jets.

It is commonly understood that the low energy and high energy peaks of the SED
are produced by synchrotron and inverse Compton processes, respectively [1].
In this synchrotron self-compton model (SSC), a population of electrons emit
synchrotron radiation and then the same electron population emit gamma rays
through Inverse Compton scattering the synchrotron photons.

At a redshift of $z=0.2365$, RBS 1366 is an AGN in giant elliptical host
galaxy with a central black hole mass of $9.31\pm 0.32$ M⊙ [2]. RBS 1366 has
been a prime candidate for detection and classification as an EHSP [3, 4, 5].
It was selected for VERITAS observation based on the high synchrotron peak
frequency in the Swift XRT X-ray band ($>10^{17}$ Hz), along with detection in
the 20 MeV to 300 GeV Fermi-LAT band displaying a hard spectrum (and placement
in the 3FHL catalog [6]). RBS 1366 was expected to behave very much like the
EHBL 1ES 0229+200 (see Table 1 of [5]) with a measured a redshift of z = 0.237
based on Ca II, G band, Fe I, Mg I and Na absorption [7]. In the 1-300 GeV
band, Toomey (2020) [3] found a power law fit with normalization
$k=$7.20[165]\text{\times}{10}^{-11}\text{\,}\mathrm{c}\mathrm{m}^{-2}\mathrm{s}^{-1}\mathrm{G}\mathrm{e}\mathrm{V}^{-1}$$
and index $\gamma=$1.63[8]$$.

It is notable that a VERITAS upper limit has been reported [8, 5],
anticipating that a detection of this source with an accompanying spectrum
would enable further EHSP characterization. A differential upper limit of
$1.7\text{\times}{10}^{-11}\text{\,}\mathrm{c}\mathrm{m}^{-1}\mathrm{s}^{-1}\mathrm{e}\mathrm{r}\mathrm{g}^{-1}$
at $327\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ was calculated for this source
by VERITAS based on 10 hours of observations in 2016 [8].

##  2 Methods

###  2.1 Fermi-LAT

The Large Area Telescope (LAT) on board the Fermi satellite has operated since
2008 [9]. It is sensitive to photons between ${\sim}20\text{\,}\mathrm{MeV}$
and ${\sim}1\text{\,}\mathrm{TeV}$ and has an ${\sim}60\degree$ field of view,
enabling it to survey the entire sky in about 3 hours.

RBS 1366 was observed from 2008 to 2023, where the data was fit with a power
law as the base spectral model. The Fermi-LAT 4FGL catalog defines this source
as J1417.9+2543 located at $14^{h}17^{m}58.6^{s}$, $None$.

The publicly available Fermi-LAT data were analyzed using with the Fermitools
suite of tools provided by the Fermi Science Support Center (FSSC). Using the
Fermipy analysis package [10]111https://fermipy.readthedocs.io/en/latest/ ;
v1.2, the data were prepared for a binned likelihood analysis in which a
spatial spectral model is fit over the energy bins. The data were selected
using the SOURCE class of events, which are optimized for point-source
analysis, within an angle of $15\degree$ from the analysis target position. A
$90\degree$ zenith angle cut was applied to remove any external background
events due to the effect of the Earth. The standard background models were
applied to the test model, incorporating an isotropic background and a
galactic diffuse emission model without any modifications ($gll\\_iem\\_v07$
and $iso\\_P8R3\\_SOURCE\\_V2\\_v1$). The standard 4FGL catalog was then
queried for sources within the field of view and their default model
parameters [11]. With the improvements to Fermi-LAT low-energy sensitivity in
PASS8 reconstruction, the energy range was expanded to $100$ MeV - $1$ TeV.

###  2.2 VERITAS

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an
Imaging Atmospheric Cherenkov Telescope (IACT) array consisting of four 12 m
telescopes separated by approximately $100\text{\,}\mathrm{m}$, at the Fred
Lawrence Whipple Observatory (FLWO) in southern Arizona, USA [12, 13]. The
observatory is sensitive to photons within the energy range
$\backsim$$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ to
$\backsim$$30\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, with the ability to
detect 1% of the emission of the Crab Nebula in 25 hours (at $5\sigma$). The
instrument has an angular resolution (68% containment) of $\backsim$0.1° at
$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, an energy resolution of
$\backsim$15% at $1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, and 3.5° field of
view.

VERITAS observed RBS 1366 for $56\text{\,}\mathrm{h}\mathrm{r}$ between 2008
and 2021, under dark sky conditions. The data in this paper were taken using
“wobble" pointing mode, where the source is offset from the center of the
camera by $0.5\degree$. This mode creates space for a radially symmetric "off"
region to be used for background estimation in the same field of view, saving
time from targeted background observations that contain the same data
observing conditions.

The data were processed with standard VERITAS calibration and reconstruction
pipelines, and then cross-checked with a separate analysis chain [14, 15].
Specifically, we used an Image Template Method (ITM) to improve event angular
and energy reconstruction [16], where analysis cuts are determined with a set
of a priori data-selection cuts optimized on sources with a soft power law
index (from 2.5 to 3).

##  3 Results

###  3.1 Fermi-LAT

Using the Fermi-LAT data from this work, RBS 1366 (as 4FGL J1417+2543) was
detected with a test statistic TS=$353.9$ ($\sim 18\sigma$) over an
observation period of 14.5 years in the energy range from 100 GeV to 1 TeV.
The source was fit to a power law model $dn/dE=(4.13\pm 0.76\times
10^{-15})(E/10.4\mathrm{~{}GeV})^{-1.6\pm 0.1}$ MeV-1 cm-2 s-1. The spectrum
is shown in the left of Figure 1.

RBS 1366 was also tested for time variability by constructing a light curve
over the entire observing period in 6 month bins, integrating over the same
energy range of 100 GeV to 1 TeV. The light curve is binned into periods
determined by the Bayesian Blocks algorithm where a false alarm rate of
$p0=0.0027$ (equivalent to 3 sigma) is used to compute the prior of the number
of bins, yielding 3 bins of relative stability, see Figure 2[17]. This is
consistent with the 4FGL catalog entry for RBS 1366, where the variability
index is measured to be 17.4 (an index above 24.72 over 12 intervals indicates
<1% chance of being a steady source [18]).

###  3.2 VERITAS

VERITAS has detected RBS 1366 at a significance of 6.5 sigma or standard
deviation, with an exposure of 56.8 hours. A summary of the detection is shown
in Table 1. The observed spectrum was fit to a power law model, shown in
Figure 1.

Table 1: Summary of VERITAS detection. The quality-selected live time, number of gamma-ray-like events in the on- and off-source regions, the normalization for the larger off-source region, the observed excess of gamma-rays and the corresponding statistical significance are shown. For each observation epoch, the integral flux corresponding to the observed excess is given. The flux is reported above the observation threshold of 200 GeV, and is also given in percentage of Crab Nebula flux above the same threshold. Epoch | T | On | Off | Norm | Excess | Significance | Flux (>200 GeV) | Crab
---|---|---|---|---|---|---|---|---
| [hr] |  |  |  |  | [$\sigma$] | [${10}^{-12}\text{\,}\mathrm{c}\mathrm{m}^{-2}\mathrm{s}^{-1}$] | %
2008-2022 | 57 | 3006 | 29180 | 0.0909 | 353 | 6.53 | $1.7\pm 0.5$ | 0.5

An effective energy threshold of $\sim
200\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ is calculated. The spectrum is fit
between $200\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and
$2\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ where the photon count in each bin
is above 10 excess counts and the significance is above $1\sigma$, since this
source is weakly detected. Upper limit points were calculated to a confidence
limit of 99%. The spectrum is fit with a power law, where upper limits are not
included. The power law fit gives ($2.39\times 10^{-12}\pm 5.2\times
10^{-13})\times(E/400\mathrm{~{}GeV})^{(-2.7\pm 0.5)}$ TeV^-1 cm^-1 s^-1.

Figure 1: Left) Differential spectrum for Fermi-LAT source 4FGL J1417+2543,
fit using a power law. Right) VERITAS spectrum. A TS threshold was set to
$1\sigma$. Grey open square: observed differential spectrum with power law fit
to data. Red circle: EBL deabsorbed differential spectrum using the Dominguez
model [19], for redshift z=0.237 [7]. For low significance bins (TS<1), the
Rolke upper limit is derived for a power law of index 2.5, and is not included
in the overall fit of the larger energy range [20].  Figure 2: Light curve for
VERITAS (red) and Fermi-LAT (blue) observations. VERITAS data is binned into
30 day bins, and the flux is calculated above$350\text{\,}\mathrm{GeV}$. The
Fermi-LAT data is binned into 6 month bins and the flux is calculated above
$100\text{\,}\mathrm{MeV}$. Dotted vertical lines mark the VERITAS observation
on the Fermi-LAT light curve. A significance threshold of $3\sigma$ was used
to determine upper limits, and the Bayesian Block algorithm was used with a
false positive probability equivalent to $3\sigma$ to determine bin edges in
the Fermi-LAT light curve.

VHE photons are absorbed by the extragalactic background light (EBL)
throughout the universe; the flux must be corrected to account for this
effect. This absorption is energy and redshift dependent. Deabsorption is
applied to the flux using the model of [19]. After deabsorbing the spectrum
from the EBL, the power law fit is $(7.49\times 10^{-12}\pm 1.7\times
10^{-12})\times(E/400\mathrm{~{}GeV})^{(-1.1\pm 0.6)}$ TeV^-1 cm^-1 s^-1.
Although these flux points are harder than the overlapping Fermi-LAT data,
they are consistent within error. Both observed and EBL deabsorbed spectral
points are shown in Figure 1.

The combined HE and VHE spectral fit is shown in Figure 3. Although the
VERITAS flux point are hardened by the lowest energy point at 200 GeV, the
combined power law fit is consistent within error with a $\chi^{2}$ statistic
of 13.49 (ndof=12). These fit results are summarized in Table 2.

Figure 3: Combined Fermi-LAT and VERITAS spectrum with a power law fit (dotted line). The peak of the Compton component is not clear, but is above 0.1 TeV and possibly higher. Table 2: Power law spectral fit to HE and VHE data. Parameter | Value | Uncertainty | Unit
---|---|---|---
Index | 1.62 | 0.05 |
Amplitude | $9.0$ | $1.4$ |  $10^{-12}$ cm-2 s-1 TeV-1
Reference | 0.4 | - | TeV

In addition to the VERITAS SED, a light curve was also computed over 30 day
bins, integrating above 300 GeV. This is shown in Figure 2 along with the
corresponding Fermi-LAT light curve.

##  4 Discussion & Conclusion

The source exhibited a hard HE+VHE spectrum (index < 2) after correcting for
EBL absorption. With a synchrotron peak at log frequency of $17.7\pm 0.2$
($\sim 2\text{\,}\mathrm{keV}$), a redshift of $0.237$, and the known hardness
of the Fermi-LAT power law, the evidence already points toward an EHSP
classification. The combined HE+VHE fit in this work supports this
classification and merits investigation into the subclassification posited by
[5]. The power law index $1.62\pm 0.05$ on the combined EBL-deabsorbed HE+VHE
spectrum in Figure 3 gives a VHE gamma-ray slope ($S=2-\gamma$ where $\gamma$
is power law index) of $S=0.38$, which places RBS 1366 similar to other
blazars such as 1ES 1101-232, 1ES 0347-121 and 3FGL J0710+5808, which were
labeled as Hard-TeV EHBLs in [5]]. A full multi-wavelength-SED will enable a
complete classification in a future study.

The HE-VHE luminosity of RBS 1366 appears to peak above an energy of $\sim 1$
TeV, where it reaches an intrinsic luminosity of $3\times 10^{44}$ erg/s
(assumed to be isotropic). However, since the turnover of the Compton
component is not yet evident, we are not able to claim that the Compton
component is not dominant over the synchrotron component, which peaks at a log
frequency of $17.7\pm 0.2$ ($\sim 2\text{\,}\mathrm{keV}$) with a luminosity
of $\sim 10^{45}$ erg/s [5]. A more sensitive observation from the Cherenkov
Telescope Array (CTA), a next generation IACT observatory, in the future may
yield a more comprehensive estimation of the Compton peak with improved
sensitivity above $1$ TeV [21, 22].

A basic estimate of the relativistic boosting in the SSC model fit to
available published data is very high, between $\delta\sim 50-100$, implying a
significantly efficient acceleration that is common to EHSPs. This indicates
the emission region is probably in the jet where the emission region is in the
range of $10^{17}$ cm and the magnetic field is very weak. While variability
has been noted in the optical band [8], the lack of variability in the HE-VHE
band is expected for the rising end of the Compton peak seen here.

The observational circumstances above all support the conclusion that RBS 1366
is an EHSP now detected in the TeV energy range. Figure 4 is a reproduction of
the index versus synchrotron peak frequency plot sampling the other TeV
detection EHSPs in [5], which cluster into Hard-TeV and HBL-like blazars.
Using the results presented in this work, RBS 1366 may be classified as a
Hard-TeV EHSP. A complete study using a full multi-wavelength SED to confirm
the Hard-TeV EHSP classification is underway.

Figure 4: Using a selection of TeV detected EHSPs, the HE-VHE power law index
in the range of 0.1 to 10 TeV vs the log frequency of the synchrotron peak.
RBS 1366 is shown in green. At these energies, EHSPs diverge into Hard-TeV
blazars with high synchrotron peaks (clustered to the right), and HBL-like
variable blazars with low synchrotron peak and softer spectrum (clustered to
the left). This is reproduced from [5].

## Acknowledgments

This work was partially supported by NSF award PHY 2110737. This research is
supported by grants from the U.S. Department of Energy Office of Science, the
U.S. National Science Foundation and the Smithsonian Institution, by NSERC in
Canada, and by the Helmholtz Association in Germany. This research used
resources provided by the Open Science Grid, which is supported by the
National Science Foundation and the U.S. Department of Energy’s Office of
Science, and resources of the National Energy Research Scientific Computing
Center (NERSC), a U.S. Department of Energy Office of Science User Facility
operated under Contract No. DE-AC02-05CH11231. We acknowledge the excellent
work of the technical support staff at the Fred Lawrence Whipple Observatory
and at the collaborating institutions in the construction and operation of the
instrument.

## References

  * [1] Tavecchio, F., Maraschi, L. and Ghisellini, G., APJ 509 (Dec., 1998) 608–619.
  * [2] Keenan, M. et al, MNRAS 505 (June, 2021) 4726–4745.
  * [3] Toomey, M.W., Oikonomou, F. and Murase, K., MNRAS 497 (Sept., 2020) 2455–2468. arXiv:2005.11581 [astro-ph].
  * [4] Costamante, L. and Ghisellini, G., A & A 384 (Mar., 2002) 56–71.
  * [5] Foffano, L. et al, MNRAS null (2019) null. tex.arxivid: 1903.07972.
  * [6] Ajello, M. et al, APJ Supplement Series 232 (Sept., 2017) 18.
  * [7] Halpern, J.P. et al, APJ 302 (Mar., 1986) 711. ADS Bibcode: 1986ApJ…302..711H.
  * [8] Archambault, S. et al, The Astronomical Journal 151 (May, 2016) 142\.
  * [9] Atwood, W.B. et al, Astrophysical Journal 697 (2009) 1071–1102. Publisher: Institute of Physics Publishing.
  * [10] Wood, M. et al in PoS, vol. 301 of International Cosmic Ray Conference, p. 824, Jan., 2017.
  * [11] Abdollahi, S. et al, APJ Supplement Series 247 (Mar., 2020) 33.
  * [12] Weekes, T. et al, Astroparticle Physics 17 (May, 2002) 221–243.
  * [13] Holder, J. et al, Astroparticle Physics 25 (July, 2006) 391–401.
  * [14] Cogan, P., VEGAS, the VERITAS Gamma-ray Analysis Suite, in Proceedings of the 30th International Cosmic Ray Conference, ICRC 2007, vol. 3, pp. 1385–1388, 2008.
  * [15] Maier, G. and Holder, J. in PoS, vol. 301 of 35th International Cosmic Ray Conference (ICRC2017), p. 747, Sissa Medialab Srl, Jan., 2017\.
  * [16] Christiansen, J. in PoS, (Trieste, Italy), p. 789, Sissa Medialab, Aug., 2017\.  ISSN: 18248039.
  * [17] Scargle, J.D. et al, APJ 764 (Feb., 2013) 167. ADS Bibcode: 2013ApJ…764..167S.
  * [18] Abdollahi, S. et al, APJ Supplement Series 260 (June, 2022) 53.
  * [19] Domínguez, A. et al, MNRAS 410 (Feb., 2011) 2556–2578. Publisher: Oxford Academic.
  * [20] Rolke, W.A., López, A.M. and Conrad, J., Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 551 (Mar., 2005) 493–503.
  * [21] Funk, S. and Hinton, J.A., Astroparticle Physics 43 (Mar., 2013) 348–355.
  * [22] The CTA Consortium, Science with the Cherenkov Telescope Array.  World Scientific, Mar., 2019.

## Full Author List: VERITAS Collaboration

A. Acharyya1, C. B. Adams2, A. Archer3, P. Bangale4, J. T. Bartkoske5, P.
Batista6, W. Benbow7, J. L. Christiansen8, A. J. Chromey7, A. Duerr5, M.
Errando9, Q. Feng7, G. M. Foote4, L. Fortson10, A. Furniss11,12, W. Hanlon7,
O. Hervet12, C. E. Hinrichs7,13, J. Hoang12, J. Holder4, Z. Hughes9, T. B.
Humensky14,15, W. Jin1, M. N. Johnson12, M. Kertzman3, M. Kherlakian6, D.
Kieda5, T. K. Kleiner6, N. Korzoun4, S. Kumar14, M. J. Lang16, M. Lundy17, G.
Maier6, C. E McGrath18, M. J. Millard19, C. L. Mooney4, P. Moriarty16, R.
Mukherjee20, S. O’Brien17,21, R. A. Ong22, N. Park23, C. Poggemann8, M.
Pohl24,6, E. Pueschel6, J. Quinn18, P. L. Rabinowitz9, K. Ragan17, P. T.
Reynolds25, D. Ribeiro10, E. Roache7, J. L. Ryan22, I. Sadeh6, L. Saha7, M.
Santander1, G. H. Sembroski26, R. Shang20, M. Splettstoesser12, A. K.
Talluri10, J. V. Tucci27, V. V. Vassiliev22, A. Weinstein28, D. A. Williams12,
S. L. Wong17, and J. Woo29

1Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL
35487, USA

2Physics Department, Columbia University, New York, NY 10027, USA

3Department of Physics and Astronomy, DePauw University, Greencastle, IN
46135-0037, USA

4Department of Physics and Astronomy and the Bartol Research Institute,
University of Delaware, Newark, DE 19716, USA

5Department of Physics and Astronomy, University of Utah, Salt Lake City, UT
84112, USA

6DESY, Platanenallee 6, 15738 Zeuthen, Germany

7Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA

8Physics Department, California Polytechnic State University, San Luis Obispo,
CA 94307, USA

9Department of Physics, Washington University, St. Louis, MO 63130, USA

10School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA

11Department of Physics, California State University - East Bay, Hayward, CA
94542, USA

12Santa Cruz Institute for Particle Physics and Department of Physics,
University of California, Santa Cruz, CA 95064, USA

13Department of Physics and Astronomy, Dartmouth College, 6127 Wilder
Laboratory, Hanover, NH 03755 USA

14Department of Physics, University of Maryland, College Park, MD, USA

15NASA GSFC, Greenbelt, MD 20771, USA

16School of Natural Sciences, University of Galway, University Road, Galway,
H91 TK33, Ireland

17Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

18School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

19Department of Physics and Astronomy, University of Iowa, Van Allen Hall,
Iowa City, IA 52242, USA

20Department of Physics and Astronomy, Barnard College, Columbia University,
NY 10027, USA

21 Arthur B. McDonald Canadian Astroparticle Physics Research Institute, 64
Bader Lane, Queen’s University, Kingston, ON Canada, K7L 3N6

22Department of Physics and Astronomy, University of California, Los Angeles,
CA 90095, USA

23Department of Physics, Engineering Physics and Astronomy, Queen’s
University, Kingston, ON K7L 3N6, Canada

24Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam-
Golm, Germany

25Department of Physical Sciences, Munster Technological University,
Bishopstown, Cork, T12 P928, Ireland

26Department of Physics and Astronomy, Purdue University, West Lafayette, IN
47907, USA

27Department of Physics, Indiana University-Purdue University Indianapolis,
Indianapolis, IN 46202, USA

28Department of Physics and Astronomy, Iowa State University, Ames, IA 50011,
USA

29Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027,
USA


# Tailoring accidental double bound states in the continuum in all-dielectric
metasurfaces

###### Abstract

Bound states in the continuum (BICs) have been thoroughly investigated due to
their formally divergent Q-factor, especially those emerging in all-
dielectric, nanostructured metasurfaces from symmetry protection at the
$\Gamma$ point (in-plane wavevector $k_{||}=0$). Less attention has been paid
to accidental BICs that may appear at any other $k_{||}\not=0$ in the band
structure of supported modes, being in turn difficult to predict. Here we make
use of a coupled electric/magnetic dipole model to determine analytical
conditions for the emergence of accidental BICs, valid for any planar array of
meta-atoms that can be described by dipolar resonances, which is the case of
many nanostructures in the optical domain. This is explored for all-dielectric
nanospheres through explicit analytical conditions that allow us in turn to
predict accidental BIC positions in the parameter space
$(\omega,\bf{k_{||}}$). Finally, such conditions are exploited to determine
not only single, but also double (for both linear polarizations) accidental
BICs occurring at the same position in the dispersion relation
$\omega-\bf{k_{||}}$ for realistic semiconductor nanodisk meta-atoms. This
might pave the way to a variety of BIC-enhanced light-matter interaction
phenomena at the nanoscale such as lasing or non-linear conversion, that
benefit from emerging at wavevectors away from the $\Gamma$ point (off-normal
incidence) overlapping for both linear polarizations.

###### keywords:

accidental bound states in the continuum, all-dielectric metasurfaces, Mie-
resonant nanostructures

Diego R. Abujetas* Jorge Olmos-Trigo José A. Sánchez-Gil*

Dr. Diego R. Abujetas
Physics Department, Fribourg University, Chemin de Musée 3, 1700 Fribourg
Switzerland
Email Address<EMAIL_ADDRESS>
Dr. Jorge Olmos-Trigo
Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastian,
Spain

Dr. José A. Sánchez-Gil
Instituto de Estructura de la Materia (IEM-CSIC), Consejo Superior de
Investigaciones Científicas, Serrano 121, 28006 Madrid, Spain
Email Address<EMAIL_ADDRESS>
##  1 Introduction

Planar arrays of particles, often called metasurfaces, are versatile thin
platforms designed to control the properties of light at both the far- and the
near- field, that have opened up new opportunities to tailor light-matter
interaction by customizing the electromagnetic field [1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16]. The shape, position, and material properties
of each particle in the array can be freely designed in order to produce
desired wavefronts (known as Huygens metasurfaces) or can in turn be
periodically ordered to construct periodic metasurfaces, where the particles
that form the repeating unit cell are call meta-atoms. In this context, the
term metasurface is coined for arrays with sub-wavelength lattice constant
where only the non-diffractive orders can propagate to the far field, whereas
metagrating is used to indicate that diffraction can be relevant.

The properties of a system can be fully described in terms of their resonant
modes, where a resonant mode is a self-consistent solution of the
electromagnetic field in the absence of external sources. For infinitely
extended systems, as periodic planar arrays, the nature of the resonant modes
can be divided in two kinds: leaky and confined modes. Typically, the resonant
frequencies are complex inside the continuum of radiation (for the in-plane
component of the wavevector smaller than the free propagating wavevector),
with an imaginary part proportional to the timescale at which the energy is
leaked out to the far field. However, it is also possible to find confined
modes, characterized by real frequencies, inside the continuum of radiation.
Notably, despite there being available radiation channels in which the energy
can decay, these modes cannot couple into them, receiving the name of bound
states in the continuum (BICs). Thus, BICs are resonant states with infinitely
high Q-factors that cannot be excited by far field radiation [17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29]. Nonetheless, under perturbation in the
parameter space around the BIC condition, it is possible (from the far field)
to excite resonances (quasi-BICs) with arbitrary large Q-factors and huge
enhancement of the electromagnetic field at the near field, holding promise
(both BICs and quasi-BICs) of unprecedented planar devices in Nanophotonics.
These interesting properties have been extensively investigated for diverse
photonic applications such as enhanced sensing [30, 31, 32, 33, 34], filtering
[35], lasing [19, 36, 37, 38, 39, 40, 41], electromagnetically-induced
transparency [42], chirality [43, 44], and non-linear conversion [45, 46, 47,
48].

Depending on the mechanism that prevents the coupling of BIC with the
continuum of radiation, they are classified as symmetry-protected or
accidental BICs [20, 29, 49]. Symmetry-protected BICs have been largely
explored and they can only be found at the $\Gamma$ point, where there is a
mismatch between the symmetry of the field mode and the available radiation
channel imposed by the in-phase oscillation of the particles in the array. By
contrast, accidental BICs emerge when different modes are coupled and they
interfere destructively at a specific angle, typically arising in high-order
resonant bands in photonic crystal slabs [19, 34, 36, 41, 49, 50, 51]. Despite
both photonic crystals and metasurfaces can be considered analogous in the
non-diffracting regime as quasi-planar periodic arrangements, there are subtle
theoretical differences in the manner that connected domains (slab with holes
in photonic crystals) polarize light, as compared with isolated domains (meta-
atoms in the case of metasurfaces) that allow for resonant electromagnetic
confinement inside, with non-trivial consequences in the resulting
phenomenology; and Babinet-like principles are not at all evident as long as
holes/meta-atoms are thick and penetrable [52]. In the case of metasurfaces,
since particles in the array suffer from depolarization effects from the own
lattice, the conditions for the formation of accidental BICs are not evident,
and there are no guidelines to design them in the wavevector space; in fact,
it has been assumed that multipolar contributions are needed to yield
accidental BICs in metasurfaces [49], condition that we will demonstrate not
strictly necessary, as in the case of one-dimensional grating [23], with
periodicity only along one in-plane direction, invariant along the other.

In this work, the emergence of accidental BICs in a rectangular metasurface of
dielectric particles is investigated on the basis of a coupled electric and
magnetic dipole (CEMD) formulation [53, 54]. First, generic conditions are
established, demonstrating that, for arrays of axially symmetric particles,
accidental BIC can only arise along the symmetry lines $\Gamma X$ (and $\Gamma
Y$ for rectangular arrays). Analytic expressions are then derived to determine
the condition for accidental BICs along these symmetry lines, showing that the
individual polarizability of the particles must be related to the lattice
depolarization Green function. These conditions only depend on the geometry of
the lattice and meta-atom polarizibilities, enabling us to establish a general
guideline for the formation of accidental BICs. Later, the accidental BIC
condition for rectangular arrays of all-dielectric spheres is studied, giving
examples of systems supporting TE and TM accidental BICs. Interestingly, it is
shown that rectangular arrays can support more than one accidental BIC for the
same metasurface at a given polarization. Finally, a feasible example with a
square array of semiconductor nanodisks is presented, where the aspect ratio
of the disk offers us a new degree of freedom for engineering the
polarizability of the meta-atom. In addition, rectangular metasurfaces also
allow to tailor the emergence of, not only single, but also double (for both
linear polarizations) accidental BICs at the same position in $\omega,k$-space
(any off-normal angle of incidence along $\Gamma X$ or $\Gamma Y$ for a given
frequency). This in turn confirmed through full numerical simulations, in good
agreement with our quasi-analytical CEMD calculations. Although all the
phenomenology is analyzed for rectangular arrays and one meta-atom per unit
cell, the principles used in this study can be extended to metasurfaces with
any geometry as long as the CEMD approach remains valid.

##  2 CEMD formulation

First, let us consider an infinite rectangular array of identical particles,
labeled as $(n,m)$ and placed at

${\bm{r}}_{nm}=x_{n}\hat{x}+y_{n}\hat{y}=na\hat{x}+mb\hat{y},$ (1)

where $a$ and $b$ are the lattice constants along the $x$ and $y$ axis,
respectively. For the homogeneous problem in the absence of external
illumination, each particle in the array is excited by the waves emitted from
the rest of the array. The self-consistent incident field on the $(n,m)=(0,0)$
particle (${\bm{r}}_{00}=\bm{0}$), $\bm{\Psi}_{\text{inc}}(\bm{0})$, is then
given by the solution of

$\displaystyle\bm{\Psi}_{\text{inc}}(\bm{0})=\sum_{nm}{{}^{\prime}}k^{2}{\overleftrightarrow{{\bf
G}}}(-{\bm{r}}_{nm}){\overleftrightarrow{{\bm{\alpha}}}}\bm{\Psi}_{\text{inc}}({\bm{r}}_{nm}),$
(2)

where $\sum_{nm}{{}^{\prime}}$ means that the sum runs for all indices except
for $(n,m)=(0,0)$. ${\overleftrightarrow{{\bf G}}}({\bm{r}})$ and
${\overleftrightarrow{{\bm{\alpha}}}}$ are matrices representing the dyadic
Green function and the dipolar polarizability of the particles, respectively,
and their representations depend on the basis chosen to describe the
electromagnetic fields, $\bm{\Psi}({\bm{r}})$. For the fields, the time
dependence $\exp\left(-i\omega t\right)$ is assumed; $\omega$ is the angular
frequency, related to the modulus of the wavevector through $k=\omega/c$, $c$
being the speed of light. The dyadic Green function is obtained from the
scalar Green function, $g\left({\bm{r}}\right)$, by applying a linear
differential operator, $\mathbf{\mathcal{L}}$, that also depends on the chosen
basis [53]

For periodic arrays the Bloch’s theorem holds,
$\bm{\Psi}_{\text{inc}}({\bm{r}}_{nm})=\bm{\Psi}_{\text{inc}}(\bm{0})\exp({\mathrm{i}}k_{x}na)\exp({\mathrm{i}}k_{y}mb)=\bm{\Psi}_{\text{inc}}(\bm{0})e^{{\mathrm{i}}\phi_{nm}}$,
where $k_{x}$ and $k_{y}$ are the in-plane components of the wavevectors.
Thus, the self-consistent incident field can be written as

$\displaystyle\bm{\Psi}_{\text{inc}}(\bm{0})=k^{2}\left[\sum_{nm}{{}^{\prime}}{\overleftrightarrow{{\bf
G}}}(-{\bm{r}}_{nm})e^{{\mathrm{i}}\phi_{nm}}\right]{\overleftrightarrow{{\bm{\alpha}}}}\bm{\Psi}_{\text{inc}}(\bm{0})\equiv
k^{2}{\overleftrightarrow{{\mathbb{G}}}}_{b}{\overleftrightarrow{{\bm{\alpha}}}}\bm{\Psi}_{\text{inc}}(\bm{0}).$
(3)

We have defined ${\overleftrightarrow{{\mathbb{G}}}}_{b}$, the lattice
depolarization dyadic (or return Green function), as

$\displaystyle{\overleftrightarrow{{\mathbb{G}}}}_{b}(k,k_{x},k_{y})\equiv\sum_{nm}{{}^{\prime}}{\overleftrightarrow{{\bf
G}}}(-{\bm{r}}_{nm},k,k_{x},k_{y})e^{{\mathrm{i}}\phi_{nm}(k,k_{x},k_{y})},$
(4)

and we have explicitly shown the dependence on both the modulus wavevector and
the in-plane components of the wavevector, $\mathbf{k_{||}}=(k_{x},k_{y})$.
${\overleftrightarrow{{\mathbb{G}}}}_{b}$ tells us about the coupling strength
between particles, and is crucial to determine all the lattice properties and
the nature of the modes supported by the metasurface. Finally, we rearrange
Equation (3) as follows:

$\left[{\overleftrightarrow{{\bf
I}}}-k^{2}{\overleftrightarrow{{\mathbb{G}}}}_{b}{\overleftrightarrow{{\bm{\alpha}}}}\right]\bm{\Psi}_{inc}(\bm{0})=\bm{0},$
(5)

where ${\overleftrightarrow{{\bf I}}}$ is the unit dyadic.

In order to find the resonant states of the metasurfaces we need to find a
solution to the homogeneous linear system of equations Equation (5), appearing
only when the determinant is equal to zero:

$\displaystyle\left|{\overleftrightarrow{{\bf
I}}}-k^{2}{\overleftrightarrow{{\mathbb{G}}}}_{b}{\overleftrightarrow{{\bm{\alpha}}}}\right|=\left|\dfrac{1}{k^{2}{\overleftrightarrow{{\bm{\alpha}}}}}-{\overleftrightarrow{{\mathbb{G}}}}_{b}\right|=0.$
(6)

The complex frequencies at which Equation (6) is satisfied are the
eigenfrequencies, denoted by $\nu=\nu^{\prime}-i\nu^{\prime\prime}$, where
$\nu^{\prime}$ and $\nu^{\prime\prime}$ are real numbers that denote the real
and imaginary parts, respectively. In the latter eigenmode equation, it is
more convenient to use the second expression rather than the first one because
the imaginary part is well defined. Since
${\overleftrightarrow{{\mathbb{G}}}}_{b}$ can be expressed in the reciprocal
space [53], Equation (6) can thus be employed to determine the dispersion
relation of resonant modes in metasurfaces. For the sake of simplicity, we
will assume in what follows that the array is embedded in a uniform medium
(vacuum, indeed), so that ${\overleftrightarrow{{\bf G}}}$ in Equation (3) is
the free space Green dyadic. This is equivalent in an experimental
configuration with a metasurface on a substrate to adding an index-matching
layer [27]. Nonetheless, the impact of a substrate could be straightforwardly
incorporated in our CEMD formulation by including the Green dyadic for two
semi-infinite media.

##  3 Generic condition for accidental BICs

In the more generic case (with any restriction for the values of
$\bf{k_{||}}$), the determinant can be only factorized into two terms (for
diagonal polarizabilities), one for each polarization:

$\displaystyle\left|\dfrac{1}{k^{2}{\overleftrightarrow{{\bm{\alpha}}}}}-{\overleftrightarrow{{\mathbb{G}}}}_{b}\right|=\left|\eta^{(TE)}_{xyz}\right|\left|\eta^{(TM)}_{xyz}\right|,$
(7)

with

$\displaystyle\eta_{xyz}^{(TE)}$ $\displaystyle=$
$\displaystyle\left(\dfrac{1}{k^{2}\alpha_{x}^{(e)}}-G_{bxx}\right)\left(\dfrac{1}{k^{2}\alpha_{y}^{(e)}}-G_{byy}\right)\left(\dfrac{1}{k^{2}\alpha_{z}^{(m)}}-G_{bzz}\right)$
$\displaystyle-$
$\displaystyle\dfrac{1}{k^{2}\alpha_{x}^{(e)}}G_{byz}^{2}-\dfrac{1}{k^{2}\alpha_{y}^{(e)}}G_{bzx}^{2}-\dfrac{1}{k^{2}\alpha_{z}^{(m)}}G_{bxy}^{2}+2G_{bxy}G_{byz}G_{bzx},$
$\displaystyle\eta_{xyz}^{(TM)}$ $\displaystyle=$
$\displaystyle\left(\dfrac{1}{k^{2}\alpha_{x}^{(m)}}-G_{bxx}\right)\left(\dfrac{1}{k^{2}\alpha_{y}^{(m)}}-G_{byy}\right)\left(\dfrac{1}{k^{2}\alpha_{z}^{(e)}}-G_{bzz}\right)$
(8) $\displaystyle-$
$\displaystyle\dfrac{1}{k^{2}\alpha_{x}^{(m)}}G_{byz}^{2}-\dfrac{1}{k^{2}\alpha_{y}^{(m)}}G_{bzx}^{2}-\dfrac{1}{k^{2}\alpha_{z}^{(e)}}G_{bxy}^{2}+2G_{bxy}G_{byz}G_{bzx}.$

In these expressions, $G_{bij}$ are the matrix elements of
${\overleftrightarrow{{\mathbb{G}}}}_{b}$, $\alpha_{i}^{(e),(m)}$ are the
diagonal matrix elements of the electric $(e)$ and magnetic $(m)$
polarizability, and $TM$, $TE$ stand for transverse electric and magnetic
modes, respectively. To simplify, the upper index (EM) used in ref. [53]
(indicating that they connect electric and magnetic dipoles), is omitted in
$G_{byz}$ and $G_{bzx}$, understood for physical reasons. Note that $G_{bxy}$
couples either electric only or magnetic only dipoles, while $G_{byz}$ and
$G_{bzx}$ connect electric with magnetic dipoles. Also, for these expressions,
the sign of $G_{byz}$ and $G_{bzx}$ is taken in such a way that:

$\displaystyle\mathrm{Im}\,\left[G_{byz}\right]=-\dfrac{1}{2abk_{z}}\dfrac{k_{x}}{k},\quad\mathrm{Im}\,\left[G_{bzx}\right]=-\dfrac{1}{2abk_{z}}\dfrac{k_{y}}{k},$
(9)

since their signs must be consistent with their definitions.

The resonant surface modes arising from the zeroes of $\eta^{(TE)}_{xyz}$
(respectively, $\eta^{(TM)}_{xyz}$) represent hybrid modes with different
polarization, where the electric (respectively, magnetic) in-plane dipoles are
coupled among them (by $G_{bxy}$) and with the magnetic (respectively,
electric) dipoles along the $z$ axis (by $G_{byz}$ and $G_{bzx}$). The
interference between the different dipolar modes leads to interesting
phenomenology, with the open question yet on whether accidental BICs can be
supported: recall that the imaginary part of the above mentioned resonant
modes must vanish ($\nu^{\prime\prime}=0$) to be considered BICs.

To this end, let us analyze Equations (8) in a general manner. First, note
that it is possible to solve for the zeroes of each of Equations (8) as a
system of two equations for the real and imaginary parts (separately) with
three unknowns: the real parts of the inverse of the polarizabilities along
the three axis. Then the accidental BIC conditions can be solved for either
mode as a function of a parameter (the real part of the inverse of one of the
polarizabilities), leading to plausible solutions (zeroes) of Equations (8)
along any direction in $k$-space by tuning the remaining polarizibilities.
This is true for non axially symmetric meta-atoms (as an ellipsoidal disk or a
regular disk with its axis in the $xy$ plane) with three different
polarizibilities along the main axes.

Nonetheless, if the particles possess axial symmetry around the $z$ axis (as
in a disk), in-plane polarizibilities are equal, $\alpha_{x}=\alpha_{y}$, and
less degrees of freedom are available. In fact, it can be shown that, if both
$k_{x}\not=0$ and $k_{y}\not=0$, no real values for the inverse of the
polarizability can be found for which Equation (7) is fulfilled. Therefore, it
is not possible to design accidental BIC away from the symmetry lines $\Gamma
X$ and $\Gamma Y$ (for rectangular arrays) for axially symmetric arrays of
particles (in the dipolar regime). Moreover, no accidental BIC can be formed
either along the symmetry lines $\Gamma M$ and $\Gamma S$.

By contrast, $G_{bxy}$ and $G_{bzx}$ (or $G_{byz}$) become zero along the
symmetry line $\Gamma X$ (respectively, $\Gamma Y$), and only the interference
between the electric dipole along the $y$ axis (respectively, $x$ axis) and
the magnetic dipole along the $z$ axis is relevant (and vice versa exchanging
electric/magnetic). Therefore, since the degrees of freedom remain the same,
polarizibilities satisfying Equation (7) can be found, so that accidental BICs
along the symmetry line $\Gamma X$ (or $\Gamma Y$) can be supported by
metasurfaces with axially symmetric meta-atoms, as we reveal in the following
sections.

Figure 1: Sketch of BICs supported by metasurfaces with one (dipolar) meta-
atom per unit cell. (a) Symmetry-protected BIC formed by the in-phase
oscillation of dipoles along the $z$ axis. (b) The destructive interference
between in-plane electric (magnetic) dipoles and out-plane magnetic (electric)
dipoles allows for the formation of accidental BICs. The arrows indicate the
directions along which the emission pattern vanishes.

##  4 Accidental BICs along $\Gamma X$

If we thus consider modes propagating along the $x$ axis ($k_{x}\not=0$ and
$k_{y}=0$) and diagonal polarizabilities, Equation (6) can be factorized in
four terms [53]:

$\displaystyle\left|\dfrac{1}{k^{2}{\overleftrightarrow{{\bm{\alpha}}}}}-{\overleftrightarrow{{\mathbb{G}}}}_{b}\right|=\left|\eta^{(TE)}_{x}\right|\times\left|\eta^{(TM)}_{x}\right|\times\left|\eta^{(TE)}_{yz}\right|\times\left|\eta^{(TM)}_{yz}\right|,$
(10)

with

$\displaystyle\eta^{(TE)}_{x}=\dfrac{1}{k^{2}\alpha_{x}^{(e)}}-G_{bxx},\quad\eta^{(TE)}_{yz}=\left(\dfrac{1}{k^{2}\alpha_{y}^{(e)}}-G_{byy}\right)\left(\dfrac{1}{k^{2}\alpha_{z}^{(m)}}-G_{bzz}\right)-G_{byz}^{2},$
$\displaystyle\eta^{(TM)}_{x}=\dfrac{1}{k^{2}\alpha_{x}^{(m)}}-G_{bxx},\quad\eta^{(TM)}_{yz}=\left(\dfrac{1}{k^{2}\alpha_{y}^{(m)}}-G_{byy}\right)\left(\dfrac{1}{k^{2}\alpha_{z}^{(e)}}-G_{bzz}\right)-G_{byz}^{2}.$
(11)

For waves propagating along the $y$ axis ($k_{x}=0$ and $k_{y}\not=0$), the
phenomenology is the same by replacing the index $y$ by $x$ in Equation (11).
Among the different resonant surface modes represented by each term in Eq.
(10), we exclude $\eta^{(TE)}_{x}$ and $\eta^{(TM)}_{x}$ since they have been
shown to yield no BIC whatsoever (all their poles inside the continuum of
radiation have non-negligible widths, $\nu^{\prime\prime}\not=0$).

We thus focus on the hybrid modes, $\eta^{(TE)}_{yz}$ and $\eta^{(TM)}_{yz}$,
emerging from electric (magnetic) dipoles along the $y$ axis coupled to
magnetic (electric) dipoles along the $z$ axis. At the $\Gamma$ point,
$k_{x}=k_{y}=0$, well known symmetry-protected BICs arise [53]. For the sake
of completeness, we explicitly show the conditions, not shown in ref. [53],
yielding symmetry-protected BICs given by the in-phase oscillation of lossless
dipoles along the $z$ axis:

$\displaystyle\begin{array}[]{cc}\mathrm{TE\ symmetry-protected\
BIC}&\quad\mathrm{TM\ symmetry-protected\ BIC}\\\
\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{z}^{(m)}}\right]_{sBIC}=\mathrm{Re}\,\left[G_{bzz}\right],&\quad\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{z}^{(e)}}\right]_{sBIC}=\mathrm{Re}\,\left[G_{bzz}\right].\end{array}$
(14)

Therefore, for a metasurface to support symmetry-protected BICs, the real part
of the inverse of the polarizability of the meta-atoms must be consistent with
the lattice properties.

Searching for accidental BICs at $k_{x}\not=0$ (with $k_{y}=0$), also for
lossless meta-atoms, the condition $\eta_{yz}=0$ yields a system of two
equations (real and imaginary parts) with two unknowns: the real parts of the
inverse of the polarizabilities along the $y$ and $z$ axes. After a
straightforward derivation, the accidental BIC conditions are obtained:

$\displaystyle\begin{array}[]{cc}\mathrm{TE\ accidental\ BIC}&\mathrm{TM\
accidental\ BIC}\\\
\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{y}^{(e)}}\right]_{aBIC}=\mathrm{Re}\,\left[G_{byy}+\dfrac{k}{k_{x}}G_{byz}\right],&\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{y}^{(m)}}\right]_{aBIC}=\mathrm{Re}\,\left[G_{byy}+\dfrac{k}{k_{x}}G_{byz}\right],\\\
\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{z}^{(m)}}\right]_{aBIC}=\mathrm{Re}\,\left[G_{bzz}+\dfrac{k_{x}}{k}G_{byz}\right],&\mathrm{Re}\,\left[\dfrac{1}{k^{2}\alpha_{z}^{(e)}}\right]_{aBIC}=\mathrm{Re}\,\left[G_{bzz}+\dfrac{k_{x}}{k}G_{byz}\right];\end{array}$
(18)

note that two conditions, for $1/\alpha_{y}$ and $1/\alpha_{z}$ separately,
must be satisfied simultaneously.

Equations (14) and (18) give us an interesting perspective about the
differences between symmetry-protected and accidental BICs. Symmetry-protected
BICs only need to fulfil one condition at a specific frequency for the
polarizability (with in-plane modulus wavevector $k_{||}=0$), whereas
accidental BICs need to fulfil two conditions at the same frequency and
longitudinal component of the wavevector. Then, accidental BICs obey more
restrictive conditions that are not always reachable, depending on the precise
cancellation of radiation from different dipolar modes. These characteristics
are illustrated in Figure 1: a symmetry-protected BIC may simply arise from
vertical dipolar resonances, whereas an accidental BIC requires the
interference between different in-plane and out-of-plane dipolar modes to
cancel radiation to the far-field at specific off-normal angles. Quite
importantly, Equation (18) is a guideline to engineer accidental BICs for a
wide variety of all-dielectric (and plasmonic) photonic metasurfaces with a
common (rectangular) lattice symmetry.

##  5 Accidental BICs on all-dielectric nanosphere metasurfaces

Figure 2:  Accidental BIC condition in parameter space $(a/\lambda,\theta)$
for: (a) a square array ($a=b$) in TE polarization at $n_{s}=$ 2.8, 3, 3.5 and
4; (d) a rectangular array ($a/b=$ 0.66) in TM polarization at $n_{s}=$ 1.8,
2.5, 3 and 3.5. Color lines indicate fixed $r/a$ ratios and the shadowed
sections correspond to diffracting spectral regions. Dispersion relation of
modes that evolve to accidental BICs, marked with geometric symbols in (a, d),
for (b) TE polarization and (e) TM polarization. (c, f) Quality factors of the
resonant modes defined as $Q_{f}=\nu^{\prime}/\nu^{\prime\prime}$, revealing
the divergence associated to accidental BIC emergence.

To shed light on the phenomenology resulting from the previous theoretical
formulation, we represent in Figure 2 the accidental BIC position as a
function both of the normalized frequency ($a/\lambda$) and of the angle of
incidence ($\theta$) for rectangular arrays (lattice parameters $a,b$) of
dielectric spheres with different values of their refractive index. The
polarizability of the spheres is calculated from Mie theory, taking into
account only the first two (dipolar) terms of the harmonic expansion. The
color of each line is correlated with the normalized sphere radius $r/a$ that
the metasurface must fulfil in order to support an accidental BIC, where
$r/a=0.5$ represents a metasurface of touching spheres. The results for TE
accidental BICs in square arrays ($a,b$) are analyzed in Figure 2a. It can be
seen that the general tendency is that a higher $r/a$ ratio is needed to
support accidental BICs as the refractive index of the particles becomes
smaller and as the longitudinal component of the wavevector increases. Also,
note that the curves corresponding to $n_{s}=2.8$ and $n_{s}=3$ end before
achieving the condition $k=k_{x}$ (being $k_{x}$ related to $\theta$ through
the expression $k_{x}=k\sin\theta$). Indeed, it can be found that no
accidental BIC can be formed for approximately $n_{s}<2.5$. On the other hand,
as the value of $n_{s}$ increases, accidental BICs become only accessible at
higher values of $k_{x}$, reaching the diffraction limit, as illustrated in
Figure 2a for $n_{s}=4$.

By contrast, we have found no accidental BICs in TM polarization for square
arrays of all-dielectric spheres. For dielectric spheres, the dipolar mode
$\alpha_{y}^{(m)}$ resonates at frequencies lower than that of
$\alpha_{z}^{(e)}$, precluding its formation according to Equations (18).
Alternatively, we study TM accidental BICs in rectangular arrays that indeed
support them, as revealed in Figure 2d for $a/b=0.66$. As in the case of TE
polarization, larger spheres are needed for lower refractive indices.
Interestingly, these rectangular metasurfaces support accidental BICs for
dielectric spheres with relative small values of $n_{s}$; namely, high-
refractive indices are not strictly necessary, which widens the range of
dielectric materials that could be exploited in this respect. For example,
accidental BICs can be engineered around $\theta=20^{\circ}-50^{\circ}$ for
spheres with $n_{s}=1.8$ and meaningful radii. For higher values of $n_{s}$,
accidental BICs arise at higher frequencies, being possible to engineer them
for $n_{s}\approx 5$, but they are not accessible for all values of $\theta$.
Also, it is remarkable that two accidental BICs can be found in parameter
space for the same $n_{s}$ and $r/a$ parameters (i.e., for the same
metasurface). For example, there are accidental BICs at
$(a/\lambda,\theta)\approx(0.68,19)$ and $(a/\lambda,\theta)\approx(0.59,42)$
for $n_{s}=3.5$ and $r/a=0.2$.

By way of example, we show in Figure 2 the dispersion relations and associated
quality factors for modes approaching the accidental BIC condition in four
specific metasurfaces. These magnitudes are calculated by solving Equation (6)
and using the formalism described in Ref. [53]. The dispersion relations
$\nu^{\prime}(k_{x})$ are shown in the parameter space given by angle $\theta$
(recall that $k_{x}=k\sin\theta$) and normalized frequency $a/\lambda$ (with
$\lambda=c/\nu^{\prime}$); the quality factor is defined as
$Q_{f}=\nu^{\prime}/\nu^{\prime\prime}$ and calculated along the mode
dispersion. TE modes in square arrays are shown in Figure 2b,c for two
different metasurfaces: high-refractive index particles $n_{s}=4$ and
$r/a=0.2$ (blue curves); moderate high-refractive index particles $n_{s}=2.8$
and $r/a=0.3$ (red curves). As expected, the quality factor of the modes
diverges at the predicted position in parameter space of the accidental BICs,
around $\theta=43^{\circ}$ and $\theta=18^{\circ}$ for the two cases shown in
Figure 2b,c. The dispersion relations of TM modes in rectangular arrays
($a/b=0.66$) are also represented in Figure 2e,f. The blue curves correspond
to spheres with low index of refraction, $n_{s}=1.8$ and $r/a=0.37$, whereas
red curves are for high refraction index particles, $n_{s}=3.5$ and $r/a=0.2$.
Interestingly, both metasurfaces support two accidental BICs for the same
polarization, and their relative distance in parameter space can be tailored
by the $r/a$ ratio. Incidentally, if $r/a$ is decreased, a threshold value is
found at which both accidental BICs collapse into a single one at the same
angle and frequency, below which no accidental BIC is found (although a narrow
resonance is still present).

Figure 3: Color maps of the reflectance as a function of normalized frequency
and angle of incidence for four particular cases, marked with geometric
symbols in Figure 2b,e, with the following parameters: TE polarization for a
square array with (a) $n_{s}=4$ and $r/a=0.2$; (b) $n_{s}=2.8$ and $r/a=0.3$;
TM polarization for a rectangular array ($a/b=0.66$) with (c) $n_{s}=1.8$ and
$r/a=0.37$; (d) $n_{s}=3.5$ and $r/a=0.2$.

For the sake of illustration, the reflectance is shown in Figure 3 for
different configurations supporting accidental BICs, marked as black symbols
in Figure 2a,e. The reflectance maps have been calculated using the CEMD
formulation [53]. For TE polarization, two square arrays have been considered
in Figure 3a,b. First, an array of spheres with high refractive index
$n_{s}=4$ and radius $r/a=0.2$ is examined in Figure 3a. The expected behavior
in parameter space (formally equivalent to the $\omega,k$) is observed,
associated with a narrow (accessible) quasi-BIC band with high Q-factors
turning into an (inaccessible) BIC, with a vanishing reflectance corresponding
to the BIC at the predicted angle around $\theta=43^{\circ}$. The other
example is shown in Figure 3b for a lower (but moderately high) refractive
index of $n_{s}=2.8$ at $r/a=0.3$, for which the accidental BIC appears around
$\theta=18^{\circ}$ over a broad high reflectance background.

For TM polarization, we have considered the two examples as marked as black
symbols in Figure 2e. First, the interesting case mentioned above for a low
index of refraction ($n_{s}=1.8$) and $r/a=0.37$: as shown in Figure 3c, two
accidental BICs are found in a rectangular array ($a/b=0.66$) around
$\theta\approx 25^{\circ}$ and $38^{\circ}$ . Note that there are very narrow
high-reflectance bands between both BIC conditions, barely visible. Indeed,
the accidental BIC occurring at large angles is zoomed in in the inset of
Figure 3c.

Finally, Figure 3d presents the reflectance for the rectangular array of high
refractive index spheres with $n_{s}=3.5$ and $r/a=0.2$. This metasurface also
supports two accidental BICs with an interesting phenomenology. The accidental
BIC around $\theta\approx 42^{\circ}$ shows the typical vanishing feature,
i.e., there is no background. Conversely, the BIC around $\theta\approx
19^{\circ}$ occurs over a broad reflectance background with a complex pattern
stemming from the combination of a broad resonance; then the quasi-BIC to BIC
transition appears as a narrow dip within the broad band, stemming from the
interference that fades away within the high reflectance background, as a BIC-
induced transparency band [42].

Although all preceding cases have been discussed in normalized frequencies
(which could be easily extrapolated to any spectral regime), we would like to
emphasize that all of them can be associated to realistic nanosphere
metasurfaces supporting accidental BICs in the optical domain, e.g., with
lattice parameters of the order of $a=300$ nm, and dielectric spheres made of
low (oxides) or high (semiconductors) refractive indices. To further stress
the ability of Equations (18) to design all-dielectric metasurfaces with
accidental BICs, we apply it next to a realistic scenario similar to that
already used for symmetry-protected BICs [39].

##  6 Accidental BICs on semiconductor nanodisk metasurfaces

Figure 4: Accidental BIC condition in parameter space for square arrays of
semiconductor nanodisks with refractive index $n_{d}=3.5$, as a function of
the diameter, $D$, and height, $H$. The lattice constants of the array are
$a=b=$ 300 nm: (a) TE polarization; (b) TM polarization. The dashed curves are
calculated by interpolating the polarizabilities for different heights at a
constant diameter, as a guide to the eye. The shadow sections correspond to
diffracting spectral regions.

In order to give a broader view of accidental BICs in all-dielectric
metasurfaces, the phenomenology is also investigated for a more practical
case: a square array of dielectric nanocylinders of diameter, $D$, and height,
$H$. Moreover, the additional degree of freedom, given by the aspect ratio of
the particle, offers a new parameter to engineer the properties of accidental
BICs. Unlike for spherical particles, the relative spectral position of the
resonances can be tuned with the aspect ratio of cylindrical particles. To
this end, the polarizabilities of the disks are numerically calculated through
SCUFF [55, 56], rigorously accounting for the different components [57]; a
refractive index of $n=3.5$ is assumed, as a typical value for dielectric
semiconductors at optical and near-IR frequencies [39]. Specifically, the
polarizabilities of disks with diameters $D=120,150$, $180$, and $210$ nm, and
heights $H=80,110,140,170,$ and $200$ nm (with all different combinations) are
calculated.

The accidental BIC position for an array of semiconductor disks with lattice
constant $a=b=300$ nm (square array) is presented in Figure 4. Therein the
position of the symbols marks the point in the parameter space at which the
metasurface supports an accidental BIC for the different disks calculated
numerically through Equations (18). In addition, the dashed lines indicate the
evolution of the spectral position as a function of the diameter, where the
polarizability is calculated by interpolating those for different heights at
constant diameter. As a general trend, the angular position of the accidental
BIC increases with the diameter in both polarizations, up to a point for
larger disks that the conditions are no longer fulfilled. In addition, it
should be emphasized that accidental BICs can be supported in square arrays of
nanodisks for TM polarization, in contrast to the case of a square array of
nanospheres.

Figure 5: (a) Condition for the emergence of double accidental BICs for a
rectangular array of semiconductor ($n_{d}=3.5$) nanodisks of diameter $D=210$
nm and height $h=200$ nm, as a function of $\theta$ and $b$ for different
values of $a$. The dashed curves correspond to the cases at which the relation
dispersion and the reflectance is calculated. Dispersion relation at both
polarization for (b) $a=450$ nm and $b=845$ nm and (d) $a=550$ nm and $b=985$
nm. The insets show the detuned between the resonances, $\delta$. The
associated quality factors are shown in (c, e).

The appearance of accidental BICs for TM polarization is possible due to the
tunability of the relative spectral position of the dipolar resonances of the
individual nanodisks with different aspect ratios. Moreover, exploiting such
versatility, it is also possible to find configurations at which the
metasurface supports accidental BICs for both polarizations at the same
position in the parameter space. The requirements are: (i) equal electric and
magnetic in-plane polarizabilities, (ii) equal electric and magnetic out-of-
plane polarizabilities; and (iii) fulfilment of Equations (18). We have
checked that a cylinder of $D=210$ nm and $H=200$ nm closely meets (i) and
(ii) around the wavelength $\lambda\sim 910$ nm, where the parameters of the
lattice can be tuned to fulfil the third requirement. Therefore, based on the
properties of this specific disk and fixing the operation wavelength at
$\lambda=910$ nm, the double accidental BIC condition is found as a function
of the lattice parameter along the $y$ axis, $b$, and the in-plane wavevector
along the $x$ axis, $k_{x}$ (where $k_{x}=k\sin\theta$), by changing the
lattice parameter along the $x$ axis, $a$. The results are illustrated in
Figure 5a. As the lattice parameter $a$ becomes smaller, the angle at which
the BIC is found moves to higher angles of incidence, up to the point below
$a\sim 420$ nm such that no accidental BIC emerges (not shown in the graph).
Also, it can be seen that the lattice parameter $b$ must be always bigger than
$a$.

Now let us calculate the dispersion relation of the resonant modes in two
specific configurations, proceeding as in the case of sphere arrays shown in
Figure 2b,e. Note that, for a proper calculation of the dispersion relation,
the polarizability at complex frequencies is needed, trivial for the
analytical polarizabilities of spheres, but not for those of disks.
Nonetheless, since we are interested in the dispersion relation around the
accidental BIC, the imaginary part of the mode frequency, $\nu$”, is
vanishingly small. Therefore, the polarizability at complex frequencies can be
safely approximated to the polarizability at the real part of the complex
frequency. This approximation is very accurate for modes with $Q_{f}>10^{3}$.

First, Figure 5b,c shows the dispersion relations and quality factors for both
polatizations for a metasurface with lattice parameters $a=450$ nm and $b=845$
nm. Both dispersion relations nearly overlap in parameter space. For the shake
of clarity, the inset shows the detuning parameter between resonances,
$\delta$, defined as
$\delta=|\nu^{\prime}_{TM}-\nu^{\prime}_{TE}|/(\nu^{\prime}_{TM}+\nu^{\prime}_{TE})$:
note that $\delta<2\times 10^{-5}$ in the range
$\theta=71^{\circ}-78^{\circ}$. Also, it is remarkable that the quality factor
of both modes diverges around $\theta=72^{\circ}$, as predicted. A similar
phenomenology is observed in a rectangular array with lattice constant $a=550$
nm and $b=985$ nm. As illustrated in Figure 5d,e, both modes lie down again
close to each other in parameter space, although in this configuration the
detuning is higher in the vicinity of the minimum. Nonetheless, around
$\theta=37^{\circ}$, where the quality factor becomes very large, the detuning
is near zero. Significantly, the quality factor of the TM mode exhibits two
divergences, associated to two accidental BICs, a phenomenology that was also
shown above for dielectric sphere arrays in Figure 2f.

By choosing $a=450$ nm, the specular reflectance is represented in Figures
6a,b at the double accidental BIC condition. For both polarizations, the
accidental BIC is found at $\theta\sim 72^{\circ}$ (and at $\lambda\sim 910$
nm) when the lattice constant along the $y$ axis is $b=845$ nm. Moreover,
Figure 6c,d shows the reflectance at both polarizations for a rectangular
array with lattice constants $a=550$ nm and $b=985$ nm. As predicted by Figure
5a, a double accidental BIC is supported at $\theta\sim 37^{\circ}$.

Figure 6: Reflectance color maps for (a,c) TE and (b,d) TM polarization for
rectangular arrays of Si nanodisks ($D=210$ nm and $H=200$ nm) with lattice
parameters: (a,b) $a=450$ nm and $b=845$ nm, supporting a double accidental
BIC at $\theta\sim 72^{\circ}$ and $\lambda\sim 910$ nm; (c,d) $a=550$ nm and
$b=985$ nm, supports a double accidental BIC at $\theta\sim 37^{\circ}$ and
$\lambda\sim 910$ nm. Insets zoomed in the spectral region of the BIC
conditions. Black or white horizontal/vertical dashed lines mark the BIC
angular/spectral position. Blue dashed lines delimit the diffractive region.

Finally, we verify the accuracy of our theoretical model by carrying out
numerical calculations through SCUFF [55] (free software implementation based
on the method of moments), for a specific case corresponding to Figure 6c,d: a
rectangular array (lattice constants $a=550$ nm and $b=985$ nm) of
semiconductor nanodisks of diameter $D=210$ nm and height $H=200$ nm. The
resulting reflectances are shown in Figure 7 confirming that accidental BICs
are observed in both polarizations with similar dispersion relations and Q
factors. Nevertheless, there is a slight spectral/angular shift with respect
to the predictions of our CEMD model, especially in TM polarization: the TE-
BIC emerges at $\theta\sim 37^{\circ}$ and $\lambda\sim 902$ nm, whereas the
TM-BIC appears at $\theta\sim 41^{\circ}$ and $\lambda\sim 930$ nm. This is
somewhat expected since our model only retains the dipolar responses. Still,
the agreement is remarkable by all means (other resonant bands and related
features are reproduced too), bearing in mind in turn that for high resolution
reflectance maps, CEMD calculations take a few hours, about two orders of
magnitude less than typical numerical simulations. Our CEMD can thus predict
in a fast manner the range of parameters wherein double accidental BICs emerge
in real configurations, to be fine tuned later on through full numerical
simulations in a limited range of disk/lattice parameters.

Bear in mind that double accidental BICs can be very useful to enhance light-
matter interactions within all-dielectric nanodisk metasurfaces: this implies
that two resonances extending all over the metasurface with extremely large
(in principle, diverging) Q-factors can be simultaneously excited with mixed
polarized light, with the advantage that the arrays can be designed to obtain
ad-hoc positions in the $(\omega,k)$ space (namely, wavelength and angle
nearly at will).

Figure 7:  Numerically calculated (SCUFF) reflectance color maps for (a) TE
and (b) TM polarization, for the same case as that theoretically calculated
through our CEMD in Figure 6c,d: namely, a rectangular array of Si nanodisks
($D=210$ nm and $H=200$ nm) with lattice parameters $a=550$ nm and $b=985$ nm
that supports a double accidental BIC. White horizontal/vertical dashed lines
mark the BIC angular/spectral position. Blue dashed lines delimit the
diffractive region.

##  7 Conclusions

To summarize, with the help of our coupled electric/magnetic dipole model for
infinite planar arrays, valid for meta-atoms properly described by dipolar
resonances, we have determined analytical conditions for the emergence of
accidental BICs in all-dielectric metasurfaces. This is demonstrated first for
all-dielectric spheres in rectangular arrays through explicit conditions that
incorporate their dipolar responses through analytical Mie polarizabilities.
We have explored the analytical conditions in normalized dimensions that can
be easily extrapolated to any spectral regime as long as the refractive index
of the spheres remains constant within the studied spectral domain; in
particular, realistic nanosphere metasurfaces supporting accidental BICs in
the optical domain can be inferred from them, using either oxides and
semiconductors as sphere material. Moreover, such conditions are also
exploited to design semiconductor metasurfaces supporting accidental BICs in
the optical domain, typically fabricated through standard lithographic means,
such as nanodisk arrays, by properly accounting for the polarizabilities of
the cylinder-shaped meta-atoms. In this regard, thanks to the additional
degree of freedom (aspect ratio) allowed by the cylindrical shape as compared
to spheres, we are able to determine conditions for the emergence of, not only
single, but also double (both linear polarizations) accidental BICs at the
same $\omega$ and $k_{||}$: this in turn implies that two resonances extending
all over the metasurface with extremely large (in principle, diverging)
Q-factors can be simultaneously excited with mixed polarized light, with the
advantage that the arrays can be designed to obtain ad-hoc positions in the
$(\omega,k)$ space (namely, wavelength and angle nearly at will). Such
phenomenology opens a new venue for BIC-enhanced light-matter interaction in
all-dielectric metasurfaces in the optical domain (and throughout the
electromagnetic spectrum), paving the way to a variety of phenomena that may
benefit for the emergence of two accidental BICs away from the $\Gamma$-point
for both polarizations, such as non-linear processes, lasing, chirality, etc.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the
author.

Acknowledgements
Financial support is acknowledged from the Swiss National Science Foundation
through the project 197146 and from the Spanish Ministerio de Ciencia e
Innovación (MICIU/ AEI/ FEDER, UE) through the grants PGC2018-095777-B-C21
(MELODIA) and PID2019-109905GA-C22.

Conflict of Interest
The authors declare no conflict of interest.

Data Availability Statement
The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O’Hara, J. Booth, D. R. Smith,  _IEEE Antennas Propag. Mag._ 2012, _54_ , 2 10.
  * [2] S. B. Glybovski, S. A. Tretyakov, P. A. Belov, Y. S. Kivshar, C. R. Simovski,  _Phys. Rep._ 2016, _634_ 1\.
  * [3] A. M. Shaltout, A. V. Kildishev, V. M. Shalaev,  _J. Opt. Soc. Am. B_ 2016, _33_ , 3 501.
  * [4] G. Li, S. Zhang, T. Zentgraf,  _Nat. Rev. Mater._ 2017, _2_ , 5 1.
  * [5] P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, R. Devlin,  _Optica_ 2017, _4_ , 1 139.
  * [6] D. N. Neshev, I. Aharonovich,  _Light Sci. Appl._ 2018, _7_ , 1 58.
  * [7] P. Qiao, W. Yang, C. J. Chang-Hasnain,  _Adv. Opt. Photonics_ 2018, _10_ , 1 180.
  * [8] S. Baur, S. Sanders, A. Manjavacas,  _ACS Nano_ 2018, _12_ , 2 1618.
  * [9] S. Sun, Q. He, J. Hao, S. Xiao, L. Zhou,  _Adv. Opt. Photonics_ 2019, _11_ , 2 380.
  * [10] A. S. Kupriianov, Y. Xu, A. Sayanskiy, V. Dmitriev, Y. S. Kivshar, V. R. Tuz,  _Phys. Rev. Appl._ 2019, _12_ , 1 014024.
  * [11] R. Paniagua-Domínguez, B. Luk’yanchuk, A. Miroshnichenko, J. A. Sánchez-Gil,  _J. Appl. Phys._ 2019, _126_ , 15 150401.
  * [12] I. Staude, T. Pertsch, Y. S. Kivshar,  _ACS Photonics_ 2019, _6_ , 4 802.
  * [13] B. Sain, C. Meier, T. Zentgraf,  _Adv. Photonics_ 2019, _1_ , 02 1.
  * [14] T. Pertsch, Y. Kivshar,  _MRS Bull._ 2020, _45_ , 3 210.
  * [15] O. Tsilipakos, A. C. Tasolamprou, A. Pitilakis, F. Liu, X. Wang, M. S. Mirmoosa, D. C. Tzarouchis, S. Abadal, H. Taghvaee, C. Liaskos, A. Tsioliaridou, J. Georgiou, A. Cabellos‐Aparicio, E. Alarcón, S. Ioannidis, A. Pitsillides, I. F. Akyildiz, N. V. Kantartzis, E. N. Economou, C. M. Soukoulis, M. Kafesaki, S. Tretyakov,  _Adv. Opt. Mater._ 2020, _8_ , 17 2000783.
  * [16] A. Cuartero-González, S. Sanders, L. Zundel, A. I. Fernández-Domínguez, A. Manjavacas,  _ACS Nano_ 2020, _14_ , 9 11876.
  * [17] D. C. Marinica, A. G. Borisov, S. V. Shabanov,  _Phys. Rev. Lett._ 2008, _100_ , 18 183902.
  * [18] J. Lee, B. Zhen, S.-L. Chua, W. Qiu, J. D. Joannopoulos, M. Soljačić, O. Shapira,  _Phys. Rev. Lett._ 2012, _109_ , 6 67401.
  * [19] C. W. Hsu, B. Zhen, J. Lee, S.-L. Chua, S. G. Johnson, J. D. Joannopoulos, M. Soljačić,  _Nature_ 2013, _499_ , 7457 188.
  * [20] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, M. Soljačić,  _Nat. Rev. Mater._ 2016, _1_ , 9 16048.
  * [21] K. L. Koshelev, S. Lepeshov, M. Liu, A. A. Bogdanov, Y. Kivshar,  _Phys. Rev. Lett._ 2018, _121_ , 19 193903.
  * [22] M. Minkov, I. A. D. Williamson, M. Xiao, S. Fan,  _Phys. Rev. Lett._ 2018, _121_ , 26 263901.
  * [23] H. M. Doeleman, F. Monticone, W. den Hollander, A. Alù, A. F. Koenderink,  _Nat. Photonics_ 2018, _12_ , 7 397.
  * [24] K. L. Koshelev, A. Bogdanov, Y. Kivshar,  _Sci. Bull._ 2019, _64_ , 12 836.
  * [25] D. R. Abujetas, Á. Barreda, F. Moreno, J. J. Sáenz, A. Litman, J.-M. Geffrin, J. A. Sánchez-Gil,  _Sci. Rep._ 2019, _9_ , 1 16048.
  * [26] D. R. Abujetas, N. van Hoof, S. ter Huurne, J. Gómez Rivas, J. A. Sánchez-Gil,  _Optica_ 2019, _6_ , 8 996.
  * [27] S. Murai, D. R. Abujetas, G. W. Castellanos, J. A. Sánchez-Gil, F. Zhang, J. G. Rivas,  _ACS Photonics_ 2020, _7_ , 8 2204.
  * [28] X. Chen, W. Fan,  _Nanomaterials_ 2020, _10_ , 4 623.
  * [29] S. I. Azzam, A. V. Kildishev,  _Adv. Opt. Mater._ 2021, _9_ , 1 2001469.
  * [30] A. A. Yanik, A. E. Cetin, M. Huang, A. Artar, S. H. Mousavi, A. Khanikaev, J. H. Connor, G. Shvets, H. Altug,  _Proc. Natl. Acad. Sci._ 2011, _108_ , 29 11784.
  * [31] F. Yesilkoy, E. R. Arvelo, Y. Jahani, M. Liu, A. Tittl, V. Cevher, Y. Kivshar, H. Altug,  _Nat. Photonics_ 2019, _13_ , 6 390.
  * [32] D. R. Abujetas, J. J. Sáenz, J. A. Sánchez-Gil,  _J. Appl. Phys._ 2019, _125_ , 18 183103.
  * [33] M. L. Tseng, Y. Jahani, A. Leitis, H. Altug,  _ACS Photonics_ 2021, _8_ , 1 47.
  * [34] S. Romano, M. Mangini, E. Penzo, S. Cabrini, A. C. De Luca, I. Rendina, V. Mocella, G. Zito,  _ACS Nano_ 2020, _14_ , 11 15417.
  * [35] J. M. Foley, S. M. Young, J. D. Phillips,  _Phys. Rev. B - Condens. Matter Mater. Phys._ 2014, _89_ , 16 165111.
  * [36] A. Kodigala, T. Lepetit, Q. Gu, B. Bahari, Y. Fainman, B. Kanté,  _Nature_ 2017, _541_ , 7636 196.
  * [37] S. T. Ha, Y. H. Fu, N. K. Emani, Z. Pan, R. M. Bakker, R. Paniagua-Domínguez, A. I. Kuznetsov,  _Nat. Nanotechnol._ 2018, _13_ , 11 1042.
  * [38] E. Khaidarov, Z. Liu, R. Paniagua-Domínguez, S. T. Ha, V. Valuckas, X. Liang, Y. Akimov, P. Bai, C. E. Png, H. V. Demir, A. I. Kuznetsov,  _Laser Photonics Rev._ 2020, _14_ , 1 1900235.
  * [39] S. Murai, G. W. Castellanos, T. V. Raziman, A. G. Curto, J. G. Rivas,  _Adv. Opt. Mater._ 2020, _8_ , 16 1902024.
  * [40] M. Wu, S. T. Ha, S. Shendre, E. G. Durmusoglu, W.-K. Koh, D. R. Abujetas, J. A. Sánchez-Gil, R. Paniagua-Domínguez, H. V. Demir, A. I. Kuznetsov,  _Nano Lett._ 2020, _20_ , 8 6005.
  * [41] M.-S. Hwang, H.-C. Lee, K.-H. Kim, K.-Y. Jeong, S.-H. Kwon, K. Koshelev, Y. Kivshar, H.-G. Park,  _Nat. Commun._ 2021, _12_ , 1 4135.
  * [42] D. R. Abujetas, Á. Barreda, F. Moreno, A. Litman, J. M. Geffrin, J. A. Sánchez-Gil,  _Laser Photonics Rev._ 2021, _15_ , 1.
  * [43] M. V. Gorkunov, A. A. Antonov, V. R. Tuz, A. S. Kupriianov, Y. S. Kivshar,  _Adv. Opt. Mater._ 2021, _2100797_ 2100797\.
  * [44] K. Kim, J. Kim,  _Adv. Opt. Mater._ 2021, _9_ , 22 2101162.
  * [45] L. Carletti, K. L. Koshelev, C. De Angelis, Y. Kivshar,  _Phys. Rev. Lett._ 2018, _121_ , 3 033903.
  * [46] A. P. Anthur, H. Zhang, R. Paniagua-Domínguez, D. A. Kalashnikov, S. T. Ha, T. W. Maß, A. I. Kuznetsov, L. Krivitsky,  _Nano Lett._ 2020.
  * [47] N. Bernhardt, K. L. Koshelev, S. J. White, K. W. C. Meng, J. E. Fröch, S. Kim, T. T. Tran, D.-Y. Choi, Y. Kivshar, A. S. Solntsev,  _Nano Lett._ 2020, _20_ , 7 5309.
  * [48] L. Kang, Y. Wu, X. Ma, S. Lan, D. H. Werner,  _Adv. Opt. Mater._ 2021, _2101497_ 2101497\.
  * [49] Z. F. Sadrieva, K. Frizyuk, M. Petrov, Y. Kivshar, A. A. Bogdanov,  _Phys. Rev. B - Condens. Matter Mater. Phys._ 2019, _100_ , 11 115303.
  * [50] M. Kang, S. Zhang, M. Xiao, H. Xu,  _Phys. Rev. Lett._ 2021, _126_ , 11 117402.
  * [51] S. Han, P. Pitchappa, W. Wang, Y. K. Srivastava, M. V. Rybin, R. Singh,  _Adv. Opt. Mater._ 2021, _2002001_ 2002001\.
  * [52] J. D. Ortiz, J. P. del Risco, J. D. Baena, R. Marqués,  _Appl. Phys. Lett._ 2021, _119_ , 16 1.
  * [53] D. R. Abujetas, J. Olmos-Trigo, J. J. Sáenz, J. A. Sánchez-Gil,  _Phys. Rev. B - Condens. Matter Mater. Phys._ 2020, _102_ , 12 125411.
  * [54] D. R. Abujetas, J. A. Sánchez-Gil,  _Nanomaterials_ 2021, _11_ , 4 998.
  * [55] M. T. H. Reid, S. G. Johnson,  _IEEE Trans. Antennas Propag._ 2015, _63_ , 8 3588.
  * [56] M. H. Reid,  _http://homerreid.github.io/scuff-em-documentation/_ .
  * [57] D. A. Bobylev, D. A. Smirnova, M. A. Gorlach,  _Phys. Rev. B - Condens. Matter Mater. Phys._ 2020, _102_ , 11 1.

Table of Contents

Accidental bound state in the continuum sketch and reflectance map.