G-0K8J8ZR168
검색
검색 팝업 닫기

Ex) Article Title, Author, Keywords

## Article

Curr. Opt. Photon. 2021; 5(6): 597-605

Published online December 25, 2021 https://doi.org/10.3807/COPP.2021.5.6.597

## Fourier Modal Method for Optical Dipole Radiation in Photonic Structures

Sungjae Park1, Joonku Hahn2, Hwi Kim1

1Department of Electronics and Information Engineering, Korea University Sejong Campus, Sejong 30019, Korea
2School of Electronic and Electrical Engineering, Kyungpook National University, Daegu 41566, Korea

Corresponding author: *hwikim@korea.ac.kr, ORCID 0000-0002-4283-8982

Received: June 17, 2021; Revised: July 27, 2021; Accepted: September 23, 2021

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

An extended Fourier modal method (FMM) for optical dipole radiation in three-dimensional photonic structures is proposed. The core elements of the proposed FMM are the stable bidirectional scattering matrix algorithm for modeling internal optical emission, and a novel optical-dipole-source model that prevents numerical errors induced by the Gibbs phenomenon. Through the proposed scheme, the FMM is extended to model a wide range of source-embedded photonic structures.

Keywords: Fourier modal method, Numerical modeling

OCIS codes: (000.3860) Mathematical methods in physics; (050.1755) Computational electromagnetic methods; (050.1960) Diffraction theory

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].

### II. Bidirectional scattering-matrix formulation

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 Lk and Bk respectively. The bidirectional multiblock structure is composed of three separate parts. The active-source block denoted by L0 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 L0 with a permittivity of εr. The left and right parts of the active-source block are multiple stacked single blocks Lk, with the indices −(NL + 1) ≤ k ≤ −1 and 1 ≤ k ≤ (NR + 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 L0 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 L0 should be εr = 1 for a dipole embedded in free space, whereas we need to set the permittivity of L0 as εr(r0)when the material permittivity at the point r = r0 is εr(r0). 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)+,andH(k)(g) represented by the pseudo-Fourier series are prepared for all single blocks [1]. E(k)(g)±andH(k)(g)±
are the electric field and magnetic field of the gth positive (negative) eigenmode in the kth block Lk. 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

EH=g=1M+C (k),g+ E (k) (g)+ H (k) (g)++g=1MC (k),g E (k) (g) H (k) (g),forBkzB(k1),1kNL,

EH=g=1M+C (k),g+ E (k) (g)+ H (k) (g)++g=1MC (k),g E (k) (g) H (k) (g),forBk1zBk,1kNR,

where C(k),g±andC(k),g± 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,

EH=g=1MT (NL1),g E (NL1) (g) H (NL1) (g),forzBNL,

EH=g=1M+C (NR+1),g+ E (NR+1) (g)+ H (NR+1) (g)+,forBNRz,

where T(NL1),gandT(NR+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 T(NL1),g are denoted by T, whereas those of the right block T(NR+1),g are represented by T. The coupling coefficients of the left and right sides, C(k),g±andC(k),g±, are represented by the collective coupling-coefficient operators C(L)andC(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

S(L)= T (L) R (L) R (L) T (L),

S(R)= T (R) R (R) R (R) T (R).

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

Tp=T(R)(I R (L) R (R) )1,

Tp=T(L)(I R (R) R (L) )1R(R),

Tn=T(R)(I R (L) R (R) )1R(L),

Tn=T(L)(I R (R) 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 UandU respectively. Considering the first component TpU, we can expand it as

TpU=T(R) I R (L) R (R) 1U=T(R)U+T(R)R(L)R(R)U+T(R) R (L) R (R) 2U+T(R) R (L) R (R) 3U.

The first term T(R)U is the direct transmission of U through the right block, and the second term T(R)R(L)R(R)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 UandU as

T= T pU+ T nU,

T= T pU+ T nU.

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

CP=C(R) IR (L) R (R) 1,

Cp=C(L) IR (R) R (L) 1R(R),

Cn=C(R) IR (L) R (R) 1R(L),

Cn=C(L) IR (R) R (L) 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,

C= C pU+ C nU,

C= C pU+ C nU.

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, UandU, as explained in the next section.

### III. Aperiodic dipole-radiation model for the Fourier modal method

To represent the dipole-radiation field, a proper source operator model, UandU , 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 UandU 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 L0.

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

×E=jωμ0H+jωμ0Mδ(r),

×H=jωε0εrEjωPδ(r),

where P and M are the external polarization and magnetization vectors respectively. ε0 and μ0 denote respectively the electric permittivity and magnetic permeability in free space. εr is the relative permittivity of L0 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 L0 by setting the magnetization vector M = 0. The dyadic Green function of the optical dipole is expressed by the Weyl representation [8],

G±(r)=jεr2π Q ± kz ej(kxx+kyy+kzz)dkxdky,

where kz2=k02εrkx2ky2. Q+andQ 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,

Q±=1p2pq±pmpq1q2±qm±mp±mq1m2for±z0,

where p, q, and m are given by

(p,q,m)=(kx/(k0εr),ky/(k0εr),1p2q2).

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

E=μ0ω24πεrG±(r)P= μ0ω24πεrG +(r)Pforz0 μ0ω24πεrG (r)Pforz<0,

H=jωμ01×E

The radiation power of a dipole, P = (0, 0, pez), in a homogeneous medium with refractive index n=εr is n times larger than that in free space, as

Wn=ω2ImEP*r=0=n μ0 4π ω4 3cpez2=nWfree.

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

H=ω02ε0μ04πG±(r)Mfor±z0,

E=jωε0εr1×H.

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

Wn=ωμ02Im HM* r=0=n3μ04πω4 3c3 mez2.

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

G±(r)=ΔkxΔkyjεr2π m=MMn=NNQ m,n±e j( kx,m,nx+ ky,m,ny+ kz,m,nz),

where Qm,n± is the discretized form of Q± for kx,m,n, ky,m,n, and kz,m,n. The electric fields are represented, from Eq. (9a), as

U±=(Ux±,Uy±,Uz±)=U0m=MMn=NN Qm,n ±Pej(kx,m,nx+ky,m,ny+kz,m,nz)for±z0,

where U0 is jμ0ω2ΔkxΔky/(8π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, UandU, as

U=U0 Q m,n+Py Q m,n+Px,

U=U0 Q m,nPy Q m,nPx,

where [Qm,n±p]y(x) is the raw-major-order-column-vector representation of the y(x)-directional component of Qm,n±p. The (m, n)th component of [Qm,n±p]y(x) is located in (m + M)(2N + 1) + n + N + 1, where m and n range as −MmM and −NnN 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, TandT, and the modal spectra of the internal optical field, CandC, in the photonic structure induced by a dipole polarization source, as follows:

T= T pU+ T nU,

T= T pU+ T nU,

C= C pU+ C nU,

C= T pU+ T nU,

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 Ex 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) Ex 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

σmn=sinc(m/M)sinc(n/N),

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

U±=(Ux±,Uy±,Uz±)=U0m=MMn=NN σmnQm,n ±Pej(kx,m,nx+ky,m,ny+kz,m,nz)for±z0,

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

U_=U0 σ mnQ m,n+Py σ mnQ m,n+Px,

U_=U0 σ mnQ m,nPy σ mnQ m,nPx.

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 Ez, 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 Uz (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 B0 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) Ex 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 Ex field distributions, (b) without, and (c) with the sigma filter. 2D dipole radiation near a 40-nm-thick Ag film: (d) schematic and Ex 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).

1. H. Kim, J. Park and B. Lee in Fourier Modal Method and Its Applications in Computational Nanophotonics, (CRC Press, FL, 2012).
2. I. M. Fradkin, S. A. Dyakov and N. A. Gippius, “Nanoparticle lattices with bases: Fourier modal method and dipole approximation,” Phys. Rev. B 102, 045432 (2020).
3. A. B. Petrin, “An elementary dipole emitter located either on the boundary of a layered structure or inside it,” Opt. Spectrosc. 128, 1809-1827 (2020).
4. K. Postava, T. Fördös, H. Jaffrès, L. Halagačka, H. J. Drouhin and J. Pištora, “Modeling of anisotropic grating structures with active dipole layers,” Proc. SPIE 9516, 95160O (2015).
5. D. Delbeke, P. Bienstman, R. Bockstaele and R. Baets, “Rigorous electromagnetic analysis of dipole emission in periodically corrugated layers: the grating-assisted resonant-cavity light-emitting diode,” J. Opt. Soc. Am. A 19, 871-880 (2002).
6. J. M. Jarem, “Rigorous coupled-wave-theory analysis of dipole scattering from a three-dimensional, inhomogeneous, spherical dielectric, and permeable system,” in IEEE Trans. Microw. Theory Tech. 45, 1193-1203 (1997).
7. H. Kim and B. Lee, “Mathematical modeling of crossed nanophotonic structures with generalized scattering-matrix method and local Fourier modal analysis,” J. Opt. Soc. Am. B 25, 518-544 (2008).
8. H. Kim, I.-M. Lee and B. Lee, “Extended scattering-matrix method for efficient full parallel implementation of rigorous coupled-wave analysis,” J. Opt. Soc. Am. A 24, 2313-2327 (2007).
9. L. Novotny and B. Hecht in Principles of Nano-Optics, 1st ed, (Cambridge University Press, Cambridge, 2006).
10. H. Kim, G.-W. Park and C.-S. Kim, “Investigation on the convergence behavior with fluctuation features in the Fourier modal analysis of a metallic grating,” J. Opt. Soc. Korea 16, 196-202 (2012).
11. C. Lanczos in Discourse on Fourier series, (Oliver and Bovd, Edinburgh, 1966).
12. E. Silberstein, P. Lalanne, J.-P. Hugonin and Q. Cao, “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A 18, 2865-2875 (2001).

### Article

#### Article

Curr. Opt. Photon. 2021; 5(6): 597-605

Published online December 25, 2021 https://doi.org/10.3807/COPP.2021.5.6.597

## Fourier Modal Method for Optical Dipole Radiation in Photonic Structures

Sungjae Park1, Joonku Hahn2, Hwi Kim1

1Department of Electronics and Information Engineering, Korea University Sejong Campus, Sejong 30019, Korea
2School of Electronic and Electrical Engineering, Kyungpook National University, Daegu 41566, Korea

Correspondence to:*hwikim@korea.ac.kr, ORCID 0000-0002-4283-8982

Received: June 17, 2021; Revised: July 27, 2021; Accepted: September 23, 2021

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

An extended Fourier modal method (FMM) for optical dipole radiation in three-dimensional photonic structures is proposed. The core elements of the proposed FMM are the stable bidirectional scattering matrix algorithm for modeling internal optical emission, and a novel optical-dipole-source model that prevents numerical errors induced by the Gibbs phenomenon. Through the proposed scheme, the FMM is extended to model a wide range of source-embedded photonic structures.

Keywords: Fourier modal method, Numerical modeling

### I. INTRODUCTION

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].

### II. Bidirectional scattering-matrix formulation

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 Lk and Bk respectively. The bidirectional multiblock structure is composed of three separate parts. The active-source block denoted by L0 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 L0 with a permittivity of εr. The left and right parts of the active-source block are multiple stacked single blocks Lk, with the indices −(NL + 1) ≤ k ≤ −1 and 1 ≤ k ≤ (NR + 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 L0 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 L0 should be εr = 1 for a dipole embedded in free space, whereas we need to set the permittivity of L0 as εr(r0)when the material permittivity at the point r = r0 is εr(r0). 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)± and H(k)(g)±$
are the electric field and magnetic field of the gth positive (negative) eigenmode in the kth block Lk. 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

$EH=∑g=1M+C← (−k),g+ E (−k) (g)+ H (−k) (g)++∑g=1M−C← (−k),g− E (−k) (g)− H (−k) (g)−, for B−k≤z≤B−(k−1), 1≤k≤NL,$

$EH=∑g=1M+C→ (k),g+ E (k) (g)+ H (k) (g)++∑g=1M−C→ (k),g− E (k) (g)− H (k) (g)−, for Bk−1≤z≤Bk, 1≤k≤NR,$

where $C←(−k),g± and C→(k),g±$ 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,

$EH=∑g=1M−T← (−NL−1),g E (−NL−1) (g)− H (−NL−1) (g)−, for z≤B−NL,$

$EH=∑g=1M+C→ (NR+1),g+ E (NR+1) (g)+ H (NR+1) (g)+, for BNR≤z,$

where $T←(−NL−1),g and T→(NR+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 $T←(−NL−1),g$ are denoted by $T←$, whereas those of the right block $T→(NR+1), g$ are represented by $T→$. The coupling coefficients of the left and right sides, $C←(−k),g± and C→(k),g±$, are represented by the collective coupling-coefficient operators $C←(L) and 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

$S(L)= T→ (L) R→ (L) R← (L) T← (L),$

$S(R)= T→ (R) R→ (R) R← (R) T← (R).$

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

$T→p=T→(R)(I− R → (L) R← (R) )−1,$

$T←p=T←(L)(I− R ← (R) R→ (L) )−1R←(R),$

$T→n=T→(R)(I− R → (L) R← (R) )−1R→(L),$

$T←n=T←(L)(I− R ← (R) 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 $U→ and U←$ respectively. Considering the first component $T→pU→$, we can expand it as

$T→pU→=T→(R) I− R→ (L) R← (R) −1U→=T→(R)U→+T→(R)R→(L)R←(R)U→+T→(R) R→ (L) R← (R) 2U→+T→(R) R→ (L) R← (R) 3U→⋯.$

The first term $T→(R)U→$ is the direct transmission of $U→$ through the right block, and the second term $T→(R)R→(L)R←(R)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 $U→ and U←$ as

$T→= T → pU→+ T → nU←,$

$T←= T ← pU→+ T ← nU←.$

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

$C→P=C→(R) I−R →(L) R← (R) −1,$

$C←p=C←(L) I−R ←(R) R→ (L) −1R←(R),$

$C→n=C→(R) I−R →(L) R← (R) −1R→(L),$

$C←n=C←(L) I−R ←(R) R→ (L) −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,

$C→= C → pU→+ C → nU←,$

$C←= C ← pU→+ C ← nU←.$

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, $U→ and U←$, as explained in the next section.

### III. Aperiodic dipole-radiation model for the Fourier modal method

To represent the dipole-radiation field, a proper source operator model, $U→ and 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 $U→ and 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 L0.

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

$∇×E=jωμ0H+jωμ0Mδ(r),$

$∇×H=−jωε0εrE−jωPδ(r),$

where P and M are the external polarization and magnetization vectors respectively. ε0 and μ0 denote respectively the electric permittivity and magnetic permeability in free space. εr is the relative permittivity of L0 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 L0 by setting the magnetization vector M = 0. The dyadic Green function of the optical dipole is expressed by the Weyl representation [8],

$G↔±(r)=jεr2π∫−∞∞ ∫ −∞ ∞ Q↔ ± kz ej(kxx+kyy+kzz)dkxdky,$

where $kz2=k02εr−kx2−ky2$. $Q↔+ and 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,

$Q↔±=1−p2−pq±pm−pq1−q2±qm±mp±mq1−m2 for ±z≥0,$

where p, q, and m are given by

$(p,q,m)=(kx/(k0εr),ky/(k0εr),1−p2−q2).$

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

$E=μ0ω24πεrG↔±(r)⋅P= μ0ω24πεrG↔ +(r)⋅P for z≥0 μ0ω24πεrG↔ −(r)⋅P for z<0,$

$H=jωμ0−1∇×E$

The radiation power of a dipole, P = (0, 0, pez), in a homogeneous medium with refractive index $n=εr$ is n times larger than that in free space, as

$Wn=ω2ImE⋅P*r=0=n μ0 4π ω4 3cpez2=nWfree.$

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

$H=ω02ε0μ04πG↔±(r)⋅M for ±z≥0,$

$E=−jωε0εr−1 ∇×H.$

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

$Wn=ωμ02Im H⋅M* r=0=n3μ04πω4 3c3 mez2.$

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

$G↔±(r)=ΔkxΔkyjεr2π∑ m=−MM∑n=−NNQ↔ m,n±e j( kx,m,nx+ ky,m,ny+ kz,m,nz),$

where $Q↔m,n±$ is the discretized form of $Q↔±$ for kx,m,n, ky,m,n, and kz,m,n. The electric fields are represented, from Eq. (9a), as

$U±=(Ux±,Uy±,Uz±)=U0∑m=−MM∑n=−NN Q↔m,n ±⋅Pej(kx,m,nx+ky,m,ny+kz,m,nz) for ±z≥0,$

where U0 is $jμ0ω2ΔkxΔky/(8π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, $U→ and U←$, as

$U→=U0 Q↔ m,n+⋅Py Q↔ m,n+⋅Px,$

$U←=U0 Q↔ m,n−⋅Py Q↔ m,n−⋅Px,$

where $[Q↔m,n±⋅p]y(x)$ is the raw-major-order-column-vector representation of the y(x)-directional component of $Q↔m,n±⋅p$. The (m, n)th component of $[Q↔m,n±⋅p]y(x)$ is located in (m + M)(2N + 1) + n + N + 1, where m and n range as −MmM and −NnN 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, $T→ and T←$, and the modal spectra of the internal optical field, $C→ and C←$, in the photonic structure induced by a dipole polarization source, as follows:

$T→= T → pU→+ T → nU←,$

$T→= T← pU→+ T← nU←,$

$C→= C → pU→+ C → nU←,$

$C←= T ← pU→+ T ← nU←,$

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 Ex 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) Ex 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

$σmn=sinc(m/M)sinc(n/N),$

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

$U±=(Ux±,Uy±,Uz±)=U0∑m=−MM∑n=−NN σmnQ↔m,n ±⋅Pej(kx,m,nx+ky,m,ny+kz,m,nz) for ±z≥0,$

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

$U→_=U0 σ mnQ↔ m,n+⋅Py σ mnQ↔ m,n+⋅Px,$

$U←_=U0 σ mnQ↔ m,n−⋅Py σ mnQ↔ m,n−⋅Px.$

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 Ez, 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 Uz (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 B0 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) Ex 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 Ex field distributions, (b) without, and (c) with the sigma filter. 2D dipole radiation near a 40-nm-thick Ag film: (d) schematic and Ex field distributions. (e) without, and (f) with the sigma filter. A surface plasmon is observed in the back side of the Ag film.

### IV. CONCLUSION

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.

### ACKNOWLEDGMENT

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

### Fig 1.

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.
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### Fig 2.

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.
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### Fig 3.

Figure 3.Gibbs-induced error in the field distribution obtained without the sigma operator. (a) Ex and (b) |E| in free space.
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### Fig 4.

Figure 4.Discrete Weyl representations of Uz (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).
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### Fig 5.

Figure 5.Reduction of the Gibbs-induced error in the field distribution with the sigma filter. (a) Ex and (b) |E| in free space.
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### Fig 6.

Figure 6.3D optical dipole radiation in free space: (a) schematic and Ex field distributions, (b) without, and (c) with the sigma filter. 2D dipole radiation near a 40-nm-thick Ag film: (d) schematic and Ex field distributions. (e) without, and (f) with the sigma filter. A surface plasmon is observed in the back side of the Ag film.
Current Optics and Photonics 2021; 5: 597-605https://doi.org/10.3807/COPP.2021.5.6.597

### References

1. H. Kim, J. Park and B. Lee in Fourier Modal Method and Its Applications in Computational Nanophotonics, (CRC Press, FL, 2012).
2. I. M. Fradkin, S. A. Dyakov and N. A. Gippius, “Nanoparticle lattices with bases: Fourier modal method and dipole approximation,” Phys. Rev. B 102, 045432 (2020).
3. A. B. Petrin, “An elementary dipole emitter located either on the boundary of a layered structure or inside it,” Opt. Spectrosc. 128, 1809-1827 (2020).
4. K. Postava, T. Fördös, H. Jaffrès, L. Halagačka, H. J. Drouhin and J. Pištora, “Modeling of anisotropic grating structures with active dipole layers,” Proc. SPIE 9516, 95160O (2015).
5. D. Delbeke, P. Bienstman, R. Bockstaele and R. Baets, “Rigorous electromagnetic analysis of dipole emission in periodically corrugated layers: the grating-assisted resonant-cavity light-emitting diode,” J. Opt. Soc. Am. A 19, 871-880 (2002).
6. J. M. Jarem, “Rigorous coupled-wave-theory analysis of dipole scattering from a three-dimensional, inhomogeneous, spherical dielectric, and permeable system,” in IEEE Trans. Microw. Theory Tech. 45, 1193-1203 (1997).
7. H. Kim and B. Lee, “Mathematical modeling of crossed nanophotonic structures with generalized scattering-matrix method and local Fourier modal analysis,” J. Opt. Soc. Am. B 25, 518-544 (2008).
8. H. Kim, I.-M. Lee and B. Lee, “Extended scattering-matrix method for efficient full parallel implementation of rigorous coupled-wave analysis,” J. Opt. Soc. Am. A 24, 2313-2327 (2007).
9. L. Novotny and B. Hecht in Principles of Nano-Optics, 1st ed, (Cambridge University Press, Cambridge, 2006).
10. H. Kim, G.-W. Park and C.-S. Kim, “Investigation on the convergence behavior with fluctuation features in the Fourier modal analysis of a metallic grating,” J. Opt. Soc. Korea 16, 196-202 (2012).
11. C. Lanczos in Discourse on Fourier series, (Oliver and Bovd, Edinburgh, 1966).
12. E. Silberstein, P. Lalanne, J.-P. Hugonin and Q. Cao, “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A 18, 2865-2875 (2001).

Wonshik Choi,
Editor-in-chief