We present an open-geometry Fourier modal method based on a new combination of open boundary conditions and an efficient -space discretization. The open boundary of the computational domain is obtained using basis functions that expand the whole space, and the integrals subsequently appearing due to the continuous nature of the radiation modes are handled using a discretization based on nonuniform sampling of the space. We apply the method to a variety of photonic structures and demonstrate that our method leads to significantly improved convergence with respect to the number of degrees of freedom, which may pave the way for more accurate and efficient modeling of open nanophotonic structures.
© 2016 Optical Society of America
Many important properties of photonic structures such as cavities  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 , factor of photonic resonators as well as the spontaneous emission (SE) factor in waveguides are important figures of merit in the analysis of nanolasers  and single-photon sources , 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 , while for numerical techniques relying on finite-sized computational domains like the finite-difference time-domain method  and the finite element method , this is achieved using artificial absorbing boundaries, typically in the form of so-called perfectly matched layers (PMLs) .
In Fourier-based modal-expansion techniques [10,11], PMLs can be implemented using complex coordinate transforms . 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 factors  and errors of dipole coupling to radiation modes  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 -space discretization is required that we, however, have the freedom to choose. Similar ideas have previously been reported for two-dimensional (2D)  and rotationally symmetric three-dimensional (3D) [18–20] structures, but without discussion of the important problem of choosing the -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  and reflection of the fundamental mode from a waveguide–metal interface . 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, factors, and modal reflection coefficients in Section 4. Finally, Section 5 concludes the work.
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 -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 dependence is treated by combining -invariant sections using the scattering matrix formalism (see, e.g.,  and  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 and in Fourier series in the lateral coordinates 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 function basis since it leads to reduced dimensionality of the problem . 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 and [written using cylindrical coordinates in Eqs. (A1)–(A6) in Appendix A], where is the permittivity of the medium, is the vacuum permeability, is the angular frequency, and and are the electric and magnetic fields, we obtainA7)–(A9). The fields in the single -invariant section can be expanded using the eigenmodes of the system as
Using the discretized eigenfunction expansion in Eqs. (2) and (3), the fields in each -invariant section can be expressed with column vectors consisting of electric and magnetic field expansion coefficients, , all denoted with the single index in the following. Thus, taking the 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 transformA15) and (A16)]. While in the analytical definition of the Fourier transform the expansion basis in Eq. (5) is infinite, for the numerical calculations the basis must be truncated as 3. This is in contrast to previous approaches [17,18], and we show later that such a nonuniform discretization is a significant improvement. The expansions in the cylindrical coordinate system are given by Eqs. (A15)–(A16). Furthermore, the elements of are given in Eqs. (A23)–(A26). Solving for the eigenvectors and eigenvalues of matrix yields the expansion coefficients and propagation constants in the -invariant structures, while the fields in the full structure are then obtained by combining the -invariant sections using the scattering matrix formalism.
B. Dipole Emission
The field emitted by a point dipole placed at inside a -invariant structure can be represented as24]. The coupling coefficient depends on the dipole position and dipole moment through a dot-product . For the sake of notational clarity, we omit these dependencies in the following. Furthermore, are the expansion coefficients for mode , and are the basis functions.
The emitted field [Eq. (7)] consists of three contributions : guided modes, radiation modes, and evanescent modes. In a waveguide surrounded by air, the modes are guided if the propagation constant obeys , where is the refractive index of the waveguide. In contrast, the modes are radiating if and evanescent if . We will apply this classification in Section 3 when we investigate discretization schemes.26] , where and are the emission rates to mode and in a bulk medium, respectively. In the following, we will use only the normalized unitless quantity for the emission rate. Thus, the spontaneous emission factor (i.e., the factor), defined as the ratio of emission to the fundamental mode (FM) and the total emission , is obtained as
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 -space discretization, which allows for a high cut-off value together with dense-sampling -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 and how to sample the 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  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 plane with equidistant angles , as shown in Fig. 1. This is also known as the Chebyshev grid . Then the different transverse wavenumber values are given by , where the equidistantly sampled angles are measured from the axis. Although the values of are selected uniformly, the values of are more densely sampled in the proximity of , 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 .
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 (see Fig 2). The radial component of the emitted electric field can be written by rearranging the terms as follows [cf. Section 2.B, Eqs. (6) and (7), and Appendix A]:A15). The sum over describes the modes of the structure, while index accounts for the expansion of the modes using the selected basis functions. The first summation inside the brackets in Eq. (11) describes the guided-mode () contribution, while the second summation describes the radiation mode () contribution to the total emission with radial wavenumber . Figure 2 shows the guided and radiation mode contributions and their sum as functions of radial wavenumber. The emitted electric field has a peak around . When the radius increases, a peak around also gradually builds up, while in the bulk limit (), the peak around disappears. These results indicate that (i) for wires with radius , the space should be densely sampled around ; (ii) for wider structures, dense sampling around is also required; while (iii) in bulk media, dense sampling is required only around . Thus, since we are mainly interested in the region where the wire radius is of the order of or smaller than the wavelength, we will use the following discretization scheme that is dense and symmetric around and where the discretization step-size gradually increases toward the cut-off value. Let be the fixed number of values on the intervals , , and , respectively. Then we can write 3.
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.
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 () 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.
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 and , respectively. In contrast to the bulk medium case in Fig. 4(a) where the cut-off was (), 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 . 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 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 , we consider a dipole emitter oriented along the wire axis in an infinitely long nanowire with , surrounded by air. Figure 5(a) presents the factor and the emission rates to the fundamental guided mode (), the second guided mode (), 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 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 , 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 to , 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.
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 . The refractive indices of the nanowire and metal are and at the wavelength .
Figure 6 shows the calculated reflection coefficient as a function of the nanowire diameter using both (a) the equidistant sampling of the discretization and (b) the nonuniform discretization with several different numbers of discretization modes. In the nonuniform discretization, the -space values are sampled more densely close to as discussed in Section 3. With small wire diameter, the reflection coefficients are essentially determined by the air–metal reflection () 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 (). Nevertheless, the figures show that faster convergence is obtained using the nonuniform discretization scheme instead of the equidistant discretization scheme.
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 . We select a wire having diameter of 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 values are chosen such that, when the cut-off is increased, extra points are added to the original grid. The results show that the calculations converge around .
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.
We have demonstrated an open-geometry Fourier modal method formalism relying on open boundary conditions and nonuniform -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].
APPENDIX A: FOURIER–BESSEL EXPANSION IN CYLINDRICAL COORDINATES
The derivation of the open BC method in a rotationally symmetric case is outlined following the approach presented in . We use cylindrical coordinates . Since the considered structures are rotationally symmetric, the angular dependence is expanded using Fourier series . The contributions for different orders are decoupled, and it is thus possible to solve for each order individually. This advantage is exploited to reduce the 2D lateral eigenvalue problem to an effective 1D problem.
Using the Fourier expansion, the time-harmonic Maxwell’s equations and can be written component-wise as
The Helmholtz equation for each Fourier component is given as , which, in component-wise form, reads as
Equations (A7) and (A8) are of the form of Bessel differential equations and couple the radial and angular components of the electric field. In order to simplify calculations, these equations are decoupled using the notation
The transverse components of the Helmholtz equation can then be written as
These Bessel-differential equations have general solutions
Equivalent equations are obtained for magnetic fields by substituting and . The components are obtained using the time-harmonic Maxwell’s equations (A3) and (A6), the above expansions, and the derivation rules for Bessel functions as
To obtain expression for , we expand Eq. (A18) using , integrate both sides of Eq. (A18) with , and use the orthogonality of the Bessel functions. We then obtain the expression for , which is substituted to the expansion of , giving
The expansion coefficients and are obtained representing the system as an eigenvalue problem by applying the differential method as follows. The dependence of the Maxwell’s equations’ expansion coefficients are written as an eigenvalue problem:
Here the dependence is of the form . The derivatives of the electric field expansion coefficients couple only to the magnetic field components and vice versa, so that the propagation constants and expansion coefficients can be solved from the eigenvalue problem
The magnetic field expansion coefficients are obtained from the electric field ones using matrix . Equivalently, the eigenvalue problem can be written for magnetic field coefficients.A25) and (A26) involving can be calculated using the direct rule [30,31] as
The expansion for electric displacement is obtained by inverting as
European Metrology Research Programme (EMRP) (EXL02); Danish Research Council for Technology and Production (DFF-4005-00370); Villum Fonden.
This work was funded by project SIQUTE (contract EXL02) of the European Metrology Research Programme (EMRP). The EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union. Furthermore, support from the Danish Research Council for Technology and Production via the Sapere Aude project LOQIT (DFF-4005-00370), and support from the Villum Foundation via the VKR Centre of Excellence NATEC are gratefully acknowledged.
1. K. J. Vahala, “Optical microcavities,” Nature 424, 839–846 (2003). [CrossRef]
2. G. Lecamp, P. Lalanne, and J. P. Hugonin, “Very large spontaneous-emission β factors in photonic-crystal waveguides,” Phys. Rev. Lett. 99, 023902 (2007). [CrossRef]
3. V. S. C. Manga Rao and S. Hughes, “Single quantum-dot Purcell factor and β factor in a photonic crystal waveguide,” Phys. Rev. B 75, 205437 (2007). [CrossRef]
4. S. Strauf and F. Jahnke, “Single quantum dot nanolaser,” Laser Photon. Rev. 5, 607–633 (2011). [CrossRef]
5. N. Gregersen, P. Kaer, and J. Mørk, “Modeling and design of high-efficiency single-photon sources,” IEEE J. Sel. Top. Quantum Electron. 19, 1–16 (2013). [CrossRef]
6. A. V. Lavrinenko, J. Lægsgaard, N. Gregersen, F. Schmidt, and T. Søndergaard, Numerical Methods in Photonics (CRC Press, 2014), Chap. 7, pp. 197–249.
7. A. Taflove, Computational Electrodynamics: The Finite-difference Time-domain Method, Antennas and Propagation Library (Artech, 1995).
8. J. Reddy, An Introduction to the Finite Element Method (McGraw-Hill, 2005).
9. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]
10. E. Noponen and J. Turunen, “Eigenmode method for electromagnetic synthesis of diffractive elements with three-dimensional profiles,” J. Opt. Soc. Am. A 11, 2494–2502 (1994). [CrossRef]
11. M. G. Moharam, E. B. Grann, D. A. Pommet, and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12, 1068–1076 (1995). [CrossRef]
12. J. P. Hugonin and P. Lalanne, “Perfectly matched layers as nonlinear coordinate transforms: a generalized formalization,” J. Opt. Soc. Am. A 22, 1844–1849 (2005). [CrossRef]
13. N. Gregersen, S. Reitzenstein, C. Kistner, M. Strauss, C. Schneider, S. Höfling, L. Worschech, A. Forchel, T. Nielsen, J. Mørk, and J.-M. Gérard, “Numerical and experimental study of the Q factor of high-Q micropillar cavities,” IEEE J. Quantum Electron. 46, 1470–1483 (2010). [CrossRef]
14. M. Pisarenco, J. Maubach, I. Setija, and R. Mattheij, “Aperiodic Fourier modal method in contrast-field formulation for simulation of scattering from finite structures,” J. Opt. Soc. Am. A 27, 2423–2431 (2010). [CrossRef]
15. P. T. Kristensen, C. V. Vlack, and S. Hughes, “Generalized effective mode volume for leaky optical cavities,” Opt. Lett. 37, 1649–1651 (2012). [CrossRef]
16. J. R. de Lasson, “Modeling and simulations of light emission and propagation in open nanophotonic systems,” Ph.D. thesis (Technical University of Denmark, 2015), available at http://orbit.dtu.dk/files/119895633/PhDThesis_Jakobrdl_Oct2015.pdf.
17. B. Guizal, D. Barchiesi, and D. Felbacq, “Electromagnetic beam diffraction by a finite lamellar structure: an aperiodic coupled-wave method,” J. Opt. Soc. Am. A 20, 2274–2280 (2003). [CrossRef]
18. N. Bonod, E. Popov, and M. Nevière, “Differential theory of diffraction by finite cylindrical objects,” J. Opt. Soc. Am. A 22, 481–490 (2005). [CrossRef]
19. G. P. Bava, P. Debernardi, and L. Fratta, “Three-dimensional model for vectorial fields in vertical-cavity surface-emitting lasers,” Phys. Rev. A 63, 023816 (2001). [CrossRef]
20. M. Dems, I.-S. Chung, P. Nyakas, S. Bischoff, and K. Panajotov, “Numerical methods for modeling photonic-crystal VCSELs,” Opt. Express 18, 16042–16054 (2010). [CrossRef]
21. J. Claudon, N. Gregersen, P. Lalanne, and J.-M. Gérard, “Harnessing light with photonic nanowires: fundamentals and applications to quantum optics,” ChemPhysChem 14, 2393–2402 (2013). [CrossRef]
22. I. Friedler, P. Lalanne, J. P. Hugonin, J. Claudon, J. M. Gérard, A. Beveratos, and I. Robert-Philip, “Efficient photonic mirrors for semiconductor nanowires,” Opt. Lett. 33, 2635–2637 (2008). [CrossRef]
23. L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13, 1024–1035 (1996). [CrossRef]
24. A. V. Lavrinenko, J. Lægsgaard, N. Gregersen, F. Schmidt, and T. Søndergaard, Numerical Methods in Photonics (CRC Press, 2014), Chap. 6, pp. 139–195.
25. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).
26. L. Novotny and B. Hecht, Principles of Nano-Optics, 2nd ed. (Cambridge University, 2012), Chap. 8, pp. 224–281.
27. J. R. de Lasson, T. Christensen, J. Mørk, and N. Gregersen, “Modeling of cavities using the analytic modal method and an open geometry formalism,” J. Opt. Soc. Am. A 29, 1237–1246 (2012). [CrossRef]
28. J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed. (Dover, 2001).
29. L. Novotny and B. Hecht, Principles of Nano-Optics, 2nd ed. (Cambridge University, 2012), Chap. 10, pp. 313–337.
30. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]
31. E. Popov, M. Nevière, and N. Bonod, “Factorization of products of discontinuous functions applied to Fourier-Bessel basis,” J. Opt. Soc. Am. A 21, 46–52 (2004). [CrossRef]