1. INTRODUCTION

Many important properties of photonic structures such as cavities [1] and waveguides [2,3] depend on the radiative losses that stem from coupling of energy into freely propagating optical modes that escape the system. The quality, or $Q$, factor of photonic resonators as well as the spontaneous emission (SE) $\beta$ factor in waveguides are important figures of merit in the analysis of nanolasers [4] and single-photon sources [5], for example, and these quantities depend sensitively on the radiative losses. In modeling such open photonic systems, the choice of boundary conditions (BCs) at the computational domain edges becomes crucial and may impact the results not just quantitatively but also qualitatively. Integral equation Green’s function formulations inherently adopt this openness [6], while for numerical techniques relying on finite-sized computational domains like the finite-difference time-domain method [7] and the finite element method [8], this is achieved using artificial absorbing boundaries, typically in the form of so-called perfectly matched layers (PMLs) [9].

In Fourier-based modal-expansion techniques [10,11], PMLs can be implemented using complex coordinate transforms [12]. The absorbing boundaries are implemented by mapping the real spatial coordinates into complex ones, which is straightforward to implement. However, it is unclear which complex coordinate transform to implement and why, and there have been no systematic studies on the influence of PML parameters and the size of the computational domain on computed quantities of interest. In addition to Fourier resolution convergence checks, the size of the computational domain should be varied to estimate the computational accuracy, but this is rarely done [13–15]. In our experience, different choices of PML parameters and domain sizes lead to results that agree qualitatively, but may vary substantially—for example, errors of $Q$ factors $\sim 20\%$ [13] and errors of dipole coupling to radiation modes $\sim 15–25\%$ [16] have been reported.

Instead of searching an extremely large PML parameter space without intuitive or clear guidelines, we propose a different technique that relies on finite-sized structures and open BCs, with fields expanded via Fourier integrals instead of Fourier series. The use of Fourier integrals, in principle, gives an exact description but, for numerical implementation, a $k$-space discretization is required that we, however, have the freedom to choose. Similar ideas have previously been reported for two-dimensional (2D) [17] and rotationally symmetric three-dimensional (3D) [18–20] structures, but without discussion of the important problem of choosing the $k$-space discretization. Furthermore, the important example of dipole emission, which depends sensitively on the proper implementation of the open BCs, was not treated in these works. In this paper, we address both these central issues. Our examples include calculations of light emission from emitters placed in rotationally symmetric waveguides [21] and reflection of the fundamental mode from a waveguide–metal interface [22]. We term this new approach an open-geometry Fourier modal method (oFMM).

This paper is organized as follows. Section 2 outlines the theory of the oFMM approach, while the details are given in Appendix A. The details of the new discretization scheme are discussed in Section 3. The method is tested for several structures by calculating dipole emission rates, $\beta$ factors, and modal reflection coefficients in Section 4. Finally, Section 5 concludes the work.

2. THEORY

In this section, we outline the derivation of the open BC formalism and introduce the theoretical concepts required to understand the results of the following sections. The detailed derivations of the open BC formalisms in rotationally symmetric geometry are given in Appendix A.

A. Open Boundary Condition Formalism

We employ the open BC formalism to describe the electromagnetic field propagation in rotationally symmetric structures. Complete vectorial description is used in connection with Fourier expansion to describe the Maxwell’s equations in a $z$-invariant material section. Using cylindrical coordinates in the rotationally symmetric case allows for simplification of the problem to 1D expansion in the radial coordinate. The $z$ dependence is treated by combining $z$-invariant sections using the scattering matrix formalism (see, e.g., [23] and [24] for details). Our task is then to determine the lateral electric and magnetic field components of the eigenmodes, which are subsequently used as an expansion basis for the optical field. In the conventional FMM, this is done by expanding the field components as well as the permittivity $\epsilon ({\mathbf{r}}_{\perp })$ and $\eta ({\mathbf{r}}_{\perp })\equiv 1/\epsilon ({\mathbf{r}}_{\perp })$ in Fourier series in the lateral coordinates ${\mathbf{r}}_{\perp }$ on a finite-sized computational domain, implying that these functions vary periodically in these coordinates. In the open boundary formalism, we instead consider an infinite-sized computational domain and employ expansions in Fourier integrals. Essentially, any expansion basis can be used, e.g., the plane wave basis in the general 3D case. However, in rotationally symmetric structures, we use the Bessel $J$ function basis since it leads to reduced dimensionality of the problem [18]. In the following, we describe the general steps and equations required to expand the field components and to solve for the expansion and propagation coefficients. The specific equations and derivations are given in Appendix A and referenced throughout this section.

Starting from the time-harmonic Maxwell’s equations $\nabla \times \mathbf{E}(\mathbf{r})=\mathrm{i}\omega {\mu }_0\mathbf{H}(\mathbf{r})$ and $\nabla \times \mathbf{H}(\mathbf{r})=-\mathrm{i}\omega \epsilon (\mathbf{r})\mathbf{E}(\mathbf{r})$ [written using cylindrical coordinates in Eqs. (A1)–(A6) in Appendix A], where $\epsilon$ is the permittivity of the medium, ${\mu }_0$ is the vacuum permeability, $\omega$ is the angular frequency, and $\mathbf{E}$ and $\mathbf{H}$ are the electric and magnetic fields, we obtain

Using the discretized eigenfunction expansion in Eqs. (2) and (3), the fields in each $z$-invariant section can be expressed with column vectors $\mathbf{a}$ consisting of electric and magnetic field expansion coefficients, ${[a_j,a_l\mathrm{\Delta }{k}_l]}^{\mathrm{T}}$, all denoted with the single index $j$ in the following. Thus, taking the $z$ derivative of Eq. (2), we can formulate an eigenvalue problem describing the fields in the system as

Since the eigenfunctions are specific to each layer, we choose a general basis and expand the eigenfunctions in each layer using the common basis. Thus, any function (the field components and the relative permittivity) can be expanded as a Fourier transform

B. Dipole Emission

The field emitted by a point dipole placed at ${\mathbf{r}}_{\mathrm{pd}}$ inside a $z$-invariant structure can be represented as

The emitted field [Eq. (7)] consists of three contributions [25]: guided modes, radiation modes, and evanescent modes. In a waveguide surrounded by air, the modes are guided if the propagation constant ${\beta }_j$ obeys $k_0^2<{\beta }_j^2\le {(n_w{k}_0)}^2$, where $n_w$ is the refractive index of the waveguide. In contrast, the modes are radiating if $0<{\beta }_j^2\le k_0^2$ and evanescent if ${\beta }_j^2<0$. We will apply this classification in Section 3 when we investigate discretization schemes.

The normalized power emitted by dipole to a selected mode can be expressed as [24,26]

3. DISCRETIZATION SCHEME

In addition to the open BCs described in the previous section, an advantage of the presented method is that it enables the use of a nonuniform $k$-space discretization, which allows for a high cut-off value together with dense-sampling $k$-space regions while still maintaining a moderate total number of modes, i.e., achieving the required accuracy with less computational power. In this section, we investigate how to select the cut-off value $k_{\mathrm{cut}\text{-}\mathrm{off}}$ and how to sample the $k$ space effectively. The numerical tests in Section 4 show that faster convergence is achieved using an appropriate mode-sampling scheme.

The transverse wavenumber values in the conventional modal expansion approach [10] are selected equidistantly:

In a bulk medium, light emission occurs with equal weights in all directions. Therefore, a natural starting point for the discretization scheme is to sample the wavevectors in the $(\beta ,{\mathbf{k}}_{\perp })$ plane with equidistant angles [27], as shown in Fig. 1. This is also known as the Chebyshev grid [28]. Then the different transverse wavenumber values are given by $k_m=n{k}_0\sin ({\theta }_m)$, where the equidistantly sampled angles $0<{\theta }_m<\pi /2$ are measured from the $\beta$ axis. Although the values of ${\theta }_m$ are selected uniformly, the values of $k_m$ are more densely sampled in the proximity of $n{k}_0$, cf. Fig. 1. If, instead of a bulk medium, we consider a structure like a nanowire consisting of several materials, it is necessary also to account for the modes beyond $n{k}_0$.

 

Fig. 1. Nonuniform discretization scheme: in a bulk medium, all propagation directions have equal weights. Therefore, the wavevector $\mathbf{k}$ is sampled in the $(\beta ,{\mathbf{k}}_{\perp })$ plane using equidistant angles, as shown by $\theta$ in the figure. Due to the uniform angle distribution, the $k_{\perp }$ discretization is more dense close to $n{k}_0$.

Download Full Size | PDF

To obtain insight into the discretization in different types of structures, we first investigate emission from a radially oriented point dipole placed on the axis of rotationally symmetric infinite semiconductor nanowires having radius from subwavelength to several wavelengths and a refractive index $n_w$ (see Fig 2). The radial component of the emitted electric field $E_r(r)=\sum _j{a}_j{E}_{r,j}(r)=\sum _j{a}_j\sum _m{c}_m{g}_{r,m}(r)\mathrm{\Delta }{k}_m$ can be written by rearranging the terms as follows [cf. Section 2.B, Eqs. (6) and (7), and Appendix A]:

 

Fig. 2. Fourier components of a point dipole emission defined in Eq. (11). The figures show the calculated radiation and guided mode contributions and the total emission as function of the radial wave number in $z$-invariant nanowires of varying radii. The nanowire has a refractive index of $n_w=3.5$ and the wavelength is $\lambda =950\mathrm{nm}$. An equidistant $k$ discretization with 1500 points and $k_{\max }=10k_0$ was used.

Download Full Size | PDF

In the next section, we use these discretization schemes in the modeling of various structures and compare the convergence and required computational power with those obtained using a conventional discretization scheme. When comparing the different discretization schemes, we use the same cut-off value and the same number of modes for both of the schemes.

 

Fig. 3. Example of discretization step sizes $\mathrm{\Delta }{k}_m$ for a nonuniform discretization with $M_1=M_2=M_3=10$ and $k_{\mathrm{cut}\text{-}\mathrm{off}}/k_0=4$. For comparison, the equidistant discretization is also shown with the corresponding number of modes and cut-off values.

Download Full Size | PDF

4. RESULTS AND DISCUSSION

Next, after introducing the principles of the open BC formalism together with the new discretization strategy, we are ready to test the method with several numerical examples. The purpose of these selected examples is to show that the calculations using oFMM formalism converge toward a well-defined open-geometry limit and that faster convergence can be achieved using the discretization schemes introduced in Section 3 compared to using the conventional equidistant discretization. We start with calculating the dipole emission rates (or emission power) in a bulk medium and close to an interface, since these results can be verified analytically. After these basic checks, we investigate the performance of our method for the cases of light emission from emitters in waveguides as well as the case of reflection at a waveguide–metal interface, all of which depend critically on a correct and accurate description of the open boundaries.

A. Dipole Emission in Bulk and Close to an Interface

As a first example, we consider dipole emission in a bulk medium and close to a bulk–bulk interface. Both of these examples can also be solved analytically [26,29], allowing easy comparison of the convergence of the results. Figure 4(a) shows the dipole emission power in a bulk material ($n_b=1$) calculated using the rotationally symmetric model and normalized with the analytical result. Numerical results are calculated using both the equidistant discretization and the nonuniform discretization presented in Section 3. The obtained results show that applying the nonuniform discretization leads to much faster convergence of the emission rates.

 

Fig. 4. (a) Calculated emission power $P_{\mathrm{num}}$ in a bulk medium ($n_b=1$) normalized with analytical result $P_{\mathrm{ana}}$ with a fixed wavelength $\lambda =950\mathrm{nm}$. Both numerical schemes have the wavenumber cut-off value $n_b{k}_0$ in the bulk medium, and the horizontal axis shows the number of modes included in the calculations. (b) Normalized dipole emission power in air in front of a glass ($\epsilon =2.25$) half-space. The dipole is parallel to the interface. (c) Normalized power emitted by point dipole placed in air close to an air–metal interface ($\epsilon =-41+2.5\mathrm{i}$). The dipole is perpendicular to the interface. Numerical results in (b) and (c) are calculated using a cut-off value of $2k_0$ and 200 modes. The powers are normalized with the bulk medium value and the distance $z_0$ from the interface with the wavelength $\lambda =950\mathrm{nm}$.

Download Full Size | PDF

In the bulk medium case, only propagating modes contribute to the light emission, and the emission rate converges provided that enough propagating modes are included in the calculation. In contrast, in the case of a dipole emitter in close proximity to an interface, the evanescent modes also contribute through evanescent mode scattering at the interface and re-excitation of the propagating modes. As a next example, we therefore investigate the interface case.

Figure 4(b) shows the power emitted by a dipole close to an air–glass interface while Fig. 4(c) shows the power emitted by a dipole close to an air–metal interface. The values of the metal and glass permittivities are $\epsilon =-41+2.5\mathrm{i}$ and $\epsilon =2.25$, respectively. In contrast to the bulk medium case in Fig. 4(a) where the cut-off was $n_b{k}_0$ ($n_b=1$), we now need to include the evanescent modes. Figures 4(b) and 4(c) show the separate contributions from propagating and evanescent modes to the emission rate. Again, the nonuniform discretization leads to faster convergence, especially for the contribution from the evanescent modes.

B. Dipole Emitter in a Rotationally Symmetric Waveguide

Next, we investigate the emission in waveguides by calculating the emission rates to selected modes and the spontaneous emission factor $\beta$. In contrast to the bulk medium and interface cases investigated in the previous section, we face an additional computational challenge, which is to compute the radiation modes accurately. The waveguides considered in nanophotonics usually support only a few guided modes. However, the total emission rate and thus the $\beta$ factor depend on emission to a continuum of radiation modes that can radiate light out of the waveguide. Calculating the radiation modes accurately requires more extensive calculations than the emission on bulk medium, as will be seen in the following examples.

Similar to the calculations represented in [21], we consider a dipole emitter oriented along the wire axis in an infinitely long nanowire with $n_w=3.45$, surrounded by air. Figure 5(a) presents the $\beta$ factor and the emission rates to the fundamental guided mode (${\mathrm{HE}}_{11}$), the second guided mode (${\mathrm{E}\mathrm{H}}_{11}$), and the radiation modes, all normalized to the bulk medium emission rate (see Section 2.B) as functions of the nanowire diameter. While the rates calculated using both the equidistant and nonuniform discretization schemes with 1200 modes and a cut-off value $25k_0$ agree well qualitatively, discrepancies are observed in the emission rate to radiation modes. Figure 5(b) shows a convergence investigation of the emission rate to radiation modes. We fix the nanowire geometry by setting the diameter as $0.3\lambda$, use both discretization schemes, and vary the cut-off value of the transverse wavenumber as well as the number of modes. The results show that only a slight improvement is achieved by increasing the cut-off from $20k_0$ to $25k_0$, while the results depend on the number of modes for small mode numbers and converge around 500. At high mode numbers and cut-off values, the results converge to the same value.

 

Fig. 5. Emission from a point dipole placed on the axis of an infinitely long rotationally symmetric nanowire of diameter $d$. (a) $\beta$ factor and normalized emission rates to the first and second guided modes ${\mathrm{HE}}_{11}$, ${\mathrm{E}\mathrm{H}}_{11}$ and radiation modes as functions of $d$. The nanowire refractive index is $n=3.45$ and the wavelength is $\lambda =950\mathrm{nm}$. In both discretization schemes, 1200 modes and cut-off value of $25k_0$ were used. (b) The emission rate to radiation modes calculated with a fixed nanowire diameter $0.3\lambda$. The horizontal axis shows the number of discretization modes, and the legend shows the cut-off value of the wavenumber in units of $k_0$.

Download Full Size | PDF

C. Reflection From Semiconductor Nanowire–Metal Interface

Finally, we investigate convergence of the method in a structure consisting of a nanowire standing on top of a metallic bulk layer. We calculate the reflection coefficient of the fundamental mode from a nanowire–metal interface similar to the setup investigated in [22]. The refractive indices of the nanowire and metal are $n_w=3.5$ and $n_{\mathrm{Ag}}=\sqrt{-41+2.5\mathrm{i}}$ at the wavelength $\lambda =950\mathrm{nm}$.

Figure 6 shows the calculated reflection coefficient as a function of the nanowire diameter using both (a) the equidistant sampling of the $k_{\perp }$ discretization and (b) the nonuniform $k_{\perp }$ discretization with several different numbers of discretization modes. In the nonuniform discretization, the $k$-space values are sampled more densely close to $k_0$ as discussed in Section 3. With small wire diameter, the reflection coefficients are essentially determined by the air–metal reflection ($R_{\mathrm{Air}\text{-}\mathrm{Ag}}\approx 0.98$) since in this limit the fundamental mode is mainly located in air. In contrast, in the limit of large nanowires, the fundamental mode is primarily located in the GaAs wire ($R_{\mathrm{GaAs}\text{-}\mathrm{Ag}}\approx 0.95$). Nevertheless, the figures show that faster convergence is obtained using the nonuniform discretization scheme instead of the equidistant $k$ discretization scheme.

 

Fig. 6. Reflection coefficient of the fundamental mode calculated using (a) an equidistant grid and (b) a nonuniform grid with varying number of modes (shown in the legend) and $k_{\mathrm{cut}\text{-}\mathrm{off}}=20k_0$ as a function of the nanowire diameter. The wire and the metal have refractive indices of $n_w=3.5$ and $n_{\mathrm{Ag}}=\sqrt{41+2.5\mathrm{i}}$, respectively, at wavelength $\lambda =950\mathrm{nm}$. (c) The reflection coefficient of the fundamental mode using equidistant (dotted lines) and nonuniform (dashed lines) discretization and varying the cut-off of $k_m$ for a nanowire having diameter of $0.22\lambda$. The values of $k_m$ are chosen such that $k_m$ is equidistantly/nonuniformly sampled up to value $2n_w{k}_0$ ($n_w=3.5$), with $M$ shown in the legend. Then extra $k_m$ values are added according to the scheme when the cut-off value is increased.

Download Full Size | PDF

The reflection coefficients in Figs. 6(a) and 6(b) are obtained for a fixed cut-off value. Next, we fix the geometry and study the effect of the cut-off value of $k_m$. We select a wire having diameter of $0.22\lambda$ since the reflection coefficients shown in Figs. 6(a) and 6(b), calculated with different discretization schemes and with varying number of modes, have large variations around this diameter. Reflection coefficients as functions of the cut-off value calculated using both discretization schemes, with several different numbers of included modes, are shown in Fig. 6(c). The $k_m$ values are chosen such that, when the cut-off is increased, extra points are added to the original $k_m$ grid. The results show that the calculations converge around $5n_w{k}_0$.

The convergence checks in the selected waveguide examples represented in Figs. 5 and 6 show convergence for the investigated waveguide sizes and structures. Although these examples do not guarantee the convergence of our method in all waveguide sizes and geometries, we expect our method to converge in various types and sizes of waveguides provided that geometry-specific modifications to the discretization scheme are implemented.

5. CONCLUSIONS

We have demonstrated an open-geometry Fourier modal method formalism relying on open boundary conditions and nonuniform $k$-space sampling. Due to the inherent open boundary conditions, we avoid the artificial absorbing boundary conditions that in some cases lead to numerical artifacts. We have tested the approach by investigating the dipole emission in a bulk medium, close to an interface, and in waveguide structures, and by calculating the reflection coefficient of the fundamental waveguide mode for a nanowire–metal interface. Our simulations show that the calculations based on the open-geometry Fourier modal method formalism indeed converge toward an open geometry limit when varying the cut-off and the number of modes, and that the use of the nonuniform discretization scheme leads to a faster convergence of the simulations compared to using the conventional equidistant discretization. We expect that our new method will prove useful in the accurate modeling of a variety of nanophotonic structures, for which the open boundaries are inherently difficult to describe. Also, extension of the formalism to the three-dimensional Fourier modal method is straightforward and could be used for accurate modeling of, for example, light emission in photonic crystal membrane waveguides [2,3].