Fourier finite element modeling of light emission in waveguides: 2.5-dimensional FEM approach

We present a Fourier finite element modeling of light emission of dipolar emitters coupled to infinitely long waveguides. Due to the translational symmetry, the three-dimensional (3D) coupled waveguide-emitter system can be decomposed into a series of independent 2D problems (2.5D), which reduces the computational cost. Moreover, the reduced 2D problems can be extremely accurate, compared to its 3D counterpart. Our method can precisely quantify the total emission rates, as well as the fraction of emission rates into different modal channels for waveguides with arbitrary cross-sections. We compare our method with dyadic Green's function for the light emission in single mode metallic nanowire, which yields an excellent agreement. This method is applied in multi-mode waveguides, as well as multi-core waveguides. We further show that our method has the full capability of including dipole orientations, as illustrated via a rotating dipole, which leads to unidirectional excitation of guide modes. The 2.5D Finite Element Method (FEM) approach proposed here can be applied for various waveguides, thus it is useful to interface single-photon single-emitter in nano-structures, as well as for other scenarios involving coupled waveguide-emitters.

In the filed of waveguides, Chen et al. [24] studied emission rates γ of waveguide modes, using 2D mode profiles by taking advantage of the translation symmetry of the waveguide. However, the presence of the dipolar emitter breaks the translation symmetry of the solution, which gives rise to difficulty in calculating the total emission rate γ total . Thus, in order to obtain γ total , or the fraction ( γ m /γ total ) of emitted photon into a specific waveguide mode, the authors in [24] further developed a 3D FEM to study the total emission rate, consisting of bound modes, radiation modes, and other non-radiative channels. Chen's 3D FEM model, though equipped with non-trivial mode matching technique, is limited by the fact that only one dominatingly guided mode is allowed to be excited. However, multi-mode excitation is needed as the crosssection of the waveguide becomes large, or when random orientated dipoles excite two different polarization modes in an elliptical waveguide [29]. We refer to Ref. [24] for the motivation of our choice of choosing FEM as our method to study the coupled waveguide-emitter system.
In this paper, we lift the restriction of the single dominating mode excitation in [24] by developing a 2.5D FEM approach. The aforementioned difficulty can be solved by fully taking advantage of the point nature of the dipolar emitter as well as the 2D profile of material indices, meaning that the coupling between the waveguide modes and the emitter, which are associated with different propagation constants (k nz ), turn out to be decoupled. This leads to a Fourier finite element implementation of light emission from point dipoles in waveguides. Essentially, the 3D coupled waveguide-emitter problem is decomposed into a series of 2D uncoupled waveguide-emitter problems, each of which is associated with a single k nz . Our method resembles a spectrum presentation of the dyadic Green's function of the waveguide in terms of k nz . The 2.5D FEM approach can handle multi-mode excitation for waveguide with arbitrary cross-section, in which the dipole orientation can be arbitrary as well.
The paper is organized as follows. Section 2 describes the foundation of a 2.5D FEM approach, and the basic theory of light emission. In Section 3, we first provide an example of a single metallic nanowire as a benchmark of the 2.5 FEM approach against the Green's function techniques. We further apply our 2.5D FEM approach to multi-mode and a multi-core waveguide. Finally, we show that the 2.5D FEM approach can be used to realize unidirectional excitation by engineering dipole orientation. Section 4 presents the main conclusions of the paper.

Theoretical foundation of 2.5D FEM approach
In this section, we shall discuss how the 2.5D FEM approach is realized. We consider the coupling between single emitter and an infinite long waveguides. With the combination of Maxwell's equations and the constitutive relations, we have the following wave equation, whereμ µ µ r andε ε ε r (r) are relative permeability and relative dielectric constant respectively, k 0 is vacuum wave number, and p(r 0 , ω) is the dipole moment of the emitter. We introduce a weighting function F(r, ω) to test Eq. (1), and integrate over the entire modeling domain, where ∂V denotes the surface that encloses volume V , and n denotes the outward normal unit vector to the surface of the modeling domain. We want to find E(r, ω) that minimizes Eq. 2 over all test functions. However, Eq. (2) is a triple integral where the length of z axis is infinite, and it is a challenging problem for the domain discretization and construction of interpolation functions. Next, we will show how to decompose the 3D problem into a sequence of 2D problems by Fourier transforming of electric field E(r) and dipole moment p(r 0 ), and by selecting appropriate test function for each mode of the waveguide. For short, p(r 0 , ω), E(r, ω) and F(r, ω) will be denoted as p(r 0 ), E(r) and F(r) onwards. We introduce Fourier expansion over the relevant quantities, i.e., E(r) = 1 2π z (x, y, k 1z ). Due to the invariance of the dielectric function along z direction, i.e., ε(r) = ε(x, y), the resulting formulation becomes a sequence of uncoupled 2D problems (2.5D problem) rather than a sequence of coupled 2D problems (amount to 3D problem). Substituting E(r) and p(r 0 ) into Eq. (2) and selecting a mono-modal test function of the form F n (r) = F n (x, y, z) = F n (x, y, k nz )e ik nz z , we can rewrite the weak formulation as where Test functions (F n ) are selected as Fourier modes, including the waveguide continuum, forming a complete and orthogonal basis in which any arbitrary field can be expanded. In Eq. (2), such typical selection of test function allows us to obtain a trivial integration along z axis without coupling Fourier modes. Subsequently, we obtain a series of decoupled planar integration in terms of E(x, y, k nz ), CurlE(x, y, k nz ), F n (x, y, k nz ) and CurlF n (x, y, k nz ), which contain the same propagation constant k nz . Notably, the 2D integration given by Eq. (3) for each k nz can be solved as a standard 2D boundary-value problem by employing traditional finite element solution procedures, including discretization, selection of basis functions and assembling of the sparse matrix [30]. Coefficients k nz correspond to the Fourier frequencies in the z-direction.
The number of such frequencies is often smaller than the number of basis funcions that one may select on the z-direction on a traditional finite element method, since higher-order methods like Fourier exhibits higher-order convergence. Additionally, the uncoupling among Fourier modes produces great computational savings with respect to traditional finite element methods, since it enables to decouple the original 3D problem into a sequence of uncoupled 2D problems, namely, one per Fourier mode. While a 3D problem can be interpreted from the algebraic point of view as a dense matrix of 2D problems, in the proposed approach, the resulting system becomes a diagonal matrix of 2D problems, which can be solved in a fraction of the time. In Eq. (3), the dipole moment is invariant with respect to the propagating constant, which amounts to a time-harmonic line current source independent of k nz . Fourier frequencies k nz are selected in such a way that the source is well approximated in the frequency domain when using its Fourier series expansion.

Rate of energy dissipation in inhomogeneous environment
According to Poynting's theorem, the radiated power of any current distribution with a harmonic time dependence in a linear medium has to be identical to the rate of energy dissipation P [31][32][33]. If we introduce the dipole's current density such as in Eq. 2, i.e., J(r, ω) = −iωpδ (r − r 0 ), we obtain the emitted power as where E(r 0 ) is evaluated at the dipole's origin r 0 . To formulate the spectrum presentation of the radiated power, we introduce differential emission rate P(k nz ), the integration of which over k nz gives the total emission rate P. We estimate the normalized differential energy dissipation rate P(k nz )/P 0 for a particular k nz , where P(k nz ) is given by and P 0 = +∞ −∞ P 0 (k nz )dk nz is the emitted power by the same dipole source into all modes of a homogeneous material where the emitter is seated.
The dyadic Green's function provides us with an alternative approach to study spontaneous decay rate γ of a two-level quantum system in an arbitrary optical environment [34,35]. The differential spontaneous decay rate γ(k nz ) into each mode with propagating constant k nz can be normalized by the free-space spontaneous decay rate γ 0 = ω 3 0 |p| 2 3πε 0h c 3 , and written in terms of spectrum expanding of Green's tensor as with n p being the unit vector in the direction of the dipole moment p , and G(r 0 , r 0 , ω; k nz ) is the Green's tensor at the dipole position. The normalized differential spontaneous decay rate γ(k nz ) γ 0 into each mode calculated by the dyadic Green's function method is equivalent to the normalized differential energy dissipation rate P(k nz ) P 0 calculated by the 2.5 FEM approach, i.e., γ(k nz ) γ 0 = P(k nz ) P 0 . We also introduce the β factor to describe the fraction of the emitted energy that is coupled to the guided mode. This factor is defined as β m = P m P total , where P m is the dissipation rate into the m guided mode obtained by intergrating P(k nz ) around the guided mode resonance, and P total is the total dissipation rate into all modes.

Single mode cylindrical metallic nanowire
To illustrate the validity of the 2.5 FEM approach, we compare our numerical results with analytical solutions obtained by using the Green's dyadic method. The coupling between waveguide and emitter is sketched in Fig. 1(a), in which we consider a dipole (frequency f 0 = 300T Hz) placed close to the surface of an infinitely long metallic nanowire (R = 20nm) coupled into the surface plasmons of the conducting waveguide. The wire with electric permittivity ε 1 = −50 + 0.6 j is surrounded by a material with ε 2 = 2, and the dipole moment is oriented along the x axis. The inset of Fig. 1(a) shows the profile of a 2D cross-section of the model. The normalized differential emission rate with modal index n e f f , i.e., P(k nz )/P 0 , is calculated by our 2.5 FEM approach. As can be seen from Fig. 1(b), our FEM result show an excellent agreement with that calculated by the dyadic Green's function method. Moreover, there are three critical features shown in Fig. 1(b), which match extremely well when employing the two independent methods. First, there exist two peaks when the effective indexes is around ±2.29, which are essentially the effective index of the fundamental mode of the waveguide (n eig = 2.29). The two methods also provide the same maximal value of the peak for waveguide with the same losses, as Fig. 1(b) shows. Second, when n e f f > √ ε 2 = √ 2, the plasmons are evanescent in the radial direction, and the dissipation rate of the dipole rapidly declines. Taking into account the medium absorption, i.e., a small value of Im(ε 1 ), the Dirac-delta function spectrum of the bound mode is approximately broadened into a Lorentz-like lineshape, as shown in the inset of Fig. 1 (b), the half-width at half maximum (HWHM) is determined by Im(ε 1 ). That is to say, in the evanescent region, the main contribution to the dipole emission rate comes from excitation of plasmonic modes on the nanowire, which leads to an enhancement of the decay rate of the emitter and channeling of its emission into a single propagating plasmonic mode. Third, it is worth noticing, in the traveling-wave part ( n e f f < √ ε 2 = √ 2), the contribution of the emission rate comes from the freely propagating photons. Thus, we conclude that our 2.5 FEM approach is verified from the comparison with the independent Green's function technique. To compare the computational efficiencies of 2.5D FEM with its 3D counterpart, we plot Fig. 1(c) in which we show the length dependence of computation time for the 3D counterpart by using mode matching method. The 2.5D FEM approach can decompose a 3D problem into a series of decoupled 2D problems, which can be computed in parallel so the computational efficiency is significantly increased. According to the discussion of [24], the variation in the total decay rate is reduced by increasing the length of nanowire. The relative error at L = 4200nm is about 0.26%, and the computing time is 825.3s, which is calculated by 3D FEM. However, with appropriate selections of k nz and its calculating ranges (we selected 228 k nz with the range of [−5k 0 , 5k 0 ]), the relative error is about 0.19%, along with a calculating time of 172s calculated by 2.5D FEM in the sequential (no parallel) version, and 57s with 3 cores. Thus, the 2.5D FEM indeed needs less computation time and provides higher precision compared with its 3D counterpart.

Multi-mode elliptical nanowire
As a known fact, it is difficult to study the emission of single emitter to elliptical waveguides or others with irregular cross section [29]. However, such problems can be easily handled by using our 2.5D FEM approach. In the following, we consider a semiconductor photonic nanowire made of GaAs (refractive index n 1 = 3.45) and surrounded by an air cladding (n 2 = 1). The wire has an elliptical cross section, with major diameter 360nm, twice the length of the minor diameter, and the emitter is placed on the long axis of a elliptical section, as the inset of Fig. 2(a) shows.
We obtain the normalized differential dissipation rate of the dipole P(k nz ) P 0 as a function of the different modes from Fig. 2(a). The four guided modes are labelled as mode A, B, C, and D respectively. We introduce a little material loss Im(ε 1 ) to the wire to better visualize the emission peaks given by the four modes in our numerical calculation. We display the fraction of the emission (β factor) coupled into each guided mode as a function of the distance between dipole position and the center of ellipse section in Fig. 2(b). Figures 2(c)-2(f) discribe the electric field distribution of the four guided modes. As shown in Fig. 2(e), the electric field of mode C (n 3 e f f ) in the center of waveguide goes along the y axis, which implies the mode C can only be stimulated by a central dipole with y-component. Similarly, mode D (n 4 e f f ) can not be stimulated by a central dipole without x-component. Hence, in order to excite all the guided modes of the waveguide, We set the emitter with two equal dipole components, i.e., p p p =xp 0 +ŷp 0 . The electric field distributions of mode C and D show that the field intensity decrease gradually as the distance D increase, as shown in Figs. 2(e)-2(f). Hence, the coupling between the dipole and the modes will decrease for large D, which is consistent with the results shown on Fig. 2(b). The field distributions of the other two guided mode B and mode C indicate that field intensity increase to a its maximum as the distance D increases and crosses an optimal distance, as illustrated in Figs. 2(c)-2(d), which is also consistent with the tendency of coupling efficiency change shown in Fig. 2(b).

Multi-core nanowire
We proceed to discuss the application of our 2.5D FEM approach to multi-core waveguide, which has the potential to expand the information capacity in our current optical communication network, as discussed by Tu et al. [36,37]. We propose a multi-core plasmonic waveguide, and place the emitter in the middle of two single-mode metallic nanowires (r 1 = 20nm) and a multimode nanowire (r 2 = 120nm). The metallic cores with ε 1 = −50 + 0.6 j are surrounded by a material with ε 2 = 2, and the dipole has two equal dipole moment components along x and y axis, as shown in the inset of Fig. 3(a). Figure. 3(a) shows the normalized differential emission rate P(k nz )/P 0 as a function of the mode effective index of the waveguide. The multi-core waveguide supports four guided modes, labeled by A, B, C, and D. Figures 3(b)-3(e) show the electric filed distribution |E E E| corresponding to the four excited plasmonic modes, and the arrows denote the intensity and direction of the electric filed at the dipole position. Obviously, the mode D with the maximum |E E E| at the dipole position indicates that the coupling efficiency of the dipole to the guided mode is highest, as confirmed by the maximal value of the β factor (β D ). On the other hand, the minimum |E E E| of mode B indicates the lowest coupling efficiency, and corresponds to the minimum value of the β factor (β B ). The β factors of Fig. 3(b)-3(e) also illustrate that the energy mainly dissipates into the four guided plasmonic modes rather than other evanescent waves or freely propagating photons. It is important to point out that the merit of the 2.5D FEM is that we can obtain the emission rate coupled into arbitrary modal channels for waveguides with arbitrary cross-sections, as well as obtain each mode coupling efficiency for arbitrary dipole position.

Unidirectional excitation of electromagnetic guided mode
In this section, we discuss the flexibility of dipole orientation in our 2.5D FEM approach and possible applications. As a concrete example, we study the unidirectional excitation of bound modes using rotating dipoles [38][39][40][41][42]. For the non-rotating dipoles in the aforementioned discussions, waveguides support guided modes that can be excited in pairs, which share the same effective index n e f f but propagate in opposite directions. Interestingly, it is found that rotating