Doubly-Resonant Photonic Crystal Cavities for Efficient Second-Harmonic Generation in III–V Semiconductors

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 via localized defects in a photonic crystal structure. We built on this approach by applying an evolutionary optimization algorithm in connection with Maxwell equations solvers; the proposed design recipe can be applied to any material platform. We explicitly calculated the second-harmonic generation efficiency for doubly-resonant photonic crystal cavity designs in typical III–V semiconductor materials, such as GaN and AlGaAs, while targeting a fundamental harmonic at telecom wavelengths and fully accounting for the tensor nature of the respective nonlinear susceptibilities. These results may stimulate the realization of small footprint photonic nanostructures in leading semiconductor material platforms to achieve unprecedented nonlinear efficiencies.

1. Introduction

Nonlinear optical processes are notoriously poorly efficient, due to small higher-order nonlinear susceptibilities $\left({\chi }^{\left(n\right)}\right)$ that mediate such processes in conventional materials, such as semiconductors or insulators, in their transparency spectral range [1]. A way out of 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 to 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, thereby enhancing the efficiency of nonlinear effects [2]. 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 [3,4], where the radiation is trapped only for the first-harmonic (FH) frequency, ${\omega }_{\mathrm{FH}}=\omega$, and doubly-resonant PhC cavities [2,5,6], 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 made to fulfill a doubly-resonant condition that could boost SH generation efficiency, until very recently [7]. The main obstacle comes from the difficulty of designing PhCs with photonic band gaps around the FH and SH frequencies, which would allow obtaining two doubly-resonant cavity modes [8]. Recently, a novel approach based on the so-called bound-states in the continuum (BICs) that overcomes the aforementioned obstacle has been proposed [6]. This strategy has later been shown to be practically effective by realizing doubly-resonant SH generation in a III–V semiconductor PhC cavity [7]. The BICs are particular states that stay confined despite the fact that they lie above the light line—i.e., their dispersion falls within the continuum of radiative modes [9]. This feature is indispensable for achieving good temporal confinement of the SH cavity mode, even in the absence of a photonic bandgap [6].

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 optimization (PSO) algorithm combined with a numerical solver for Maxwell equations [10]; in particular, we use the finite-difference time-domain (FDTD) algorithm here [11], 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}^{\left(2\right)}$ [2], 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 ($\overline{\beta }$) between the FH and SH modes, taking into account the ${\chi }^{\left(2\right)}$ tensor. The latter can be regarded as a generalization of the phase-matching condition for the confined systems. The factor $\overline{\beta }$ is directly related to the estimation of the SHG efficiency, which can be approximated to scale as $Q_1^2{Q}_2{\left|\overline{\beta }\right|}^2$, $Q_1$ and $Q_2$ being the FH and SH modes quality factors, respectively [2,5]. 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 $\mathrm{III}-\mathrm{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 }_{\mathrm{FH}}\approx$ 1550 $\mathrm{n}$$\mathrm{m}$ and ${\lambda }_{\mathrm{SH}}\approx$ 775 $\mathrm{n}$$\mathrm{m}$. It was found that the $[0,0,1]$ growth direction of the zincblende crystalline materials leads 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, leads to an isotropic generation efficiency that does not depend on the PhC’s orientation against the unitary crystal cell.

The manuscript is organized as follows. In Section 2 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 Section 3 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 Section 4 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.

2. Evolutionary Optimization of Doubly Resonant Photonic Crystals

2.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 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 parts of complex eigenmodes [8]. Here, we consider only PhC slabs with holes carved in the dielectric slabs 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, the (dimensionless) radius $r/a$ and the slab thickness $d/a$. GME starts from a decomposition of the electromagnetic field on the 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 [8]. Moreover, perturbation theory allows one 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 [12]. This approach turns out to be computationally very efficient and quantitatively accurate to within a few percent [13], albeit it has approximate and polarization-dependent accuracy. Hence, it has been applied for rapid and coarse scanning of the parameters, and for 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 [11]. 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 with 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 by Lumerical-Ansys. Moreover, a particle swarm optimization (PSO) algorithm has been exploited in connection with the 3D-FDTD solver in order to finely tune the PhC parameters [10]. Such a routine allows one to efficiently minimize 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 is provided in Appendix B. In the present case, the PSO was chosen, since it allows one to reach convergence within a few evolutionary steps when the starting point is not too far from the target objectives, as is most often the case when targeting the doubly-resonant condition.

2.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 [6]. This approach has also been recently shown to experimentally enhance SHG efficiency in GaN PhC cavities [7]. 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 summarized as follows:

• Calculate the photonic bandstructure of the PhC slab with GME.
• 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.
• Check that FH and SH modes have non-vanishing overlap integral $\overline{\beta }$, depending on the ${\chi }^{\left(2\right)}$ tensor of the given material.
• Vary the hole radius and the slab thickness in order to approximately match the doubly-resonant condition ${\omega }_2\approx 2{\omega }_1$.
• 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.
• Spatially confine the FH and SH mode by suitably designing a heterostructure PhC cavity.

2.3. Doubly-Resonant Photonic Crystal Slabs

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) [8,13]. These are shown for the specific case of a GaN PhC slab in Figure 2 and Figure 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 Figure 2b,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$, and for SH modes $n=2.31$. These values correspond to the GaN refractive indices at the target wavelengths ${\lambda }_{\mathrm{FH}}=$ 1550 $\mathrm{n}$$\mathrm{m}$ and ${\lambda }_{\mathrm{SH}}=$ 775 $\mathrm{n}$$\mathrm{m}$, 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 Figure 2a and Figure 3a, and the imaginary part is directly compared only for the target BIC mode in the SH case, Figure 3b. 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 $\mathsf{\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, and 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 [6]. 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 $\mathsf{\Gamma }$, which is a BIC in the $\mathsf{\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.

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 optimization algorithm vary the radius of the holes and the slab thickness in the search-space $A=[r_{\min },r_{\max }]\times [d_{\min },d_{\max }]$, where A is a small neighbourhood around the best GME parameters; while the figure of merit to be minimized has been defined as:

$$\mathrm{FOM}=\left|\frac{{\omega }_2-2{\omega }_1}{{\omega }_2}\right|,$$

(1)

where ${\omega }_1$ and ${\omega }_2$ are the FH and SH mode frequency, respectively. We imposed the condition $\mathrm{FOM}\le 1\%$ to be fulfilled, such that the doubly-resonant condition is achieved with a given tolerance.

To show the generality of the approach, the very same procedure has been applied to both GaN and $\mathrm{A}{\mathrm{l}}_{\mathrm{x}}\mathrm{G}{\mathrm{a}}_{1-\mathrm{x}}\mathrm{A}\mathrm{s}$ 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 summarized in dimensionless units in

Table 1. Main design parameters of the GaN and $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ photonic crystal (PhC) slabs obtained with the PSO algorithm applied to a 3D-FDTD solver.

Material	${\mathit{n}}_{\mathbf{FH}}$	${\mathit{n}}_{\mathbf{SH}}$	$\mathit{d}/\mathit{a}$	$\mathit{r}/\mathit{a}$	${\mathit{\omega }}_1\mathit{a}/\left(2\mathit{\pi }\mathit{c}\right)$	${\mathit{\omega }}_2\mathit{a}/\left(2\mathit{\pi }\mathit{c}\right)$	a (nm)	d (nm)	r (nm)
GaN	$2.28$	$2.31$	$0.337$	$0.206$	$0.414$	$0.832$	$642.3$	$216.4$	$132.2$
$\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$	$3.23$	$3.47$	$0.24$	$0.175$	$0.3356$	$0.6650$	$520.1$	$124.8$	$91.2$

. They can also easily be converted to dimensional units by imposing the target wavelength, e.g., ${\lambda }_{\mathrm{FH}}=$ 1550 $\mathrm{n}$$\mathrm{m}$ for the FH mode, yielding $a={\tilde{\omega }}_1{\lambda }_{\mathrm{FH}}$, ${\tilde{\omega }}_1$ being the FH dimensionless frequency, from which the lattice constant a, and all other quantities, can be calculated.

3. 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 [14]. 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 Figure 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 [6]. For consistency, we set $N_c=6$, $N_t=4$ and $N_o=10$, for simulations on both GaN and $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ platforms. The slab thickness and the hole radius in the core region are inherited from the PhC slabs parameters, and 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.

We studied the properties of the same heterostructure cavity design realized in both GaN and $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ 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 one to fully recover the frequencies and quality factors of both FH and SH cavity modes. The outcomes are reported in

Table 2. Main design parameters of the GaN and $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ cavities, and the related doubly-resonant features obtained from 3D FDTD simulations.

Material	a (nm)	d (nm)	${\mathit{r}}_{\mathit{c}}$ (nm)	${\mathit{r}}_{\mathit{t}}$ (nm)	${\mathit{r}}_{\mathit{o}}$ (nm)	${\mathit{\lambda }}_{\mathbf{FH}}$	${\mathit{\lambda }}_{\mathbf{SH}}$	${\mathit{Q}}_{\mathbf{FH}}$	${\mathit{Q}}_{\mathbf{SH}}$
GaN	642	200	122	145	165	1503	753	$3.6\times {10}^4$	1100
$\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$	520	124	91	105	121	1523	759	$1.1\times {10}^5$	400

.

Moreover, from the FDTD simulations it is possible to retrieve the electromagnetic field profiles inside the cavities, and 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 Figure 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 [7].

4. 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, and their confinement volume. For second-order nonlinear processes, a nonlinear overlap between the FH and SH modes can generally be defined as [2]:

$$\overline{\beta }=\frac{\int d\mathbf{r}{\sum }_{ijk}{\overline{\chi }}_{ijk}^{\left(2\right)}{E}_{2\omega i}^{*}{E}_{\omega j}{E}_{\omega k}}{\left(\int d\mathbf{r}{\epsilon }_{\omega }\left(\mathbf{r}\right){\left|{\mathbf{E}}_{\omega }\right|}^2\right){\left(\int d\mathbf{r}{\epsilon }_{2\omega }\left(\mathbf{r}\right){\left|{\mathbf{E}}_{2\omega }\right|}^2\right)}^{1/2}}\sqrt{{\lambda }_{\mathrm{FH}}^3},$$

(2)

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

$$\begin{array}{cc}\hfill & \mathrm{Zincblende}:{\chi }^{\left(2\right)}={\chi }_{zxy}^{\left(2\right)}\hfill \\ \hfill & \mathrm{Wurtzite}:{\chi }^{\left(2\right)}={\chi }_{xxz}^{\left(2\right)}.\hfill \end{array}$$

(3)

The overlap factor in Equation (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 Equation (2) through the tensor elements ${\overline{\chi }}_{ijk}^{\left(2\right)}$. Following the Kleinmann symmetry [1], the non-vanishing tensor elements of the wurtzite (e.g., GaN) and zincblende (e.g., $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$) crystal structures grown along the $[0,0,1]$ direction are:

$$\begin{array}{cc}\hfill & \mathrm{Zincblende}:{\overline{\chi }}_{xyz}^{\left(2\right)}={\overline{\chi }}_{zxy}^{\left(2\right)}={\overline{\chi }}_{yzx}^{\left(2\right)}={\overline{\chi }}_{xzy}^{\left(2\right)}={\overline{\chi }}_{yxz}^{\left(2\right)}={\overline{\chi }}_{zyx}^{\left(2\right)}=1\hfill \\ \hfill & \mathrm{Wurtzite}:{\overline{\chi }}_{xxz}^{\left(2\right)}={\overline{\chi }}_{zxx}^{\left(2\right)}={\overline{\chi }}_{xzx}^{\left(2\right)}={\overline{\chi }}_{yyz}^{\left(2\right)}={\overline{\chi }}_{zyy}^{\left(2\right)}={\overline{\chi }}_{yzy}^{\left(2\right)}=1,{\overline{\chi }}_{zzz}^{\left(2\right)}=-2.\hfill \end{array}$$

(4)

We supposed ${\overline{\chi }}_{zzz}^{\left(2\right)}=-2$, this condition is not always experimentally verified but since the FH is mostly TE-polarized, the actual value of the ${\overline{\chi }}_{zzz}^{\left(2\right)}$ tensor element has a negligible effect on the beta-factor of Equation (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}(\varphi ,\theta ,\psi )$ in which $(\varphi ,\theta ,\psi )$ are the Euler angles (see Appendix A for details). For example, the transformation that allows one to switch from the $[0,0,1]$ direction to the $[1,1,1]$ growth direction is

$$\begin{array}{cc}\hfill & {\overline{\chi }}_{\alpha ,\beta ,\gamma }^{\left(2\right)[1,1,1]}=\sum _{a,b,c}{U}_{\alpha ,a}{U}_{\beta ,b}{U}_{\gamma ,c}{\overline{\chi }}_{a,b,c}^{\left(2\right)[0,0,1]},\hfill \\ \hfill & \mathbf{U}\doteq \mathbf{R}\left(\frac{3}{4}\pi ,{cos}^{-1}\left(\frac{1}{\sqrt{3}}\right),\frac{\pi }{2}\right).\hfill \end{array}$$

(5)

Here, we calculated the overlap integral (2) in the case of $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ 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 Figure 6a for both material platforms. A definition of the rotation angle is schematically given in Figure 6b.

The overlap factor, $|\overline{\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 $|\overline{\beta }{|}^2$ goes beyond ${10}^{-4}$ for any orientation angle, and in the $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ it is modulated and can vanish at specific orientation angles for the $[0,0,1]$ growth direction. On the other hand, the $|\overline{\beta }{|}^2$ does not depend on the relative orientation between the crystallographic directions and the PhC lattice for the $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ 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 }^{\left(2\right)}$ tensor elements, which is material dependent and is not explicitly taken into account in Equation (2). This additional parameter can play a crucial role in calculating the final nonlinear efficiency, which can be estimated with the following relation [2]:

$$\frac{P_o}{P_i^2}=\frac{8}{{\omega }_1}{\left(\frac{{\chi }^{\left(2\right)}}{\sqrt{{\epsilon }_0{\lambda }_{\mathrm{FH}}^3}}\right)}^2{\left|\overline{\beta }\right|}^2{Q}_{\mathrm{FH}}^2{Q}_{\mathrm{SH}},$$

(6)

where ${\epsilon }_0$ is the vacuum permittivity. A detailed derivation of Equation (6) is reported in Appendix C for completeness. We also notice that Equation (6) is valid under the undepleted-pump approximation, and considering perfect in- and out-coupling of the radiation for both FH and SH modes [6]. Notice that similar scaling relations can be found for other nonlinear generation processes, e.g., SPDC [16]. Notice also that this expression is explicitly obtained for narrow-band harmonic resonances, at variance with SHG in broadband photonic crystal waveguide structures [17]. We evaluated the conversion efficiencies for a wide range of $Q_{\mathrm{FH}}$ and $Q_{\mathrm{SH}}$. We set ${\chi }^{\left(2\right)}=$ 100 $\mathrm{p}$$\mathrm{m}$ ${\mathrm{V}}^{\mathrm{-1}}$ for $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$, and ${\chi }^{\left(2\right)}=$ 10 $\mathrm{p}$$\mathrm{m}$ ${\mathrm{V}}^{\mathrm{-1}}$ for GaN, respectively [18,19]. Moreover, we considered the axes orientations maximizing the nonlinear overlap factor according to the results in Figure 6a. The results are reported in Figure 7, which inverts the outcome of the previous Figure: in fact, the $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ doubly-resonant cavity (grown along [1,1,1]) is theoretically more efficient than the GaN one, despite the dimensionless $|\overline{\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 $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$, which is characterized by an higher ${\chi }_{\mathrm{eff}}^{\left(2\right)}$. 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 ${\mathrm{W}}^{\mathrm{-1}}$ for the $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$, and $2.9$ ${\mathrm{W}}^{\mathrm{-1}}$ for the GaN doubly-resonant PhC cavity, as summarized in

Table 3. Calculated nonlinear figures of merit for both the GaN and $\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$ doubly-resonant PhC cavities designed in this work.

Material	${\mathit{Q}}_{\mathbf{FH}}$	${\mathit{Q}}_{\mathbf{SH}}$	$|{\overline{\mathit{\beta }}}_{\mathbf{max}}{|}^2$	${\mathit{\chi }}^{\left(2\right)}$ (pm V${}^{-1}$)	${\mathit{P}}_{\mathit{o}}/{\mathit{P}}_{\mathit{i}}^2$ (W${}^{-1}$)
GaN	$3.6\times {10}^4$	1100	$1.0\times {10}^{-4}$	10	$2.9$
$\mathrm{A}{\mathrm{l}}_{0.3}\mathrm{G}{\mathrm{a}}_{0.7}\mathrm{A}\mathrm{s}$	$1.1\times {10}^5$	400	$1.2\times {10}^{-5}$	100	119

.

5. 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.