Modeling open nanophotonic systems using the Fourier modal method: Generalization to 3D Cartesian coordinates

TL;DR

This paper extends the Fourier modal method to 3D Cartesian coordinates for open nanophotonic systems, introducing a non-uniform Fourier space sampling technique that enhances accuracy and convergence in modeling open 3D structures.

Contribution

The authors generalize the open boundary Fourier modal method to 3D Cartesian coordinates with a novel non-uniform Fourier space sampling, improving modeling accuracy for open nanophotonic systems.

Findings

Enhanced accuracy in modeling radiation modes

Improved convergence compared to conventional methods

Effective application to various 3D optical structures

Abstract

Recently, an open geometry Fourier modal method based on a new combination of an open boundary condition and a non-uniform $k$-space discretization was introduced for rotationally symmetric structures providing a more efficient approach for modeling nanowires and micropillar cavities [J. Opt. Soc. Am. A 33, 1298 (2016)]. Here, we generalize the approach to three-dimensional (3D) Cartesian coordinates allowing for the modeling of rectangular geometries in open space. The open boundary condition is a consequence of having an infinite computational domain described using basis functions that expand the whole space. The strength of the method lies in discretizing the Fourier integrals using a non-uniform circular "dartboard" sampling of the Fourier $k$ space. We show that our sampling technique leads to a more accurate description of the continuum of the radiation modes that leak out from the structure. We also compare our approach to conventional discretization with direct and inverse factorization rules commonly used in established Fourier modal methods. We apply our method to a variety of optical waveguide structures and demonstrate that the method leads to a significantly improved convergence enabling more accurate and efficient modeling of open 3D nanophotonic structures.

Equations

$$\nabla \times \mathbf{E}(x,y,z)$$

$$\nabla \times \mathbf{H}(x,y,z)$$

$$f(x,y)$$

$$f(x,y)\simeq \sum _{\xi }{c}_f(k_x^{\xi },k_y^{\xi })g(k_x^{\xi },k_y^{\xi })\Delta k_{\xi },$$

$$\left[\begin{array}{cc}{\mathbf{k}}_x{\mathbf{\epsilon }}_{\mathrm{Tot}}^{-1}{\mathbf{k}}_y& -{\mathbf{k}}_x{\mathbf{\epsilon }}_{\mathrm{Tot}}^{-1}{\mathbf{k}}_x+k_0^2\mathbf{I}\\ {\mathbf{k}}_y{\mathbf{\epsilon }}_{\mathrm{Tot}}^{-1}{\mathbf{k}}_y-k_0^2\mathbf{I}& -{\mathbf{k}}_y{\mathbf{\epsilon }}_{\mathrm{Tot}}^{-1}{\mathbf{k}}_x\end{array}\right]\left[\begin{array}{c}{\mathbf{h}}_x\\ {\mathbf{h}}_y\end{array}\right]$$

$$\left[\begin{array}{cc}-{\mathbf{k}}_x{\mathbf{k}}_y& {\mathbf{k}}_x^2-k_0^2{\mathbf{\epsilon }}_{\mathbf{y}}\\ k_0^2{\mathbf{\epsilon }}_{\mathbf{x}}-{\mathbf{k}}_y^2& {\mathbf{k}}_y{\mathbf{k}}_x\end{array}\right]\left[\begin{array}{c}{\mathbf{e}}_x\\ {\mathbf{e}}_y\end{array}\right]$$

$$\mathbf{E}(x,y,z)=\sum _j{a}_j({\mathbf{r}}_{\mathrm{pd}},\mathbf{p}){\mathbf{E}}_j(x,y,z)$$

$$\frac{P_j}{P_{\mathrm{Bulk}}}$$

$${\mathbf{T}}_{12}$$

$${\mathbf{R}}_{12}$$

$$(k_x^m,k_y^l)=(-k_{\mathrm{cut}-\mathrm{off},x}+m\Delta k_x,-k_{\mathrm{cut}-\mathrm{off},y}+l\Delta k_y),$$

$$\begin{array}[]{ll}k_{m}=k_{0}\sin(\theta_{m}),&1\leq m\leq N_{s}/2\\ k_{m}=k_{0}[2-\sin(\theta_{m})],&N_{s}/2+1\leq m\leq N_{s},\end{array}$$

$${\partial }_y{E}_z\mp \mathrm{i}\beta E_y$$

$$\pm \mathrm{i}\beta E_x-{\partial }_x{E}_z$$

$${\partial }_x{E}_y-{\partial }_y{E}_x$$

$${\partial }_y{H}_z\mp \mathrm{i}\beta H_y$$

$$\pm \mathrm{i}\beta H_x-{\partial }_x{H}_z$$

$${\partial }_x{H}_y-{\partial }_y{H}_x$$

$$f(x,y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{c}_f(k_x,k_y)g(k_x,k_y)\mathrm{d}{k}_x\mathrm{d}{k}_y.$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }g(k_x,k_y)g^{\ast }(k_x^{\prime },k_y^{\prime })\mathrm{d}x\mathrm{d}y=(2\pi {)}^2\delta (k_x-k_x^{\prime })\delta (k_y-k_y^{\prime }).$$

$$c_f(k_x,k_y)=\frac{1}{(2\pi {)}^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)g^{\ast }(k_x,k_y)\mathrm{d}x\mathrm{d}y.$$

$$\epsilon (x,y)$$

$$\eta (x,y)$$

$$c_{\epsilon }(k_x,k_y)={\epsilon }_B\delta (k_x)\delta (k_y)+c_{\Delta \epsilon }(k_x,k_y)$$

$$c_{\Delta \epsilon }(k_x,k_y)=\frac{1}{(2\pi {)}^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\Delta \epsilon (x,y)g^{\ast }(k_x,k_y)\mathrm{d}x\mathrm{d}y.$$

$$\int \int \left[k_y{c}_{E_z}(k_x,k_y)\mp \beta c_{E_y}(k_x,k_y)\right]g(k_x,k_y)\mathrm{d}{k}_x\mathrm{d}{k}_y=$$

$$\omega {\mu }_0\int \int c_{H_x}(k_x,k_y)g(k_x,k_y)\mathrm{d}{k}_x\mathrm{d}{k}_y,$$

$$(2\pi {)}^2\int \int \left[k_y{c}_{E_z}(k_x,k_y)\mp \beta c_{E_y}(k_x,k_y)\right]$$

$$k_y{c}_{E_z}(k_x,k_y)\mp \beta c_{E_y}(k_x,k_y)=\omega {\mu }_0{c}_{H_x}(k_x,k_y),$$

$${\mathbf{k}}_y{\mathbf{e}}_z\mp \beta {\mathbf{e}}_y=\omega {\mu }_0{\mathbf{h}}_x,$$

$$\pm \beta {\mathbf{e}}_x-{\mathbf{k}}_x{\mathbf{e}}_z$$

$${\mathbf{k}}_x{\mathbf{e}}_y-{\mathbf{k}}_y{\mathbf{e}}_x$$

$$\int \int \left[k_x{c}_{H_y}(k_x,k_y)-k_y{c}_{H_x}(k_x,k_y)\right]$$

$$k_x^{\prime }{c}_{H_y}(k_x^{\prime },k_y^{\prime })-k_y^{\prime }{c}_{H_x}(k_x^{\prime },k_y^{\prime })$$

$${\mathbf{k}}_x{\mathbf{h}}_y-{\mathbf{k}}_y{\mathbf{h}}_x=-\omega {\epsilon }_0\left[\mathbf{\Delta }\mathbf{\epsilon }\Delta \mathbf{k}+{\epsilon }_B\mathbf{I}\right]{\mathbf{e}}_z,$$

$${\mathbf{e}}_z=-\frac{1}{\omega {\epsilon }_0}{\left[\mathbf{\Delta }\mathbf{\epsilon }\Delta \mathbf{k}+{\epsilon }_B\mathbf{I}\right]}^{-1}\left[{\mathbf{k}}_x{\mathbf{h}}_y-{\mathbf{k}}_y{\mathbf{h}}_x\right].$$

Full text

Fourier modal method, Computational electromagnetic methods, Micro-optics, Waveguides, Mathematical methods in physics, Numerical approximation and analysis

I Introduction

Numerous nanophotonic devices including microcavity resonators Vahala2003 , slow-light waveguides Lecamp2007b ; MangaRao2007 ; Mork2009 and single-photon sources Gregersen2013 ; Aharonovich2016 are open systems with properties strongly characterized by their leakage of light into the surroundings, that in principle extend to infinity. With the exception of the Green’s function integral equation approach Lavrinenko2014 , the majority of conventional approaches for modeling photonic nanostructures including the finite-difference time-domain technique Taflove2004 ; Lavrinenko2014 and the finite-elements method Lavrinenko2014 inherently rely on a limited computational domain combined with either periodic, closed or artificially absorbing boundary conditions (BCs). That is, most conventional methods cannot fully account for the openness of a system, although this is required to correctly model radiative losses. Therefore, simulations of open systems require careful treatment of the boundaries of the computational domain to avoid artificial reflections from the domain wall Berenger1994 ; Hugonin2005a ; Gregersen2008a ; Gregersen2010 .

To circumvent the problem of selecting a proper artificial absorbing BC DeLasson2015a , we have developed a Fourier modal method based on a new combination of an open BC and an efficient discretization scheme Hayrynen2016 , called oFMM in the following. The formalism presented in Hayrynen2016 was, however, limited to rotationally symmetric structures. In this work, we apply both the open boundary formalism and the efficient sampling of the $k$ space to model general 3D structures in Cartesian coordinates, in particular the rectangular waveguide. We remark that the new oFMM formalism affects only the eigenmode calculation; when they have been computed, the oFMM formalism is otherwise identical to the well-established Fourier modal method.

The open boundary of the computational domain can be described by using basis functions, plane waves in this case, expanding the whole infinite space and by using the Fourier transformation. This formalism replaces the usual Fourier series expansion Li1997 ; Li1996b ; Lalanne1998 , which inherently assumes periodicity of the field components. While the use of the Fourier integral transformation gives an exact description of the structure in the limit of continuous $k$-space sampling, the numerical implementation does require a $k$-space discretization. An advantage of the new approach, however, is that we have the freedom to choose the $k$-space discretization in a way that leads to a more efficient mode sampling. Similar ideas have also been reported for two-dimensional (2D) Guizal2003 and rotationally symmetric three-dimensional (3D) Bonod2005 ; Bava2001 ; Dems2010 structures, but without applying efficient $k$-space discretization schemes. In contrast, in our recent work Hayrynen2016 we developed the oFMM approach based on open BCs and a Chebyshev grid Boyd2001 ; DeLasson2012 for rotationally symmetric structures, an approach which we here generalize for any 3D system.

In addition, we discuss how to use Li’s factorization rules Li1997 ; Li1996b ; Lalanne1998 in connection with the 3D oFMM method. It turns out that, while in the rotationally symmetric case Li’s factorization rules are straightforwardly adopted for any $k$-space discretization Bonod2005 ; Hayrynen2016 , for the general 3D approach we can only apply the inverse factorization rule when using the conventional discretization scheme. In spite of this subtlety, we will show that our new discretization scheme leads to a faster convergence compared to traditional discretization schemes.

The manuscript is organized as follows. Section II outlines the theory of the oFMM approach. The details of the new discretization scheme are discussed in Section III. The method is tested by calculating the dipole emission in a waveguide and the reflection of the fundamental mode from a waveguide-metal interface in Section IV. After a discussion of advantages and limitations of the method in Section V, conclusions are drawn in Section VI, and detailed derivations of our theory are provided in the Appendix.

II Theory

In this section, we follow the approach of Ref. Hayrynen2016 and generalize the results for the 3D Cartesian coordinate system. We outline the derivation of the open BC formalism and introduce the theoretical concepts required to understand the results of the following sections. As important examples we show how the oFMM approach is applied to calculate the emission from a dipole placed inside a waveguide and to compute the reflection from a waveguide-metal interface. In Appendix A, we give the detailed derivations of the open geometry formalism and discuss the applicability of the Fourier factorization rules.

II.1 Open boundary condition formalism

We use a complete vectorial description of Maxwell’s equations based on Fourier expansion and open BCs to describe the electro-magnetic (EM) fields in a $z$-invariant material section. The $z$ dependence can be treated by combining $z$-invariant sections using the scattering matrix formalism (see, e.g. Li1996a ; Lavrinenko2014 for details); this part of the calculation remains unchanged by the new oFMM formalism which only alters the way that modes of each $z$-invariant section are calculated. The task is then to compute the lateral electric and magnetic field components of the eigenmodes, which form the expansion basis for the EM field. In the conventional FMM, this is done by expanding the field components as well as the permittivity profile in Fourier series in the lateral coordinates $(x,y)$ 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. We use a plane-wave expansion as basis functions. In the following, we describe the general steps and equations required to expand the field components and to solve for the expansion coefficients and the propagation constant. The specific equations and derivations are given in Appendix A and referenced throughout this section.

We start by considering a $z$-invariant part of the space where the lateral structure is defined by the relative permittivity $\varepsilon(x,y)$ and impermittivity $\eta(x,y)\equiv 1/\varepsilon(x,y)$. For simplicity, we consider a non-magnetic material having vacuum permeability $\mu_{0}$. In such a region of space, we write the Maxwell’s equations using a harmonic time dependence of the form $\exp\left(-\mathrm{i}\omega t\right)$ as

[TABLE]

where $\omega$ is the angular frequency and $\mathbf{E}$ and $\mathbf{H}$ are the vectorial electric and magnetic fields respectively. We then write the fields in a component-wise representation as shown in Eqs. (13)-(18) in Appendix A and introduce a $z$ dependence of the form $\exp\left(\pm\mathrm{i}\beta z\right)$, where $\beta$ is the propagation constant of a particular eigenmode. The individual field components and the permittivity and impermittivity functions are then expanded on plane wave basis functions $g(k_{x},k_{y})=\exp[\mathrm{i}(k_{x}x+k_{y}y)]$ as

[TABLE]

where in the last row the integral expansions are discretized using a Riemann sum on a $(k_{x}^{m},k_{y}^{l})$ grid for numerical calculations. The double summation over the indices $m$ and $l$ (3) is valid for the conventional separable discretization scheme, where the discretization grid coordinates along the $k_{x}$ and $k_{y}$ axes in $k$ space are defined independently of each other. However, when using a non-separable representation as we will do in the following, the Riemann sum is instead written as

[TABLE]

where a single index $\xi$ is used to describe the discretization points in the 2D $k$ space and $\Delta k_{\xi}$ is the discretization area for the $\xi$’th $k$ point. In the particular case of the separable discretization in (3), we have $\Delta k_{\xi}=\Delta k_{x}^{m}\Delta k_{y}^{l}$. This discretization area $\Delta k_{\xi}$ will generally vary as function of $\xi$. Furthermore, as will be discussed in detail in Section III, the selection of the wave number values $k_{x}^{\xi}$ and $k_{y}^{\xi}$ defines the Fourier expansion basis. The computational efficiency of our approach depends crucially on the choice of this expansion basis.

When expanding the field components using the separable discretization, we treat the product of the permittivity function and the electric field components in Eq. (2) using Li’s factorization rules Li1997 ; Li1996b ; Lalanne1998 . However, as discussed in Appendix A, this is not possible when using a non-separable discretization. The details of the expansions are given in Eqs. (26)-(28), (32)-(33) and (37)-(43). After inserting the expansions into Maxwell’s equations and eliminating the $z$ components of the EM fields, we arrive at two set of equations that couple the lateral field components (Eqs. (36) and (46) in Appendix A)

[TABLE]

where $\mathbf{e}_{x}$, $\mathbf{e}_{y}$, $\mathbf{h}_{x}$ and $\mathbf{h}_{y}$ are the vectors of the expansion coefficients of $E_{x}$, $E_{y}$, $H_{x}$, and $H_{y}$, respectively, and $\mathbf{k}_{x}$ and $\mathbf{k}_{y}$ are diagonal matrices of the discretized $k_{x}^{\xi}$ and $k_{y}^{\xi}$ values. Furthermore, $\bm{\varepsilon}_{\mathrm{Tot}}=\bm{\Delta\varepsilon}\bm{\Delta k}+\varepsilon_{B}\mathbf{I}$ with $\bm{\Delta\varepsilon}$ being the Toeplitz matrix defined in Eq. (34), $\mathbf{I}$ is the identity operator and $\bm{\Delta k}$ is the diagonal matrix containing the elements $\Delta k_{\xi}$. When using a separable discretization grid, $\bm{\varepsilon}_{x}$ and $\bm{\varepsilon}_{y}$ are given by Eqs. (42) and (43) respectively. Combining Eqs. (5) and (6) allows us to compute, for example, the lateral electric field components $E_{x,j}(x,y)$ and $E_{y,j}(x,y)$ of the eigenmode $j$ and its propagation constant $\beta_{j}$, after which the lateral magnetic field components $H_{x,j}(x,y)$ and $H_{y,j}(x,y)$ and the longitudinal field components $E_{z,j}(x,y)$ and $H_{z,j}(x,y)$ can be derived. In Appendix A, we show how Li’s factorization rules are correctly used with the oFMM based on equidistant discretization. However, our non-separable ”dartboard” discretization introduced in Section III is not compatible with the inverse factorization rule, and for this reason we employ only the direct factorization rule, which means that we use $\bm{\varepsilon}_{x}=\bm{\varepsilon}_{y}=\bm{\varepsilon}_{\mathrm{Tot}}$.

II.2 Field emitted by a point dipole

In the modal expansion method, the emission from a point dipole placed in a photonic structure can be described Lavrinenko2014 as an expansion of eigenmodes with expansion coefficients proportional to the electric field strength of the corresponding eigenmode obtained from Eqs. (5)-(6) at the emitter position. The total field emitted by a point dipole p placed at $\mathbf{r}_{\mathrm{pd}}$ inside a $z$-invariant structure can be represented as

[TABLE]

where $a_{j}(\mathbf{r}_{\mathrm{pd}},\mathbf{p})$ is the dipole coupling coefficient to mode $j$, which can be calculated using the Lorentz reciprocity theorem Lavrinenko2014 . The coupling coefficient depends on the dipole position $\mathbf{r}_{\mathrm{pd}}$ and dipole moment $\mathbf{p}$ through a dot product $\mathbf{p}\cdot\mathbf{E}_{j}(\mathbf{r}_{\mathrm{pd}})$. For the sake of notational clarity, we omit these dependencies in the following. Furthermore, $c_{j,\xi}$ are the expansion coefficients for mode $j$, and $\mathbf{g}_{\xi}(x_{\mathrm{p}},y_{\mathrm{p}})$ are the vectorial plane wave basis functions.

The emitted field (7) consists of three contributions Snyder1983 : guided modes, radiating modes, and evanescent modes. In a waveguide surrounded by air, the eigenmode $j$ is guided if the propagation constant $\beta_{j}$ obeys $k_{0}^{2}<\beta_{j}^{2}\leq(n_{w}k_{0})^{2}$, where $n_{w}$ is the refractive index of the waveguide. In contrast the mode is radiating if $0<\beta_{j}^{2}\leq k_{0}^{2}$, and evanescent if $\beta_{j}^{2}<0$. We will apply this classification in Section IV when we investigate the performance of the discretization schemes.

The normalized power emitted by a dipole to a selected mode can be expressed as Novotny2012Chap8

[TABLE]

where $P_{\mathrm{Bulk}}=|\textbf{p}|^{2}n_{B}\omega^{4}/(12\pi\epsilon_{0}c^{3})$ is the emitted power in a bulk medium of refractive index $n_{B}$. The normalized power is equal to the normalized spontaneous emission rate Novotny2012Chap8 $\gamma_{j}/\gamma_{\mathrm{Bulk}}=P_{j}/P_{\mathrm{Bulk}}$, where $\gamma_{j}$ and $\gamma_{\mathrm{Bulk}}$ are the spontaneous emission rates to the mode $j$ and to a bulk material, respectively. In the following we will only use the normalized unitless quantity $\Gamma_{j}=\gamma_{j}/\gamma_{\mathrm{Bulk}}$ for the emission rates.

II.3 Reflection at an interface

While the theory above holds for a structure with uniformity along the $z$ axis, most geometries of interest consist of several $z$-invariant sections. The full structure can be described by combining the solutions of Eqs. (5)-(6) using a scattering matrix approach Li1996a ; Lavrinenko2014 . Since our oFMM is based on expanding the fields in each layer using the same basis function, the reflections and transmission of the eigenmodes can be calculated conveniently using the expansion coefficients as described in the following.

Let $\mathbf{C}^{E}_{i}$ and $\mathbf{C}^{H}_{i}$, where $i=1,2$ is the layer index, be matrices whose columns contain the vector expansion coefficients for the lateral electric and magnetic fields respectively computed using (5) and (6). Then, at the interface of material layers 1 and 2, the transmission and reflection matrices are given as Lavrinenko2014

[TABLE]

III Discretization scheme

We have now, via the modal representation in (4), developed a formalism based on a non-uniform $k$-space discretization, which is a generalization of the uniform $k$-space discretization traditionally used in the Fourier modal method. In this section, we describe the important point of how to efficiently sample the $k$ space, before proceeding to example calculations

The lateral expansion basis function $g(k_{x},k_{y})=\exp[\mathrm{i}(k_{x}x+k_{y}y)]$ are plane waves defined entirely by the discretized values of the lateral wavenumbers $k_{x}$ and $k_{y}$. To discretize the transverse expansion basis efficiently in a general 3D approach, we generalize the non-uniform strategy used in the rotational symmetric case Hayrynen2016 .

First, in the conventional equidistant mode discretization approach, the spatial grid in $k$ space is given by

[TABLE]

where $\Delta k_{\alpha}=2k_{\mathrm{cut-off},{\alpha}}/(N_{\alpha}-1)$ and $m,l=0,\dots,(N_{\alpha}-1)$, with $k_{\mathrm{cut-off},\alpha}$ being the cut-off value of the wavenumber and $N_{\alpha}$ the number of modes along the $\alpha=x,y$ axis, see Fig. 1. In the following, when we apply the equidistant discretization scheme, we will use identical cut-off values $k_{\mathrm{cut-off},x}=k_{\mathrm{cut-off},y}$ and modes $N_{x}=N_{y}$ along the $k_{x}$ and $k_{y}$ axes.

Now, the proposed non-uniform circular non-separable discretization approach, which we in the following refer to as the ”dartboard” scheme, is defined as follows. We consider the in-plane wavevector in polar coordinates and set $N_{\phi}$ rays on equidistantly placed angles, cf. Fig. 1. Along each of the rays, the wavenumber values are sampled so that we use dense sampling in the interval $[0,2k_{0}]$ symmetrically placed around $k_{0}$, and in the interval $[2k_{0},k_{\mathrm{cut-off}}]$ a fixed step-size $\Delta k_{\mathrm{tail}}$ is used. The symmetric dense mode sampling is defined using a Chebyshev grid Boyd2001 ; DeLasson2012 as

[TABLE]

where $\theta_{m}=\frac{m\pi}{N_{s}+1}$ and $N_{s}$ is the number of modes in the interval $[0,2k_{0}]$. Thus, in the dartboard discretization approach we have four parameters $N_{\phi}$, $N_{s}$, $\Delta k_{\mathrm{tail}}$, and $k_{\mathrm{cut-off}}$. The motivation of using symmetric dense sampling around $k_{0}$ is to accurately account for the radiating modes as discussed in detail in Hayrynen2016 . In the next section, we show that the dartboard discretization approach outperforms the conventional equidistant mode sampling. As pointed out in Hayrynen2016 , the dartboard mode sampling approach described here is not necessarily the universally optimal, and geometry specific variations may be adopted instead. However, with the proposed approach significant improvement is achieved in terms of the required number of modes and thus of the required computational power.

IV Results

Next, after introducing the principles of the oFMM formalism and the efficient mode sampling scheme, we test our method by investigating its performance for the two cases of light emission by a dipole in a square waveguide as well as of reflection at a waveguide-metal interface. Both examples depend critically on a correct and accurate description of the open BC. We also compare the new discretization scheme to the conventional discretization used in connection with Li’s factorization rules. As already mentioned, in Appendix A we show how Li’s factorization rules are correctly used with oFMM based on equidistant discretization, whereas the equations implementing the non-separable dartboard discretization used in this manuscript are not compatible with the inverse factorization rule. However, our results will demonstrate that, even without the inverse factorization rule the dartboard discretization approach outperforms the equidistant discretization implemented using Li’s factorization rules.

IV.1 Dipole emission in a square waveguide

We first investigate light emission in a square waveguide by calculating the emission rates to the guided modes and to the radiation modes. Additionally, we compute the spontaneous emission factor $\beta$ (not to be confused with the propagation constant $\beta_{j}$) describing the ratio of emitted light coupled to the fundamental guided mode. While typical nanophotonic waveguides support only a few guided modes, the total emission rate and thus the $\beta$ factor depend on the emission into the continuum of radiation modes leaking out of the waveguide. The strength of the oFMM method becomes apparent when determining the light emission to the radiation modes.

Similar to the investigations presented in Lecamp2006 , we consider a dipole emitter oriented along the $x$ axis placed on the axis of an infinitely long square waveguide with varying edge length $w_{x}=w_{y}$ and refractive index $n_{w}=3.5$ surrounded by air. Figure 2(a) presents the $\beta$ factor and the emission rates to the guided modes and to the radiation modes as functions of the waveguide size calculated using the dartboard discretization. The rates are normalized to the bulk emission rate (see Section IIII.2). Figure 2(b) shows the same properties of the waveguide calculated using an equidistant square sampling using either the direct or the inverse factorization rules. The emission rate to the guided modes calculated with the three approaches agree well, and in contrast, a clear difference is seen in the emission rate to the radiation modes and therefore also in the $\beta$ factor. In particular, the coupling to the radiation modes exhibit a spike around a normalized width of 1.15 with the square sampling, which gives an unphysical kink in the $\beta$ factor. The discretization parameters used in Fig. 2 are given in the figure caption and were selected as a result of the convergence investigations presented in the following.

To further investigate the performances of the three approaches we fix the waveguide geometry by setting the width to $w=1.15\lambda/n_{w}$ and vary the cut-off value of the transverse wavenumber as well as the number of modes. This waveguide size is selected for the convergence investigations since a clear difference of the results is seen for this diameter in Fig. 2. Figures 3(a,b) show the convergence investigations of the total emission rate as a function of the cut-off value with several different mode numbers, while Fig. 3(c) shows the convergence of the total emission rate as a function of the number of modes $N_{s}$ in the interval $[0,2k_{0}]$ for the dartboard discretization scheme. The dartboard approach shows clear convergence around a cut-off of $\approx 15k_{0}$. In contrast, the equidistant discretization scheme does not guarantee convergence even with cut-off value of $30k_{0}$. The maximum number of modes used in the calculations of Fig. 3(a) are on the upper limit of the performance of our HPC cluster computer. This is also the case for the highest number of modes used in Figs. 3(b,c). However, in Figs. 3(b,c) the convergence is achieved also for the cases with smaller number of modes.

When using the equidistant discretization, numerical artifacts in the form of large oscillations are observed at particular values of the number of modes and cut-off as displayed in Figs. 2(b) and 3(a) (as well as in Figs. 4(b) and 5(a)). As discussed in Appendix B, the oFMM together with the equidistant square discretization scheme mathematically corresponds to having periodic BCs and to using a Fourier series expansion, where the periodic lengths of the computational domain are inversely proportional to $\Delta k_{x}$ and $\Delta k_{y}$. For geometries with periodic BCs, destructive or constructive interference due to light emission in the neighboring periodic elements may occur leading to the observed large oscillations of the emission rates, that thus are an inherent consequence of the equidistant discretization scheme.

A common approach to circumvent these artifacts due to periodic BCs is to use artificial absorbing BCs, often in the form of the so-called perfectly matched layers (PMLs) Berenger1994 ; Hugonin2005a . However, for the modal method with a PML BC, the convergence of the emission properties with the PML parameters towards the open geometry limit Gregersen2008a ; Gregersen2010 is not well-established with errors in some cases as high as $\approx$ 20 % DeLasson2015a . In contrast, the oFMM with the efficient discretization scheme relies on a truly open computational domain, and therefore avoids using artificial or periodic BCs leading to improved accuracy and convergence towards the true open geometry limit.

IV.2 Reflection from dielectric waveguide-metal interface

As a second example, we investigate convergence of the method for a structure consisting of an infinite waveguide standing on top of a metallic mirror by computing the reflection coefficient of the fundamental guided mode from the waveguide-metal interface. The refractive indices of the waveguide and metal are $n_{w}=3.5$ and $n_{\mathrm{Ag}}=\sqrt{-41+2.5\mathrm{i}}$ at the wavelength $\lambda$ = 1 $\mu$m.

Figure 4 shows the calculated reflection coefficient as a function of the waveguide size $w_{x}=w_{y}$ using (a) the dartboard sampling and (b) the equidistant discretization employing the direct and inverse factorization rules with several different number of discretization modes. The cut-off in all cases is $k_{\mathrm{cut-off}}/k_{0}=14$. Furthermore, for the dartboard discretization fixed values of $N_{\phi}=14$ and $\Delta k_{\mathrm{tail}}/k_{0}=0.2$ were used and only $N_{s}$ was varied. These parameters were chosen to achieve convergence according to the investigations discussed in the next paragraph. In narrow waveguides, the reflection coefficients are essentially determined by the air-metal reflection ($R_{\mathrm{Air-Ag}}\approx 0.98$) since in this limit the fundamental mode is mainly localized in the air surrounding the waveguide. In contrast, in the limit of large waveguides the fundamental mode is primarily confined in the GaAs waveguide ($R_{\mathrm{GaAs-Ag}}\approx 0.95$). A dramatic difference in the results is seen in the region around $w/(\lambda/n_{w})\approx 0.6$, where the reflectivity drops due to a surface-plasmon mediated coupling predominantly to radiation modes propagating in directions perpendicular to the waveguide axis Friedler2008 . When a substantial amount of light is propagating in the $x$-$y$ plane, the performance of the open boundary condition becomes critical, and comparison of Figs. 4(a) and 4(b) clearly demonstrates that this light emission is better resolved using the dartboard discretization.

Whereas the reflection coefficients in Figs. 4(a) and (b) are obtained for a fixed cut-off value, we now fix the geometry and study the effect of the cut-off value of $k_{m}$. We select a waveguide width of $w_{x}=w_{y}=0.63\lambda/n_{w}$, since Fig. 4(b) reveals this to be a challenging computational point. The convergence investigation is shown in Fig. 5. The dartboard discretization (Figs. 5(b,c) again leads to convergence with respect to all of the four discretization parameters. In contrast, no clear convergence is seen when using the equidistant discretization, while we also in this case approach the performance limit of our HPC cluster computer. As discussed in the previous section, the peaks observed in Figs. 4(b) and 5(a) are a consequence of the periodicity of the computational domain when using the equidistant discretization scheme.

V Discussion

The convergence checks in the selected waveguide examples presented in Figs. 3 and 5 show that our method converges for the investigated waveguide sizes and structures. The non-separable nature of our discretization scheme prevents the use of Li’s factorization rules, but even when using the standard direct factorization, a clear improvement in the performance is obtained using the proposed dartboard discretization scheme compared to the conventional equidistant discretization of the basis functions. Although these examples do not guarantee the convergence of our method for all imaginable waveguide sizes and geometries, we generally expect our method to deliver improved performance for various types of waveguides, possibly with additional geometry specific modifications to the discretization scheme.

In high-index-contrast structures as the examples presented here, the FMM method, due to the difficulty of resolving large discontinuities using a plane wave expansion, generally requires a significant amount of modes to achieve convergence. Whereas this may not be a computational difficulty in a rotational symmetric case which in the lateral plane reduces to a 1D problem, the size of the eigenvalue problem in the general planar 2D case rapidly explodes when the number of modes are increased DeLasson2015a . Thus, we expect that a further improvement in terms of computational efficiency could be obtained by combining the dartboard discretization scheme with an adaptive spatial coordinate scheme Essig2010 or by introducing a semi-analytical approach for defining the eigenmodes. In the rotationally symmetric case, exact analytical descriptions of the eigenmodes exist Snyder1983 , while in the rectangular case approximate solutions Raghuwanshi2009 could be used.

VI Conclusion

We have generalized the recently reported open geometry Fourier modal method formalism relying on open boundary conditions and a non-uniform circular ”dartboard” $k$-space sampling for general 3D systems, allowing e.g. the modeling of rectangular waveguides. By applying open boundary conditions, we avoid using the artificial absorbing BCs. We have demonstrated the efficiency of the approach by investigating dipole emission in a square waveguide structure and by studying the reflection coefficient of the fundamental waveguide mode for a waveguide-metal mirror interface, that both are problems of fundamental interest when designing nanophotonic devices. We expect that our new method will prove useful in accurate modeling of a variety of nanophotonic structures, for which correct treatment of an open boundary is crucial.

Acknowledgment

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.

Appendix A Derivation of the eigenvalue problem in 3D open geometry

A.1 Fourier expansion of the field components

The vector components of Maxwell’s equations in Cartesian coordinates are Lavrinenko2014 :

[TABLE]

where the harmonic time-dependence $\exp(-\mathrm{i}\omega t)$ of the fields is assumed and the propagation along the $z$ axis is treated analytically as $\exp(\pm\mathrm{i}\beta z)$. Thus, the field components in a uniform layer only depend on the lateral coordinates $(x,y)$ and are represented as

[TABLE]

The basis function $g(k_{x},k_{y})=\exp[\mathrm{i}(k_{x}x+k_{y}y)]$ are plane waves and satisfy the following orthogonality condition

[TABLE]

The expansion coefficients in Eq. (19) are obtained by multiplying with $g^{*}(k_{x}^{\prime},k_{y}^{\prime})$, integrating over the transverse plane and using the orthogonality relation (20) leading to

[TABLE]

The material properties of the structure are described by the permittivity and impermittivity functions, which are written as a sum between a constant background value and a position dependent deviation from the background value, as

[TABLE]

where $\Delta\varepsilon(x,y)$ and $\Delta\eta(x,y)$ are functions with compact support, such that $\Delta\varepsilon=\Delta\eta=0$ outside a finite domain.

The expansion coefficients of the Fourier transform of the permittivity function can then be written as

[TABLE]

where

[TABLE]

The Fourier transforms of the position dependent deviations, $\Delta\varepsilon$ and $\Delta\eta$, are thus obtained by calculating finite integrals, whereas the constant $\varepsilon_{B}$ and $\eta_{B}$ contributions are handled analytically using Dirac delta functions.

In order to factorize Eqs. (13)-(18) by insertion of the expansion in Eq. (19), Li’s factorization rules Li1997 ; Li1996b ; Lalanne1998 should be considered. Eqs. (13)-(15) and Eq. (18) do not contain any products between two functions with concurrent jumps (discontinuities) and therefore the direct rule applies in these equations. However, in Eq. (16) and (17) the product $\varepsilon E_{x,y}$ is continuous, but both $\varepsilon(x,y)$ and $E_{x,y}(x,y)$ are discontinuous functions, thus they have concurrent jumps, and the inverse factorization rule should - ideally - be used.

A.2 Direct factorization rule

We start by factorizing Eqs. (13)-(18) and writing them in a matrix form one-by-one. Inserting the function expansion in Eq. (19) into Eq. (13) leads to

[TABLE]

where the integration limits (from $-\infty$ to $\infty$) have been omitted for notational clarity. Multiplying with $g^{*}(k_{x}^{\prime},k_{y}^{\prime})$, integrating over $x$ and $y$ and using the orthogonality relation (20) lead to

[TABLE]

Performing the integrations in Eq. (27) we arrive at

[TABLE]

which after discretization of the $k$ space is written in matrix form as

[TABLE]

where $\mathbf{e}_{y}$ is a vector with $c_{E_{y}}^{\xi}$ as elements. $\mathbf{k}_{x}$ and $\mathbf{k}_{y}$ are diagonal matrices with elements $k_{x}^{\xi}$ and $k_{y}^{\xi}$.

Using a similar approach, Eqs. (14) and (15) are written in matrix form as

[TABLE]

Next we prepare Eq. (18) in a discretized form in order to eliminate $\mathbf{e}_{z}$ from Eqs. (29) and (30), which can also be performed by applying the direct factorization rule. Expanding the field components, using (24)-(25) and performing a change of variables $\hat{k}_{x,y}=k_{x,y}+k_{x,y}^{\prime}$ lead to

[TABLE]

We then multiply with $\exp[-\mathrm{i}(k_{x}^{\prime}x+k_{y}^{\prime}y)]$, integrate over $x$ and $y$ and employ the orthogonality condition (20) and obtain

[TABLE]

In discretized form Eq. (33) is written as

[TABLE]

where $\bm{\Delta\varepsilon}$ is the Toeplitz matrix containing the elements $c_{\Delta\varepsilon}^{\xi}=c_{\Delta\varepsilon}(k_{x}^{\xi},k_{y}^{\xi})$, $\mathbf{I}$ is the identity matrix and $\Delta\mathbf{k}$ is the diagonal matrix containing the discretized area elements $\Delta k_{\xi}$ in $k$ space. Thus, $\mathbf{e}_{z}$ equals to

[TABLE]

allowing us to write Eqs. (29) and (30) in the form of an eigenvalue problem that couples the lateral electric field components to the lateral magnetic field components as

[TABLE]

where $\bm{\varepsilon}_{\mathrm{Tot}}=\bm{\Delta\varepsilon}\Delta\mathbf{k}+\varepsilon_{B}\mathbf{I}$.

From Eqs. (16) and (17) we can write similar set of equations that couples the lateral components so that Eqs. (16) and (17) together with Eq. (36) allows us to eliminate the magnetic field components and form an eigenvalue problem for the lateral electric field components (or vice versa). However, Eqs. (16) and (17) need special treatment due to the product $\varepsilon E_{x,y}$.

A.3 Inverse factorization approach

In the following the application of the inverse rule for open boundaries with a separable discretization grid in $k$ space will be presented. As discussed in Appendix B, an equidistant discretization with an open BC is mathematically equivalent to implementing a periodic BC and a Fourier series expansion. Furthermore, as will become apparent in the course of deriving the inverse factorization for the separable discretization, the inverse factorization approach is not applicable for our dartboard discretization scheme defined in Section III.

The factorization will be performed on Eq. (16) to illustrate how the inverse rule is implemented for the product $\varepsilon E_{x}$. The matrix representation for the $\varepsilon$ function used in the product $\varepsilon E_{x}$ will be denoted $\bm{\varepsilon}_{x}$, indicating that it accommodates for continuity of the product along the $x$ direction, where the inverse rule is applied as in Li1997 ; Lalanne1998 . Now, $E_{x}$ is discontinuous in the $x$ direction but continuous in the $y$ direction. $\varepsilon$ is discontinuous in both the $x$ and $y$ direction. Their product, $\varepsilon E_{x}$, is continuous in the $x$ direction and discontinuous in the $y$ direction, thus the inverse rule is used for the $x$ direction and the direct rule for the $y$ direction. The way this is done computationally is to divide the structure into sections separated by the interfaces in the $y$ direction and apply the inverse rule to each of these sections. This is illustrated in Figure 6.

In general the expansion coefficients for all $(x,y)$-dependent functions are given as in Eq. (21). The integration over the $y$ coordinate is then separated into sections where the function is uniform along the $y$ axis. Using Figure 6 as the example, the $y$ integration is separated into three parts

[TABLE]

where

[TABLE]

Here the notation $f(x,y_{S_{i}})$ means that the function is evaluated within section $S_{i}$ and is only dependent on the $x$ coordinate within that section. With this separation, it is possible to factorize $\varepsilon$ using the correct factorization rules provided that the discretized basis set features separable $k_{x}$ and $k_{y}$ dependency as in (11). If this is the case, we can index the $k_{x}$ and $k_{y}$ contributions to the basis mode $k$ vector as $(k_{x}^{m},k_{y}^{l})$ using separate indices $m$ and $l$. It is then possible to apply the inverse rule to the product $\varepsilon_{x}E_{x}$ factorized along the $x$ direction by first preparing the Fourier transform along the $x$ axis of the inverse permittivity as

[TABLE]

We then form the Toeplitz matrix for the $\eta_{x,S_{i}}$ function discretized on the $k_{x}^{m}$ grid. Since the product of the expansions of $\varepsilon$ and $E_{x}$ involves an integration over $k$ space as in (33), the Toeplitz matrix is given by

[TABLE]

where $\bm{\Delta\eta}_{x,S_{i}}$ is the Toeplitz matrix containing the elements $\Delta\eta_{x,S_{i}}^{m}=\Delta\eta_{x,S_{i}}(k_{x}^{m})$ and $\Delta\mathbf{k}_{x}$ is the diagonal matrix with $\Delta k_{x}^{m}$ as elements. According to the inverse rule, we then take the inverse of this matrix and Fourier transform the resulting elements along the $y$ axis as

[TABLE]

where $\bm{\Delta\eta}_{x,\mathrm{Tot}}^{\mathrm{Inv}}(y)=\bm{\eta}_{x,S_{i},\mathrm{Tot}}^{-1}-\varepsilon_{B}\bm{I}$, which is piece-wise constant over the various regions $S_{i}$ as discussed above.

The final Toeplitz matrix $\bm{\varepsilon}_{x}$ is then obtained by introducing the discretization on the $k_{y}^{l}$ grid, and its elements are given by

[TABLE]

Similarly for the product $\varepsilon_{y}E_{y}$ we obtain

[TABLE]

The integrals in Eqs. (42) and (43) can be carried out analytically when the matrix $\Delta\bm{\eta}_{x(y),\mathrm{Tot}}^{\mathrm{Inv}}$ has been found for each $y(x)$-independent section.

The factorization of Eqs. (16)–(17) thus become

[TABLE]

Eliminating $\mathbf{h}_{z}$ using Eq. (31) finally leads to the following eigenvalue problem

[TABLE]

The splitting of the factorization along the $x$ and $y$ axes such that the inverse rule can be used along the $x$ axis and the direct rule along the $y$ axis relies on the separability of the $k_{x}$ and $k_{y}$ dependencies of the discretization grid such that the discretization in (40) can be performed in a well-defined manner. However, for our dartboard discretization scheme, this separation is not possible, and for this reason, we simply use the direct rule for the factorization with $\bm{\varepsilon}_{x}=\bm{\varepsilon}_{y}=\bm{\varepsilon}_{\mathrm{Tot}}$.

Appendix B Relationship between open and periodic boundary conditions

To understand the equivalence between the open BC formalism with equidistant discretization and the periodic BC formalism, let us consider the representation of a function $f(x)$ with compact support such that $f(x)=0$ for $|x|>L/2$. The continuous integral expansion of this function is given by

[TABLE]

where the integration domain in (48) has been reduced from $[-\infty,\infty]$ to $[-L/2,L/2]$ since $f(x)=0$ outside this range.

We now implement the equidistant discretisation scheme with a discretization step $\Delta k$ such that (47) becomes

[TABLE]

where $k_{n}=n\Delta k$.

Let us compare this equation to the Fourier series expansion of the same function over the interval $[-L/2,L/2]$ given by

[TABLE]

where $k_{n}=n2\pi/L$. Now, the integral expansion (47)-(48) should ideally reproduce a function $f(x)$ for which $f(x)=0$ for $|x|>L/2$. However, we observe that the representation in (49) implementing the equidistant discretization is mathematically equivalent to the standard Fourier series representation (50)-(51) of a periodic function $f(x)=f(x+L)$, where the periodicity is given by

[TABLE]

When representing the optical fields using an open BC and equidistant discretization, we are thus in practice reintroducing a periodic BC with the associated numerical artifacts due to the presence of the neighboring elements. The artifacts can be suppressed by decreasing $\Delta k$, in which case the Riemann sum representation of the Fourier transform approaches the exact value of the integral. However, this occurs at the expense of significant computational cost, and a non-uniform discretization scheme is thus strongly preferred.