The Fourier modal method (FMM) is one of the most powerful electromagnetic analysis methods in the field of photonics [1]. In FMM, the linear modal theory for light-matter interactions is constructed based on Fourier analysis of Maxwell’s equations. The modal analysis is the crucial feature of FMM, which expands the vectorial optical field into a linear superposition of electromagnetic Bloch eigenmodes. In FMM, the dispersion relation and coupling rate of each eigenmode are extracted completely. On the basis of FMM, the modal spectrum of the vectorial optical field in photonic structures can be manifested [1]. For recent active-photonic applications such as quantum-dot light emitters or photoluminescence in nanowires, where light sources are embedded inside photonic structures, the modal analysis is also of prominent interest. A few previous studies dealt with the numerical modeling of a dipolar point source. The previous studies mainly investigated a periodic dipole arrangement in periodic structures. The discrete dipole approximation in a small-particle grid [2], the one-dimensional integral model of a dipole located within or at the boundary of a planar layered structure [3], periodic arrangement of dipoles [4, 5], RCWA analysis of a central electric dipole in a nonuniform spherical system [6], and more, have been reported. However, generalizing FMM to an aperiodic structure is strongly needed, to analyze recent state-of-the-art light-source-embedding photonic structures, including larger-scale aperiodic active photonic structures such as organic light-emitting diodes (OLEDs). In contrast, our method introduces a generalized FMM technique for aperiodic multiblock structures using a scattering matrix (S matrix), and presents an aperiodic-optical-dipole-source model that prevents numerical errors due to the Gibbs phenomenon. In this paper, the aperiodic model is enabled using perfect matched layers (PML), and the validity is verified through simulation of the real OLED model.

Here a novel FMM scheme for optical dipole radiation in a photonic structure is proposed, a scheme that is expected to be able to provide the requisite mode-theoretical analysis framework for light-source-embedded structures not only for the optical field distribution, but also for power flow and energy loss. For our theoretical description, we follow the conventions and symbols used in the previous FMM study [1].

Two contrasting optical field distributions, generated by two dipole emitters with orthogonal polarization embedded in the active block, are presented in Fig. 1. In Fig. 1, the example multiblock structure is schematically shown, with specifications. A two-dimensional dipole emission source of wavelength 532 nm is placed at r = (0, 0, 0) in the active block. The vectorial field distribution for the dipole with polarization P = (1, 0, 0), and that for the dipole with polarization P = (0, 0, 1), are presented in Figs. 1(a) and 1(b) respectively.

Figure 1. Photonic multiblock structure with an optical-source-embedding block, and the optical field distributions for (a) a dipole with polarization P = (1, 0, 0) and (b) a dipole with polarization P = (0, 0, 1). The thicknesses of block #1, the cathode, block #2, the active block, block #3, and the anode are 70 nm, 8 nm, 31 nm, 38 nm, 157 nm, and 114 nm respectively, and their refractive indices (n + jk) are 1.88, 0.191 + 3.24j, 1.71, 1.84 + 0.00272j, 1.94, and 0.129 + 3.19j, respectively.

The field distribution is strongly dependent on the polarization state of the emission source. These field distributions are calculated by applying the method proposed in this paper to a basic radiative-device structure, to begin the description of the concrete theory more conveniently. This calculation is carried out by FMM, which allows us to analyze the modal structure of a field distribution. The modal-structure analysis provides in-depth insight about the optical power flow and energy loss, and their structural dependency. These topics will be dealt with in our further papers in the near future, but here we focus on the development of the novel extended FMM for optical dipole radiation.

Such a light-source-embedded photonic structure is modeled by the bidirectional multiblock structure shown in Fig. 2(a). The blocks and their interfaces are denoted by L_k and B_k respectively. The bidirectional multiblock structure is composed of three separate parts. The active-source block denoted by L₀ is placed in the center of the structure. It is assumed in the proposed model that a dipole emitter is placed in the homogeneous block L₀ with a permittivity of ε_r. The left and right parts of the active-source block are multiple stacked single blocks L_k, with the indices −(N_L + 1) ≤ k ≤ −1 and 1 ≤ k ≤ (N_R + 1) respectively. The notation follows the conventional structure-modeling method for FMM [1]. On the whole, the structure takes the form of a finite multiblock structure surrounded by the left and right semi-infinite blocks L_−(NL+1) and L_(NR+1).

Figure 2. Bidirectional S-matrix representation of a photonic structure with an optical-source-embedding block: (a) schematic diagram of the source-embedded multiblock structure and (b) bidirectional S-matrix model.

In the proposed scheme, the thickness of the dipole-embedding slab L₀ is set to zero, in the sense of a limit. Even though the thickness of the block is zero, the permittivity of the block must be taken into account to define the internal impedance of the dipole emitter. For example, the permittivity of L₀ should be ε_r = 1 for a dipole embedded in free space, whereas we need to set the permittivity of L₀ as ε_r(r₀)when the material permittivity at the point r = r₀ is ε_r(r₀). This is critical for matching the dipole’s internal impedance to that of the surrounding medium.

The optical field in the multiblock structure is represented by the Fourier modal expansions in the respective blocks. Here it is assumed that the complete set of the Bloch eigenmodes $E_{(k)}^{(g)+},E_{(k)}^{(g)-},H_{(k)}^{(g)+},and{H}_{(k)}^{(g)-}$ represented by the pseudo-Fourier series are prepared for all single blocks [1]. $E_{(k)}^{(g)\pm }and{H}_{(k)}^{(g)\pm }$
are the electric field and magnetic field of the g^th positive (negative) eigenmode in the k^th block L_k. Here the time-harmonic function exp(−jωt) is assumed, where ω is the angular frequency. The field distributions in the left and right parts of the multiblock structure are represented respectively by

(1a)$$\left(\begin{array}{l}E\\ H\end{array}\right)={\displaystyle \sum _{g=1}^{M^{+}}{\overleftarrow{C}}_{(-k),g}^{+}\left(\begin{array}{l}{E}_{(-k)}^{(g)+}\\ H_{(-k)}^{(g)+}\end{array}\right)}+{\displaystyle \sum _{g=1}^{M^{-}}{\overleftarrow{C}}_{(-k),g}^{-}\left(\begin{array}{l}{E}_{(-k)}^{(g)-}\\ H_{(-k)}^{(g)-}\end{array}\right)},for{B}_{-k}\le z\le B_{-(k-1)},1\le k\le N_{L,}$$

(1b)$$\left(\begin{array}{l}E\\ H\end{array}\right)={\displaystyle \sum _{g=1}^{M^{+}}{\overrightarrow{C}}_{(k),g}^{+}\left(\begin{array}{l}{E}_{(k)}^{(g)+}\\ H_{(k)}^{(g)+}\end{array}\right)}+{\displaystyle \sum _{g=1}^{M^{-}}{\overrightarrow{C}}_{(k),g}^{-}\left(\begin{array}{l}{E}_{(k)}^{(g)-}\\ H_{(k)}^{(g)-}\end{array}\right)},for{B}_{k-1}\le z\le B_k,1\le k\le N_R,$$

where ${\overleftarrow{C}}_{(-k),g}^{\pm }and{\overrightarrow{C}}_{(k),g}^{\pm }$ are the coupling coefficients of the Bloch eigenmodes in the left and right finite blocks. The field distribution in the left and right semi-infinite blocks are explicitly given as, respectively,

(1c)$$\left(\begin{array}{l}E\\ H\end{array}\right)={\displaystyle \sum _{g=1}^{M^{-}}{\overleftarrow{T}}_{(-N_L-1),g}^{}\left(\begin{array}{l}{E}_{(-N_L-1)}^{(g)-}\\ H_{(-N_L-1)}^{(g)-}\end{array}\right)},forz\le B_{-N_L},$$

(1d)$$\left(\begin{array}{l}E\\ H\end{array}\right)={\displaystyle \sum _{g=1}^{M^{+}}{\overrightarrow{C}}_{(N_R+1),g}^{+}\left(\begin{array}{l}{E}_{(N_R+1)}^{(g)+}\\ H_{(N_R+1)}^{(g)+}\end{array}\right)},for{B}_{N_R}\le z,$$

where ${\overleftarrow{T}}_{(-N_L-1),g}and{\overrightarrow{T}}_{(N_R+1),g}$ are the transmission coefficients in the left and right semi-infinite blocks respectively. The scattering matrices (S matrices) contain complete information about the electromagnetic characteristics of the target multiblock structure [1, 7]. The coupling-coefficient operators [1, 8] indicate the excitation rates of the internal Bloch eigenmodes. Those are essential to visualizing the optical field, as well as analyzing the modal spectrum of the internal optical field. The distribution of the coupling coefficients can be straightforwardly interpreted as the modal spectrum.

In terms of the S-matrix method, the field propagation and distribution inside the structure is completely represented in the modal spectrum’s domain. For that, the transmission and coupling coefficients are treated in the form of matrix operators and associated recurrence relations [1]. The recurrence relations reflect the internal multiple reflection and transmission processes between the individual blocks. The operators of the transmission coefficients for the left semi-infinite block ${\overleftarrow{T}}_{(-N_L-1),g}$ are denoted by $\overleftarrow{T}$, whereas those of the right block ${\overrightarrow{T}}_{(N_R+1),g}$ are represented by $\overrightarrow{T}$. The coupling coefficients of the left and right sides, ${\overleftarrow{C}}_{(-k),g}^{\pm }and{\overrightarrow{C}}_{(k),g}^{\pm }$, are represented by the collective coupling-coefficient operators ${\overleftarrow{C}}^{(L)}and{\overrightarrow{C}}^{(R)}$ in FMM [1, 8].

The conventional S matrix has not been extended to represent the complicated optical wave propagation in the optical-source-embedding structure. We propose a bidirectional S-matrix algorithm for stable field propagation, for that problem. The proposed bidirectional S-matrix formulation is schematically illustrated in Fig. 2(b). The S matrices of the left and right multiblocks, S^(L) and S^(R), are given by

(2a)$$S^{(L)}=\left(\begin{array}{cc}{\overrightarrow{T}}^{(L)}& {\overrightarrow{R}}^{(L)}\\ {\overleftarrow{R}}^{(L)}& {\overleftarrow{T}}^{(L)}\end{array}\right),$$

(2b)$$S^{(R)}=\left(\begin{array}{cc}{\overrightarrow{T}}^{(R)}& {\overrightarrow{R}}^{(R)}\\ {\overleftarrow{R}}^{(R)}& {\overleftarrow{T}}^{(R)}\end{array}\right).$$

The transmission-coefficient operators in the left and right parts, induced by the positive input field, are denoted by ${\overrightarrow{T}}_{p}and{\overleftarrow{T}}_p$ (with the subscript p for ‘positive’) and ${\overrightarrow{T}}_{n}and{\overleftarrow{T}}_n$ (with the subscript n for ‘negative’) respectively. In Fig. 2(b), the thickness of the intermediate source-embedding block L₀ is indicated by d_v. We can construct the bidirectional S-matrix algorithm following the recursive ray-tracing method [1], taking the limiting process of lim d_v → 0. The four components of the bidirectional S-matrix algorithm are derived as

(3a)$${\overrightarrow{T}}_p={\overrightarrow{T}}^{(R)}{(I-{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)})}^{-1},$$

(3b)$${\overleftarrow{T}}_p={\overleftarrow{T}}^{(L)}{(I-{\overleftarrow{R}}^{(R)}{\overrightarrow{R}}^{(L)})}^{-1}{\overleftarrow{R}}^{(R)},$$

(3c)$${\overrightarrow{T}}_n={\overrightarrow{T}}^{(R)}{(I-{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)})}^{-1}{\overrightarrow{R}}^{(L)},$$

(3d)$${\overleftarrow{T}}_n={\overleftarrow{T}}^{(L)}{(I-{\overleftarrow{R}}^{(R)}{\overrightarrow{R}}^{(L)})}^{-1}.$$

More specifically, we define the internal-light-source operators along the positive z direction and negative z direction denoted by identity matrices $\overrightarrow{U}and\overleftarrow{U}$ respectively. Considering the first component ${\overrightarrow{T}}_p\overrightarrow{U}$, we can expand it as

(3e)$${\overrightarrow{T}}_p\overrightarrow{U}={\overrightarrow{T}}^{(R)}{\left(I-{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\right)}^{-1}\overrightarrow{U}={\overrightarrow{T}}^{(R)}\overrightarrow{U}+{\overrightarrow{T}}^{(R)}{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\overrightarrow{U}+{\overrightarrow{T}}^{(R)}{\left({\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\right)}^2\overrightarrow{U}+{\overrightarrow{T}}^{(R)}{\left({\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\right)}^3\overrightarrow{U}\cdots .$$

The first term ${\overrightarrow{T}}^{(R)}\overrightarrow{U}$ is the direct transmission of $\overrightarrow{U}$ through the right block, and the second term ${\overrightarrow{T}}^{(R)}{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\overrightarrow{U}$ is interpreted as the transmission after a single round-trip reflection in the zero-thickness block between the left and right parts. The transmission-coefficient operators for the internal-light-emission source are obtained in the form of the sum of the field components excited by $\overrightarrow{U}and\overleftarrow{U}$ as

(4a)$$\overrightarrow{T}={\overrightarrow{T}}_p\overrightarrow{U}+{\overrightarrow{T}}_n\overleftarrow{U},$$

(4b)$$\overleftarrow{T}={\overleftarrow{T}}_p\overrightarrow{U}+{\overleftarrow{T}}_n\overleftarrow{U}.$$

Also, taking into account the internal field distribution, the coupling-coefficient operators are given by the following algorithm:

(5a)$${\overrightarrow{C}}_P={\overrightarrow{C}}^{(R)}{\left(I-{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\right)}^{-1},$$

(5b)$${\overleftarrow{C}}_p={\overleftarrow{C}}^{(L)}{\left(I-{\overleftarrow{R}}^{(R)}{\overrightarrow{R}}^{(L)}\right)}^{-1}{\overleftarrow{R}}^{(R)},$$

(5c)$${\overrightarrow{C}}_n={\overrightarrow{C}}^{(R)}{\left(I-{\overrightarrow{R}}^{(L)}{\overleftarrow{R}}^{(R)}\right)}^{-1}{\overrightarrow{R}}^{(L)},$$

(5d)$${\overleftarrow{C}}_n={\overleftarrow{C}}^{(L)}{\left(I-{\overleftarrow{R}}^{(R)}{\overrightarrow{R}}^{(L)}\right)}^{-1}.$$

Similar to the formation of Eqs. (4a) and (4b), the coupling-coefficient operators of the internal field distributions in the left and right parts are obtained by, respectively,

(6a)$$\overrightarrow{C}={\overrightarrow{C}}_p\overrightarrow{U}+{\overrightarrow{C}}_n\overleftarrow{U},$$

(6b)$$\overleftarrow{C}={\overleftarrow{C}}_p\overrightarrow{U}+{\overleftarrow{C}}_n\overleftarrow{U}.$$

Consequently, the bidirectional S-matrix method has been established. Practical application of optical sources requires additional modeling of a specific optical source, which is given by the discrete column vectors of the bidirectional angular spectrum, $\overrightarrow{U}and\overleftarrow{U}$, as explained in the next section.

To represent the dipole-radiation field, a proper source operator model, $\overrightarrow{U}and\overleftarrow{U}$ , should be designed. For this, we investigate the discretization of the Weyl representation of the Green dyadic function [9] and propose the discrete Weyl representation of an optical dipole-emission source for FMM. In the aforementioned S-matrix formulation, the source operators $\overrightarrow{U}and\overleftarrow{U}$ need to be specified only in the source-embedding block. The emitted optical field propagating toward the other blocks is calculated by the bidirectional S-matrix algorithm established in Section 2. Thus, our goal is to construct the operators of the dipole source in L₀.

The optical emission in L₀ is described by the inhomogeneous Maxwell’s equations:

(7a)$$\nabla \times E=j\omega {\mu }_{0}H+j\omega {\mu }_{0}M\delta (r),$$

(7b)$$\nabla \times H=-j\omega {\epsilon }_0{\epsilon }_{r}E-j\omega P\delta (r),$$

where P and M are the external polarization and magnetization vectors respectively. ε₀ and μ₀ denote respectively the electric permittivity and magnetic permeability in free space. ε_r is the relative permittivity of L₀ and δ(r) is the three-dimensional Dirac delta function. The permittivity can be generally represented by an anisotropic tensor, but for convenience ε_r is set to a constant here. The electric and magnetic dipole sources are placed at the origin (r = 0).

First, we take into account the electric dipole radiation in L₀ by setting the magnetization vector M = 0. The dyadic Green function of the optical dipole is expressed by the Weyl representation [8],

(8a)$${\overleftrightarrow{G}}^{\pm }(r)=\frac{j{\epsilon }_r}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\frac{{\overleftrightarrow{Q}}^{\pm }}{k_z}}}{e}^{j(k_{x}x+k_{y}y+k_z\left|z\right|)}d{k}_{x}d{k}_y,$$

where $k_z^2=k_0^2{\epsilon }_r-k_x^2-k_y^2$. ${\overleftrightarrow{Q}}^{+}and{\overleftrightarrow{Q}}^{-}$ are the operators generating the three-dimensional dipole field propagating in the positive and negative z directions respectively, and are given by two partial operator matrices,

(8b)$${\overleftrightarrow{Q}}^{\pm }=\left(\begin{array}{ccc}1-p^2& -pq& \pm pm\\ -pq& 1-q^2& \pm qm\\ \pm mp& \pm mq& 1-m^2\end{array}\right)for\pm z\ge 0,$$

where p, q, and m are given by

(8c)$$(p,q,m)=(k_x/(k_0\sqrt{{\epsilon }_r}),k_y/(k_0\sqrt{{\epsilon }_r}),\sqrt{1-p^2-q^2}).$$

The electric and magnetic fields induced by a dipole source P are obtained respectively from

(9a)$$E=\frac{{\mu }_0{\omega }^2}{4\pi {\epsilon }_r}{\overleftrightarrow{G}}^{\pm }(r)\cdot P=\left\{\begin{array}{c}\frac{{\mu }_0{\omega }^2}{4\pi {\epsilon }_r}{\overleftrightarrow{G}}^{+}(r)\cdot Pforz\ge 0\\ \frac{{\mu }_0{\omega }^2}{4\pi {\epsilon }_r}{\overleftrightarrow{G}}^{-}(r)\cdot Pforz<0\end{array}\right.,$$

(9b)$$H={\left(j\omega {\mu }_0\right)}^{-1}\nabla \times E$$

The radiation power of a dipole, P = (0, 0, p_ez), in a homogeneous medium with refractive index $n=\sqrt{{\epsilon }_r}$ is n times larger than that in free space, as

(9c)$$W_n=\frac{\omega }{2}\mathrm{Im}{\left[E\cdot P^{*}\right]}_{r=0}=n\left(\frac{{\mu }_0}{4\pi }\right)\left(\frac{{\omega }^4}{3c}\right)p_{ez}^2=n{W}_{free}.$$

Similarly, the magnetic and electric fields induced by the magnetic dipole can be obtained as

(10a)$$H=\left(\frac{{\omega }_0^2{\epsilon }_0{\mu }_0}{4\pi }\right){\overleftrightarrow{G}}^{\pm }(r)\cdot M\left|for\pm z\ge 0,\right.$$

(10b)$$E=-{\left(j\omega {\epsilon }_0{\epsilon }_r\right)}^{-1}\nabla \times H.$$

In the case of a magnetic dipole, the radiation power is n³ times as large as that in free space,

(10c)$$W_n=\frac{\omega {\mu }_0}{2}\mathrm{Im}{\left[H\cdot M^{*}\right]}_{r=0}=n^3\left(\frac{{\mu }_0}{4\pi }\right)\left(\frac{{\omega }^4}{3c^3}\right)m_{ez}^2.$$

The power calculations shown above are rigorous, because the calculation is performed on a continuum support. Since the FMM is a numerical method in the spatial Fourier domain, we need to build up a numerical framework. The core of our discussion is the discretization of the Weyl representation of the Green dyadic function. Direct discretization results in the discrete representation

(11a)$${\overleftrightarrow{G}}^{\pm }(r)=\Delta k_x\Delta k_y\frac{j{\epsilon }_r}{2\pi }{\displaystyle \sum _{m=-M}^M{\displaystyle \sum _{n=-N}^N{\overleftrightarrow{Q}}_{m,n}^{\pm }{e}^{j(k_{x,m,n}x+k_{y,m,n}y+k_{z,m,n}\left|z\right|)}}},$$

where ${\overleftrightarrow{Q}}_{m,n}^{\pm }$ is the discretized form of ${\overleftrightarrow{Q}}^{\pm }$ for k^x,m,n, k^y,m,n, and k^z,m,n. The electric fields are represented, from Eq. (9a), as

(11b)$$U^{\pm }=(U_x^{\pm },U_y^{\pm },U_z^{\pm })=U_0{\displaystyle \sum _{m=-M}^M{\displaystyle \sum _{n=-N}^N{\overleftrightarrow{Q}}_{m,n}^{\pm }\cdot P{e}^{j(k_{x,m,n}x+k_{y,m,n}y+k_{z,m,n}\left|z\right|)}}}for\pm z\ge 0,$$

where U⁰ is $j{\mu }_0{\omega }^2\Delta k_x\Delta k_y/(8{\pi }^2)$. The Fourier coefficients of the x and y components of the positive and negative components U⁺ and U⁻ in Eq. (11b) are rearranged to column-vector forms, $\overrightarrow{U}and\overleftarrow{U}$, as

(12a)$$\overrightarrow{U}=U_0\left(\begin{array}{c}{\left[{\overleftrightarrow{Q}}_{m,n}^{+}\cdot P\right]}_y\\ {\left[{\overleftrightarrow{Q}}_{m,n}^{+}\cdot P\right]}_x\end{array}\right),$$

(12b)$$\overleftarrow{U}=U_0\left(\begin{array}{c}{\left[{\overleftrightarrow{Q}}_{m,n}^{-}\cdot P\right]}_y\\ {\left[{\overleftrightarrow{Q}}_{m,n}^{-}\cdot P\right]}_x\end{array}\right),$$

where ${[{\overleftrightarrow{Q}}_{m,n}^{\pm }\cdot p]}_{y(x)}$ is the raw-major-order-column-vector representation of the y(x)-directional component of ${\overleftrightarrow{Q}}_{m,n}^{\pm }\cdot p$. The (m, n)th component of ${[{\overleftrightarrow{Q}}_{m,n}^{\pm }\cdot p]}_{y(x)}$ is located in (m + M)(2N + 1) + n + N + 1, where m and n range as −M ≤ m ≤ M and −N ≤ n ≤ N respectively. Therefore, by substituting Eqs. (12a) and (12b) into Eqs. (4a) and (4b) and Eqs. (5a) and (5b), we can obtain the modal spectra of the radiated field, $\overrightarrow{T}and\overleftarrow{T}$, and the modal spectra of the internal optical field, $\overrightarrow{C}and\overleftarrow{C}$, in the photonic structure induced by a dipole polarization source, as follows:

(13a)$$\overrightarrow{T}={\overrightarrow{T}}_p\overrightarrow{U}+{\overrightarrow{T}}_n\overleftarrow{U},$$

(13b)$$\overrightarrow{T}={\overleftarrow{T}}_p\overrightarrow{U}+{\overleftarrow{T}}_n\overleftarrow{U},$$

(13c)$$\overrightarrow{C}={\overrightarrow{C}}_p\overrightarrow{U}+{\overrightarrow{C}}_n\overleftarrow{U},$$

(13d)$$\overleftarrow{C}={\overleftarrow{T}}_p\overrightarrow{U}+{\overleftarrow{T}}_n\overleftarrow{U},$$

With the obtained modal spectra, the optical field distribution is reconstructed according to Eqs. (1a)–(1d).

In Fig. 3, the exemplary optical field distribution is calculated by the aforementioned bidirectional S-matrix method. For simplicity, the simulation is conducted for a two-dimensional dipole source (N = 0), which presents the results of the proposed formulation more clearly. In Figs. 3(a) and 3(b), the x-directional electric field distribution E_x and its intensity profile |E| in the xz plane, which are excited by a two-dimensional dipole source of wavelength 532 nm and polarization P = (1, 0, 0), are presented.

Figure 3. Gibbs-induced error in the field distribution obtained without the sigma operator. (a) E_x and (b) |E| in free space.

In both Figs. 3(a) and 3(b), the field profiles near z = 0 show an abrupt high amplitude, which is considered to be an unphysical numerical error originated from the Gibbs phenomenon. The Gibbs phenomenon, associated with convergence in field and structure in FMM, is a well-known issue [10]. However, in this case the Gibbs phenomenon is related to the numerical representation of the optical source.

With respect to numerical analysis, when a dipole source of subwavelength dimension is represented by the discrete angular spectrum (i.e. the discrete Weyl representation), the Gibbs phenomenon becomes problematic and induces significant errors, as shown in Fig. 3. Therefore, it is important to develop an effective method to reduce or eliminate those Gibbs-induced errors in the discrete Weyl representation, Eq. (11b).

In terms of signal processing, to reduce the Gibbs phenomenon without losing physical meaning, the specific window is wrapped to the angular-spectrum profile of the discrete Weyl representation. We choose the sigma filter [11], σ_mn, which is defined by

(14a)$${\sigma }_{mn}=\sin c(m/M)\sin c(n/N),$$

and try to apply σ_mn to the discrete Weyl representation of Eq. (11b) as

(14b)$$U^{\pm }=(U_x^{\pm },U_y^{\pm },U_z^{\pm })=U_0{\displaystyle \sum _{m=-M}^M{\displaystyle \sum _{n=-N}^N{\sigma }_{mn}{\overleftrightarrow{Q}}_{m,n}^{\pm }\cdot P{e}^{j(k_{x,m,n}x+k_{y,m,n}y+k_{z,m,n}\left|z\right|)}}}for\pm z\ge 0,$$

The optical field spectra of the discrete Weyl representation refined with the sigma filter, which are comparable to Eqs. (12a)–(12b), read as

(15a)$$\underset{\_}{\overrightarrow{U}}=U_0\left(\begin{array}{c}{\left[{\sigma }_{mn}{\overleftrightarrow{Q}}_{m,n}^{+}\cdot P\right]}_y\\ {\left[{\sigma }_{mn}{\overleftrightarrow{Q}}_{m,n}^{+}\cdot P\right]}_x\end{array}\right),$$

(15b)$$\underset{\_}{\overleftarrow{U}}=U_0\left(\begin{array}{c}{\left[{\sigma }_{mn}{\overleftrightarrow{Q}}_{m,n}^{-}\cdot P\right]}_y\\ {\left[{\sigma }_{mn}{\overleftrightarrow{Q}}_{m,n}^{-}\cdot P\right]}_x\end{array}\right).$$

Figure 4 compares the field profile of the z-directional electric field at r = (x, 0, 0) calculated from the discrete Weyl representation, Eqs. (12a) and (12b), to that from the refined Weyl representation, Eqs. (15a) and (15b). The unphysical, highly oscillatory features in E_z, ascribed to the Gibbs phenomenon, are observed in Fig. 4(a), while the effect of the sigma filter is clearly seen in the field profile in Fig. 4(b), where the use of the sigma filter dampens the Gibbs-induced numerical errors effectively.

Figure 4. Discrete Weyl representations of U_z (a) without the sigma filter, and (b) with the sigma filter. The sigma filter’s reduction of Gibbs-induced error is clearly observed in (b).

As presented in Fig. 5, the optical field distributions obtained with the sigma filter demonstrate that the Gibbs-induced errors have been removed, and continuity in the field distribution is confirmed at r = (x, 0, 0), which corresponds to the interface B₀ in Fig. 2(a). We compare three-dimensional (3D) field distributions of a dipole source at r = (0, 0, 0) in free space.

Figure 5. Reduction of the Gibbs-induced error in the field distribution with the sigma filter. (a) E_x and (b) |E| in free space.

A three-dimensional perfectly matched layer (PML) [12] is used to model this structure, as shown in in Fig. 6(a). Without the sigma filter [Fig. 6(b)], the Gibbs phenomenon occurs across the whole xy plane (z = 0). In contrast, the field distribution with the sigma filter exhibits significantly reduced Gibbs error [Fig. 6(c)], indicating that the proposed discrete Weyl representation is an appropriate dipole-source model for the FMM. As an additional example, the optical dipole radiation near a 40-nm-thick silver film placed at z = −200 nm [Fig. 6(d)]. is simulated with and without the sigma filter, and the results are compared in Figs. 6(e) and 6(f), verifying the proposed FMM scheme with the proposed discrete Weyl representation and bidirectional S-matrix method.

Figure 6. 3D optical dipole radiation in free space: (a) schematic and E_x field distributions, (b) without, and (c) with the sigma filter. 2D dipole radiation near a 40-nm-thick Ag film: (d) schematic and E_x field distributions. (e) without, and (f) with the sigma filter. A surface plasmon is observed in the back side of the Ag film.

We have formulated the basis for FMM analysis of general aperiodic source-embedding photonic structures: the bidirectional S-matrix method with the discrete Weyl representation and perfectly matched layer. The proposed scheme can provide an FMM-based modal-analysis framework for the design and analysis of advanced photonic structures with embedded optical sources.

Funding was provided by the National Research Foundation of Korea (NRF) (NRF-2019R1A2C1010243).