1. Introduction

High-quality cavities are useful in nanophotonics for a wide range of applications: sensors, filters, lasers, switches, etc. Recently, a highly efficient cavity mechanism was introduced [1], and experimentally demonstrated [2]. The simple geometry involved creates a host of possibilities for wavelength-scale, high-quality cavity devices. We therefore employed this cavity structure as the basis for a practical modeling exercise, with the added value of studying a hybrid material device.

Hybrid material systems, which join two or more different materials, aim to combine functions that are difficult in one or the other material. Silicon photonics e.g. has already been established as a very fruitful platform, as it combines strong confinement, due to a large refractive index, and low fabrication costs, when large scale lithography systems are employed [3]. However, active functionalities in Si remain a challenging problem, as the indirect bandgap lowers electron transition probabilities. On the other hand, III–V materials, such as GaAs and InP, are very effective for amplification and lasing purposes. Therefore much attention has been paid to find effective ways to combine Si and III–V material systems [4].

In this paper we introduce a design which applies the gradual photonic crystal cavity concept to a hybrid material context. We study the cavity modes of this design with multiple modeling approaches. The simulation results elucidate the mechanisms at play in the cavity. In addition they provide a challenging 3D modeling tool comparison, which was created in the context of COST Action MP 0702.

Experimental results and numerical results reported in the literature often vary by several orders of magnitude. Typically model errors (i.e. neglected roughness, material parameters or geometrical imperfections) are blamed for the experimental-numerical mismatch. However, in this study we show that the mismatch is likely due to insufficient numerical convergence. This situation is especially common when high-Q cavities are involved.

The paper is structured as follows. In section 2 we introduce the device structure. Section 3 discusses the simulation methods. The main part is section 4 which presents and interprets the simulation results in various situations.

2. Hybrid cavity structure definition

The design consists of a gradual cavity in InP, coupled to an (input or output) waveguide in Si (Figs. 1(a)–1(d)). The cavity is in a III–V material, as this is the natural choice for an active component. The waveguide is in silicon, as this could constitute the backbone of a larger, multifunctional photonic chip with relatively long silicon waveguides, bends, splitters etc. The material around the III–V and Si is a benzocyclobuten-based polymer (BCB), because this material has been proven effective in bonding Si wafers with InP dies [5]. The inclusion of BCB all around the structure also renders the cavity more symmetric, which often improves the quality factor of embedded cavities.

 

Fig. 1 Geometry of the simulated device, with parameters and axes indicated. (a) xz-side view. (b) yz-side view. (c) Top view. (d) Perspective view with only InP and Si shown.

Download Full Size | PDF

The InP stack cavity is created by gradually decreasing the size of a number of InP veins towards the center. We keep the central positions of the veins constant (always one period apart), but we vary their longitudinal thickness (along x), so the thinnest veins are at the center. An important parameter is N_cav: the number of modulated veins on each side of the center, so in total there are 2N_cav veins in the cavity with slightly smaller widths. The center of the cavity is in BCB, in between the first two veins, which have the same width. The modulated cavity vein widths are

As index values we employ n = 3.46 for Si, n = 1.45 for silicon oxide (the substrate material for the Si waveguide), n = 3.17 for InP and n = 1.54 for BCB.

The mechanism of confinement is explained as follows [1]. In the center of the InP structure a Bloch mode is guided. However, when this mode propagates from the center to the side, the veins become thicker. This increase of high-index material generally means that the dispersion of the mode lowers in frequency. Eventually, at the operating frequency of the cavity, this mode becomes cut-off in the mirror sections on the side. The mode, which was propagating in the center, thus encounters a band-gap, energy can be reflected, or scattered into radiation. However, because the vein thickness change is gradual, the adjustment of the Bloch mode is very slow, leading to a substantially high probability that the forward-propagating Bloch mode is reflected into the backward-propagating mode. The latter leads to very high reflections, and, when this process happens on both sides of the center, to high-quality cavity modes. More gradual cavities will generally lead to higher quality confinement, but the mode volume will increase.

3. Modeling methods

This section introduces the four employed simulation methods. All methods simulate the devices in full 3D, and can determine the quality factor of cavity modes, however the numerical effort varies strongly.

3.1. Finite-difference time-domain (FDTD)

The finite-difference time-domain method offers a large flexibility and is widely used for electromagnetic calculations [6]. Fields are calculated on a rectangular spatial grid, in combination with a time-stepping procedure. We employ the freely available MEEP software package, which includes a subpixel averaging technique to decrease the fineness of the required calculation grid [7]. Even with this technique we employed 20 pixels for a period (grid size about 18nm), leading to large 3D simulations for which we employed the supercomputer infrastructure at Ghent University.

FDTD is popular because of its large flexibility and extendability, many types of calculations can be performed through fairly robust algorithms. To determine quality factors we employ the ‘harminv’ feature available in MEEP, which employs the filter diagonalization method to extract decay patterns in time series data [8].

3.2. Bidirectional eigenmode propagation (BEP)

This technique, which belongs to the class of frequency-domain modal methods, assumes that the analysed structure is composed of longitudinally uniform waveguide sections. The refractive index profile in any section is a function of transversal coordinates only, i.e., n = n(y, z) in Fig. 1. Then the solution of the Maxwell equations is achieved in two steps [9]: (1) the total field in each section is expanded into the orthogonal set of waveguide modes, and (2) the field continuity conditions are used at the section interfaces.

We have used our own implementation of full-vector BEP for 3D structures. The waveguide modes are searched by means of finite-element commercially available software COMSOL Multiphysics. Modes of homogeneous sections are calculated analytically. The propagation technique follows the standard formulation which can be found, e.g., in [10, 11] and therefore is not repeated here. Note that the waveguide sections and their interfaces are described by means of scattering matrices [12] where the interface matrices (the reflection and transmission matrices) are determined from overlap integrals of modal fields.

In order to find resonance frequencies and quality factors (Q), we adopted the algorithm which uses eigenvalues of the cavity reflectivity matrix (reflectivity matrix of symmetric half of the whole cavity) similar to that described in [13]. Note that the algorithm has been found to provide more reliable results than the standard approach based on locating the maximum and bandwidth of a resonance curve.

3.3. Finite-element method (FEM)

The finite-element method is a numerical method for solving Maxwell’s equations in various formulations, or other partial differential equations. Main advantages of finite-element methods are the exact representation of complex geometrical features present in real-world applications, the very good convergence properties which allow for very accurate results in relatively short computation times, and the stability and general applicability of the method to different simulation settings. We have used the FEM software JCMsuite developed by JCMwave and Zuse Institute Berlin. This program package contains finite-element method based solvers for time-harmonic Maxwell eigenvalue and for Maxwell scattering problems. The method is based on higher order vectorial elements, adaptive unstructured grids, and on a rigorous treatment of transparent boundaries [14].

Here we have used the resonance mode solver: given the geometrical setup of the nanocavity, we compute an electric field distribution E and a complex eigenfrequency ω which satisfy Maxwell’s time-harmonic wave equation without sources. Transparent boundary conditions take into account the specific geometry of the problem where waveguides are modelled to extend to infinity in the exterior domain. When the eigenmode (E, ω) is computed, the respective Q-factor is deduced from the real and imaginary parts of the complex eigenfrequency, Q = −ℜ(ω)/2ℑ(ω), the resonance wavelength λ₀ is given by λ₀ = 2πc/ℜ(ω), with c the speed of light.

In previous studies on high-Q cavities the very good convergence properties of the method have been demonstrated, and the obtained values (λ₀, Q) have been confirmed by simulations of transmission spectra of light incident to the cavities [15, 16]. Accurate results on high-Q cavity properties are obtained in computation times between few seconds and few minutes on standard workstations [16]. Further, good agreement with experimental data has been shown [17].

3.4. Aperiodic rigorous coupled wave analysis (aRCWA)

This is a Fourier expansion scheme implemented in our in-house robust 3D tool which effectively combines both 2D mode solver with advanced schemes of scattering matrix formalism [18]. The aRCWA method is based on the 2D periodic grating-based RCWA method [19], currently with all efficient up-to-date techniques already included (as the efficient Fourier factorization techniques). In order to proceed towards the 2D mode solver, we followed the same idea as for the 1D solver case [20], i.e. the isolating boundary conditions applied on the boundaries of a 2D period have been utilized again. For these artificial absorbing (isolating) boundaries, both complex coordinate transform [21] as well as uniaxial perfectly matched layers [20] have been implemented. Additionally, in order to achieve faster convergence, we have incorporated the 2D adaptive spatial resolution technique [22, 23].

However, for the full 3D case, we have to cope with the problem that the number of expansion terms in the Fourier modal method rapidly increases; the total number of expansion terms is namely twice the product of the particular expansion numbers (originally diffraction orders) in each of the two transversal directions. To reduce this problem at least partially, we have fully utilized structural symmetries of the modeled objects as well [24].

For the application in this paper, resonance wavelengths and Q factors are calculated from the eigenvalues of the cavity reflectivity matrix [13], basically with the same procedure as in the case of the BEP technique described above.

4. Simulation results

A central parameter under investigation is the quality factor Q of the cavity. In the following sections we examine how the total Q of the structure changes when cavity and coupling properties are adjusted.

4.1. Cavity without waveguide

The most basic situation for the cavity is when there is no waveguide and no substrate nearby, we examine the InP stack structure completely surrounded by BCB. This provides us with what can be called the ‘intrinsic’ Q_int of the device. In later sections the ‘total’ or coupled quality factor will be adjusted mainly by coupling to the waveguide.

We simulated Q_int and the normalized resonance frequency f = period/λ = a/λ (Fig. 2). The trend in Fig. 2(a) clearly shows a strongly increasing Q_int in function of a larger N_cav. Thus a more gentle and spatially distributed cavity section leads to diminished radiation losses, and to a more perfect reflection into the central Bloch mode. As N_cav increases the rate of increase of Q diminishes gradually. We believe that for small N_cav the structure is quite sensitive to increasing N_cav, leading to a rapid decline of radiation losses and increasing Q. Therefore, at higher N_cav a slightly slower exponential increase is observed. Moreover, at very high Q values secondary loss mechanisms start to influence the cavity, losses through the N_mir mirror veins being the most probable candidate. However, leakage in perpendicular directions is only slightly affected by additional cavity blocks, therefore when leakage in perpendicular directions gets significant (for high N_cav), addition of further cavity blocks can only slightly increase Q.

 

Fig. 2 Simulation results for the cavity without nearby waveguide. (a) Quality factor Q and (b) normalized resonance frequency a/λ versus N_cav.

Download Full Size | PDF

A picture of the cavity mode is shown in Fig. 3. Confinement by the bandgap effect in the horizontal x-direction is clear by a gradual decay of the Bloch mode amplitude towards the mirroring edges. The quality factor of the cavity is mainly limited by the reflection losses, which decrease as the mirror becomes more gradual.

 

Fig. 3 Electric field of the cavity mode (lower part). Shown is the E_y-component in a horizontal plane (parallel to xy), through the cavity center. The upper part depicts the cavity geometry.

Download Full Size | PDF

The normalized resonance frequency (Fig. 2(b)) increases slightly in function of N_cav and levels off towards the upper limit of the Bloch mode frequency in the center of the cavity. The increasing and saturating trend can be understood: A larger N_cav means more veins in the cavity with a smaller thickness close to the minimum central thickness. So the resulting frequency will end up closer to the Bloch mode frequency of this thinnest central section.

The agreement between the four numerical methods in Figs. 2(a) and 2(b) is reasonably good. The correct trend for Fig. 2(a) is expected to be smooth, so deviations around the central trend are likely resulting from numerical artifacts, such as imperfect boundary conditions. Note that quantitative differences for Q are much bigger than for the resonance frequency. This is because the calculated Q is strongly affected by numerical inaccuracies (which are inevitable in any numerical technique) and by proper choice of numerical parameters. We will mention this problem in more details in the next paragraph.

4.2. Cavity with waveguide

We found that the change of Q and resonance frequency is negligible when the SiO_x (silicon oxide) substrate is included at a distance of 1 μm from the cavity (not shown). Indeed, the similar index (1.54 for BCB vs 1.45 for SiO_x) and the mode confinement indicate a very slight disturbance by the substrate. However, the situation changes drastically when a Si waveguide is included, but the change depends strongly on the waveguide properties.

We first examine the influence of the lateral Si waveguide width Si_y on the total Q, for various constant values of N_cav, see Fig. 4. All curves present a similar trend: Q is mostly constant, but shows a dip around a certain value of Si_y = 0.35μm. This pattern is explained by phase matching. At all calculated values of Si_y there is a guided mode available in the waveguide, but only at the particular width of Si_y = 0.35μm one obtains phase matching between the waveguide mode and the (central) Bloch mode in the cavity.

 

Fig. 4 Quality factor Q versus Si waveguide width Si_y, at various constant values of N_cav.

Download Full Size | PDF

The phase matching argument is proven more quantitatively by Bloch mode calculations using MPB (Fig. 5). With the software MPB [25] we can determine the propagation constants of relevant features in the device: the Si waveguide and specific periods of the InP cavity. Figure 5 shows the dispersion for Si waveguide widths just below, near and above the minimum Q-point. In addition, we plot the Bloch mode dispersion of the central cavity period (w_cav = 0.15μm), and the mirror periods (w_mir = 0.2μm). We see that near the minimum Q point (Fig. 5, middle) the Si waveguide propagation constant is exactly in the zone where we expect the (major part of) cavity Bloch section dispersions. Indeed, the cavity mode will have a large proportion of k_x-components at the edge of the Brillouin zone, as reflection is provided by the bandgap-effect in the mirror sections.

 

Fig. 5 Bloch and waveguide dispersion relations. The horizontal dotted line indicates the cavity resonance frequency. The long-dashed line shows the (folded) Si waveguide dispersion, for the indicated widths. The solid lines are identical in the three graphs, and provide the Bloch mode propagation constant with period as in the center of the cavity (upper mode), and with period as in the mirror sections (lower mode), respectively.

Download Full Size | PDF

In this way the Si waveguide propagation constant matches the Bloch propagation constant of the central cavity section around Si_y = 0.35μm, leading to a strong coupling or phase matching, just as in conventional directional couplers. Note that the crossing with the folded dispersion at Si_y = 0.4μm does not lead to coupling: indeed for the Si waveguide the folding effect is artificial (it has no periodic corrugation).

This minimum Q point provides an interesting operating point for active devices [26]. The two loss channels (radiation and coupling) provide a total Q_tot given by: 1/Q_tot = 1/Q_coup + 1/Q_int. For lasers good coupling to a nearby waveguide is achieved when Q_int/Q_coup is large [27], which is indeed achieved at the dip.

Next we show Q versus N_cav, at a constant value of Si_y (Fig. 6). This provides another image of the previous phase matching process. If Si_y is different than the value 0.35μm, Q_tot is nearly identical to Q_int, as coupling is negligible. For Si_y = 0.35μm, we see a saturation of Q_tot versus N_cav, because Q_tot cannot exceed Q_coup. Q_coup is indeed determined by Si_y, and is independent of N_cav.

 

Fig. 6 Cavity with waveguide. Q versus N_cav for different values of Si_y.

Download Full Size | PDF

The qualitative agreement between the modeling methods is reasonably good in Figs. 4 and 6, however the exact value of Q_tot varies, especially at the minimum (i.e. Si_y = 0.35μm in Figs. 4 and 6). This may be surprising because the used numerical techniques are reliable and provide good agreement in many other situations. Without offering details, we checked for all techniques that numerical parameters were chosen judiciously. Differences in Q_tot are attributed to the resonant behavior of the structure: any inevitable numerical errors are amplified in the resonance, namely for high-Q cases. The worst disagreement observed in the minimum can be explained by the high sensitivity of the phase matching process, where slight differences, e.g. in spatial resolution of the techniques, can result in large mismatches. In addition, the effectiveness of transparent boundary conditions (represented with perfectly matched layers) can influence the results strongly, namely in the FDTD and FEM methods [16]. This holds true in both cases, for strong coupling when a guided mode in the silicon waveguide should be crossing the computational domain boundaries to the exterior domain without artificial reflections, and for the case of weak coupling, when the low leakage of light from the cavity to the various material regions in the exterior domain has to be computed at high accuracy.

We highlight another facet of the phase matching process in Fig. 7, where we vary the BCB buffer thickness between the cavity and the waveguide. Again we see for all curves a dip of Q_cav around the same value of Si_y = 0.35μm. However, as BCB_z increases, the dip becomes smaller. This is in accordance with mode coupling processes: a larger distance between waveguide and cavity results in a reduced mode overlap, so that the maximum amount of coupling diminishes (Q_coup increases). Again the main discrepancy between the simulation methods occurs at the maximum coupling point.

 

Fig. 7 Cavity with waveguide. Quality factor Q versus waveguide thickness Si_y for different values of the buffer thickness BCB_z. Here, N_cav = 5, other parameters are as in Fig. 4. To be compared with Fig. 4 (BCB_z = 1.0μm), case N_cav = 5.

Download Full Size | PDF

4.3. Offset influence

Finally we show results when the waveguide is not directly underneath the cavity, we include the possibility of an offset in the y-direction, see Fig. 8. Note that still the waveguide and the cavity remain parallel. Increasing the offset between waveguide and cavity clearly increases the Q-factor, until it saturates to the intrinsic value Q_int. This is consistent again with the phase matching picture: a larger offset leads to a diminished coupling, so that eventually the separate cavity quality factor is recovered.

 

Fig. 8 Cavity Q as a function of waveguide offset in the y-direction, for N_cav = 5.

Download Full Size | PDF

Remark that the minimum Q-point is initially fairly robust for small offsets, up to 200 nm variation provides only small differences. This reduces the constraints on fabrication tolerances.

The simulation methods provide qualitatively similar trends in function of the waveguide offset, thus showing the same physical process. The precise minimum and maximum Q values differ, which are the limiting cases as already discussed before. Detailed comparison of convergence behavior of the various used methods and comparison to results from further methods will be helpful to clarify the causes of the quantitative disagreement for these situations.

5. Conclusions

Using multiple rigorous simulation methods we examine a hybrid high-Q cavity device, combining a possibly active III–V material with a waveguide in silicon. The cavity operation stems from a gradual decrease of the near-periodic material sections, leading to a highly efficient 3D confinement.

Our results show that the coupling between the cavity and the waveguide stems from a sensitive phase-matching process, which needs to be controlled to achieve an optimum design. All simulations methods provide the same qualitative trends with respect to the quality factor in various situations. However, the main discrepancies show up when strong coupling to the waveguide is involved, as this involves a delicate phase-mismatch and the need for more stringent boundary conditions. It is seen that rigorous and accurate calculation of such 3D resonant devices still remains a challenging problem which should be further investigated. In addition, if larger quality factors are required, the tapering of other parameters could be considered. For now however, the fabrication of the proposed device seems possible with current technologies, and the design can have direct applications for sensing and lasing.