# Doubly resonant photonic crystal cavities for efficient
second-harmonic generation in III-V semiconductors

Simone Zanotti  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
27100 Pavia (Italy)  Momchil Minkov  Department of Electrical Engineering, and
Ginzton Laboratory, Stanford University, Stanford, California 94305 (USA)
Shanhui Fan  Department of Electrical Engineering, and Ginzton Laboratory,
Stanford University, Stanford, California 94305 (USA)  Lucio C. Andreani
Dipartimento di Fisica, Università di Pavia, via Bassi 6, 27100 Pavia (Italy)
Institute for Photonics and Nanotechnologies (IFN)-CNR, 20133 Milano (Italy)
Dario Gerace<EMAIL_ADDRESS>Dipartimento di Fisica, Università di
Pavia, via Bassi 6, 27100 Pavia (Italy)

###### Abstract

Second-order nonlinear effects, such as second-harmonic generation, can be
strongly enhanced in nanofabricated photonic materials when both fundamental
and harmonic frequencies are spatially and temporally confined. Practically
designing low-volume and doubly resonant nanoresonators in conventional
semiconductor compounds is challenging owing to their intrinsic refractive
index dispersion. In this work we review a recently developed strategy to
design doubly resonant nanocavities with low mode volume and large quality
factor by localized defects in a photonic crystal structure. We build on this
approach by applying an evolutionary optimisation algorithm in connection with
Maxwell equations solvers, showing that the proposed design recipe can be
applied to any material platform. We explicitly calculate the second-harmonic
generation efficiency for doubly resonant photonic crystal cavity designs in
typical III-V semiconductor materials, such as GaN and AlGaAs, targeting a
fundamental harmonic at telecom wavelengths, and fully accounting for the
tensor nature of the respective nonlinear susceptibilities. These results may
stimulate the realisation of small footprint photonic nanostructures in
leading semiconductor material platforms to achieve unprecedented nonlinear
efficiencies.

##  I Introduction

Nonlinear optical processes are notoriously poorly efficient, due to small
higher-order nonlinear susceptibilities $(\chi^{(n)})$ that mediate such
processes in conventional materials such as semiconductors or insulators in
their transparency spectral range Boyd (2008). A way out this limiting factor
can be tackled by strongly confining the electromagnetic radiation in
dielectric resonators made of such nonlinear materials. The confinement of the
radiation in a small volume ($V\sim\lambda^{3}$) and for a long time (high-
quality factor $Q$) allows for a strong interaction between the
electromagnetic field and the nonlinear dielectric, thus enhancing the
efficiency of nonlinear effects Lin  _et al._ (2016). Among the possible
nanoresonators, photonic crystal (PhC) cavities offer one of the most
promising platforms to achieve such extreme confinement conditions. In the
last few years, many efforts have been made to achieve enhanced second-order
nonlinearities in both singly-resonant PhC cavities Rivoire  _et al._ (2009);
Buckley  _et al._ (2014), where the radiation is trapped only for the first-
harmonic (FH) frequency, $\omega_{\mathrm{FH}}=\omega$, and doubly-resonant
PhC cavities Lin  _et al._ (2016); Liscidini and Claudio Andreani (2006);
Minkov  _et al._ (2019), in which the confinement simultaneously occurs for
both FH and second-harmonic (SH) frequency, $\omega_{\mathrm{SH}}=2\omega$.
Hence, doubly-resonant PhC cavities are also among the most promising choices
for enhancing second-order nonlinear effects, e.g. second harmonic-generation
(SHG) (see Figure 1) and spontaneous parametric down-conversion (SPDC), since
their efficiency benefits from the simultaneous confinement of both FH and SH
modes, provided they are coupled by the nonlinear tensor elements. While PhC
cavities possess several degrees of freedom to be used in order to tailor
their optical properties, very few of such devices have been studied to fulfil
a doubly resonant condition that could boost SH generation efficiency, until
very recently Wang  _et al._ (2020). The main obstacle comes from the
difficulty in designing PhCs with photonic band gaps around the FH and SH
frequencies, which would allow obtaining two doubly-resonant cavity modes
Andreani and Gerace (2006). Recently, a novel approach based on the so-called
bound-states in the continuum (BIC) has been proposed Minkov  _et al._ (2019),
overcoming the aforementioned obstacle. This strategy has later been shown to
be practically effective by realizing doubly resonant SH generation in a III-V
semiconductor PhC cavity Wang  _et al._ (2020). The BICs are particular states
that stay confined despite the fact they lie above the light line, i.e., their
dispersion falls within the continuum of radiative modes Hsu  _et al._ (2016).
This feature is indispensable to achieve good temporal confinement of the SH
cavity mode, even in the absence of a photonic bandgap Minkov  _et al._
(2019).

Figure 1: Schematic representation of the second-harmonic generation process,
where the electromagnetic radiation is converted from the frequency $\omega$
to its double frequency $2\omega$.

In the present article, we review the design strategy based on matching the FH
of a PhC lattice with a suitable BIC mode at SH frequency, and we generalize
the procedure to show that doubly-resonant PhC cavities based on BICs can be
properly designed by applying a particle swarm optimisation (PSO) algorithm
combined with a numerical solver for Maxwell equations Robinson and Rahmat-
Samii (2004); in particular, we use the finite-difference time-domain (FDTD)
algorithm here Taflove and Hagness (2005), which allows a precise
determination of the higher order mode resonances and is amenable to be
extended to large cavity structures. The latter are shown to allow for both FH
and SH mode confinement in the plane through a heterostructure engineering of
the PhC lattice. Moreover, in determining the nonlinear properties we directly
take into account the tensor nature of the second-order nonlinear
susceptibility, $\chi_{i,j,k}^{(2)}$ Lin  _et al._ (2016), a key feature to
consider when estimating their efficiency, which is sometimes overlooked in
the literature. This allowed us to numerically calculate the SHG efficiency
and its dependence on the crystal growth direction with respect to the PhC
cavity axes, through a dimensionless overlap integral ($\bar{\beta}$) between
the FH and SH modes taking into account the $\chi^{(2)}$ tensor. The latter
can be regarded as a generalisation of the phase-matching condition for the
confined systems. The factor $\bar{\beta}$ is directly related to the
estimation of the SHG efficiency, which can be approximated to scale as
$Q_{1}^{2}Q_{2}|\bar{\beta}|^{2}$, $Q_{1}$ and $Q_{2}$ being the FH and SH
modes quality factors, respectively Lin  _et al._ (2016); Liscidini and
Claudio Andreani (2006). Our numerical calculations show that different
crystal structures lead to different dependencies of the SHG efficiency on the
crystal axes orientation. In particular, we investigated the two most common
crystal structures in $\textup{III}\\!-\\!\textup{V}$ semiconductors, namely
the zincblende and the wurtzite whose typical examples are AlGaAs and GaN,
respectively. These compounds have gained increasing interest in nonlinear
applications due to their strong nonlinear properties, and for being
transparent at both FH/SH frequencies targeted at telecom/near infrared
wavelengths, i.e., $\lambda_{\text{FH}}\approx$ $1550\text{\,}\mathrm{nm}$ and
$\lambda_{\text{SH}}\approx$ $775\text{\,}\mathrm{nm}$. It is found that the
$[0,0,1]$ growth direction of the zincblende crystalline materials lead to a
vanishing overlap integral for specific orientations of the crystal axes with
respect to the PhC lattice. In these unfavourable cases, the conversion
efficiency is theoretically suppressed, despite the strong temporal and
spatial confinements. However, we show in the present work that a suitable
growth direction or crystal structure, such as $[1,1,1]$ grown zincblende, or
the wurtzite structure, lead to an isotropic generation efficiency that does
not depend on the PhC orientation against the crystal unitary cell.

The manuscript is organized as follows. In Sec. II we review the design
procedure to obtain a doubly resonant photonic crystal structure in which the
SH mode is a BIC, generalize it by combining Maxwell solvers with evolutionary
optimization, and explicitly apply it to the most promising III-V platforms in
nonlinear optics. In Sec. III we show that a simultaneous spatial confinement
of both FH and SH modes can be achieved by applying a varying hole radii
photonic heterostructure. In Sec. IV we calculate the nonlinear conversion
efficiency in such doubly resonant photonic crystal cavities, taking into
account the crystal structure of the III-V material and the orientation of the
photonic lattice with respect to the semiconductor unit cell.

##  II Evolutionary optimization of doubly resonant photonic crystals

###  II.1 Numerical methods

The present design strategy to obtain doubly-resonant PhC cavities relies on
two numerical methods to solve Maxwell equations in nanofabricated dielectric
materials, namely the guided-mode expansion (GME) and the FDTD methods. The
former method is well suited to the analysis of layered systems with in-plane
periodicity of the dielectric function, such as PhC slabs, providing a
solution for the frequency-wavevector dispersion of the real and imaginary
part of complex eigenmodes Andreani and Gerace (2006). Here, we consider only
PhC slabs with holes carved in the dielectric slab with a triangular lattice,
but these are not intrinsic constraints and other lattices could be considered
as well. In this case, the main parameters are the lattice constant $a$, and
the (dimensionless) radius $r/a$ and slab thickness $d/a$. The GME starts from
a decomposition of the electromagnetic field in a basis of plane waves and
guided modes of an effective planar waveguide, which leads to an eigenvalue
problem whose eigenvalues correspond to the real part of the complex
eigenfrequencies of the PhC modes Andreani and Gerace (2006). Moreover,
perturbation theory allows to calculate the coupling of the PhC slab modes to
the out-of-plane propagating radiation modes of the effective planar waveguide
with average refractive index, which can be considered the loss rate of the
PhC slab modes defined as the imaginary part of the complex eigenmode
frequencies Ochiai and Sakoda (2001). This approach turns out to be
computationally very efficient and quantitatively accurate to within a few
percent error Minkov  _et al._ (2020), albeit approximate and polarization-
dependent accuracy. Hence, it has been applied for a rapid and coarse scanning
of the parameters, and for a fast identification of higher order modes in the
SH frequency range.

In parallel, the FDTD method provides the electrodynamics solution of Maxwell
equations for arbitrary complex systems, by discretizing and solving the
finite difference formulation of the coupled differential equations Taflove
and Hagness (2005). The FDTD algorithm provides a numerically exact method,
i.e. with sufficiently fine discretization its solution should converge to the
exact results. In this case, the imaginary part of the frequency of a given
electromagnetic mode of a structure is obtained by fitting the exponential
decay of the fields associated to that mode. This is best achieved by
selectively exciting the mode using symmetry considerations, and waiting for
low-Q modes to decay first. Since the FDTD is accurate but computationally
costly, it has been applied only after the target modes have been identified
with the GME solution. Throughout this work, we have employed a commercial
implementation of the three dimensional (3D) FDTD solver provided from
Lumerical-Ansys. Moreover, a particle swarm optimisation (PSO) algorithm has
been exploited in connection with the 3D-FDTD solver in order to finely tune
the PhC parameters Robinson and Rahmat-Samii (2004). Such a routine allows to
efficiently minimise a given function $f:A\to\mathbb{R}$, also named figure of
merit, by varying the parameters in the search-space $A\subset\mathbb{R}^{n}$,
further information on this algorithm are provided in Appendix B. In the
present case, the PSO was chosen since it allows to reach convergence within a
few evolutionary steps when the starting point is not too far from the target
objectives, as it is most often the case when targeting the doubly resonant
condition.

###  II.2 Design procedure

The present design strategy to obtain doubly-resonant PhC cavities builds on a
previously presented band engineering recipe, aimed at finding a SH band with
a mode displaying infinite lifetime for emission along the vertical direction
with respect to the PhC plane, i.e. a BIC Minkov  _et al._ (2019). This
approach has also been recently shown to experimentally enhance SHG efficiency
in GaN PhC cavities Wang  _et al._ (2020). Here we show that the above
strategy can be further generalized and automated, and it can be applied to an
arbitrary material platform. The procedure can be summarised as follows:

  1. 1.

Calculate the photonic bandstructure of the PhC slab with GME.

  2. 2.

Identify a band-edge mode below the light line as FH mode at frequency
$\omega_{1}$, and a non-degenerate BIC mode above the light line at frequency
$\omega_{2}$ as SH mode, whose photonic mode dispersion depends on the
refractive indices of the target material.

  3. 3.

Check that FH and SH modes have non-vanishing overlap integral $\bar{\beta}$,
depending on the $\chi^{(2)}$ tensor of the given material.

  4. 4.

Vary the hole radius and the slab thickness in order to approximately match
the doubly resonant condition $\omega_{2}\approx 2\,\omega_{1}$.

  5. 5.

Starting from the parameters calculated with GME, apply the PSO in connection
with the FDTD Maxwell solver to finely tune the doubly resonant condition and
reach $\omega_{2}=2\,\omega_{1}$ within the accepted tolerance.

  6. 6.

Spatially confine the FH and SH mode by suitably designing a heterostructure
PhC cavity.

###  II.3 Doubly-resonant photonic crystal slabs

Figure 2: (a) Real part of the GaN even modes (TE-like) bandstructure
calculated with the GME and the FDTD methods, respectively. The parameters
chosen for these simulations are: $n=2.28$, $d/a=0.337$ and $r/a=0.206$. The
$k-$space region surrounding of the FH mode at the M point is highlighted,
while the light lines are dashed. (b) Representation of the thin dielectric
slab with a triangular photonic lattice, with an indication of the relevant
symmetry plane (at $z=0$). (c) Main symmetry points considered within the
Brillouin zone.

With the target of optimizing the doubly resonant condition for III-V
semiconductors as the relevant nonlinear materials, we consider that FH mode
with dominant in-plane polarization might be coupled by the second-order
tensor elements to SH modes with dominant vertical polarization. For this
reason, we consider the photonic mode dispersion calculated for membrane PhC
slabs (i.e., with symmetric air claddings above and below the patterned
semiconductor), with FH modes having even symmetry (transverse-electric, TE-
like) with respect to the horizontal ($x,y$) plane bisecting the planar
waveguide, and SH modes possessing odd symmetry (transverse magnetic, TM-like)
Andreani and Gerace (2006); Minkov  _et al._ (2020). These are shown for the
specific case of a GaN PhC slab in Figs. 2 and 3, respectively. The photonic
lattice with the definition of the relevant real space axes and the Brillouin
zone in reciprocal space are also shown, in Fig. 2(b)-(c). In order to take
into account the material dispersion, the two solutions have been obtained
with different refractive indices: for the FH modes, we set $n=2.28$, while
for SH modes $n=2.31$. These values correspond to the GaN refractive indices
at the target wavelengths $\lambda_{\text{FH}}=$ $1550\text{\,}\mathrm{nm}$
and $\lambda_{\text{SH}}=$ $775\text{\,}\mathrm{nm}$, respectively. The
numerical solutions are explicitly shown for both GME and 3D-FDTD solvers, and
directly compared to each other. For the FDTD solution, broad-band, point-like
electric dipole sources are randomly placed in the whole elementary cell in
order to excite all the possible modes, and proper phases are given to each
dipole in order to satisfy the Bloch boundary conditions for given in-plane
wave vectors; thus, each wave vector in the first Brillouin zone of the
corresponding photonic lattice requires a separate FDTD simulation, and by
running over all the k-vectors, the whole band structure is reconstructed.
The dispersion of photonic bands (real part of complex eigenfrequencies) is
shown for both FH and SH modes in Figs. 2(a) and 3(a), while the imaginary
part is directly compared only for the target BIC mode in the SH case, Fig.
3(b). It is worth stressing that the results obtained with the two methods are
in overall good agreement, quite remarkably considering that they rely on
conceptually different methods, and moreover the GME method is an approximate
method. These results suggest that the GME method is quite effective in
simulating such photonic crystal devices. Notice that the imaginary part of
the BIC mode goes to zero at the $\Gamma$ point in both solutions. The
agreement is slightly worse for the odd modes, which may also expected owing
to the approximations inherent to the GME approach. Hence, this comparison
justifies employing GME for a fast parameter scan, while the FDTD method is
particularly suited to give accurate results in the SH case with TM-like
symmetry.

In view of defining confined modes, we identify the FH mode with the second
band of the TE-like bandstructure at the M point, which is located below the
air light line and is likely to give rise to a high quality factor cavity mode
upon confinement Minkov  _et al._ (2019). On the other hand, the SH mode
inevitably lies above the air light line, as it is clear from the TM-like
bandstructure; for this reason, one of the most promising modes to consider
for a potentially high quality factor confined mode at SH is the seventh band
at $\Gamma$, which is a BIC in the $\Gamma$ point as already stressed. These
modes are selected as a starting point for the evolutionary optimization aimed
at searching the photonic structure parameters giving doubly resonant
condition.

Figure 3: (a) Real part of the GaN odd modes (TM-like) bandstructure
calculated with the GME and the FDTD methods. The parameters chosen for these
simulations are: $n=2.31$, $d/a=0.337$ and $r/a=0.206$. The $k-$space region
surrounding the SH mode at the $\Gamma$ point is highlighted. (b) Imaginary
parts of the eigenmodes highlighted in the bandstructure (a). The vanishing
Im$(\omega)$ indicates the presence of a BIC at the $\Gamma$ point, both from
GME and FDTD solutions.

Once the candidate FH and SH modes are identified with the GME solution, a
coarse span of the parameters is necessary to ”move close” to the doubly
resonant condition. At this point, the PSO can be implemented in the FDTD
algorithm to accurately tune the PhC slab parameters. In this case, we let the
optimisation algorithm vary the radius of the holes and the slab thickness in
the search-space
$A=[r_{\text{min}},r_{\text{max}}]\times[d_{\text{min}},d_{\text{max}}]$,
where $A$ is a small neighbourhood around the best GME parameters; while the
figure of merit to be minimised has been defined as:

$\text{FOM}=\bigg{\lvert}\frac{\omega_{2}-2\omega_{1}}{\omega_{2}}\bigg{\rvert},$
(1)

where $\omega_{1}$ and $\omega_{2}$ are the FH and SH mode frequency,
respectively. We imposed the condition $\text{FOM}\leq 1\%$ to be fulfilled,
such that the doubly resonant condition is achieved with a given tolerance.

Table 1: Main design parameters of the GaN and Al_0.3Ga_0.7As PhC slabs obtained with the PSO algorithm applied to a 3D-FDTD solver. Material | $n_{\text{FH}}$ | $n_{\text{SH}}$ | $d/a$ | $r/a$ | $\omega_{1}a/(2\pi c)$ | $\omega_{2}a/(2\pi c)$ |  $a$ ($\mathrm{nm}$) |  $d$ ($\mathrm{nm}$) |  $r$ ($\mathrm{nm}$)
---|---|---|---|---|---|---|---|---|---
GaN | $2.28$ | $2.31$ | $0.337$ | $0.206$ | $0.414$ | $0.832$ | $642.3$ | $216.4$ | $132.2$
Al_0.3Ga_0.7As | $3.23$ | $3.47$ | $0.24$ | $0.175$ | $0.3356$ | $0.6650$ | $520.1$ | $124.8$ | $91.2$

To show the generality of the approach, the very same procedure has been
applied to both GaN and Al_x Ga_1-x As materials (assuming $x=0.3$ for the
choice of room temperature refractive indices), targeting the same bands in
the same high symmetry points within the Brillouin zone to optimize the doubly
resonant condition. All the quantities are summarised in dimensionless units
in table 1. They can also easily be converted to dimensional units by imposing
the target wavelength, e.g. $\lambda_{\text{FH}}=$ $1550\text{\,}\mathrm{nm}$
for the FH mode, yielding $a=\tilde{\omega}_{1}\lambda_{\text{FH}}$,
$\tilde{\omega}_{1}$ being the FH dimensionless frequency, from which the
lattice constant $a$, and all other quantities, can be calculated.

Figure 4: (a) The heterostructure cavity is defined by core, transition, and
outer regions. The hole radius is gradually increased on going outwards from
the core region. (b) Schematic representation of the SHG process in the doubly
resonant PhC cavity.

##  III Doubly-resonant photonic crystal cavities

In the previous section we have outlined how to achieve the doubly resonant
condition in PhC slabs. While in this situation the electromagnetic field is
bound in the vertical direction, $z$, it is still fully delocalized in the
$x,y$ plane. In view of enhancing the nonlinear efficiency, it is necessary to
introduce a heterostructure design to achieve the field confinement also along
the other two spatial directions. Such a heterostructure ultimately aims at
producing a confining potential for the photons, under the very same principle
allowing for charged carriers to be confined within a semiconductor
heterostructure. As it is commonly observed, for the photonic crystal slab
system considered here, the resonators with highest quality factors tend to
confine photons in the regions of higher dielectric constant Ge  _et al._
(2018). Hence, the cavity is designed by assuming three concentric hexagonal
regions with gradually increasing radii, such that an effective potential well
is created for photons. The structure geometry is shown in Fig. 4, and it is
characterized by the side length of the hexagons ($N_{c},N_{t},N_{o}$) in
terms of lattice constant and the holes radii ($r_{c},r_{t},r_{o}$), in which
the subscripts ($c,t,o$) refer to core, transition and outer regions,
similarly to Ref. Minkov  _et al._ (2019). For consistency, we set $N_{c}=6$,
$N_{t}=4$ and $N_{o}=10$, for simulations on both GaN and Al_0.3Ga_0.7As
platforms. The slab thickness and the hole radius in the core region are
inherited from the PhC slabs parameters, while the radii of holes within the
transition and outer regions are slightly increased to allow for a gentle
confinement of the electromagnetic field. Due to small frequency shifts, a
slight tuning of the parameters may be necessary in order to restore the
doubly resonance condition in the cavity design.

Table 2: Main design parameters of the GaN and Al_0.3Ga_0.7As cavities, and the related doubly resonant features obtained from 3D FDTD simulations. Material |  $a$ ($\mathrm{nm}$) |  $d$ ($\mathrm{nm}$) |  $r_{c}$ ($\mathrm{nm}$) |  $r_{t}$ ($\mathrm{nm}$) |  $r_{o}$ ($\mathrm{nm}$) | $\lambda_{\text{FH}}$ | $\lambda_{\text{SH}}$ | $Q_{\text{FH}}$ | $Q_{\text{SH}}$
---|---|---|---|---|---|---|---|---|---
GaN | $642$ | $200$ | $122$ | $145$ | $165$ | $1503$ | $753$ | $3.6\cdot 10^{4}$ | $1100$
Al_0.3Ga_0.7As | $520$ | $124$ | $91$ | $105$ | $121$ | $1523$ | $759$ | $1.1\cdot 10^{5}$ | $400$
Figure 5: (a)-(c) FH and SH cavity mode intensity profiles at the centre of
the slab $z=0$. (b)-(d) Colour maps representing the emission intensity
profile of the FH and SH cavity modes. Every green circle identifies
$10\text{\,}\mathrm{\SIUnitSymbolDegree}$ from the vertical direction of
emission with respect the slab plane. The white arrows show the vector field
defined by ($\Re{\left(E_{x}\right)},\Re{\left(E_{y}\right)}$).

We studied the properties of the same heterostructure cavity design realized
in both GaN and Al_0.3Ga_0.7As materials, by using 3D FDTD simulations. We
employed the same real space meshing for both simulations at fundamental and
second-harmonic frequencies, by checking that the second-harmonic one,
corresponding to about 30 mesh steps per wavelength (i.e., to a minimum step
mesh of 25 nm for both material platforms), was at convergence (i.e.,
resonance frequencies stable to within percent accuracy). The cavity modes are
excited by point-like sources of electromagnetic radiation, once the emission
of the sources is elapsed, the decaying fields are recorded in space and time
domains. Such information allows to fully recover the frequencies and quality
factors of both FH and SH cavity modes. The outcomes are reported in table 2.

Moreover, from the FDTD simulations it is possible to retrieve the
electromagnetic field profiles inside the cavities, as well as the emission
properties of the cavity modes calculated from a near-to-far field projection
of the 2D cavity modes recorded at the surface of the PhC slab. From the
intensity profiles of the FH and SH modes that are shown in Fig. 5, it is
clear how the radiation is well confined in the core region in both cases. The
FH farfield intensity patterns are highly focused around the vertical
direction of emission, which is crucial to facilitate the vertical in-coupling
of the electromagnetic radiation coming from a pump beam with a gaussian
intensity profile; while the SH mode exhibits a donut-shaped farfield
intensity pattern, which can be directly attributed to its quasi-BIC nature
Wang  _et al._ (2020).

Figure 6: (a) Nonlinear overlap factor between FH and SH modes in the doubly
resonant PhC cavity design, $|\bar{\beta}|^{2}$, for GaN and Al_0.3Ga_0.7As
grown along different directions, as a function of the rotation angle in the
PhC lattice plane where $None$ indicates the [1,0,0] direction. A schematic of
the unit cells and the definition of the rotation angles for the three cases
considered is reported in (b); these crystallographic pictures were obtained
through the VESTA software Momma and Izumi (2011).

##  IV Overlap factor and nonlinear efficiency

Having both spatial and temporal confinements of FH and SH modes is not a
sufficient condition to guarantee an efficient second-order nonlinear
efficiency. In fact, the spatial overlap between the FH and SH fields has to
be suitably considered, as well as their confinement volume. For second-order
nonlinear processes, a nonlinear overlap between the FH and SH modes can
generally be defined as Lin  _et al._ (2016):

$\bar{\beta}=\frac{\int
d\mathbf{r}\,\sum_{ijk}\bar{\chi}^{(2)}_{ijk}E_{2\omega i}^{*}E_{\omega
j}E_{\omega k}}{\left(\int
d\mathbf{r}\,\varepsilon_{\omega}(\mathbf{r})|\mathbf{E}_{\omega}|^{2}\right)\left(\int
d\mathbf{r}\,\varepsilon_{2\omega}(\mathbf{r})|\mathbf{E}_{2\omega}|^{2}\right)^{1/2}}\sqrt{\lambda_{\text{FH}}^{3}},$
(2)

where $\mathbf{E}_{\omega}$ and $\mathbf{E}_{2\omega}$ are the FH and SH
electric field profiles, $\varepsilon_{\omega}$ and $\varepsilon_{2\omega}$
are the dielectric functions at FH and SH frequencies, while
$\bar{\chi}^{(2)}_{ijk}$ are dimensionless nonlinear tensor elements defined
through $\bar{\chi}^{(2)}_{ijk}=\chi^{(2)}_{ijk}/\chi^{(2)}$ with $\chi^{(2)}$
chosen to be the main tensor element of the relative compound grown in the
$[0,0,1]$ direction:

$\displaystyle\text{Zincblende}:\qquad\chi^{(2)}=\chi^{(2)}_{zxy}$ (3)
$\displaystyle\text{Wurtzite}:\qquad\quad\chi^{(2)}=\chi^{(2)}_{xxz}.$

The overlap factor in Eq. (2) is expressed in dimensionless units, and it is
particularly useful to compare different doubly resonant structures, also made
of different materials. In fact, the constitutive material of the cavity and
its crystallographic orientation are directly involved in the calculation of
the integral in Eq. (2) through the tensor elements $\bar{\chi}^{(2)}_{ijk}$.
Following the Kleinmann symmetry Boyd (2008), the non-vanishing tensor
elements of the wurtzite (e.g., GaN) and zincblende (e.g., Al_0.3Ga_0.7As)
crystal structures grown along the $[0,0,1]$ direction are:

$\displaystyle\text{Zincblende:
}{\scriptstyle\bar{\chi}^{(2)}_{xyz}=\bar{\chi}^{(2)}_{zxy}=\bar{\chi}^{(2)}_{yzx}=\bar{\chi}^{(2)}_{xzy}=\bar{\chi}^{(2)}_{yxz}=\bar{\chi}^{(2)}_{zyx}=1}$
(4) $\displaystyle\text{Wurtzite:
}{\scriptstyle\,\bar{\chi}^{(2)}_{xxz}=\bar{\chi}^{(2)}_{zxx}=\bar{\chi}^{(2)}_{xzx}=\bar{\chi}^{(2)}_{yyz}=\bar{\chi}^{(2)}_{zyy}=\bar{\chi}^{(2)}_{yzy}=1,\;\bar{\chi}^{(2)}_{zzz}=-2.}$

We supposed $\bar{\chi}^{(2)}_{zzz}=-2$, this condition is not always
experimentally verified but since the FH is mostly TE-polarized, the actual
value of the $\bar{\chi}^{(2)}_{zzz}$ tensor element has a negligible effect
on the beta-factor of Eq. (2), as we have verified numerically. Starting from
(4), we can recover the tensor elements for an arbitrary orientation of the
crystal axes with respect to the PhC lattice. This can be accomplished by
applying a proper rotation in the three-dimensional space, described by a
rotation matrix $\mathbf{R}(\phi,\theta,\psi)$ in which $(\phi,\theta,\psi)$
are the Euler angles (see App. A for details). For example, the transformation
that allows to switch from the $[0,0,1]$ direction to the $[1,1,1]$ growth
direction is

$\displaystyle\bar{\chi}^{(2)\,[1,1,1]}_{\alpha,\beta,\gamma}=\sum_{a,b,c}U_{\alpha,a}U_{\beta,b}U_{\gamma,c}\,\bar{\chi}^{(2)\,[0,0,1]}_{a,b,c},$
(5)
$\displaystyle\mathbf{U}\doteq\mathbf{R}\left(\frac{3}{4}\pi,\cos^{-1}{\left(\frac{1}{\sqrt{3}}\right)},\frac{\pi}{2}\right).$

Here, we calculated the overlap integral (2) in the case of Al_0.3Ga_0.7As
grown both along the $[0,0,1]$ and $[1,1,1]$ directions. The GaN is only
considered to be grown along the $[0,0,1]$. Then, we also considered its
dependence as a function of a rotation in the PhC lattice plane, i.e. the
relative orientation between the PhC lattice and the crystallographic
directions of the underlying nonlinear material. The squared modulus of the
overlap integral is reported in Fig. 6(a) for both material platforms. A
definition of the rotation angle is schematically given in Fig. 6(b).

Figure 7: (a)-(b) SHG conversion efficiency for the Al_0.3Ga_0.7As and GaN doubly resonant PhC cavities, calculated for varying FH and SH cavity factors. The green circles identify the expected efficiency for the cavities designed in the present work. Table 3: Calculated nonlinear figures of merit for both the GaN and Al_0.3Ga_0.7As doubly resonant PhC cavities designed in this work. Material | $Q_{\text{FH}}$ | $Q_{\text{SH}}$ | $|\bar{\beta}_{\text{max}}|^{2}$ |  $\chi^{(2)}$ ($\mathrm{pm}\text{\,}{\mathrm{V}}^{-1}$) |  $P_{o}/P_{i}^{2}$ (${\mathrm{W}}^{-1}$)
---|---|---|---|---|---
GaN | $3.6\cdot 10^{4}$ | $1100$ | $1.0\cdot 10^{-4}$ | $10$ | $2.9$
Al_0.3Ga_0.7As | $1.1\cdot 10^{5}$ | $400$ | $1.2\cdot 10^{-5}$ | $100$ | $119$

The overlap factor, $|\bar{\beta}|^{2}$, displays different qualitative trends
depending on the crystal structure and its relative orientation. In
particular, the largest nonlinear overlap is obtained for the GaN where
$|\bar{\beta}|^{2}$ goes beyond $10^{-4}$ for any orientation angle, while in
the Al_0.3Ga_0.7As it is modulated and can vanish at specific orientation
angles for the $[0,0,1]$ growth direction. On the other hand, the
$|\bar{\beta}|^{2}$ does not depend on the relative orientation between the
crystallographic directions and the PhC lattice for the Al_0.3Ga_0.7As grown
along the $[1,1,1]$, which also displays an overall larger values as compared
to the same material grown along $[0,0,1]$.

While the dimensionless overlap factor gives the dependence of the conversion
efficiency on the crystallographic axes orientation, it only gives partial
information since the physical process of nonlinear conversion involves the
actual magnitude of the $\chi^{(2)}$ tensor elements, which is material
dependent and is not explicitly taken into account in Eq. 2. This additional
parameter can play a crucial role in calculating the final nonlinear
efficiency, which can be estimated with the following relation Lin  _et al._
(2016):

$\frac{P_{o}}{P_{i}^{2}}=\frac{8}{\omega_{1}}\left(\frac{\chi^{(2)}}{\sqrt{\varepsilon_{0}\lambda^{3}_{\text{FH}}}}\right)^{2}|\bar{\beta}|^{2}\,Q^{2}_{\text{FH}}Q_{\text{SH}},$
(6)

where $\varepsilon_{0}$ is the vacuum permittivity. A detailed derivation of
Eq. (6) is reported in App. C for completeness. We also notice that Eq. (6) is
valid under the undepleted-pump approximation, and considering perfect in- and
out-coupling of the radiation for both FH and SH modes Minkov  _et al._
(2019). Notice that similar scaling relations can be found for other nonlinear
generation processes, e.g., SPDC Helt  _et al._ (2012). Notice also that this
expression is explicitly obtained for narrow-band harmonic resonances, at
variance with SHG in broadband photonic crystal waveguide structures Rivoire
_et al._ (2011). We evaluated the conversion efficiencies for a wide range of
$Q_{\text{FH}}$ and $Q_{\text{SH}}$. We set $\chi^{(2)}=$
$100\text{\,}\mathrm{pm}\text{\,}{\mathrm{V}}^{-1}$ for Al_0.3Ga_0.7As, and
$\chi^{(2)}=$ $10\text{\,}\mathrm{pm}\text{\,}{\mathrm{V}}^{-1}$ for GaN,
respectively Duchesne  _et al._ (2011); Sanford  _et al._ (2005). Moreover, we
considered the axes orientations maximizing the nonlinear overlap factor
according to the results in Fig. 6(a). The results are reported in Fig. 7,
which inverts the outcome of the previous Figure: in fact, the Al_0.3Ga_0.7As
doubly resonant cavity (grown along [1,1,1]) is theoretically more efficient
than the GaN one, despite the dimensionless $|\bar{\beta}|^{2}$ factor being
one order of magnitude larger in favor of the latter. This can mostly be
attributed to the larger nonlinear properties of Al_0.3Ga_0.7As, which is
characterised by an higher $\chi^{(2)}_{\text{eff}}$ . If we restrict
ourselves to the quality factors calculated in the present work with perfectly
in- and out-coupled beams, we can estimate an efficiency of
$119\text{\,}{\mathrm{W}}^{-1}$ for the Al_0.3Ga_0.7As, and
$2.9\text{\,}{\mathrm{W}}^{-1}$ for the GaN doubly resonant PhC cavity, as
summarized in table 3.

##  V Conclusions

We have presented a numerical procedure to design doubly-resonant PhC cavities
by using evolutionary optimization algorithms, which can be applied to a vast
selection of cavity geometries and material platforms. The main figures of
merit, such as the quality factors of the cavity modes and their nonlinear
field overlap, were optimized targeting the main III-V semiconductors employed
in nonlinear photonics. These results could help the practical realization of
such devices for applications in high-efficiency second-order nonlinear
processes.

###### Acknowledgements.

We acknowledge useful discussions with M. Bollani, M. Clementi, C. De Angelis,
M. Galli, M. Liscidini.

## Funding

This research was partly supported from the Italian Ministry of Education and
Research (MIUR) through 2017 PRIN project 2017MP7F8F-004 “NOMEN”. M. Minkov
and S. Fan acknowledge the support of a MURI project from the U.S. Air Force
Office of Scientific Research (Grant No. FA9550-17-1-0002).

##  Appendix A Rotation of the nonlinear tensor

###  A.1 Euler angles and rotation matrices

Every rotation in the 3D space can be expressed in terms of the Euler angles
$(\phi,\theta,\psi)$, which define the rotation matrix
$\mathbf{R}(\phi,\theta,\psi)=\mathbf{A}(\psi)\mathbf{B}(\theta)\mathbf{C}(\phi)$
Goldstein  _et al._ (2002), where

$\displaystyle\scriptstyle\mathbf{A}\left(\psi\right)=\begin{pmatrix}\scriptstyle\cos\psi&\scriptstyle\sin\psi&\scriptstyle
0\\\ \scriptstyle\scriptstyle-\sin\psi&\scriptstyle\cos\psi&\scriptstyle 0\\\
\scriptstyle 0&\scriptstyle 0&\scriptstyle 1\end{pmatrix},$
$\displaystyle\scriptstyle\mathbf{B}\left(\theta\right)=\begin{pmatrix}\scriptstyle
1&\scriptstyle 0&\scriptstyle 0\scriptstyle\scriptstyle\\\ \scriptstyle
0&\scriptstyle\cos\theta&\scriptstyle\sin\theta\\\ \scriptstyle
0&\scriptstyle-\sin\theta&\scriptstyle\cos\theta\end{pmatrix},$
$\displaystyle\scriptstyle\mathbf{C}\left(\phi\right)=\begin{pmatrix}\scriptstyle\cos\phi&\scriptstyle\sin\phi&\scriptstyle
0\\\ \scriptstyle-\sin\phi&\scriptstyle\cos\phi&\scriptstyle 0\\\ \scriptstyle
0&\scriptstyle 0&\scriptstyle 1\end{pmatrix}.$

The coordinate transformation $(x,y,z)\to(x^{\prime},y^{\prime},z^{\prime})$
obtained as a result of the generic rotation $\mathbf{R}(\phi,\theta,\psi)$ is
schematically illustrated in Fig. 8.

Figure 8: Definition of the Euler angles describing an arbitrary 3D rotation,
for which the reference frame $(x,y,z)$ is transformed into
$(x^{\prime},y^{\prime},z^{\prime})$ by applying the rotation matrix
$\mathbf{R}(\phi,\theta,\psi)$.

###  A.2 Rotated nonlinear tensor

The reduced nonlinear tensor, $\bar{d}$ (see Ref. Boyd (2008) for a definition
in terms of $\chi^{(2)}$ tensor elements), is hereby reported for the
zincblende (wurtzite) in the reference frame corresponding to the
$[0,0,1]_{z}$ and $[1,1,1]_{z}$ ($[0,0,1]_{w}$) crystal growth directions.

$\displaystyle 2\bar{d}=\begin{cases}[0,0,1]_{z}:&\begin{pmatrix}\scriptstyle
0&\scriptstyle 0&\scriptstyle 0&\scriptstyle 1&\scriptstyle 0&\scriptstyle
0\\\ \scriptstyle 0&\scriptstyle 0&\scriptstyle 0&\scriptstyle 0&\scriptstyle
1&\scriptstyle 0\\\ \scriptstyle 0&\scriptstyle 0&\scriptstyle 0&\scriptstyle
0&\scriptstyle 0&\scriptstyle 1\end{pmatrix},\\\ \\\\[7.0pt]
[1,1,1]_{z}:&\begin{pmatrix}\scriptstyle\sqrt{\frac{2}{3}}&\scriptstyle-\sqrt{\frac{2}{3}}&\scriptstyle
0&\scriptstyle 0&\scriptstyle-\sqrt{\frac{1}{3}}&\scriptstyle 0\\\
\scriptstyle 0&\scriptstyle 0&\scriptstyle
0&\scriptstyle-\sqrt{\frac{1}{3}}&\scriptstyle
0&\scriptstyle-\sqrt{\frac{2}{3}}\\\
\scriptstyle-\sqrt{\frac{1}{3}}&\scriptstyle-\sqrt{\frac{1}{3}}&\scriptstyle\sqrt{\frac{4}{3}}&\scriptstyle
0&\scriptstyle 0&\scriptstyle 0\end{pmatrix},\\\ \end{cases}$ $\displaystyle
2\bar{d}=\begin{cases}[0,0,1]_{w}:&\begin{pmatrix}\scriptstyle 0&\scriptstyle
0&\scriptstyle 0&\scriptstyle 0&\scriptstyle 1&\scriptstyle 0\\\ \scriptstyle
0&\scriptstyle 0&\scriptstyle 0&\scriptstyle 1&\scriptstyle 0&\scriptstyle
0\\\ \scriptstyle 1&\scriptstyle 1&\scriptstyle-2&\scriptstyle 0&\scriptstyle
0&\scriptstyle 0\end{pmatrix}.\end{cases}$

##  Appendix B Particle swarm optimization

The particle swarm optimization (PSO) is an evolutionary algorithm employed
for the search of the minimum (or the maximum) of a given function
$f:A\to\mathbb{R}$ in the search space $A\subset\mathbb{R}^{n}$. Firstly, the
PSO generates a certain number, D, of trial solutions $\mathbf{x}\in A$, these
define the D-dimensional first generation. The subsequent generations are
calculated by evolving the previous solutions with the formula
$\mathbf{x}_{t}=\mathbf{x}_{t-1}+\mathbf{v}_{t}$, where $t$ labels the
generation, while $\mathbf{v}_{t}$ is evaluated as follows:

$\displaystyle\mathbf{v}_{t}=\omega\mathbf{v}_{t-1}+c_{1}\eta_{1}(\mathbf{p}_{t-1}-\mathbf{x}_{t-1})+c_{2}\eta_{2}(\mathbf{g}_{t-1}-\mathbf{x}_{t-1}),$

where $\omega$ is an inertial weight, $\eta_{1,2}$ are two random numbers
between $0$ and $1$, $\mathbf{p}$ is the best position explored by
$\mathbf{x}$, $\mathbf{g}$ is the global best position, while $c_{1}$ and
$c_{2}$ are the cognitive and social factors, respectively.

In the present work, the Lumerical-Ansys implementation of the PSO has been
used, in which standard values of $\omega$, $c_{1}$ and $c_{2}$ guarantee the
fast algorithm convergence in most cases.

##  Appendix C Second-harmonic conversion efficiency

The temporal coupled rate equations of two cavity modes, $a_{1}$ and $a_{2}$,
which mutually exchange energy through the coupling terms $\beta_{1}$ and
$\beta_{2}$ can be formally written as Lin  _et al._ (2016):

$\displaystyle\frac{da_{1}}{dt}=i\omega_{1}\left(1+\frac{i}{2Q_{1}}\right)a_{1}-i\omega_{1}\beta_{1}a_{1}^{*}a_{2}+\sqrt{\frac{2}{\tau_{1}}}s_{1+},$
(7)
$\displaystyle\frac{da_{2}}{dt}=i\omega_{1}\left(1+\frac{i}{2Q_{2}}\right)a_{2}-i\omega_{2}\beta_{2}a_{1}^{2},$
(8) $\displaystyle s_{1-}=\sqrt{\frac{\omega_{1}}{Q_{1}}}a_{1}-s_{1+},$ (9)
$\displaystyle s_{2-}=\sqrt{\frac{\omega_{2}}{Q_{2}}}a_{2},$ (10)

where $Q_{i}=\omega_{i}\tau_{i}/2$ ($i=1,2$) are the quality factors of the
cavity modes, $s_{1+}$ is the pump term for the first mode, while $s_{i-}$
describe the output rates. By assuming oscillating solutions, we set
$a_{1}=\tilde{a}_{1}e^{i\omega t}$, $a_{2}=\tilde{a}_{2}e^{i2\omega t}$,
$s_{1}=\tilde{s}_{1+}e^{i\omega t}$, and neglecting the down-conversion term
$i\omega_{1}\beta_{1}a_{1}^{*}a_{2}$ from Eq. 7, we get

$\frac{d\left(\tilde{a}_{1}e^{i\omega
t}\right)}{dt}=i\omega_{1}\left(1+\frac{i}{2Q_{1}}\right)\tilde{a}_{1}e^{i\omega
t}+\sqrt{\frac{2}{\tau_{1}}}\tilde{s}_{1+}e^{i\omega t}\,.$ (11)

In steady state the condition $d\tilde{a}_{1}/dt=0$ applies, hence Eq. 11
simplifies to

$i\omega\tilde{a}_{1}=i\omega_{1}\tilde{a}_{1}-\frac{1}{\tau_{1}}\tilde{a}_{1}+\sqrt{\frac{2}{\tau_{1}}}\tilde{s}_{1+}\,,$
(12)

leading to the following relation between $|\tilde{a}_{1}|^{2}$ and
$|\tilde{s}_{1+}|^{2}$

$|\tilde{a}_{1}|^{2}=\frac{2/\tau_{1}}{(\omega-\omega_{1})^{2}+1/\tau_{1}^{2}}\,|\tilde{s}_{1+}|^{2}\,.$
(13)

Following the same procedure with Eq. 8, we obtain

$\frac{d\left(\tilde{a}_{2}e^{i2\omega
t}\right)}{dt}=i\omega_{2}\left(1+\frac{i}{2Q_{2}}\right)\tilde{a}_{2}e^{i2\omega
t}-i\omega_{2}\beta_{2}\left(\tilde{a}_{1}e^{i\omega t}\right)^{2}\,.$ (14)

As in the previous case, steady state solutions can be found by imposing the
condition $d\tilde{a}_{2}/dt=0$, thus obtaining

$i2\omega\tilde{a}_{2}=i\omega_{2}\tilde{a}_{2}-\frac{1}{\tau_{2}}\tilde{a}_{2}-i\omega_{2}\beta_{2}\tilde{a}_{1}^{2}\,,$
(15)

which leads to the expression that links $|\tilde{a}_{2}|^{2}$ with
$|\tilde{a}_{1}|^{4}$:

$|\tilde{a}_{2}|^{2}=\frac{\omega_{2}^{2}\,|\beta_{2}|^{2}}{(2\omega-\omega_{2})^{2}+1/\tau_{2}^{2}}\,|\tilde{a}_{1}|^{4}\,.$
(16)

Finally, we can calculate the conversion efficiency as
$P_{o}/P_{i}^{2}=|\tilde{s}_{2-}|^{2}/|\tilde{s}_{1+}|^{4}$, in which
$|\tilde{s}_{2-}|^{2}$ is straightforwardly calculated from Eq. 10 as

$|\tilde{s}_{2-}|^{2}=\frac{\omega_{2}}{Q_{2}}|\tilde{a}_{2}|^{2}=\frac{2}{\tau_{2}}|\tilde{a}_{2}|^{2}\,,$
(17)

in which we can insert Eqs. 13 and 16 after imposing the resonant conditions,
i.e., $2\omega=\omega_{2}$ and $\omega=\omega_{1}$, from which we get

$|\tilde{s}_{2-}|^{2}=8\omega_{2}^{2}|\beta_{2}|^{2}\tau_{2}\tau_{1}^{2}||\tilde{s}_{1+}|^{4}=\frac{128}{\omega_{1}}|\beta_{2}|^{2}Q_{1}^{2}Q_{2}|\tilde{s}_{1+}|^{4}\,.$
(18)

Finally, given the expression of $\beta_{2}$

$\beta_{2}=\frac{1}{4}\frac{\int
d\mathbf{r}\;\chi^{(2)}\bar{\chi}_{ijk}E_{2\omega i}^{*}E_{\omega j}E_{\omega
k}}{\left(\int
d\mathbf{r}\,\varepsilon_{0}\varepsilon_{1}|\mathbf{E}_{\omega}|^{2}\right)\left(\int
d\mathbf{r}\,\varepsilon_{0}\varepsilon_{2}|\mathbf{E}_{2\omega}|^{2}\right)^{\frac{1}{2}}}\,,$
(19)

a dimensionless $\bar{\beta}$ can be defined as

$\bar{\beta}=\frac{\int d\mathbf{r}\;\bar{\chi}_{ijk}E_{2\omega
i}^{*}E_{\omega j}E_{\omega k}}{\left(\int
d\mathbf{r}\,\varepsilon_{1}|\mathbf{E}_{\omega}|^{2}\right)\left(\int
d\mathbf{r}\,\varepsilon_{2}|\mathbf{E}_{2\omega}|^{2}\right)^{\frac{1}{2}}}\,,$
(20)

where $\varepsilon_{1}=\varepsilon(\omega_{1})$ and
$\varepsilon_{2}=\varepsilon(\omega_{2})$. The conversion efficiency can then
be obtained from Eq. 18, in terms of $\bar{\beta}$, $Q_{1}$ and $Q_{2}$, as

$\frac{P_{o}}{P_{i}^{2}}=\frac{8}{\omega_{1}}\left(\frac{\chi^{(2)}}{\sqrt{\varepsilon_{0}\lambda_{1}^{3}}}\right)^{2}|\bar{\beta}|^{2}Q_{1}^{2}Q_{2}\,,$
(21)

which reproduces Eq. (6) in the main text, and it is consistent with the
expression already reported in the literature Lin  _et al._ (2016).

## References

  * Boyd (2008) R. W. Boyd,  _Nonlinear Optics, Third Edition_ , 3rd ed. (Academic Press, Inc., USA, 2008).
  * Lin  _et al._ (2016) Z. Lin, X. Liang, M. Lončar, S. G. Johnson, and A. W. Rodriguez, Optica 3, 233 (2016).
  * Rivoire  _et al._ (2009) K. Rivoire, Z. Lin, F. Hatami, W. T. Masselink, and J. Vučković, Opt. Express 17, 22609 (2009).
  * Buckley  _et al._ (2014) S. Buckley, M. Radulaski, J. Petykiewicz, K. G. Lagoudakis, J.-H. Kang, M. Brongersma, K. Biermann, and J. Vučković, ACS Photonics 1, 516 (2014), https://doi.org/10.1021/ph500054u .
  * Liscidini and Claudio Andreani (2006) M. Liscidini and L. Claudio Andreani, Phys. Rev. E 73, 016613 (22006).
  * Minkov  _et al._ (2019) M. Minkov, D. Gerace, and S. Fan, Optica 6, 1039 (2019).
  * Wang  _et al._ (2020) J. Wang, M. Clementi, M. Minkov, A. Barone, J.-F. Carlin, N. Grandjean, D. Gerace, S. Fan, M. Galli, and R. Houdré, Optica 7, 1126 (2020).
  * Andreani and Gerace (2006) L. C. Andreani and D. Gerace, Phys. Rev. B 73, 235114 (2006).
  * Hsu  _et al._ (2016) C. W. Hsu, B. Zhen, A. Stone, J. Joannopoulos, and M. Soljac̆ić, Nature Reviews Materials 1, 16048 (2016).
  * Robinson and Rahmat-Samii (2004) J. Robinson and Y. Rahmat-Samii, IEEE Transactions on Antennas and Propagation 52, 397 (2004).
  * Taflove and Hagness (2005) A. Taflove and S. C. Hagness,  _Computational electrodynamics: the finite-difference time-domain method_ , 3rd ed. (Artech House, Norwood, 2005).
  * Ochiai and Sakoda (2001) T. Ochiai and K. Sakoda, Phys. Rev. B 64, 045108 (2001).
  * Minkov  _et al._ (2020) M. Minkov, I. A. D. Williamson, L. C. Andreani, D. Gerace, B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, ACS Photonics 7, 1729 (2020), https://doi.org/10.1021/acsphotonics.0c00327 .
  * Ge  _et al._ (2018) X. Ge, M. Minkov, S. Fan, X. Li, and W. Zhou, Applied Physics Letters 112, 141105 (2018), https://doi.org/10.1063/1.5026433 .
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Helt  _et al._ (2012) L. G. Helt, M. Liscidini, and J. E. Sipe, J. Opt. Soc. Am. B 29, 2199 (2012).
  * Rivoire  _et al._ (2011) K. Rivoire, S. Buckley, F. Hatami, and J. Vučković, Applied Physics Letters 98, 263113 (2011).
  * Duchesne  _et al._ (2011) D. Duchesne, K. A. Rutkowska, M. Volatier, F. Légaré, S. Delprat, M. Chaker, D. Modotto, A. Locatelli, C. D. Angelis, M. Sorel, D. N. Christodoulides, G. Salamo, R. Arès, V. Aimez, and R. Morandotti, Opt. Express 19, 12408 (2011).
  * Sanford  _et al._ (2005) N. A. Sanford, A. V. Davydov, D. V. Tsvetkov, A. V. Dmitriev, S. Keller, U. K. Mishra, S. P. DenBaars, S. S. Park, J. Y. Han, and R. J. Molnar, Journal of Applied Physics 97, 053512 (2005), https://doi.org/10.1063/1.1852695 .
  * Goldstein  _et al._ (2002) H. Goldstein, C. Poole, and J. Safko,  _Classical Mechanics_ (Addison Wesley, 2002).