Generalized mode solver for plasmonic transmission lines embedded in layered media based on the Method of Moments

Accepted Manuscript Generalized mode solver for plasmonic transmission lines embedded in layered media based on the method of moments Mai O. Sallam, Guy A.E. Vandenbosch, Georges Gielen, Ezzeldin A. Soliman

PII: DOI: Reference:

S0010-4655(18)30207-8 https://doi.org/10.1016/j.cpc.2018.06.001 COMPHY 6531

To appear in:

Computer Physics Communications

Received date : 10 July 2017 Revised date : 6 May 2018 Accepted date : 1 June 2018 Please cite this article as: M.O. Sallam, G.A.E. Vandenbosch, G. Gielen, E.A. Soliman, Generalized mode solver for plasmonic transmission lines embedded in layered media based on the method of moments, Computer Physics Communications (2018), https://doi.org/10.1016/j.cpc.2018.06.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Generalized Mode Solver for Plasmonic Transmission Lines Embedded in Layered Media based on the Method of Moments Mai O. Sallama,b,c,*, Guy A. E. Vandenboschb, Georges Gielenb, and Ezzeldin A. Solimana,* a

Physics Department, The American University in Cairo, P.O. Box 74, New Cairo 11835, Egypt b

c

ESAT-TELEMIC, Katholieke Universteit Leuven, B-3001-Leuven, Belgium

Electrical Engineering Department, Ahram Canadian University, Giza, Egypt *[email protected], * [email protected]

This paper presents an integral equation formulation for the calculation of the propagation characteristics of plasmonic transmission lines embedded within a multi-layered structure. The Method of Moments (MoM) technique is adopted in this paper due to its superior advantages over other techniques including the finite difference and finite element methods. Plasmonic transmission lines consist of a number metallic strips of arbitrary shapes immersed within a stack of planar dielectric or metallic layers. These transmission lines can support one or more mode, each of which has its characteristic mode profile and it propagates with a certain propagation and attenuation constants. The developed solver is tested for different plasmonic transmission line topologies surrounded by various layered media. The obtained results have are compared to CST commercial software for verification. Very good agreement between the proposed solver and CST has been observed. The developed MoM solver requires much smaller number of unknowns if compared with those based on Finite Difference Time Domain (FD-TD) or Finite Element Method (FEM) such as Lumerical and CST.

1. INTRODUCTION Over the last few decades, the field of plasmonics received considerable attention due to the unique properties of plasmonic materials at optical and visible wavelengths making them attractive for several applications [1]. The term “plasmonics” refers to the properties of metals at very high frequencies where light interacts with their free electrons resulting in the excitation of a surface wave, called Surface Plasmon Polariton (SPP). This wave propagates with high confinement along the interface between metallic/dielectric layers and exponentially decays away from it [2]. Such strong field confinement enables the design of optical devices of sub-wavelength dimensions [3]. Thanks to the advances in the fabrication technology, the manufacture of such extremely small devices is now feasible. The high field localization for plasmonic devices opened the gate for lots of applications based on light-matter interaction [4] including sensing [5, 6] and spectroscopy [7,8]. For sensing application, the spectral behavior of the transmitted or received power features either peaks or dips which are strongly affected by any variations in the medium surrounding the plasmonic sensor [9, 10]. Plasmonic devices could also be used for other applications including energy harvesting [11, 12] and optical telecommunications [13 – 15]. For energy harvesting applications, nano-antennas convert the electromagnetic waves in the visible/infra-red range into alternating voltage across the gap of the antenna. A rectifier is placed across this gap to convert the alternating voltage into

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

direct voltage. For optical telecommunication applications, nano-antennas are used for inter/intra-chip data transfer instead of using the lossy wires for connections. As mentioned earlier, plasmonic materials is concerned with the behavior of metals at very high frequencies (THz range and above). At lower frequencies (GHz range and below) metals are good conductors and are referred to as “classical metals”. Regardless the fact that plasmonic metals have dimensions down to the nano-scale, Maxwell’s equations are still valid and no quantum effects should be taken into consideration [16]. Nevertheless, plasmonic metals have some other differences compared to classical ones. First, they have very dispersive dielectric properties. Their dielectric permittivities could be approximated using Drude’s model [17] as follows:

ˆr       j     

 p2    jc 

(1)

where   is the epsilon at infinite frequency,  p is the plasma frequency, and c is the collision frequency. This equation shows that the dielectric permittivities of plasmonic materials consist of real and imaginary components which are highly frequency dependent. The real component could be either positive or negative and it describes the polarization of the plasmonic medium. On the other hand, the imaginary component is a measure for the losses inside the plasmonic material. These properties are different from classical metals, which are characterized by their very high conductivity. It is worth noting that Drude’s model is considered a good approximation for describing the dielectric permittivity of plasmonic materials at optical frequencies. It gives the correct prediction that at low optical frequency range, 180˚ out-of-phase polarization is achieved when an external field is applied [18]. This is interpreted by the negative sign of the real component of the dielectric constant that indicates high reflectivity [18]. The second difference between plasmonic and classical metals is related to the modeling of current flowing inside them. For classical metals, conduction current is flowing very close to the surface. As for plasmonic metals, polarization current is flowing along the entire volume. Therefore, plasmonic metals cannot be treated as 2D sheets like the case of classical metals. Instead, they should be dealt with as 3D structures [19]. Accordingly, these characteristics of plasmonic metals should be modeled appropriately in order to take into account the deviations from classical metals behavior. The aforementioned problems could be solved using different numerical techniques including differential (like finite difference [20, 21] and finite element methods [22]) or integral equation based techniques (like the method of moments, MoM). For differential techniques, the unknowns to be solved are the electric and magnetic fields within a volume surrounding and including the metallic objects. On the other hand, the integral techniques, the unknowns are the conduction/polarization current along the metallic objects only [23]. As a result, the integral methods are characterized by

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

their lower memory requirements and their ability to provide faster solutions when compared to the differential techniques [23]. This is due to the fact that integral techniques require only the discretization of domains carrying currents while differential techniques necessitate the discretization of the whole domain. Moreover, integral techniques don’t require erroneous truncation of the computation domain, which make them more suited to open structures such as nano antennas. On the other hand, the analytical efforts required to set up a system of equations for the integral techniques is more complicated compared to differential techniques [19] In this paper, an integral equations formulation solved using the MoM technique is used to evaluate the propagation characteristics of plasmonic Transmission Lines (TLs) immersed inside stratified layers of dielectric/metals. Specifically, it is used to calculate the propagation and attenuation constants of each propagating mode and their corresponding current profiles. Plasmonic TLs are essential components for optical nano-antennas (nantennas) intended for telecommunication applications where nantennas are usually formed from radiating wires, i.e. transmission lines. Moreover, TLs act as the interface between the electronic/optical circuits and nantennas. The developed solver has been used to analyze different plasmonic TLs topologies embedded within various layered media. The TLs investigated in this paper include: rectangular strip, horizontally coupled strips, and U-shaped strip placed on top of silicon dioxide substrate. Furthermore, we consider the nano-strip transmission line placed on top of finite thickness silicon dioxide layer which is backed by a gold layer. The results obtained from the developed solver are compared to those obtained from the commercial tool “CST” [24]. Generally, a very good agreement between both solvers is observed. The developed solver is further used to understand the effect of the transmission line dimensions on their propagation characteristics. In previous work [25], plasmonic TLs in homogeneous medium are considered. Compared to [25], this paper provides a generic solver which calculates the propagation characteristics of plasmonic TLs embedded in layered media. Unlike the work presented in [25], where the Green’s functions in homogeneous medium have a closed form, in this work, the Green’s functions have no closed form. Thus, they are calculated within the solver since they are dependent on the layered structure surrounding the TLs. In addition, in this paper, half roof-top basis functions at the edges of metallic strips are used to represent transversal current. This eliminates the approximation of zero transversal current at the edges in [25]. Furthermore, the testing technique adopted in this paper is the Razor-blade testing compared to the Galerkin testing that was used in [25]. The usage of Razor-blade testing technique reduces the calculation time when half roof-tops basis functions are considered. So in conclusion, compared to [25], this manuscript provides a more practical and accurate solver where in reality plasmonic TLs are immersed within a supporting structure. Nevertheless, breaking the structure symmetry is associated with

3

an increase of losses leading to smaller propagation length for the propagating wave [26]. The paper is arranged as follows. Section II presents the detailed formulation of the problem under consideration. Section III presents the results obtained from the developed solver with comparison to CST. The important conclusions are summarized in section IV. 2. Problem Formulation The structure of the generic plasmonic transmission line (TL) under investigation is presented in figure 1. It consists of one or more metallic strips of any topology. Each metallic strip has finite dimensions along the xz-plane while it has infinite extension along the y-axis. These strips are located within stratified layers of dielectrics and/or metals which have finite thickness along the z-direction and extend infinitely along the xy-plane. The upper and lower layers are considered as halfspaces. These TLs are capable of supporting one or more propagating modes according to both the dimensions and the number of the metallic strips comprising the TL.

Open Half-Space

Dielctric or Metal Layers

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

ˆ r 1

z

ˆ r 2

y

x

ˆ r 3

ˆ r 4

Open Half-Space

FIG. 1. General structure of a plasmonic transmission line embedded within stacked layers of dielectric/metal.

In this paper, the Method of Moments (MoM) technique is used to calculate the propagation and attenuation constants of the propagating mode(s) and their corresponding mode profiles(s). For a finite current distribution, the electric field can be obtained as follows [25]:

Exi  x, z    dz   dxGijE , J  x, z , x , z   J xj  x , z     dz   dx  z j

x j

z j

x j

2 E ,J  Gij  x, z , x , z   J xj  x , z   x 2

   E , J z  jk y  dz   dx GijE , J   x, z , x, z   J yj  x, z     dz   dx  Gij  x, z , x , z   J zj  x , z   x x z z j x j z j x j

4

(2a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

E yi ( x, z )   jk y  dz   dx z j

x j

 E ,J  Gij  x, z , x, z   J xj  x , z     dz   dx GijE , J  x, z , x , z   J yj  x , z   x z j x j

 k y2  dz   dxGijE , J   x, z , x, z   J yj  x, z    jk y  dz   dx z j

x j

Ezi ( x, z )   dz   dx z j

x j

z j

x j

(2b)

  Ez , J   Gij  x, z, x, z   J xj  x, z    jk y  dz   dx GijEz , J   x, z , x , z   J yj  x , z   z x z z j x j

  dz   dx GijEz , J z  x, z , x, z   J zj  x, z     dz   dx z j

x j

 E , J z Gij  x, z , x , z   J zj  x , z  z

z j

x j

 2 Ez , J z Gij  x, z , x , z   J zj  x , z   z 2

(2c)

where  x j , z j  is the domain of the current distribution in the xz-plane along the cross-section of the jth strip, GijE , J , GijE , J  , GijE , J z , GijEz , J  , and GijEz , J z are the spatial domain Green’s functions. As shown, each Green’s function is described by two

superscripts. The first one denotes the orientation of the electric field which could be lateral, E , or z-directed, Ez . The second superscript represents the type and the orientation of the current density: J is a lateral electric current density, J  is the derivative of the lateral current density, J z is a z-directed electric current density, and J z is the derivative of the zdirected electric current density. In the MoM technique, the metal strips comprising the plasmonic TL are discretized with rectangular segments with dimensions  x along the x-axis, and  z along the z-axis. The current along each strip is modeled as a summation of pre-assumed basis functions weighted by unknown coefficients, as described by the following equations: Nx

J x ( x, z )   An T  x  an  R( z  bn )

(3a)

n 1

Ny

J y ( x, z )   Bn R  x  cn  R ( z  d n )

(3b)

n 1

Nz

J z ( x, z )   Cn R  x  en  T ( z  f n )

(3c)

n 1

where N x , N y , and N z are the number of basis functions required to approximate the x- directed, y-directed, and z-directed current densities respectively. An , Bn , and Cn are the expansion coefficients (weights). T  u  , and R  u  are triangular and rectangular functions, respectively, forming the current basis functions.  an , bn  ,  cn , dn  , and  en , f n  are the  x, z  centers of the x-directed, y-directed, and z-directed current basis functions, respectively, as shown in figure 2. It can be noticed from (3) that the longitudinal current J y is represented by a rectangular prism basis function, while the transversal currents J x and J z

5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

are represented by triangular prism basis functions. The triangular functions defined in (3a) and (3c) could be an even triangular function, ET  u  , a right-sided triangular function, RT  u  , or a left-sided triangular function, LT u  , as shown in figure 3. The definitions of the triangular and rectangular functions are given by:  u   u   u2 ,   u  u  0  ET (u )    u  u   u2 , 0  u   u 0, otherwise 

(4a)

   u  u2 , 0  u   u RT (u )   u otherwise 0,

(4b)

 u  u  u2 ,  u  u  0 LT (u )   otherwise 0,

(4c)

1  ,   u 2   u   u 2  R(u )   u otherwise 0,

(4d)

It is worth noting that using ET  u  function together with R  u  function gives rise to a full “roof-top” basis function, while using either RT  u  or LT u  together with R results in a “half roof-top” basis function [27]. It is essential to highlight that the half roof-top basis functions are important to accurately describe J x and J z current densities at the edges of the strip. This is because of the fact that these polarization current components represent parallel electric field components, which are intensive at the edges of the strip. Using full roof-top at the edges forces the transversal electric field to be zero, which is physically wrong. On the other hand, half roof-tops allow the transverse field to be non-zero at the edges leading to more accurate representation of the current/field distribution. Figure 2 shows that J x  J z  at the left and right (top and bottom) edges of a metallic strip are represented by green arrows denoting half roof-top basis functions. On the other hand, full rooftops indicated by blue arrows are used to express J x and J z elsewhere.

(a)

(b)

(c)

FIG. 2. Basis functions used to expand the unknown modal current: (a) x-component, (b) y-component, and (c) z-component. Blue arrows of Jx and Jz denote full roof-top basis functions while green arrows represent half roof-tops basis functions.

6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

(c)

(d)

FIG. 3. Basis functions of the current density: (a) rectangular function, R  u  , (b) even triangular function, ET  u  , (c) left-sided triangular function, LT  u  , and (d) right-sided triangular function, RT  u  .

In spectral domain, the Green’s functions can be analytically evaluated by solving Maxwell’s equations along with the boundary equations at the source-carrying and source-free interfaces of the layer structure [28]. For a strip located on top and adjacent to the interface between one of the two half-spaces and its adjacent dielectric layer, the spectral domain Green’s functions are expressed by: GijE , J 

GijE , J   

0 2

 TE e jkz1  z  z1  e jkz1 z  z1    2  k z1 k z1  

 TM e  jkz1  z  z1  e  jkz1 z  z1    2  k z1 k z1    jk z  z k z21  TM e  jk z1  z  z1  e  jkz1 z  z1  e  jkz1  z  z1  e z1 1       2  k z1 k z1  21k 2  k z1 k z1 

(5a)

G11E , J k z21  k 2 21k 2 0 2k 2

 TE  2 

(5b)

GijE , J z 

 jk z  z 1  TM e jkz1  z  z1  e z1 1   2 21  k z1 k z1

  

(5c)

GijEz , J  

 jk z  z 1  TM e jkz1  z  z1  e z1 1    2  21  k z1 k z1 

(5d)

7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

GijEz , J z 

GijEz , J z 

0  TM e jkz1  z  z1  e jkz1 z  z1    2  2  k z1 k z1 

(5e)

 jk z  z 1  TM e jkz1  z  z1  e z1 1     2  21  k z1 k z1 

(5f)

where  is the angular frequency, 0 is the permeability of free-space,  1 is the dielectric constant of the first layer, k z1 is the propagation constant along the direction of stratification: k z1  k12  k2 and k2  kx2  k y2 , where k x and k y are the spectral TM domain counterparts of the spatial variables x and y.  TE 2   2  is the reflection coefficient in the second layer of the TE

(TM) system [28]. The value of this reflection coefficient terms depends on the layered structure, as will be shown in section 3. The spatial domain Green’s functions for filament sources can be obtained from their spectral domain counterparts using the following Sommerfeld identity [29]: Gij  x, z  

1 2



 Gij  k x , z  e

jk x  x  x  

dk x 



1





 G  k , z  cos  k  x  x  dk ij

x

x

x

(6)

0

Nevertheless, calculating the integration in (6) is numerically inefficient and time consuming. This is due to the slowly decaying and highly oscillatory nature of the integrand [30]. However, if the spectral Green’s function is expanded in the form of exponentials divided by k z1 , the Sommerfeld integration can be evaluated analytically [31] as follows: 

 0

e

 jk zi  z  z1 

k zi



cos  k x  x  x   dk x  jK0 

 x  x    z  z   2

2



Re z  z  0, Re   0, x  x  0

(7)

where   k y2  ki2 , K0 is the modified Bessel function of the second kind and zero order. The Discrete Complex Image Method (DCIM) [32] can perform this task. It begins by approximating the function (the spectral domain Green’s functions in our case) into a series of complex exponentials. In this work, the Generalized Pencil of Function (GPOF) [33] method is adopted for this approximation. Compared to Prony’s method, GPOF is more stable, more efficient, and less sensitive to numerical noise. Afterwards, the inverse Fourier transform of the obtained series of exponentials can be calculated analytically. It is worth mentioning that the GPOF method requires uniform sampling of the function to be approximated. Doing this, any function f  t  can be expressed as: N

f  t    wn ent n 1

8

(8)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where N is the number of terms required to approximate the function f  t  and it is determined by the GPOF method,

t is a

variable which is sampled uniformly from 0 to tm. wn and  n represent the coefficient and exponent terms of the nth term. Due to the fast variations of the spectral domain Green’s functions at low spectral values, an accurate approximation of these functions requires large number of samples. On the contrary, at high spectral values, the spectral domain Green’s functions vary slowly. Consequently, it is more efficient to perform the DCIM along two levels 32. The first level corresponds to the high spectrum part of the function where a low sampling rate is enough for an accurate approximation. In the second level, the sampling rate must be high enough in order to capture the fast-varying behavior of the Green’s functions in the low spectrum part. Since the Sommerfeld integration along k x starts from zero to infinity, according to (6), the sampling should cover this range. Figure 4 shows the sampling paths along the k  -plane for the first level (C1) and second level (C2) of the DCIM. The function is first discretized along the high spectral values (i.e. along C1) starting from k    2 until k   1 . The optimum choice of  2 is to be slightly higher than the maximum propagation constant along the layered medium. The end point of the first path is taken at 1  50k0 , which is enough for a high accuracy. It is worth noting that along this level, the Green’s functions are sampled along the real axis of k x due to the absence of poles and branch cuts along the high spectrum. The second path starts at  3  k y , which corresponds to kx  0 , and it ends at the starting point of the first path (i.e. at  2 ). It is worth noting that sampling along the real axis is not possible along this level due to the presence of poles and branch cuts resulting in singularities. Thus, path C2 is performed along the complex plane of k  , as shown in figure 4. In terms of k zi the first path starts at  j 2 , and ends at  j 1 whereas the second path starts at  j 3 and ends at  j 2 , as shown in figure 5 where:

kzi2  ki2  k2  ki2  k y2  kx2

(9)

 1  12  ki2

(10a)

 2  22  ki2

(10b)

 3  k y2  ki2

(10c)

For the GPOF method, we assume that sampling is uniform along the variable t, which starts here at t = 0 until t = 1. Therefore, the equations relating k zi with t along the two paths C1 and C2 are given by:

9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

  j  2    1   2  t  k zi     j  3    2   3  t     t   

2



along C1 (11) along C2

jk zi  2  1

along C1

 2  1 3 jk zi  3 2 3 2

(12) along C2

FIG. 4. Sampling path across the k  - plane using the two-levk2el DCIM.

FIG. 5. Sampling path across the k zi - plane using the two-level DCIM.

After sampling the function f (t ) with uniform steps along t, the GPOF method is applied where it calculates the coefficients and exponents of the exponentials defined in (8). In terms of k zi , the function along path C1, could be represented as follows:

10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

N1

f1 (t )   w1n e

1n t

n 1

1n

N1

  w1n e

jk zi 2 1n  2  1  2  1

 f1  k zi 

e

n 1

(13)

Now, the low spectral values are not considered along path C1 and the approximation of the function along this path does not provide an accurate approximation for the low spectrum. Thus, the differences between the exact function f (t ) and

f1 (t ) must be sampled along C2 with uniform steps of t

32

. Doing so, the following equation is obtained after applying the

GPOF technique:

N2

2 n

N2

f 2 (t )  f (t )  f1 (t )   w2 n e2 nt   w2 n e n 1

3  3  2

e

 2 n

jk zi

 3  2

n 1

 f 2  k zi 

(14)

Consequently, each of the spectral domain Green’s functions along the entire spectrum can be finally expressed as: N1

1n

f (k zi )  f1 (k zi )  f 2 (k zi )   w1n e

jk zi 2 1n  2  1  2  1

e

n 1

N2

2 n

  w2n e

3  3  2

e

 2 n

jk zi

 3  2

(15)

n 1

It is important to note that the DCIM process is performed on the spectral functions appearing in (5) without including the terms e jkz1  z  z1  k z1 and e jkz1 z  z1 k z1 . In other words, the coefficients of the exponentials will only be subjected to the DCIM. The reason for that is to keep the DCIM outcome independent on the spatial distances between the test and source locations. Now, after the two-levels DCIM is applied on the spectral Green’s functions as specified, and the omitted terms are put back, and the first function in (5b) is divided by k z1 , any spectral Green’s function in (5) can take either one of the following forms:         jk z 1  z  z1  2 n   jk z 1  z  z1  1n  3 2  3  2   2  1  N1 N2     2n 1n  e  2  1 e   w2 n e  3  2  w1n e k k z1 n 1  n 1 z1 G  k zi , z        jk z 1  z  z1  1n   N1    3   2  1  N2   jk z 1  z  z1  2 n  2 n 1n 2 e  3  2   3  2  w e  2  1    w2 n e e 1n  k z1 n 1  n 1

(16)

Applying Sommerfield’s identity in (7) on the spectral Green’s function in (16), their spatial domain counterparts can be expressed as follows:

11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2  2  N1    w e 1n  2 1 K    x  x 2   z  z   1n     1n 0    2  1   n 1      2    3 2 n  N2   2n   2  3  2    K0   x  x    z  z     w2 n e   3 2     n 1    G  x, z     2    1n 2  N1  1n   2    w1n e 2 1 K 0    x  x    z  z     2  1     n 1    2    N2  3   w e 2 n  3  2 K    x  x 2   z  z    2 n   2n 0    3 2   n 1    

(17)

At this stage, all needed spatial domain Green’s functions are prepared. The consequent step is to insert these expressions in (2), and apply a suitable testing technique. In this work, razor-blade testing [34] is adopted where the test function is in the form of a rectangular prism. Appling razor-blade testing, the integrally tested electric field components can be expressed in a matrix form as follows:

 Z xy    E x    Z xx1    Z xx 2      Z yx   Z yy1    Z yy 2    E y         Z zy    Ez     Z zx 

where

 Z xx1  ,  Z xx 2  ,  Z

xy

 Z xz 

    Z yz    Z zz1    Z zz 2 

 A    B  C 

 ,  Z xz  ,  Z yx  ,  Z yy1  ,  Z yy 2  ,  Z yz  ,  Z zx  ,  Z zy  ,  Z zz1  , and

 Z zz 2  ,

(18)

are the impedance sub-

matrices. The full expressions for the element  m, n  in these matrices, where m and n are the orders of the test and basis functions, respectively, are given below:

Z xx1 (m, n) 

 j0 2

z

x

z

x

z

x

z

x

 dz  dxR  x  am  R  z  bm   dz   dxG  x, x, z, z  T  x  an  R  z   bn 

12

(19a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

j

Z xx 2 (m, n)   

z

x

z

x

z

x

z

x

x

 dz  dx

21 j

z

x

z

x

z

x

 dz  dx

21

 2 G  x, x, z , z  

 dz  dxR  x  am R  z  bm   dz   dx

21

j

z

R  x  am  x

z

x

z

x

x 2

R  z  bm   dz   dx G  x, x , z , z  

  x  x     x  x  x

T  x   an R  z   bn 

T  x   an  R  z   bn  x

R  z  bm 

z x   x  an    *  dz   dx G  x, x, z , z    T   x  an   sign *  R  z   bn  x z x   



j

 

k y 21 21

z

x

ky 21 z

x

z

x

z





z

x

z

x

z

x

z

x

 dz  dxR  x  am  R  z  bm   dz   dx R  x  am  x

 dz R  z  bm   dx

z

z

x

z

z

x

z

x

z

x



x

z

x

x  x

 G  x, x , z , z   x  x



x

z

x

21

z

x

z

x

z

x

21 j 21

 dz  dx R  x  am  R  z  bm   dz   dx

z

x

z

x

z

x

 dz  dx

R  x  am  R  z  bm  z x

R  x   cn  R  z   d n 

 dz   dx G  x, x, z, z  R  x  cn  R  z   d n 

x

 dz  dz   dx  G  x, x, z, z  

 dz  dx

x

(19b)

R  z  bm   dz   dx  G  x, x , z , z  R  x   cn  R  z   d n 

z

j

G  x, x, z , z  

  x  x     x  x  z 

x

z

j



z

x

1 21  x

Z xz (m, n) 



  z dz z dz   G  x, x, z, z  xx xan,  G  x, x, z, z  xx xan,  R  z  bm  R  z   bn   

 dz  dx  1

ky



z

 dz   dx G  x, x, z, z  x  x  G  x, x, z, z   x  x R  z  bm  T   x  an  R  z   bn 

z

ky



x

sign 21 2x

Z xy (m, n)  

z

 dz

21 x

j

z





x

z

x

 R  z  b  R  x  c  R  z   d 

 2 G  x, x, z , z   xz

m

n

n

(19c)

R  x   en  T  z   f n 

 dz   dxG  x, x, z, z R  x  en  T  z   f n 

z

  x  x     x  x    z  z     z  z  x

*

z z

x

z

x

*  dz   dxG  x, x, z , z  R  x   en  T  z   f n  

j 21 x  z

z

x



 dz   dx  G  x, x, z, z 

z

x

x  x , z  z

  G  x, x, z , z   x  x ,  G  x, x, z , z   x  x ,  G  x, x , z , z   x  x ,  z  z z  z z  z 

* R  x  en  T  z   f n 

13

(19d)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

k y

Z yx (m, n)    

k y 21 ky 21

z

x x

z

 

ky 21

 

21 j 21

x

R  x  cm  x

x

z

z

x

 dz  dz   dx  G  x, x, z, z 

21 z

x

z

x

x  x



x

z

x

z

x

z

x

z

x

z

x

z

x

z

x

z

x

z

x

z

x

 dz  dx R  x  cm  R  z  d m   dz   dx R  z  d m  z 

x

z

x

z



x

z

z

x

x

x

z

x

j

z

x

z

x

z

x

z

x

 dx  dz   dx  G  x, x, z, z 

x

z

x

z

x

z

x

z  z

 dz  dx

R  x  em    z  f m  x

(20c)

R  x   en  T  z   f n 

n

z



 G  x, x , z , z   z  z R  x  cm  R  x   en T  z   f n  

 dz  dx R  x  em  R  z  f m   dz   dx

 dz  dx

(20b)

 dz   dx G  x, x, z, z   R  x  en T  z   f n 

z

z

z

n

  z  z     z  z  z 

x

G  x, x , z , z  

(20a)

 dz   dx G  x, x, z, z  R  x  e  T  z   f 

x

21



 G  x, x , z , z   x  x R  z  d m  T  x   an R  z   bn  z

 dz  dx R  x  cm 

z

x

x

m

z

z

 dz  dx R  z  d n  R  x  cn   dz  dx G  x, x , z , z  R  z  d m  R  x  cm 

x

x

x

z

 dz  dx R  x  c 

z

z

T  x   an R  z   bn 

 dz  dx R  z  dn  R  x  cn   dz  dx G  x, x, z, z  R  z  d m  R  x  cm 

21 k y

x

x

R  z  d m   dz   dx G  x, x, z , z   T  x  an R  z   bn 

x

z

z

G  x, x, z , z  

R  z  d m   dz   dx G  x, x, z , z   T  x  an R  z   bn 

  x  x     x  x 

z

jk y2

21 z

j

z

x

ky

Z zx (m, n) 

x

 j0 2

Z yy 2 (m, n) 

21

z

 dz  dx

21 x

k y

x

x

z

Z yy1 (m, n) 



z

 dz  dx

ky

Z yz (m, n) 

x

 dz  dxR  x  cm R  z  d m   dz   dx

21 z

z

z

x

z

x

 2 G  x, x, z , z   zx

(21d) T  x   an  R  z   bn 

R  dz   dx G  x, x , z , z   T  x   an  R  z   bn 

  x  x     x  x    z  z     z  z  *

x

z z

x

z

x

*  dz   dx G  x, x , z , z   T  x   an  R  z   bn  

j 21 x  z

z

x



 dz   dx  G  x, x, z, z 

z

x

x  x z  z

  G  x, x , z , z   x  x  G  x, x , z , z   x  x  G  x, x , z , z   x  x  z  z z  z z  z 

* T  x   an  R  z   bn 

14

(21a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Z zy (m, n)    

k y

k y 21 z

21 k y 21

z

x

z

x

z

x

z

x

x

 dz  dxR  x  e 

R  z  f m  z

x

x

z

x

z

 dx  dz

k y 21 z

Z zz1 (m, n) 

x

x

z

x

 

z

x

R  x  em 

 j0 2

j 21

 j 21 j 21

z



z 

 dx  dz   dx  G  x, x, z, z 

z  z

n

x

dz   dx G  x, x , z , z  R  x   cn  R  z   d n  x 

 G  x, x, z , z   z  z

z

x

z

x

z

x

z

x

Z zz2 (m, n)  

z

z

x

R  x   cn  R  z   d n 

n

  z  z     z  z 

x



z

 dz   dx G  x, x, z, z R  x  c  R  z   d 

z

m

z

G  x, x, z , z  

 dz  dxR  x  em  R  z  f m   dz   dx



 R  x  e  R  x  c  R  z   d  m

n

(21b)

n

 dz  dx R  x  em  R  z  f m   dz   dxG  x, x, z, z R  x  en T  z   f n  j 21

z

z

x

z

x

z

x

z

x

 dz  dx R  x  em R  z  f m   dz   dx

x

 dz  dx R  x  em 

z

x

z

x

 dz  dx R  x  e 

R  z  f m  z 

z

z

x

z

x

z

x

 dz  dx

z

R  z  f m  z

m

x

z

z 2

R  x   en  T  z   f n 

G  x, x , z , z   R  x  en  T  z   f n  z 

x

 dz   dxG  x, x, z, z  R  x  e  n

z

  z  z     z  z  z

x

 dz   dx  1

 2 G  x, x, z , z  

(21c)

x

T  z   f n  z 

R  x  em 

z x   z  fn    *  dz   dxG  x, x, z , z   R  x  en   T   z   f n   sign *  z z x   



j 21 z

x

z

x

x

z

x



 dx  dz   dx G  x, x, z, z  z  z  G  x, x, z, z  z  z

sign j 21 2z

where



x

x



 dx  dx  G  x, x, z, z 

x

x



z  z z   fn



 R  x  e R  x  e T   z  f  m

  G  x, x, z , z   z  z  R  x  em  R  x  en  z   fn 

n

n

(21d)

x ( z ) and x ( z ) are the maximum and minimum, respectively, limits of the test function along the x-axis (z-

axis), while

x ( z ) and x ( z ) are the maximum and minimum, respectively, limits of the basis function along the x-axis

(z-axis). As mentioned earlier, the triangular basis function can be an even triangle (ET), right-sided triangle (RT), or left-sided (LT) triangle. The choice between the three cases depends on the location of the current sources. The derivative of the triangular functions denoted by ET   u  , RT   u  , and LT   u  , as well as the delta function   u  , which appeared in (19)(21) are defined as follows:

15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1  u2 ,   u  u  0  ET (u )  1  u2 , 0  u   u 0, otherwise 

(22a)

1 u2 , 0  u  u RT (u )   otherwise 0,

(22b)

1  2 ,  u  u  0 LT (u )   u otherwise 0,

(22c)

1  u , 0,

u0 otherwise

 (u )  

It is worth noting that in (21d) G( x, x, z, z ) z is replaced by

(22d)

G( x, x, z, z ) z  . The negative (positive) sign





corresponds to the case when the spectral domain Green’s functions are proportional to e jkzi z  z kzi e jkzi  z  z  k zi . The value 



of the variable “sign” is determined according to the type of its corresponding basis function as follows:

Function

ET u 

RT u 

LT u 

Sign

0

+1

–1

Now, the next step is to solve (18) in order to obtain the propagation constant k y . Unlike the perfect conductors where current cannot flow inside their volumes, plasmonic materials allow current to flow inside them. The electric field inside any plasmonic strip is given by [25]: E( x , y , z ) i 

1 J ( x, y, z )i j 0 (ˆri  ˆrb )

(23)

where ˆri and ˆrb are the relative dielectric constants of the ith plasmonic strip and the background material surrounding it, respectively. Combining (18) and (23), we obtain:

  Ex     Z xx       E y      Z yx       Ez     Z zx 

 Z xy   Z xz     Z yy   Z yz     Z zy   Z zz   

 A  Z xx   0     B     0  Z yy  C    0  0 

where  Z xx    Z xx1    Z xx 2  ,  Z yy    Z yy1    Z yy 2  , and

0   A  0   B    Z zz  C 

 Z zz    Z zz1    Z zz 2  .The

(24)

majority of elements in the

  ,  Z yy  , and  Z zz  are zero. The non-zero elements in these matrices are given by: matrices  Z xx  (m, n)  Z xx

  z x 1  dz  dx R  x  am  R  z  bm  T  x  an  R  z  bn ,  j 0  ˆri  ˆrb  z x

16

m n 1

(25a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

 (m, n)  Z yy

  z x 1  dz  dx R  x  cm  R  z  d m  R  x  cn  R  z  d n , m  n  j 0  ˆri  ˆrb  z x

(25b)

Z zz (m, n) 

  z x 1  dz  dx R  x  em  R  z  f m  R  x  en  T  z  f n ,  j 0  ˆri  ˆrb  z x

(25c)

m n 1

These non-zero values occur when the basis and test functions have overlapping domains. For  Z yy  , the basis and test

   Z zz  , functions overlap only if their centers coincide. Therefore, only the diagonal terms of  Z yy  are non-zero. For  Z xx the basis and test functions overlap if their x-(z-) centers coincide or are shifted by a value of  x   z  . In a compact form, (24) reduces to:  Z    Z  W    0

(26)

where  Z  is the impedance matrix relating the electric fields to the electric current densities,  Z  is the impedance matrix after satisfying the boundary conditions along the metallic strips, and W  is a vector that represents the weights of the modal current. The non-trivial solution of equation (3.34) implies that the determinant of the matrix  Z    Z  , should vanish. The dimension of this matrix is M  M where M represents the total number of the current basis functions (i.e. It is worth noting that each element of the matrix  Z    Z  is k y dependent. Hence, the “characteristic equation” can be written as follows:

M  N x  N y  N z ).

 Z    Z 

 f ky   0

(27)

The value of the complex propagation constant, ky, can be obtained using iterative methods such as Muller’s method [35] which is adopted in our solver. The last step is to obtain the mode profile of the propagating modes. This is achieved by substituting with the value of k y in (25) and applying the Single Value Decomposition (SVD) technique [36] to obtain the weights of all current components. 3. Numerical examples In this section, the developed MoM-based solver is applied to different plasmonic TL topologies in order to calculate their propagation characteristics. The TLs considered in this paper have various cross-sections, number of metallic strips, and/or are arranged in different layer structures, as shown in figure 6. The diversity of these examples demonstrates that the proposed solver is generic. For all examples considered in this paper, the metallic strips are made of gold (Au) material which has the following Drude model parameters [37]    9.069 ,  p  1.354 1016 rad/s, and f c  1.2  1014 s-12. The frequency considered in this paper ranges from (150 THz – 400 THz) corresponding to a wavelength (2 µm – 0.75 µm). 1 and  2

17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

appearing in (10) are considered to be equal to 50 k0 , and 5 k0 respectively. These values are found suitable for the layered media under investigation. In all the considered examples, the dispersion curves obtained using the developed solver are compared with the ones obtained with the frequency domain solver of CST. The effects of varying the main geometrical parameters of the TLs on their propagation characteristics are also studied.

(a)

(b)

(c)

(d)

FIG. 6. The various plasmonic TLs under investigation: (a) single gold strip on top of SiO 2 substrate, (b) horizontally coupled gold strips on top of SiO2 substrate, (c) U-shaped plasmonic gold strip placed on top of SiO2 substrate, and (d) plasmonic nanostrip placed on top of SiO2 substrate backed by gold ground layer.

a. Single strip on top of SiO2 substrate In this section, the propagation characteristics of the fundamental mode of a rectangular gold strip on top of a silicon dioxide substrate (  r  2.15 ), as shown in figure 6(a), are investigated. At optical frequencies, the guided wavelength is in the range of a few micrometers. Therefore, the thickness of the substrates is much larger than the wavelength. As a result, the assumption of a gold strip placed on top of an infinite substrate is considered a valid approximation. For the case of a current TM source located at the interface of two half spaces, the reflection coefficients  TE are expressed as: 2 , and  2

TE TE 2  RSiO2 , Air e

j 2 k z0 z2

TE , RSiO  2 , air

SiO kz  0 kz 2

0

SiO2

SiO kz  0 kz 2

18

0

SiO2



k z0  k zSiO

2

k z0  k zSiO

2

(28a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

TM TM  Rair 2 , SiO2 e

j 2 k z0 z2

TM , Rair , SiO2 

 SiO kz   0 kz 2

2

where z2 is the z-location of the air/SiO2 interface, k z2

SiO2

0

SiO2

 SiO kz   0 kz 0

(28b)

SiO2

2  k SiO  k 2 , kSiO2 is the propagation constant in the SiO2 medium and 2

 SiO is the dielectric permittivity of SiO2 . 2

The dimensions of the gold strip considered in this example are: W = 50 nm and T = 20 nm. Figure 7 shows the dispersion curves as simulated using our solver and CST. Both solvers have closely matched results. They show that as the frequency increases, the effective refractive index (neff) and the insertion loss (IL) of the strip increases. This is due to the fact that gold has higher skin depth values with increasing frequency25. Consequently, the current penetration inside the plasmonic strip increases which leads to a higher effective refractive index and high losses. At 193.55 THz (1.55 µm), the propagation constant calculated using our solver is: k y  8.23 106  3.14 105 j . At this frequency, the x-, y-, and z- current components are shown in Figure 8. As expected, the current penetrates the gold strip due to the dielectric nature of gold at optical frequencies. Furthermore, figure 8 shows that the y-component of the current density J y is dominant, but the x- and zcurrent components cannot be ignored. For J y , it is noticed that the current is maximum along the strip edges compared to its center and it is symmetric along the x-axis. On the other hand, J y is not symmetric across the z-axis where it has higher values at the lower edge, i.e. along the interface of the gold strip with the silicon dioxide layer. This is expected due to the relatively high dielectric permittivity of SiO2 compared to free-space. The x- and z- current components shown in figure 8(a) and (c) have an odd (even) symmetry along the z-axis. The maximum current for J x  J z  is along the right and left (upper and bottom) edges of the strip. It is worth mentioning that the longitudinal current J y lags the transversal currents ( J x , and

J z ) with 90 degrees in phase. The effect of varying the geometrical parameters W and T of the gold strip on its propagation characteristics is shown in figure 9, as simulated using our solver at 193.55 THz. The width (W) is varied from 30 nm till 70 nm while keeping the strip thickness fixed at T = 20 µm. On the other hand, figure 10 shows the effect of varying the thickness (T) from 10 nm till 30 nm at constant width of W = 50 µm. As shown in these figures, decreasing W or T leads to increasing both the effective refractive index and the insertion loss. This is due to the higher confinement of the wave inside the plasmonic gold strip, which is characterized by its high dielectric constant and high losses.

19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

FIG. 7. Propagation characteristics of a single gold strip of dimensions W = 50 nm and T = 20 nm on top of an infinite silicon dioxide substrate: (a) effective refractive index, and (b) insertion loss.

(a)

(b)

(c)

FIG. 8. Current distribution along the plasmonic gold strip above the SiO2 substrate at 193.55 THz as obtained using our solver: (a) xcomponent

 J x  , (b) y-component  J y  , and (c) z-component  J z  .

The gold plasmonic strip has dimensions: W = 50 nm, and T =

20nm.

(a)

(b)

FIG. 9. The effect of varying the strip width (W) of the plasmonic gold strip on top of the SiO2 substrate on the propagation characteristics as simulated using our solver when the strip thickness T = 20 nm: (a) the effective refractive index, and (b) the insertion loss.

20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

FIG. 10. The effect of varying the strip thickness (T) of the plasmonic gold strip on top of the SiO2 substrate on the propagation characteristics as simulated using our solver when the strip width W = 20 nm: (a) the effective refractive index, and (b) the insertion loss.

It is worth noting that the proposed MoM-based solver has the significant advantage of reducing the number of unknowns compared to CST. For instance, in the example of the single plasmonic strip, the proposed solver requires from 160 (208) to 460 (4016) mesh cells for the 50 nm  20 nm case. The number of mesh cells is determined according to the operating frequency. On the other hand, CST requires (18,000 – 35,000) cells to obtain the same results. The significantly reduced number of unknowns (by three orders of magnitude) results in a lower intrinsic need of computational resources when compared to CST. b. Horizontally coupled strips on top of SiO2 substrate The horizontally coupled plasmonic TL considered in this section consists of two symmetric gold strips of width W and thickness T, separated by a distance S, as shown in figure 6(b). The strips are placed on top of a silicon dioxide substrate like in the previous example. This type of TL is capable of supporting two main modes namely the even (symmetric) and the odd (asymmetric) mode. Figure 11 shows the propagation characteristics of the two modes for gold strips of the following dimensions: W = 50 µm, T = 20 µm, and S = 20 µm. The behavior of neff and IL versus frequency is similar to that of the single strip. Increasing the frequency results in higher neff and IL, as explained in the previous example. Figure 11 also shows that the even mode has lower values of the propagation constant, and consequently lower propagation losses, compared to the odd mode. The current mode profiles for the propagating modes along the horizontally coupled strips are shown in Figures 12, and 13 at 193.55 THz. The even mode ( k y  7.31106  2.26 105 j ) has symmetric J y (longitudinal) and J z (transversal) currents along the two strips and anti-symmetric J x (transversal) currents. On the contrary, for the odd mode

21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

( k y  1.16 107  4.79 105 j ), the longitudinal current J y and the transversal current J z have odd symmetry, while J x has odd symmetry.

(a)

(b)

FIG. 11. Propagation characteristics of the horizontally coupled gold strips of dimensions W = 50 nm, T = 20 nm, and S = 20 nm on top of an infinite silicon dioxide substrate: (a) effective refractive index, and (b) insertion loss.

(a)

(b)

(c)

FIG. 12. Current distribution of the even-mode along the cross-section of the horizontally coupled gold strips (W = 50 nm, T = 20 nm, and

 

S = 20 nm) on top of SiO2 substrate at 193.55 THz as obtained using the MoM technique: (a) x-component  J x  (b) y-component J y , and (c) z-component  J z  .

(a)

(b)

(c)

FIG. 13. Current distribution of the odd-mode along the cross-section of the horizontally coupled gold strips (W = 50 nm, T = 20 nm, and

 

S = 20 nm ) on top of SiO2 substrate at 193.55 THz as obtained using the MoM technique: (a) x-component  J x  (b) y-component J y , and (c) z-component  J z  .

c. U-shaped strip on top of SiO2 substrate

The developed solver has also been examined for plasmonic TLs of various topologies including the U-shaped structure, shown in figure 6(c). Similar to the previous examples, this TL is placed on top of a SiO2 substrate. The dimensions of this strip are given by W = 80 nm, T = 20 nm, Ws = 40 nm, and Ts = 10 nm. Figure 14 shows the variation of neff and IL

22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

versus frequency. Simulations using our solver and CST are in a very good agreement. As can be concluded from the previous examples, the agreement between both solvers demonstrates that the developed solver is capable of solving plasmonic transmission lines of any topology constructed from one or more plasmonic strips.

(a)

(b)

FIG. 14. Propagation characteristics of a U-shaped gold strip of dimensions: W = 80 nm, T = 20 nm, Ws = 20 nm, and Ts = 10 nm on top of an infinite silicon dioxide substrate: (a) effective refractive index, and (b) insertion loss.

d. Nano-strip on top of SiO2 substrate backed by gold ground layer The last example considered in this paper, is a gold nano-strip placed on a finite SiO2 substrate backed by a ground layer of gold, as shown in figure 6(d). The assumption made in our solver is that the bottom gold layer is extended to infinity in the lower half-space. In this case, the reflection coefficient terms  TE and TM are expressed as follows: 2 2



TE 2







TM 2

TE RSiO  TE 3 e 2 , air

1 R

TE SiO2 , air

1 R

2

TE  j 2 k zSiO2 z2 3

 e

TM TM Rair , SiO2  3 e TM air , SiO2

 j 2 k zSiO z2 j 2 k z0 z2

e

j 2 k z0 z2

(29a)

 j 2 k zSiO z2 2

 e TM 3

e

 j 2 k zSiO z2

(29b)

2

where, TE TE 3  RAu , SiO2 e

j 2 k zSiO z3 2

TE , RAu , SiO2 

 Au k z  Au k z

23

SiO2

SiO2

  SiO k z Au 2

  SiO k z Au 2



kz kz

SiO2

SiO2

 k z Au  k z Au

(30a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

TM TM  RSiO e 3 2 , Au

j 2 k zSiO z3 2

TM , RSiO  2 , Au

 Au k z  Au k z

SiO2

SiO2

  SiO k z Au 2

(30b)

  SiO k z Au 2

2 where k z2  k Au  k 2 , while z2, and z3 are the z-locations of the air/SiO2 and SiO2/Au interfaces respectively. Au

The dimensions chosen for the nano-strip are as follows: W = 50 nm, T = 20 nm, and H = 100 nm. The propagation characteristics of this TL are shown in figure 15. They have quite similar behavior compared to the single strip TL, discussed in section (3-a). The effective refractive index and the insertion losses increase with increasing frequency. Comparing with the single strip TL, the nano-strip TL has higher neff and IL. This is due to the presence of the gold layer, which has higher dielectric permittivity and losses, compared to SiO2. As shown in figure 16, by increasing the SiO 2 thickness (H), neff and IL decrease. This is due to the fact that the gold layer is located further away from the nano-strip and has weaker effect on the TL propagation characteristics. Further increasing of the SiO 2 thickness tends to make the propagation characteristics of this TL resemble that of the singular rectangular strip of the same dimensions.

(a)

(b)

FIG. 15. Propagation characteristics of a single gold nano-strip of dimensions W = 50 nm and T = 20 nm on top of a 100 nm thick SiO2 substrate backed with infinite ground gold layer: (a) effective refractive index, and (b) insertion loss.

24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

FIG. 16. The effect of varying the SiO2 substrate thickness on the propagation characteristics of the plasmonic gold nano-strip (W = 50 nm and T = 20 nm) at 193.55 THz as simulated using the MoM technique: (a) the effective refractive index, and (b) the insertion loss.

It is worth mentioning that the MoM is considered a powerful tool for the calculation of the propagation characteristics of transmission lines. The strength of the MoM arises from the fact that it requires meshing for only the surface of the metallic strips. On the other hand, CST, which is based on the finite difference method, requires substantial additional meshing for the space surrounding the strips. This additional space should be made wide enough to ensure that the fields are not strongly disturbed by the truncation boundaries. Consequently, the number of unknowns is way larger than that of the MoM. Table 1 shows a comparison of the number of unknowns used in calculating the propagation characteristics of the single rectangular strip, horizontally coupled strips, and vertically coupled strips in free-space using both solvers at 193.55 THz (1.55 μm). It is clear that our solver requires much less number of unknowns compared to CST software, illustrating the intrinsic supremacy of this technique for this type of problem. Table 1: Number of unknowns required to calculate the propagation characteristics of the plasmonic transmission lines at 1.55 μm. Transmission Line Topology

MoM Solver

CST

Single Metallic Strip

1016

39192

Horizontally Coupled Strips

2032

68040

Vertically Coupled Strips

2032

74580

4. CONCLUSION This paper presented an integral equation formulation for plasmonic transmission lines placed inside stratified layers of dielectrics/metals. The developed solver counts for the properties of plasmonic materials, which are different from classical

25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

metals. The proposed solver is characterized by its greatly reduced number of unknowns compared to differential equation based solvers. A number of plasmonic transmission line topologies were studied in this paper including the single strip, horizontally coupled strips, and the U-shaped strip placed on top of a SiO2 layer. Furthermore, a plasmonic nano-strip placed on a finite SiO2 substrate backed by a ground gold layer has also been studied. The effect of varying the main parameters of the transmission lines is also presented. A comparison between the results obtained from the developed solver and CST is made. The proposed solver produces results very close to those of CST.

REFERENCES [1] U. Guler, V. M. Shalaev, A. Boltasseva, Materialstoday 18 (2015), 227–237. [2] A. Hosseini, H. Nejati, Y. Massoud, Optics Express 15 (2007), 15280–15286. [3] W. L. Barnes, A. Dereux, T. W. Ebbesen, Nature 424 (2003), 824–830. [4] M. Chamanzar, Z. Xia, S. Yegnanarayanan , A. Adibi, Optics Express 21(2013), 32086–32098. [5] T. Wu, Y. Liu, Z. Yu, Y. Peng, C. Shu, H. Ye, Optics Express 22 (2014), 7669–7677. [6] N. Liu, M. L. Tang, M. Hentschel, H. Giessen, A. P. Alivisatos, Nature Materials 10 (2011), 631–636. [7] M. Kaniber, K. Schraml, A. Regler, J. Bartl, G. Glashagen, F. Flassig, J. Wierzbowski, J. J. Finley, Scientific Reports 6, (2016), 23203. [8] D. A. Kalashnikov, Z. Pan, A. I. Kuznetsov, and L. A. Krivitsky, Physical Review X 4 (2014), 011049. [9] B. You, J.-Y. Lu, and J.-L. Peng, Optics Express 21 (2013), 210897–21096. [10]S. M. Kandil, T. A. Ali, S. Sedky, E. A. Soliman, European Conference on Antennas and Propagation (EuCAP) (2016), 1–5.k [11] M. Hussein, N. F. F. Areed, M. F. O. Hameed, S. S. A. Obaya, IET Optoelectronics 8 (2014), 167–173. [12] I. E. Hashem, N. H. Rafat, E. A. Soliman, IEEE Transaction on Nanotechnology 13(2014), 767–778. [13] G. N. Malheiros-Silveira, G. S. Wiederhecker, H. E. Hernández-Figueroa, Optics Express 21 (2013), 1234–1239. [14] L. Yousefi, A. C. Foster, Optics Express 20 (2012), 18326–18335. [15] E. A. Soliman, M. O. Sallam, G. A. E. Vandenbosch, IEEE Journal of Lightwave Technology 32 (2014), 4296–4302. [16] Di. M. Solis, J. M. Taboada, L. Landesa, J. L. Rodriguez, and F. Obelleiro, Progress in Electromagnetics Research 154 (2015), 35–50. [17] S. A. Maier, Plasmonics: Fundamentals and Applications, Springer, 2007. [18] W. Park, Nano Convergence 1 (2014), 1–27. [19] G. A. E. Vandenbosch, V. Volski, N. Verellen, V. V. Moshchalkov, Radio Science 46 (2011), RS0E02. [20] G. Veronis, S. Fan, Journal of Lightwave Technology 25 (2007), 2511–2521. [21] D. F. P. Pile, D. K. Gramotnev, M. Haraguchi, T. Okamoto, M. Fukui, Journal of Applied Physics 100 (2006), 013101. [22] T. Holmgaard, Physics Review B 75, 245405 (2007). [23] E. A. Soliman, G. A. E. Vandenbosch, E. Beyne R. P. Mertens, IEEE Transaction on Microwave Theory and Techniques 51 (2003), 874–886. [24] CST Microwave Studio, 2012. www.cst.com [25] M. O. Sallam, G. A. E. Vandenbosch, G. Gielen, E. A. Soliman, Optics Express 22 (2014), 22388–22402. [26]W. Wei, X. Zhang, Y. Huang, X. Ren, Nanoscale Research Letters 9 (2014), 1–6. [27] E. A. Soliman, A. K. Abdelmageed, M. A. El-Gamal, International Journal of Electronics and Communications 56 (2002), 155–262. [28] E. A. Soliman, G. A. E. Vandenbosch, Progress In Electromagnetics Research 62 (2006), 21–40. [29] A. Sommerfeld, Partial Differential Equations in Physics, Academic, New York, 1949. [30] R. Golubovic, A. G. Polimeridis, J. R. Mosig, IEEE Transactions on Antennas and Propagation 60 (2012), 2409–2417.

26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[31] E. A. Soliman, P. Pieter, E. Beyne, G. A. E. Vandenbosch, IEEE Transactions on Microwave Theory and Techniques 47 (1999), 1782–1787. [32] M. I. Aksun, IEEE Transaction on Microwave Theory Technology 44 (1996), 651–658. [33] Y. Hua, T. Sarkar, IEEE Transaction on Antennas and Propagation 37 (1989), 229–234. [34] E. A. Soliman, G. A. E. Vandenbosch, E. Beyne, International Journal of RF and Microwave computer-aided engineering 10 (2000), 132–138. [35] S. D. Conte and C. de Boor, Elementary Numerical Analysis: An algorithmic Approach, McGraw-Hill, 1990. [36]MATLAB, version 7.9.0. Natick, Massachusetts: The MathWorks Inc., 2009. [37] E. A. Soliman, Microwave and Optical Technology Letters 54 (2012), 2209–2214.

27