# Continuous wave second harmonic generation enabled by quasi-bound-states in
the continuum on gallium phosphide metasurfaces

Aravind P. Anthur  Haizhong Zhang  Ramon Paniagua-Dominguez  [  Dmitry
Kalashnikov  Son Tung Ha  Tobias Wilhelm Wolfgang Mass  Arseniy I. Kuznetsov
<EMAIL_ADDRESS>Leonid Krivitsky  [
<EMAIL_ADDRESS>
###### Abstract

Resonant metasurfaces are an attractive platform for enhancing the non-linear
optical processes, such as second harmonic generation (SHG), since they can
generate very large local electromagnetic fields while relaxing the phase-
matching requirements. Here, we take this platform a step closer to the
practical applications by demonstrating visible range, continuous wave (CW)
SHG. We do so by combining the attractive material properties of gallium
phosphide with engineered, high quality-factor photonic modes enabled by bound
states in the continuum. For the optimum case, we obtain efficiencies around
5e-5 $\%$ W-1 when the system is pumped at 1200 nm wavelength with CW
intensities of 1 kW/cm2. Moreover, we measure external efficiencies as high as
0.1 $\%$ W-1 with pump intensities of only 10 MW/cm2 for pulsed irradiation.
This efficiency is higher than the values previously reported for dielectric
metasurfaces, but achieved here with pump intensities that are two orders of
magnitude lower.

###### keywords:

Second harmonic generation, bound states in the continuum, dielectric
metasurfaces, gallium phosphide.

Equal contribution. [Version submitted to Nano Letters.
Revised version: https://doi.org/10.1021/acs.nanolett.0c03601] Equal
contribution. [Version submitted to Nano Letters.
Revised version: https://doi.org/10.1021/acs.nanolett.0c03601] Equal
contribution. [Version submitted to Nano Letters.
Revised version: https://doi.org/10.1021/acs.nanolett.0c03601] IMRE,
A*STAR]Institute of Materials Research and Engineering, A*STAR (Agency for
Science, Technology and Research) Research Entities, 2 Fusionopolis Way,
$\\#$08-03 Innovis, Singapore 138634. IMRE, A*STAR]Institute of Materials
Research and Engineering, A*STAR (Agency for Science, Technology and Research)
Research Entities, 2 Fusionopolis Way, $\\#$08-03 Innovis, Singapore 138634.
SHG, BIC

Nonlinear optical processes like second harmonic generation (SHG) has
attracted the interests of researchers for more than sixty years 1, 2. SHG is
utilized for a broad range of applications including lasers, quantum optics,
spectroscopy, imaging and pulse-width measurement, to name a few 3, 4, 5, 6,
7. One of the fundamental goals for any practical use of SHG is to achieve
efficient conversion of the fundamental frequency pump to the second harmonic
signal. Several materials and device designs have been proposed and
demonstrated to enhance the efficiency. These, however, come at the price of
several constraints related to the operational wavelength ranges, necessary
pump intensities, phase-matching requirements and/or the size of these
devices. Towards reducing the device dimensions and overcoming the phase
matching constraints, there are ongoing efforts in utilizing metasurfaces,
arrays of nanostructures that exploit resonant enhancement of the pump 8, 9,
10, 11, 12. Using this concept, high efficiency has been achieved, for
example, in plasmonic metasurfaces 13, 14, 15. While leading to very high
local field, plasmonic metasurfaces have the disadvantage of having intrinsic
dissipative losses, which result in relatively low damage thresholds for this
kind of devices. As an alternative, dielectric and semiconductor metasurfaces
have recently emerged as a promising platform 16, 17, 18. They offer lower
intrinsic material losses and higher damage thresholds, as well as,
potentially, high non-linear coefficients. Using materials that lack the
inversion symmetry of the crystal lattice allows SHG in the volume, rather
than at the surface (as is the case of metals). In this regard, semiconductor
metasurfaces with multiple quantum wells are among those with the largest non-
linear coefficients, and have been used to achieve high conversion
efficiencies 19, 20. The main limitation of this approach, however, is that
the SHG signal in these structures saturates at relatively low pump
intensities, limiting the total amount of power that can be converted in the
device. Resonant dielectric and semiconductor metasurfaces using volumetric
SHG, and having higher damage and saturation thresholds, have been used to
demonstrate very high SHG efficiencies utilizing cleverly engineered optical
modes 21, 22, 23, 24, 25, 26.

In this regard, a novel design strategy of dielectric nanoantennae and
metasurfaces based on the concept of bound-states in the continuum (BIC) have
recently attracted interest in the community 26, 27, 28, 29, 30, 31, 32, 33.
The main reason is that it allows to achieve very high quality factor (Q)
modes, which can be easily tuned by varying the device geometry, enabling
highly efficient nonlinear optical processes 34, 35, 26, 36. Despite the
advances in the optical mode and material engineering, all nanoantenna- and
metasurface-based SHG devices demonstrated so far have been realized in the
high-intensity, pulsed regime (typically involving femtosecond lasers), which
hinders their use in many practical applications that work in the low
intensity CW regime. This also creates an additional limitation due to a
limited overlap of the short-pulse laser spectrum with the narrow-band
resonances of high-Q metasurface, which further limits the conversion
efficiency 37. In this work, we combine an emerging material platform for non-
linear applications, namely gallium phosphide (GaP) 38, 39, and a metasurface
design supporting a high-Q quasi-BIC mode that exploits the symmetries of the
non-linear tensor, to generate continuous wave (CW) SHG at low pumping
intensities with reasonable conversion efficiencies. We also show that, in the
pulsed regime, we achieve external efficiencies that are higher than the
results reported in the literature for dielectric metasurfaces to the best of
our knowledge, despite using two orders of magnitude lower pumping
intensities.

Figure 1: (a) Schematic of the metasurface, comprising a square lattice
(periods, $\Lambda$ = 700 nm) of dimers formed by two elliptical cylinders
(with minor and major axes D1 and D2, respectively, height of 150 nm and
center-to-center separation of d = 340 nm) with a relative tilting by an angle
$\theta$ between their major axes. (b) Schematic representation of origin of
the quasi-BIC. A single in-plane electric dipole is induced in each particle
in the dimer, containing x- and y-components ($p_{(x,y)}^{i}$). These, when
summed up, give origin to an out-of-plane magnetic dipole ($m_{z}$) and in-
plane electric quadrupole (Q${}_{xy}^{e}$), and an x-oriented electric dipole
($p_{x}$), whose amplitude can be controlled by $\theta$. Of these, only
$p_{x}$ can couple to normally incident, x-polarized plane waves. (c)
Simulated transmission spectra of x-polarized normally incident plane waves
traveling through the metasurface with $\theta=0^{\circ}$ (red) and
$\theta=5^{\circ}$ (black). (d) Amplitude of the electric near-field
distribution at the quasi-BIC resonance in the xy-plane passing through the
center of the particles. (e) Measured transmission spectra of x-polarized,
normally incident light traveling through the fabricated metasurfaces with
$\theta=0^{\circ}$ (red) and $\theta=5^{\circ}$ (black). (f) SEM image of the
fabricated sample with $\theta=5^{\circ}$ with scale bar of 1 um.

Our metasurface consists of a square array of nano-dimers comprising two GaP
cylinders with elliptical cross section, as shown in Fig. 1(a). They have
period of $\Lambda$=700 nm, minor axis length of D1=200 nm, major axis length
of D2=460 nm, height of 150 nm and center-to-center distance of d=340 nm. When
the major axes of these cylinders are perfectly parallel to each other,the
system supports a perfect BIC at the gamma point for the selected pump
wavelength ($\lambda$ $\approx$ 1225 nm) 40, 41. This BIC stems from the
excitation of an in-plane electric quadrupole mode, Q${}_{xy}^{e}$, and an
out-of-plane magnetic dipole mode, mz, in each dimer, which in turn result
from individual electric dipoles excited in each nanocylider, as shown in Fig.
1(b). Due to the symmetry of the waves radiated from these modes, they cannot
couple to plane waves emerging from the system at normal incidence and the
system becomes decoupled from the radiation continuum. A simple way to
visualize the emergence of this BIC is the following. Each unit cell,
comprising a dimer, radiates according to the sum of the in-plane electric
quadrupole and the out-of-plane magnetic dipole. Since each of these two modes
has vanishing radiation intensity in the direction normal to the metasurface,
their sum obviously also has so. Now, when all the unit-cells oscillate in
phase, a situation corresponding to the gamma point, and when the lattice
period is sub-diffractive, this normal direction is the only direction in
which the system is allowed to radiate. However, when radiation from each unit
cell vanishes in the normal direction, as is the case of these dimers, overall
radiation is forbidden from the system and the modes become bounded to the
system (hence the name BIC). By reciprocity, this also implies that no energy
can be coupled into these modes from normally incident plane waves. In order
to be able to do so, a certain radiation leakage channel is necessary. In this
system, an easy way to open such a channel is by symmetry breaking, realized
via tilting of the ellipses in opposite directions around the z-axis by a
certain angle $\theta$, as shown in Fig. 1(a). In such situation, the system
might support an additional in-plane, x-oriented electric dipole moment, px
(see Fig. 1(b)) that opens the possibility to couple power in and out of the
system, when illuminated by a normally-incident x-polarized plane wave (see
Fig. 1(a)). This effectively transforms the perfect BIC, characterized in the
lossless case for having Q $\to\infty$, into the quasi-BIC, with large but
finite Q. The amplitude of px, which can be controlled by the angle $\theta$,
determines the amount of leakage and, thus, the Q of the mode 41.

We use this quasi-BIC mode to couple power into the system and achieve
enhancement of the pump intensity, thereby achieving highly efficient SHG.
Figure 1(c) shows the simulated transmission spectra for two different
metasurfaces, one supporting a perfect BIC ($\theta=0^{\circ}$) and one
supporting a quasi-BIC ($\theta=5^{\circ}$). As can be seen, a spectrally
narrow resonance at $\sim$ 1225 nm can be observed for the case of
$\theta=5^{\circ}$, corresponding to the excitation of the quasi-BIC. The Q
factor is estimated to be more than 4000, obtained through a fitting of the
curve to the Fano formula 42. As expected, this dip is not present for the
case in which the perfect BIC is formed ($\theta=0^{\circ}$), as the mode
becomes uncoupled from external radiation and, thus, cannot be excited by the
pump. In Fig. S1, we plot the evolution of this mode as a function of
$\theta$, from which the narrowing and amplitude decrease of the resonance as
the angle reduces becomes clearly apparent, as expected for a mode evolving
towards a perfect BIC. Note that this behaviour is not observed for the broad
resonance excited at shorter wavelengths, which has its origin in a simple, in
plane electric dipole ($p_{x}$) that is almost insensitive to the tilting
angle. Figure 1(d) shows the simulated enhancement of the electric near-field
amplitude at the quasi-BIC excitation condition. The result indicates that
inside the structures, field enhancements above 50 can be obtained,
corresponding to intensity enhancement of more than 2500. While,
theoretically, arbitrarily high values of Q factors are possible by simply
reducing $\theta$, in practice they are limited by the size of the array, the
material losses and any fabrication imperfections. In our experiments, the
minimum angle for which we observed the quasi-BIC mode is $\theta=5^{\circ}$,
and is therefore the case presented here. A detailed analysis of the
multipolar moments supported by this system is presented in Fig. S2,
corroborating the interpretation of the quasi-BIC origin given above.

We choose GaP for fabricating our devices because it combines all key features
that make an attractive platform for non-linear processes at the nanoscale.
First, it has a large second-order nonlinear coefficient (d36 $\approx$ 70
pm/V for $\lambda\approx$ 1000 nm) 43. Second, it has a relatively large band
gap that pushes the cut-off wavelength down to $\lambda\approx$ 550 nm,
together with a high refractive index, n $\approx$ 3 (for $\lambda\approx$
1200 nm) 44. And last but not the least, nanofabrication of GaP is well
established in the semiconductor industry 45. Since the GaP crystal has zinc-
blende structure, it is critical to know and utilize the nonlinear tensor
appropriately to achieve high SHG efficiency. For that, we first analyze the
crystal axis orientation of our wafer, which turns out to be tilted with
respect to the surface normal 46. Then we fabricate our metasurfaces with an
orientation that maximizes the SHG. We do so using electron beam lithography
(EBL) followed by reactive ion etching (RIE). Both crystal orientation and
fabrication details are provided in the supporting information. The fabricated
samples are of high quality, as indicated by the good agreement between the
measured transmission spectra (Fig. 1(e)) with the simulated ones (Fig. 1(c)),
as well as the scanning electron microscopy (SEM) images of the sample (see
Fig. 1(f) for the $\theta=5^{\circ}$ case and Fig. S3 for all the cases
measured). Despite the slight blueshift in the measured position of the
transmission dip with respect to the simulated one, both spectra are in good
agreement, obtaining a quality factor of approximately 2000 for the
experimental quasi-BIC structure. For these transmission measurements, the
sample was illuminated using an x-polarized broad band halogen light source
(see supporting information).

The simulated and measured transmission spectra of metasurfaces with different
angle $\theta$, spanning from $5^{\circ}$ to $30^{\circ}$, are given in Figs.
2(a) and 2(b), respectively, with the corresponding SEM images given in Fig.
2(c). As expected, the Q factor experiences a dramatic increase as $\theta$
reduces, while the resonance wavelength redshifts. Simulations and experiments
are in a very good agreement with respect to the relative spectral position of
the resonance for different cases and the Q factors (see Fig. 2(f) for a
direct comparison). The overall slight blueshift ($\sim$20 nm) of the
experimental spectra with respect to the theoretical ones can be attributed to
a slight deviation in the height of the structures. The quasi-BIC resonant
wavelengths are designed to be around 1200 nm so that the SHG wavelength is
above the cut-off wavelength of the GaP. It should be noted that all the
spectra are measured under x-polarized illumination, and that no resonances
are observed for y-polarized light.

Figure 2: (a) Simulated and (b) measured transmission spectra near the quasi-
BIC resonance of x-polarized, normally incident light traveling through
metasurfaces with varying $\theta$, from 30∘ to 5∘. (c) SEM images of the
fabricated samples. (d) Simulated and (e) measured SHG spectra when the
metasurfaces in (a) and (b) are illuminated by a pulsed pump (FWHM around 5
nm) centered at the wavelength giving maximum SHG (solid lines). In (d), we
also show as shaded area the SHG spectra computed assuming an infinitely
narrow pump, whose frequency is scanned in the range shown. (f) Quality factor
(Q) retrieved from simulated and measured transmission spectra through a
fitting to the Fano formula.

We now use the near field enhancement inside the structures enabled by the
quasi-BIC resonances to obtain SHG and study its efficiency as a function of
the angle $\theta$. We first pump the metasurfaces using a pulsed laser from
an optical parametric oscillator (APE) with a pulse duration of around 200 fs
and wavelength that is tunable from 1100 nm to 1400 nm. The light is focused
onto the sample with a plano-convex lens with a focal length of 25.4 mm. The
forward generated second harmonic light is collected using an objective lens
with a numerical aperture (NA) of 0.95 in air. The SHG signal is filtered
using a bandpass filter and analyzed (see supporting information and Fig. S4
for more details on the experimental system).

In Figs. 2(d) and 2(e), we show as solid lines, respectively, the simulated
and measured SHG spectra obtained when the central wavelength of the pulsed
pump is tuned to give the maximum SH signal. In our simulations, we assume
that the pump pulse has a Gaussian shape, with 5 nm full-width at half maximum
(FWHM). In Fig. 2(d), we also show as a shaded area the simulated SHG assuming
a single frequency pump, scanned over the range of wavelengths. All simulated
results are computed in two steps. First, the electric field distribution
(E($\omega_{FF}$)) inside the nanostructures at the fundamental frequency
($\omega_{FF}$) is computed. Then, it is used to calculate the non-linear
polarization density, which serves as the source for a second simulation at
the second harmonic frequency, using the relation
P(2$\omega_{FF}$)=2$\epsilon_{0}\overline{\overline{\chi}}^{(2)}$
E($\omega_{FF}$). Here, $\epsilon_{0}$ is the vacuum permittivity and
$\overline{\overline{\chi}}^{(2)}$ is the appropriate non-linear
susceptibility tensor in the laboratory frame, which takes into account the
crystal axis orientation in the experiment, as explained in the supporting
information. Figure 2(e) is measured using a spectrometer (Ocean Optics,
USB4000) with a spectral resolution of $\sim$1 nm which explains why the SH
spectra of the 10∘ and 5∘ have similar spectral widths. From Figs. 2(d) and
(e) one can see that the wavelength of maximum SHG signal exactly corresponds
to the quasi-BIC condition. Moreover, a noticeable narrowing of the SHG
spectra is observed when the pump is tuned to the resonant wavelength of the
metasurface; its FWHM evolving from that of the pump to being determined by
the resonance itself. In the experimental data, this narrowing saturates when
it hits the resolution limit of our spectrometer.

Figure 3: (a) Measured SHG power as a function of the average pulsed pump
power. For each angle, the pulsed central wavelength is tuned to the
wavelength giving the maximum SHG. Fits to the equation
$P_{SHG}=aP_{pump}^{b}$, are shown as dotted lines, while experimental points
are shown as circles. The exponent values obtained from the fits are 2, 1.97,
1.93, and 2.1 for 30∘, 20∘, 10∘ and 5∘ respectively. The R-squared values are
$\sim$ 1\. (b) Simulated (red solid line) and measured (black circles) SHG
power as a function of the pump polarization (0∘ corresponding to the
polarization parallel to the x-axis, i.e. with the axis of the dimer) for the
metasurface with $\theta=5^{\circ}$. (c) Simulated and (d) measured SHG back-
focal plane Figure 4: images at the maximum of SHG for the metasurface with
$\theta=5^{\circ}$. (e) Simulated and (f) measured SHG spectra when the
central frequency of the pulsed pump is scanned around the quasi-BIC
condition.

Figure 4(a) shows the measured SHG power in the forward direction as a
function of the input pump power for the four cases shown in Fig. 2(e). The
input pump power is varied using a variable neutral density filter at the
input. The SHG power is measured using a power meter for the pulsed pump and
using a lock-in amplifier system with an amplified photo detector (APD) for
the CW pump (see supporting information and Fig. S4 for details). Fitting the
curve to a function with the shape $P_{SHG}=aP_{pump}^{b}$, we found the
nearly quadratic dependence expected for SHG. As can be seen from Fig. 4(a),
the average collected SHG power for $\theta=5^{\circ}$ is around 1.44 $\mu$W
for 37 mW of average input pump power. This corresponds to an external SHG
efficiency of 4e-5 ($P_{SHG}$/$P_{pump}$) and 0.1 $\%$ W-1
($P_{SHG}$/$P_{pump}^{2}$), for a pump intensity of approximately 10 MW/cm2.
This level of efficiency, to the best of our knowledge, is better than
previously reported results in the literature, despite being obtained here
using 100 times lower intensities (see Table S1). It is worth noting that, for
this input power, the SHG power is not saturated. We also note the evident
impact of the narrowing of the quasi-BIC resonance on the SHG as the angle
varies from $\theta=30^{\circ}$ to $\theta=5^{\circ}$. As expected, the
increase in Q factor has an associated increase of the SHG, reaching values
for the $\theta=5^{\circ}$ that is 2.5 times larger than the
$\theta=30^{\circ}$ case. This value is lower than what would be expected if
one simply considers the Q factor increase. This is attributed to the fact
that the width of resonance becomes significantly narrower than the pump at
small angles ($\theta$), and thus only a part of the femtosecond pump spectra
is used for the generation of the SH signal. This finally leads to a limited
enhancement in the SHG for the 5 degree compared to the 30 degree and inspires
us to pursue CW SHG in this system.

Before that, we corroborate the quasi-BIC origin of the observed SHG by
analyzing its dependence on the polarization of the pump. As mentioned before,
the quasi-BIC resonance can only be excited when the normally incident wave
has its electric field polarized along the axis of the dimer (x-axis). The
perpendicular, y-oriented polarization, cannot couple to the electric
quadrupole Q${}_{xy}^{e}$, nor to the out-of-plane magnetic dipole mode, mz,
as it follows from the symmetry considerations. Therefore, we expect
significant SHG only when the incident electric field is oriented along the
x-axis. In the experiment, the polarization of the pump is rotated using a
half wave plate and the power of the SHG at the output is measured without any
polarizer (either at the input or the output of the sample). The measurement
results are plotted as circles in Fig. 4(b), together with the simulation
predictions, represented as solid lines. As expected, the SHG power is highest
for x-polarized pump, while for the fully y-polarized pump it is below the
detection limit of our measurement system, involving amplified photodetector
and lock-in amplifier. For completeness, the simulated and measured back-focal
plane images of the far-field SHG radiation pattern in the forward direction
are given in Figs. 4(c) and 4(d), respectively, for the $\theta=5^{\circ}$
sample. The measured back-focal plane image uses an objective with an NA of
0.95 to collect the SHG signal. As can be seen from Fig. 4(d), five spots are
observed, corresponding to the direct forward emission and the SHG coupling
into the four diffraction orders opened at the SH frequency, 2$\omega_{FF}$.
Two of the five spots are more intense than the other three, in good agreement
with the simulations (Fig. 4(c)). Since the crystal axis is rotated by an
angle of 15∘ with respect to the normal to the wafer, we observe the beam
along the direction of propagation at 0∘ angle, which would otherwise be zero
due to the nonlinear tensor for GaP if the axis of the crystal would be
parallel to the normal.

Finally, we sweep the central frequency of the pulsed pump in a range of
wavelengths around the (quasi-)BIC for all samples, from $\theta=30^{\circ}$
to $\theta=0^{\circ}$, and plot the simulated and measured normalized SHG
intensity spectra in Figs. 4(e) and 4(f), respectively. As can be seen,
maximum SHG is obtained when the central wavelength corresponds to the
resonant dip observed in the transmission spectra of Figs. 2(a) and 2(b).
Also, as expected, the FWHM of the curves gets narrower for decreasing angle,
but becomes limited by the convolution with pump pulse width. In fact, as
mentioned before, the mismatch between the (narrower) spectral width of the
quasi-BIC resonance for small angles and the (broader) pulsed pump implies
that part of the incident power is not coupled to the mode and, therefore, to
the SHG process (the quasi-BIC acting effectively as a filter). Thus, we move
on to study the samples under narrow-band CW pump instead, for which the whole
input power can be efficiently used in the SHG.

In Fig. 5(a), we show the optical spectra of the maximum SHG obtained when the
different samples are excited with a tunable CW laser (see supporting
information for details). As an inset in that figure, we show the spectra of
the pump used to obtain the SHG. We note that maximum SHG is observed when the
CW pump is tuned to the corresponding resonant transmission dip of each of the
metasurface arrays. Notably, we observed and measured SHG using $\sim$ kW/cm2
levels of intensity in all the arrays, with angles ranging from
$\theta=30^{\circ}$ to $\theta=5^{\circ}$. As before, SHG is not observed for
the case of parallel ellipses, i.e. for $\theta=0^{\circ}$.

The variation of the SHG power as a function of the input power of the CW pump
is shown in Fig. 5(b). As expected, one can observe an increase in the SHG
power (and thus the efficiency) as $\theta$ decreases from $30^{\circ}$ to
$5^{\circ}$. In our experiments, the maximum SHG power is observed for the
$\theta=5^{\circ}$ case, reaching values around 70 nW for 360 mW of input pump
power. This translates into an external SHG efficiency of 2e-7
($P_{SHG}$/$P_{pump}$) and 5e-5 $\%$ W-1 ($P_{SHG}$/$P_{pump}^{2}$), for a
pump intensity as low as 1 kW/cm2. As compared to the pulsed case, since the
CW linewidth is less than 1 MHz, the narrowing of the resonance (or increase
of the Q factor) does not cut any incident power, and thus the enhancement is
not limited when $\theta$ decreases. Therefore, in Fig. 5, the impact of the
increase in Q factor is more clearly observed in the SHG, which is almost two
orders of magnitude larger for $\theta$=5∘ than it is for $\theta=30^{\circ}$,
compared to only 2.5 times increase in the pulsed case.

Figure 5: (a) Spectra of the SHG when the CW pump is tuned to the resonance of
30∘, 20∘, 10∘ and 5∘ samples, and the corresponding pump spectra given in the
inset. (b) The SHG power as a function of input CW pump power for various
angles of the dimer. Fits to the equation $P_{SHG}=aP_{pump}^{b}$, are shown
as dotted lines, while experimental points are shown as circles. The exponent
values obtained from the fits are 1.9, 2.3, 1.91, and 2.1 for 30∘, 20∘, 10∘
and 5∘ respectively. The R-squared values for the fit are $\sim$1.

The generation of the SH using CW and pulsed pump also shows that the
metasurfaces utilizing quasi-BIC based resonance enhancement of the pump
fields can be used for applications with different requirements by tuning the
symmetry-breaking of these structures thereby achieving precise tuning of the
Q factors. For example, for pulse width measurements in auto-correlators, the
spectral width of the resonance of these metasurfaces need to be larger than
the spectral width of the laser. For applications where CW pumps are used, the
asymmetry can be tuned for highest Q factor and conversion efficiency. Thus
these structures have the flexibility to be utilized in various applications
by controlling the device geometry.

## Conclusion

In conclusion, we have demonstrated high efficiency second harmonic generation
(SHG) from GaP metasurfaces with sub-wavelength thicknesses ($\sim\lambda$/8).
This is achieved by combining the emerging material platform of GaP, offering
high non-linearity, high refractive index, low loss and ease of fabrication,
with high-Q optical resonances enabled by quasi bound-states in the continuum
(BIC). We achieve these quasi-BIC resonances by breaking the metasurface unit-
cell symmetry, that controllably opens a leaky channel in the normal
direction. Opening of leaky channel allows the coupling of radiation in-and-
out of the otherwise perfect BIC, formed by in-plane electric quadrupole and
out-of-plane magnetic dipole modes. We observe a dramatic increase of the Q
factor of the resonance as the asymmetry parameter decreases, which translates
into increased SHG efficiency, $P_{SHG}$/$P_{pump}^{2}$
($P_{SHG}$/$P_{pump}$), reaching values of 0.1 $\%$ W-1 (4e-5) for pulsed
pumps with intensities of approximately 10 MW/cm2. This level of external
efficiency is obtained without reaching saturation and is higher than the
values previously reported in the literature, to the best of our knowledge,
but achieved here using 100 times lower intensities. We also demonstrated SHG
in these metasurfaces using a continuous wave (CW) pump. We generated CW SHG
using pump with intensities of 1 kW/cm2, and achieved conversion efficiencies
of 5e-5 $\%$ W-1 (2e-7) in the metasurface with the highest Q factor. We
believe that our work shows that the combination of suitable material
platforms and properly engineered optical resonances might push SHG in
metasurfaces to levels closer to those required for practical applications,
including those that demand continuous wave, second harmonic generation.

We acknowledge the support of the Quantum Technologies for Engineering (QTE)
program of A*STAR, IET A F Harvey Engineering Research Prize 2016 and A*STAR
SERC Pharos program, Grant No. 152 73 00025 (Singapore).

### Author contributions

R. P.-D. proposed and simulated the structures. H. Z. Z. fabricated the
samples. A. P. A. constructed and performed the optical characterization
measurements. S. T. H. performed the transmission measurements. T. W. W. M.
simulated the back-focal plane images. D. K. helped in the optical
characterization. A. P. A wrote the first draft of the manuscript, with inputs
from R. P.-D. and H. Z. Z. A. I. K. and L. K. conceived the idea and
supervised the project. All authors contributed to the final version of the
manuscript

### Competing interests

The authors declare no competing financial interest.

The following files are available free of charge.

  * •

Supplementary material contain the details of simulation, fabrication and
experimental measurement system. Supplementary material also contains the
following figures and table. Colour map representing the metasurface
transmission as a function of wavelength and angle of tiling ($\theta$) of the
elliptical cylinders forming the dimers in the metasurface. Multipolar
decomposition of one unit-cell of the metasurface supporting a quasi-BIC with
$\theta=5^{\circ}$. Scanning electron microscope images of the different
fabricated metasurfaces. Schematic of the experimental setup. Comparison of
the state-of-the-art performances of SHG in dielectric metasurfaces.

## Supplementary material: continuous wave second harmonic generation enabled
by quasi-bound-states in the continuum on gallium phosphide metasurfaces

### Simulations

Simulations of SHG were carried out using the Finite Element Method in a
commercial software package (COMSOL Multiphysics). The simulations consist of
two steps. In the first one, the transmission spectra and the near-fields
inside the nanostructures comprising the metasurface, E($\omega_{FF}$), are
computed at the fundamental frequency, $\omega_{FF}$. For that, we simulate
the infinite periodic system by using Periodic Boundary Conditions in the x-
and y-directions (transverse to the incident wave propagation). Periodic ports
are used at the top and bottom to excite the system and record the transmitted
and reflected fields. The GaP material parameters are taken from the software
database, which in turn are taken from the literature 44. A homogeneous
surrounding environment with refractive index 1.46 is taken. In the second
step, the near-fields recorded in the first one are used to compute the non-
linear dipole moment density, which is used as the excitation source in the
simulations. This non-linear dipole moment density is computed through the
expression:

$\mbox{{P}}(2\omega_{FF})=2\epsilon_{0}\overline{\overline{\chi}}^{(2)}\mbox{{E}}(\omega_{FF}),$
(1)

where $\epsilon_{0}$ is the vacuum permittivity and
$\overline{\overline{\chi}}^{(2)}$ is the non-linear susceptibility tensor of
the crystal in the laboratory frame. To obtain the correct results, one should
note that, in our GaP wafers, the material is grown such that the crystal axis
subtends an angle of 15∘ with respect to the normal to the wafer surface.
Thus, we first need to express the near-fields, as obtained in the first
simulation, in the crystal frame, in which the non-linear susceptibility
tensor has only off-diagonal components and they are all equal (d14 = d25 =
d36 $\approx$ 70 pm/V) 43. This can be achieved simply using the appropriate
rotation matrix. Once the fields in the crystal frame are known, we compute
the non-linear dipole moment in this frame using expression in Eq. 1. This is
brought back again to the laboratory frame by the inverse rotation matrix, to
perform the simulations. In the simulation, we use again Periodic Boundary
Conditions in the transverse directions but substitute the Periodic Ports by
Perfectly Matched Layers in the z-directions (top and bottom of the
simulation). For the computation of the far-field pattern given in Fig. 3(c),
a homemade implementation of Stratton-Chu formulas is used 47. For that, the
complex electromagnetic fields in one unit-cell are recorded in a plane at
distance of 500 nm from the structures. These are used to generate a finite
array of 50 x 50 periods, which is then rigorously propagated to the far-
field.

### Fabrication

The GaP active layer ($\sim$400 nm) is first grown on a gallium arsenide
(GaAs) substrate with an AlGaInP buffer layer by metal-organic chemical vapor
deposition (MOCVD) to reduce the lattice mismatch between the GaAs substrate
and GaP layer. Then this structure is directly bonded to the sapphire
substrate after depositing $\sim$SiO2 layer on top of both the surfaces. The
AlGaInP/GaAs substrate is then removed by wet etching. The fabrication of the
GaP nanostructures start with a standard wafer cleaning procedure (acetone,
iso-propyl alcohol and deionized water in that sequence under sonication),
followed by O2 and hexamethyl disilizane (HMDS) priming in order to increase
the adhesion between GaP and the subsequent spin-coated electron-beam
lithography (EBL) resist of hydrogen silsesquioxane (HSQ). After spin-coating
of HSQ layer with a thickness of $\sim$ 540 nm, EBL and development in 25 $\%$
tetra-methyl ammonium hydroxide (TMAH) are carried out to define the
nanostructures in HSQ resist. Inductively-coupled plasma reactive ion etching
(ICP-RIE) with N2 and Cl2 is then used to transfer the HSQ patterns to the GaP
layer. Finally, $\sim$3.2 $\mu$ m SiO2 cladding layer is deposited on top of
the structures by ICP-CVD. For optical measurements, an additional layer of
polydimethylsiloxane (PDMS) is placed on top of the cladding layer to avoid
Fabry-Perot resonances. The structures fabricated are 100 $\mu$m $\times$ 100
$\mu$m in size.

### Experimental setup

The schematic of the experimental setup is given in the Fig. S4. The laser
light from an optical parameteric oscillator (OPO) is passed through a
dichroic beam splitter cutting-off shorter wavelength components. The pulsed
pump is obtained from a Coherent APE OPO tunable from 1100 nm to 1300 nm,
pumped by a Ti:Sapphire laser at a wavelength of 830 nm. The CW OPO is a
Hubner CWave laser pumped by a 532 nm laser, that is tunable from 900 nm to
1330 nm. The pump light then passes through a variable optical neutral density
filter, quarter-wave plate (QWP) and half-wave plate (HWP). The variable
neutral density filter is used for varying the pump power levels. The QWP and
HWP are jointly used for varying the pump polarization at the input. A
polarization extinction of approximately 100 is achieved for the pump light at
the input of the sample, as given in Fig. S4(b). A polarizer is used for the
power measurements and polarization analysis of the input pump. For Fig. 3(b),
no polarizer is used at the input. All the other measurement results given in
Figs. 1, 2, 3 and 4 use a polarizer at the input of the sample. Lens, L1, is a
plano-convex lens with a focal length of 25.4 mm, to achieve a beam waist
diameter of approximately 100 um for the pulsed and CW pump. The beam width
measurement is carried out using knife edge and the results are given in Fig.
S4(c). The lens, L2, is an objective with an NA of 0.95. At the output, a
flipping mirror is used to image the far-field beam pattern and the sample. A
4f imaging system (L3-L6) with a charged coupled device (CCD) camera is used
to image the back-focal plane of the objective, to obtain Fig. 3(d). After
this imaging system, there is a filter centered at 650 nm and having a
spectral width of 150 nm (Semrock). The output of the filter is directed
through a chopper to an amplified photo-detector (APD) to measure the SHG
power for CW pump and a power meter (PM) calibrated to the APD is used for the
pulsed pump. The $P_{SHG}$ is the average measured SHG power after the filter
and $P_{pump}$ is the average measured pump power before L1. The transmission
spectrum in Fig. 2(b) is obtained using the same experimental system pumped by
a broad band super-continuum source (Leukos) and coupling the light output
through the sample, immediately after L2 using a flip mirror, to a high
resolution optical spectrum analyzer (OSA, Yokogawa) through a single mode
fiber, and taking the ratio of the transmitted signal through the sample with
and without the nanostructured array. This fiber coupling system is also used
for capturing the spectra of the SHG generated from the metasurface using the
spectrometer (Ocean Optics, USB4000).

Figure S1: Colour map representing the metasurface transmission as a function
of wavelength and angle of tilting ($\theta$) of the elliptical cylinders
forming the dimers in the metasurface, as shown in Figure 1 in the main text.
As seen, the resonance corresponding to the excitation of the quasi-BIC mode
experiences a strong narrowing of the linewidth, as well as a decrease in the
amplitude, when $\theta$ decreases and the mode evolves towards the perfect
BIC. The broad resonance observed at shorter wavelengths, corresponds to the
excitation of an in-plane, electric dipole directed along the incident
polarization ($p_{x}$), and is almost insensitive to the tilting angle. Figure
S2: Multipolar decomposition of one unit-cell of the metasurface supporting a
quasi-BIC with $\theta=5^{\circ}$. (a) Non-negligible multipole moments when
applying the multipole decomposition, only, to the fields of one of the
cylinders comprising the dimer. The origin of coordinates is taken at the
center of the cylinder. As can be seen, the only significant multipole moment
is a tilted electric dipole moment (as indicated by the presence by both x-
and y- components). (b) Non-negligible multipole moments when applying the
multipole decomposition to the whole dimer. The origin of coordinates is taken
at the center of the dimer. As can be seen, there are three significant
multipole moments: an out-of-plane magnetic dipole moment (mz), an in-plane
electric quadrupole moment (Q${}_{xy}^{e}$=Q${}_{yx}^{e}$) and, with much
lower amplitude, an in-plane electric dipole (with only x-component, double of
that of the single cylinder in (a)). The origin of these multipolar components
is detailed in the main text, in connection with Fig. 1(b). Figure S3:
Scanning electron microscope images of the different fabricated metasurfaces
studied in this work, comprising elliptical cylinder dimers with relative
tilting: (a) $\theta=0^{\circ}$, (b) $\theta=2^{\circ}$ (for which we were not
able to observe the quasi-BIC), (c) $\theta=5^{\circ}$, (d)
$\theta=10^{\circ}$, (e) $\theta=20^{\circ}$ and (f) $\theta=30^{\circ}$. The
scale bar is the same for all the images. Figure S4: (a) Schematic of the
experimental setup used to capture the back focal plane image of the
objective, measure the SHG power levels as a function of pump power and pump
polarization, and collect the transmission characteristics of the devices. (b)
Measured pump power at the output of the polarizer as a function of the HWP
rotation angle before L1. (c) Beam radius with respect to the sample position
measured using a knife edge.

Table S1: Comparison of performances of SHG in dielectric metasurfaces. Efficiency is defined as $P_{SHG}/P_{pump}$, where the power levels are average values. Platform |  Pump intensity (GW/cm2) |  Pump wavelength (nm) |  Efficiency |  Reference
---|---|---|---|---
GaP |  10e-3 |  1200 |  4e-5 |  This work (Pulse)
GaP |  1e-6 |  1200 |  2e-7 |  This work (CW)
AlGaAs |  0.22 |  1570 |  1.2e-5 |  26
AlGaAs |  7 |  1550 |  8e-6 |  48
AlGaAs |  1.72 |  1556 |  9e-6 |  22, 49
GaP |  100 |  910 - 1300 |  2e-6 |  24
GaAs |  0.32e-3 |  1000 |  6e-6 |  50
GaAs |  3.4 |  1020 |  2e-5 |  51
ZnO |  15 |  394 |  1e-8 |  52

## References

  * Boyd 2019 Boyd, R. W. _Nonlinear optics_ , 4th ed.; Academic Press, 2019
  * Rottwitt and Tidemand-Lichtenberg 2014 Rottwitt, K.; Tidemand-Lichtenberg, P. _Nonlinear optics: Principles and Applications_ , 1st ed.; CRC Press, 2014
  * Armstrong 1967 Armstrong, J. A. Measurement of picosecond laser pulse widths. _Appl. Phys. Lett._ 1967, _10_ , 16–18
  * Heinz et al. 1982 Heinz, T. F.; Chen, C. K.; Ricard, D.; Shen, Y. R. Spectroscopy of Molecular Monolayers by Resonant Second-Harmonic Generation. _Phys. Rev. Lett._ 1982, _48_ , 478–481
  * Galli et al. 2019 Galli, M. et al. Generation of deep ultraviolet sub-2-fs pulses. _Opt. Lett._ 2019, _44_ , 1308–1311
  * Stolzenburg et al. 2010 Stolzenburg, C.; Schüle, W.; Zawischa, I.; Killi, A.; Sutter, D. 700W intracavity-frequency doubled Yb:YAG thin-disk laser at 100 kHz repetition rate. _Proc. SPIE 7578, Solid State Lasers XIX: Technology and Devices_ 2010, _75780A_
  * Freund and Deutsch 1986 Freund, I.; Deutsch, M. Second-harmonic microscopy of biological tissue. _Opt. Lett._ 1986, _11_ , 94–96
  * Kildishev et al. 2013 Kildishev, A. V.; Boltasseva, A.; Shalaev, V. M. Planar Photonics with Metasurfaces. _Science_ 2013, _339_ , 1232009
  * Lapine et al. 2014 Lapine, M.; Shadrivov, I. V.; Kivshar, Y. S. Colloquium: Nonlinear metamaterials. _Rev. Mod. Phys._ 2014, _86_ , 1093–1123
  * Genevet et al. 2017 Genevet, P.; Capasso, F.; Aieta, F.; Khorasaninejad, M.; Devlin, R. Recent advances in planar optics: from plasmonic to dielectric metasurfaces. _Optica_ 2017, _4_ , 139–152
  * Krasnok et al. 2017 Krasnok, A.; Tymchenko, M.; Alu, A. Nonlinear metasurfaces: a paradigm shift in nonlinear optics. _Mater. Today Commun._ 2017, _21_ , 8–21
  * Staude et al. 2019 Staude, I.; Pertsch, T.; Kivshar, Y. S. All-Dielectric Resonant Meta-Optics Lightens up. _ACS Photonics_ 2019, _6_ , 802–814
  * Ceglia et al. 2015 Ceglia, D. d.; Vincenti, M. A.; De Angelis, C.; Locatelli, A.; Haus, J. W.; Scalora, M. Role of antenna modes and field enhancement in second harmonic generation from dipole nanoantennas. _Opt. Express_ 2015, _23_ , 1715–1729
  * Celebrano et al. 2015 Celebrano, M.; Wu, X.; Baselli, M.; Großmann, S.; Biagioni, P.; Locatelli, A.; De Angelis, C.; Cerullo, G.; Osellame, R.; Hecht, B.; Duò, L.; Ciccacci, F.; Finazzi, M. Mode matching in multiresonant plasmonic nanoantennas for enhanced second harmonic generation. _Nat. Nanotechnol._ 2015, _10_ , 412–417
  * Canfield et al. 2007 Canfield, B. K.; Husu, H.; Laukkanen, J.; Bai, B.; Kuittinen, M.; Turunen, J.; Kauranen, M. Local Field Asymmetry Drives Second-Harmonic Generation in Noncentrosymmetric Nanodimers. _Nano Lett._ 2007, _7_ , 1251–1255
  * Kuznetsov et al. 2016 Kuznetsov, A. I.; Miroshnichenko, A. E.; Brongersma, M. L.; Kivshar, Y. S.; Luk’yanchuk, B. Optically resonant dielectric nanostructures. _Science_ 2016, _354_ , aag2472
  * Kruk and Kivshar 2017 Kruk, S.; Kivshar, Y. Functional Meta-Optics and Nanophotonics Govern by Mie Resonances. _ACS Photonics_ 2017, _4_ , 2638–2649
  * Paniagua-Dominguez et al. 2019 Paniagua-Dominguez, R.; Ha, S. T.; Kuznetsov, A. I. Active and Tunable Nanophotonics With Dielectric Nanoantennas. _Proc. IEEE_ 2019, _108_ , 749–771
  * Lee et al. 2016 Lee, J.; Nookala, N.; Gomez-Diaz, J. S.; Tymchenko, M.; Demmerle, F.; Boehm, G.; Amann, M. C.; Alù, A.; Belkin, M. A. Ultrathin Second-Harmonic Metasurfaces with Record-High Nonlinear Optical Response. _Adv. Opt. Mater._ 2016, _4_ , 664–670
  * Sarma et al. 2019 Sarma, R.; De Ceglia, D.; Nookala, N.; Vincenti, M. A.; Campione, S.; Wolf, O.; Scalora, M.; Sinclair, M. B.; Belkin, M. A.; Brener, I. Broadband and Efficient Second-Harmonic Generation from a Hybrid Dielectric Metasurface/Semiconductor Quantum-Well Structure. _ACS Photonics_ 2019, _6_ , 1458–1465
  * Carletti et al. 2015 Carletti, L.; Locatelli, A.; Stepanenko, O.; Leo, G.; De Angelis, C. Enhanced second-harmonic generation from magnetic resonance in AlGaAs nanoantennas. _Opt. Express_ 2015, _23_ , 26544–26550
  * Gili et al. 2016 Gili, V. F.; Carletti, L.; Locatelli, A.; Rocco, D.; Finazzi, M.; Ghirardini, L.; Favero, I.; Gomez, C.; Lemaître, A.; Celebrano, M.; De Angelis, C.; Leo, G. Monolithic AlGaAs second-harmonic nanoantennas. _Opt. Express_ 2016, _24_ , 15965–15971
  * Camacho-Morales et al. 2016 Camacho-Morales, R. et al. Nonlinear Generation of Vector Beams from AlGaAs Nanoantennas. _Nano Lett._ 2016, _16_ , 7191–7197
  * Cambiasso et al. 2017 Cambiasso, J.; Grinblat, G.; Li, Y.; Rakovich, A.; Cortés, E.; Maier, S. A. Bridging the Gap between Dielectric Nanophotonics and the Visible Regime with Effectively Lossless Gallium Phosphide Antennas. _Nano Lett._ 2017, _17_ , 1219–1225
  * Li et al. 2020 Li, Y.; Huang, Z.; Sui, Z.; Chen, H.; Zhang, X.; Guan, H.; Qiu, W.; Dong, J.; Zhu, W.; Yu, J.; Lu, H.; Chen, Z. Optical anapole mode in nanostructured lithium niobate for enhancing second harmonic generation. _Arxiv_ 2020, _2003.09793_
  * Koshelev et al. 2020 Koshelev, K.; Kruk, S.; Melik-Gaykazyan, E.; Choi, J. H.; Bogdanov, A.; Park, H. G.; Kivshar, Y. Subwavelength dielectric resonators for nonlinear nanophotonics. _Science_ 2020, _367_ , 288–292
  * Hsu et al. 2016 Hsu, C. W.; Zhen, B.; Stone, A. D.; Joannopoulos, J. D.; Soljacic, M. Bound states in the continuum. _Nat. Rev. Mater._ 2016, _1_ , 1–13
  * Kodigala et al. 2017 Kodigala, A.; Lepetit, T.; Gu, Q.; Bahari, B.; Fainman, Y.; Kanté, B. Lasing action from photonic bound states in continuum. _Nature_ 2017, _541_ , 196–199
  * Ha et al. 2018 Ha, S. T.; Fu, Y. H.; Emani, N. K.; Pan, Z.; Bakker, R. M.; Paniagua-Domínguez, R.; Kuznetsov, A. I. Directional lasing in resonant semiconductor nanoantenna arrays. _Nature Nanotechnol._ 2018, _13_ , 1042–1047
  * Kupriianov et al. 2019 Kupriianov, A. S.; Xu, Y.; Sayanskiy, A.; Dmitriev, V.; Kivshar, Y. S.; Tuz, V. R. Metasurface Engineering through Bound States in the Continuum. _Phys. Rev. Appl._ 2019, _12_ , 1–8
  * Abujetas et al. 2018 Abujetas, D. R.; Barreda, A.; Moreno, F.; Sáenz, J. J.; Litman, A.; Geffrin, J.-M.; Sánchez-Gil, J. A. Brewster quasi bound states in the continuum in all-dielectric metasurfaces from single magnetic-dipole resonance meta-atoms. _Sci. Rep._ 2018, _9_ , 16048
  * Wu et al. 2020 Wu, M.; Ha, S. T.; Shendre, S.; Durmusoglu, E. G.; Koh, W.-K.; Abujetas, D. R.; Sánchez-Gil, J. A.; Paniagua-Domínguez, R.; Demir, H. V.; Kuznetsov, A. I. Room-Temperature Lasing in Colloidal Nanoplatelets via Mie-Resonant Bound States in the Continuum. _Nano Lett._ 2020, _20_ , 6005–6011
  * Milnikov et al. 2020 Milnikov, V.; Ha, S. T.; Pan, Z.; Valuckas, V.; Paniagua-Domínguez, R.; Demir, H. V.; Kuznetsov, A. I. Lasing action in single subwavelength particles supporting supercavity modes. _ACS Nano_ 2020, _14_ , 7338–7346
  * Xu et al. 2019 Xu, L.; Kamali, K. Z.; Huang, L.; Rahmani, M.; Smirnov, A.; Camacho-Morales, R.; Ma, Y.; Zhang, G.; Woolley, M.; Neshev, D.; Miroshnichenko, A. E. Dynamic nonlinear image tuning through magnetic dipole quasi-BIC ultrathin resonators. _Adv. Sci_ 2019, _6_ , 1802119
  * Zhou et al. 2020 Zhou, C.; Li, S.; Gong, C.; Wang, Y.; Liu, X. Resonant Asymmetric All-Dielectric Metasurface for Boosting Third-Harmonic Generation. _Arxiv_ 2020, _2004.01088_
  * Bernhardt et al. 2020 Bernhardt, N.; Koshelev, K.; White, S. J. U.; Wong Choon Meng, K.; E. Fröch, J.; Kim, S.; Tran, T. T.; Choi, D.-Y.; Kivshar, Y.; Solntsev, A. S. Quasi-BIC resonant enhancement of second-harmonic generation in WS2 monolayers. _Nano Lett._ 2020, _20_ , 5309–5314
  * Shcherbakov et al. 2019 Shcherbakov, M. R.; Shafirin, P.; Shvets, G. Overcoming the efficiency-bandwidth tradeoff for optical harmonics generation using nonlinear time-variant resonators. _Phys. Rev. A_ 2019, _100_ , 063847
  * Grinblat et al. 2020 Grinblat, G.; Zhang, H.; Nielsen, M. P.; Krivitsky, L.; Berté, R.; Li, Y.; Tilmann, B.; Cortés, E.; Oulton, R. F.; Kuznetsov, A. I.; Maier, S. A. Efficient ultrafast all-optical modulation in a nonlinear crystalline gallium phosphide nanodisk at the anapole excitation. _Science Advances_ 2020, _6_
  * Shcherbakov et al. 2020 Shcherbakov, M. R.; Zhang, H.; Tripepi, M.; Sartorello, G.; Talisa, N.; AlShafey, A.; Fan, Z.; Twardowski, J.; Krivitsky, L. A.; Kuznetsov, A. I.; Chowdhury, E.; Shvets, G. Generation of even and odd high harmonics in resonant metasurfaces using single and multiple ultra-intense laser pulses. 2020
  * Joannopoulos et al. 2008 Joannopoulos, J. D.; Johnson, S. G.; Winn, J. N.; Meade, R. D. _Photonic Crystals: Molding the Flow of Light_ , 2nd ed.; Princeton University Press, 2008
  * Koshelev et al. 2018 Koshelev, K.; Lepeshov, S.; Liu, M.; Bogdanov, A.; Kivshar, Y. Asymmetric Metasurfaces with High- Q Resonances Governed by Bound States in the Continuum. _Phys. Rev. Lett._ 2018, _121_ , 193903
  * Fano 1961 Fano, U. Effects of configuration interaction on intensities and phase shifts. _Phys. Rev._ 1961, _124_ , 1866–1878
  * Shoji et al. 1997 Shoji, I.; Kondo, T.; Kitamoto, A.; Shirane, M.; Ito, R. Absolute scale of second-order nonlinear-optical coefficients. _J. Opt. Soc. Am. B_ 1997, _14_ , 2268–2294
  * Bond 1965 Bond, W. L. Measurement of the Refractive Indices of Several Crystals. _J. Appl. Phys._ 1965, _36_ , 1674–1677
  * Wilson et al. 2020 Wilson, D. J.; Schneider, K.; Hoenl, S.; Anderson, M.; Kippenberg, T. J.; Seidler, P. Integrated gallium phosphide nonlinear photonics. _Nat. Photon._ 2020, 57–62
  * Padath Anthur et al. 2020 Padath Anthur, A.; Zhang, H.; Akimov, Y.; Ong, J. R.; Kalashnikov, D.; I. Kuznetsov, A.; Krivitsky, L. Demonstration of second harmonic generation in gallium phosphide nano-waveguides. _Arxiv_ 2020, _2001.06142v2_
  * Stratton 1941 Stratton, J. A. _Electromagnetic Theory_ , 1st ed.; McGraw-Hill Book Company, Inc., 1941
  * Kruk et al. 2017 Kruk, S. S.; Camacho-Morales, R.; Xu, L.; Rahmani, M.; Smirnova, D. A.; Wang, L.; Tan, H. H.; Jagadish, C.; Neshev, D. N.; Kivshar, Y. S. Nonlinear optical magnetism revealed by second-harmonic generation in nanoantennas. _Nano Lett._ 2017, _17_ , 3914–3918
  * Ghirardini et al. 2017 Ghirardini, L.; Carletti, L.; Gili, V.; Pellegrini, G.; Duò, L.; Finazzi, M.; Rocco, D.; Locatelli, A.; De Angelis, C.; Favero, I.; Ravaro, M.; Leo, G.; Lemaître, A.; Celebrano, M. Polarization properties of second-harmonic generation in AlGaAs optical nanoantennas. _Opt. Lett._ 2017, _42_ , 559–562
  * Vabishchevich et al. 2018 Vabishchevich, P. P.; Liu, S.; Sinclair, M. B.; Keeler, G. A.; Peake, G. M.; Brener, I. Enhanced Second-Harmonic Generation Using Broken Symmetry III-V Semiconductor Fano Metasurfaces. _ACS Photonics_ 2018, _5_ , 1685–1690
  * Liu et al. 2016 Liu, S.; Sinclair, M. B.; Saravi, S.; Keeler, G. A.; Yang, Y.; Reno, J.; Peake, G. M.; Setzpfandt, F.; Staude, I.; Pertsch, T.; Brener, I. Resonantly Enhanced Second-Harmonic Generation Using III-V Semiconductor All-Dielectric Metasurfaces. _Nano Lett._ 2016, _16_ , 5426–5432
  * Semmlinger et al. 2018 Semmlinger, M.; Tseng, M. L.; Yang, J.; Zhang, M.; Zhang, C.; Tsai, W. Y.; Tsai, D. P.; Nordlander, P.; Halas, N. J. Vacuum Ultraviolet Light-Generating Metasurface. _Nano Lett._ 2018, _18_ , 5738–5743