Numerical study
on
surface plasmon polariton
behaviors in periodic metal-dielectric structures
using a
plane-wave-assisted boundary
integral-equation method
Jyh-Yang Wang, C. C. Yang, and Yean-Woei Kiang
Graduate Institute of Communication Engineering, Graduate Institute of Electro-Optical Engineering, and Department of Electrical Engineering, National Taiwan University, Taipei 10617, Taiwan
ywkiang@ntu.edu.tw
Abstract: A novel hybrid technique based on the boundary
integral-equation method is proposed for studying the surface plasmon polariton behaviors in two-dimensional periodic structures. Considering the periodicity property of the problem, we use the plane-wave expansion concept and the periodic boundary condition instead of using the periodic Green’s function. The diffraction efficiency can then be readily calculated once the equivalent electric and magnetic currents are solved that avoids invoking the numerical calculation of the radiation integral. The numerical validity is verified with the cases of highly conducting materials and practical metals. Numerical convergence can be easily achieved even in the case of a large incident angle as 80o. Based on the numerical scheme, a metal-dielectric wavy structure is designed for enhancing the transmittance of optical signal through the structure. The excitation of the coupled surface plasmon polaritons for the high transmission is demonstrated.
©2007 Optical Society of America
OCIS codes: (050. 2770) Gratings; (240.6680) Surface plasmons; (050.1960) Diffraction theory.
References and links
1. H. Raether, Surface Plasmons (Springer-Verlag, Berlin, 1988).
2. T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and P. A. Wolff, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature 391, 667-669 (1998).
3. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302-307 (1966).
4. S. D. Gedney and R. Mittra, “Analysis of the electromagnetic scattering by thick gratings using a combined FEM/MM solution,” IEEE Trans. Antennas Propagat. 39, 1605-1614 (1991).
5. K. Yashiro and S. Ohkawa, “Boundary element method for electromagnetic scattering from cylinders,” IEEE Trans. Antennas Propagat. 33, 383-389 (1985).
6. L. C. Trintinalia and H. Ling, “Integral equation modeling of multilayered doubly-periodic lossy structures using periodic boundary condition and a connection scheme,” IEEE Trans. Antennas Propagat. 52, 2253-2261 (2004).
7. T. Sondergaard, S. I. Bozhevolnyi, and A. Boltasseva, “Theoretical analysis of ridge grating for long-range surface plasmon polaritons,” Phys. Rev. B 73, 045320 (2006).
8. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of metallic surface-relief gratings,” J. Opt. Soc. Am. A 3, 1780-1796 (1986).
9. E. Popov, B. Chernov, M. Nevière, and N. Bonod, “Differential theory: Application to highly conducting gratings,” J. Opt. Soc. Am. A 21, 199-206 (2004).
10. K. Watanabe, “Study of the differential theory of lamellar gratings made of highly conducting materials,” J. Opt. Soc. Am. A 23, 69-72 (2006).
11. E. Popov, M. Nevie`re, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A 19, 33-42 (2002).
12. K. Yasumoto and K. Yoshitomi, “Efficient calculation of lattice sums for free-space periodic Green’s function,” IEEE Trans. Antennas Propagat. 47, 1050-1055 (1999).
13. H. Rogier and D. De Zutter, “A fast converging series expansion for the 2-d periodic Green’s function based on perfectly matched layers,” IEEE Trans. Microwave Theory Tech. 52, 1199-1206 (2004). 14. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941).
15. K.-M. Chen, “A mathematical formulation of the equivalence principle,” IEEE Trans. Microwave Theory Tech. 37, 1576-1581 (1989).
16. C. A. Balanis, Advanced Engineering Electromagnetics (John Wiley & Sons, New York, 1989). 17. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13,
1870-1876 (1996).
18. D. K. Gifford and D. G. Hall, “Extraordinary transmission of organic photoluminescence through an otherwise opaque metal layer via surface plasmon cross coupling,” Appl. Phys. Lett. 80, 3679-3681 (2002). 19. S. Wedge and W. L. Barnes, “Surface plasmon-polariton mediated light emission through thin metal films,”
Opt. Express 12, 3673-3685 (2004).
20. C. Bonnand, J. Bellessa, C. Symonds, and J. C. Plenet, “Polaritonic emission via surface plasmon cross coupling,” Appl. Phys. Lett. 89, 231119 (2006).
21. U. Schroter and D. Heitmann, “Grating couplers for surface plasmons excited on thin metal films in the Kretschmann-Raether configuration,” Phys. Rev. B 60, 4992-4999 (1999).
22. I. R. Hooper and J. R. Sambles, “Coupled surface plasmon polaritons on thin metal slabs corrugated on both surfaces,” Phys. Rev. B 70, 045421 (2004).
23. D. Crouse and P. Keshavareddy, “Role of optical and surface plasmon modes in enhanced transmission and applications,” Opt. Express 13, 7760-7771 (2005).
1. Introduction
In recent years, the research of surface plasmon polaritons (SPPs) [1] in the near-infrared and visible ranges has received much attention since the phenomenon of SPP-mediated extraordinary transmission was discovered [2]. For theoretically investigating the SPP-related phenomena and designing the nanostructures for relevant applications, effective numerical methods are needed. Several numerical methods have been widely used for these purposes, including the finite-difference time-domain (FDTD) method [3], the finite-element method (FEM) [4], the integral equation method (IEM) [5-7], and the coupled-wave method (CWM) [8]. Among them, because of its computation efficiency for certain simple structures and simplicities in formulation and programming, the CWM is commonly used for calculating the diffraction efficiency of a periodic structure. However, the CWM suffers from inducing spurious features of diffraction efficiency in dealing with the structures of highly conducting materials [9, 10]. Another problem of slow convergence related to the staircase approximation for arbitrary-shaped metallic gratings is also encountered by the CWM [11].
In many SPP applications, the fine geometries of the structures are not in simple rectangular shapes and are usually quite complicated. In such a situation, the unstructured-mesh method is needed to precisely describe the problem geometry. The FEM and the IEM are the possible choices. However, in this research, we choose the boundary integral-equation method (BIEM) because it is based on a surface formulation, which automatically takes care of the discontinuities between different materials [5]. Besides, all the unknowns in the formulation are designated only on the structure boundaries and interfaces. These features result in a smaller number of unknowns when compared with other unstructured-mesh numerical techniques based on volumetric formulations, such as the FEM or the volume IEM, in which the unknowns are designated in the bulky domain. To deal with the periodicity property in a periodic problem based on the BIEM, one can usually apply the periodic Green’s function to the integral equation. However, the periodic Green’s function is usually in a form of infinite series and thus may suffer from the problem of slow convergence. Although certain techniques have been developed for speeding up the convergence, the mathematical treatment and programming implementation in applying these techniques to the IEM are quite complicated and case-dependent [12, 13]. Recently, a scheme of the BIEM was proposed, in which the periodic boundary condition was used and the free-space Green’s function was utilized in the major part of the unit cells [6]. However, the concept of the periodic Green’s function was still needed for the rest part of unit cells of homogeneous materials.
plane-BIEM in solving a problem of complex geometry. It also takes the advantage of using the plane-wave expansion method to express the fields in the homogeneous regions such as those outside the wave coupling area of major concern. In addition, the expansion coefficients can be connected to the equivalent electric and magnetic surface currents, which are the unknowns to be solved in the BIEM. Thus, we can apply the periodic boundary condition in the direction of periodicity and truncate the computation domain in the perpendicular direction. Such a choice can help in avoiding using the periodic Green’s function. The formulation can be applied to both metal and dielectric materials. Another advantage of this method is that the diffraction efficiency of any order can be directly calculated once the equivalent surface currents are solved such that the numerical calculation of the radiation integral is unnecessary. This paper is organized as follows. In section 2, the formulation of the plane-wave-assisted BIEM is presented. The numerical verifications for the proposed method and the simulations on a metal-dielectric wavy structure are shown in section 3. Finally, conclusions are drawn in section 4.
2. Plane-wave-assisted boundary integral-equation method
In this research, we consider only the two-dimensional (2D) transverse electric (TE or p-wave) case, in which all quantities are independent of z. We assume that the electric field is polarized in the x-y plane and the magnetic field is along the z direction. The transverse magnetic (TM) case can be similarly treated. Before discussing the plane-wave-assisted method, let us recall the basic concept of the BIEM. The BIEM is based on the Stratton-Chu formulas [14, 15], which state that the total field at an observation point, ρ, in an enclosed homogeneous, (ε, μ), region can be calculated, in addition to the incident field, (Einc, Hinc), by integrating all the contributions from the equivalent electric and magnetic surface currents,
( , )J M , on all boundaries of this region. If the observation point is located on a boundary (approaching from the interior of the region), the electric and magnetic fields in the 2D TE case can be expressed as
( )
( )
( ) (
)
( )
(
)
( )
(
)
( )
( )
( ) (
) ( )
(
)
2 +2 , , , 2 +2 , , inc inc E E i i J M J dl H H i M J dl ρ ρ ωμ ρ ϕ ρ ρ ρ ϕ ρ ρ ρ ϕ ρ ρ ωε ρ ρ ωε ρ ϕ ρ ρ ρ ϕ ρ ρ ⎧ = ⎪ ⎧ ⎡ ⎤ ⎫ ⎪ ⎨ ′ ′ − ′ × ∇′ ′ − ∇ ⋅′ ′ ∇′ ′ ⎬ ′ ⎢ ⎥ ⎪ ⎩ ⎣ ⎦ ⎭ ⎪⎪ ⎨ ⎪ = ⎪ ⎪ ⎡ ′ ′ + ′ × ∇′ ′ ⎤ ′ ⎣ ⎦ ⎪ ⎪⎩∫
∫
. (1)Here, the line integrals are evaluated as the Cauchy principal values, and ϕ ρ ρ( , ′) is the 2D Green’s function of the homogeneous region given by:
(
)
( )1(
)
0 , 4 i H k ϕ ρ ρ′ = ρ ρ− ′ . (2)Here, H0( )1 is the first-kind Hankel function of order zero, ω is the angular frequency, and k =
ω (με)1/2. Once we have the field expressions on each interface approaching from different regions expressed by Eq. (1), we can match them to obtain a set of boundary integral equations.
Based on the BIEM, the concept of the plane-wave-assisted method is proposed and the formulation details are discussed in the following: As shown schematically in Fig. 1, a unit cell of a periodic structure repeats itself in the x direction and extends infinitely its background materials (ε1, μ1) and (ε3, μ3) in the -y and +y directions, respectively. Region 2
solid lines describing the real structure boundaries, we introduce virtual boundaries as plotted in red dashed lines (C1, C2,..., C6) for enclosing the major part of the unit cell. We define the
region enclosed by these red dashed lines as the interior region. In contrast, the exterior region is defined as the rest part of the unit cell.
Fig. 1. A unit cell of a 2D periodic structure. The unit cell repeats itself in the x direction, with the solid lines for its true boundaries or interfaces. The red dashed lines represent virtual boundaries for applying the periodic boundary conditions and the interior-exterior connection.
On the virtual boundaries C3-C4 and C5-C6, the periodic boundary conditions are used. The
fields on these boundaries in the interior region are related to each other through the Bloch condition. As to the field matching on the virtual boundaries C1 and C2, we develop a novel
method based on the plane-wave expansion technique, which describes the field characteristics in the y-direction and provides a connection between the unknown coefficients of the expanded plane waves and those of the equivalent surface currents on the interior-exterior interfaces (C1 and C2).
In Fig. 2, we show an alternative version of Fig. 1 to focus on the discussions of the interior-exterior interfaces.
Fig. 2. An alternative version of Fig. 1 to focus on the discussions of the interior-exterior interfaces.
In region 0, a plane wave with transverse wavevector kx is incident onto the interior region,
producing the reflected wave. The total electromagnetic field includes both the incident and reflected parts, which can be expanded with a set of plane waves in the Bloch form as
( )
( ) ( )( )
( ) ( )(
)
0 0 2 2 ˆ , ˆ , ˆ , Re( ) 0, Im( ) 0 y x y x L ik y y i k G xinc ref inc
L L
ik y y i k G x
inc ref inc
L y x y y E E E E x y t h e e H H H zH x y zh e e k k k G k k η + − − =− − − + =− ⎧ = + = + ⎪ ⎪ ⎪ = + = + ⎨ ⎪ ⎪ = − + ≥ ≥ ⎪ ⎩
∑
∑
. (3)Here, η is the intrinsic impedance of the background material, h is the coefficient to be solved, and G is the primitive translation wavevector of the 1D grating with period a. The incident wave can be in any form satisfying the Bloch condition. Here we define a vector
ˆ ˆ ˆ
t = ×z k for describing any order of the reflected electric field with
(
)
ˆ ˆ kx G ˆky k x y k k + = − .As to the connection of the plane-wave expansion coefficients, h , in Eq. (3) and the unknown equivalent surface current, J, at the interface y = y0 between regions 0 and 1, we
consider the following equation for J as
( )
( )
( ) ( )
0(
0)
ˆ , 2 1 ˆ , N n n n J x xJ f x N L J x y H x y = ⎧ = = + ⎪ ⎨ ⎪− = − × ⎩∑
. (4) Here, at the interface y = y0, the tangential component of the total magnetic field is related tothe equivalent electric surface current through a set of local linear bases fn= fn,1+ fn,2 given by 1 1 1 ,1 , 0, otherwise n n n n n n x x x x x x x f − − − − ⎧ ≤ ≤ ⎪ − =⎨ ⎪⎩ , 1 1 1 ,2 , 0, otherwise n n n n n n x x x x x x x f + + + − ⎧ ≤ ≤ ⎪ − =⎨ ⎪⎩ . (5)
At a discretization point xn, we have
( )
(
)
( ) 0 0 ˆ ˆ ˆ , ˆ x n N L i k G x n n n inc n n L xJ f′ ′ x y zH x y zh e + ′= =− ⎡ ⎤ = ×⎢ + ⎥ ⎣ ⎦∑
∑
. (6)Here, we consider the matching condition only at the sampled points x0, x1, …, xN-1, because
they can provide the variations of fields or currents in one period. Thus, the transformation between the coefficients of the local linear bases and those of the global plane waves can be easily obtained as
(
)
(
)
0 1 0 0 , , L n inc n n L N n n n inc n n J H x y P h h Q J Q H x y =− − = ⎧ = + ⎪⎪ ⎨ ⎪ = ⎡⎣ − ⎤⎦ ⎪⎩∑
∑
, (7) where i k(x G x)n nP =e + , and Q are the elements of the inverse of the matrix n
[ ]
Pn . In our simulations, the condition number of[ ]
Pn is close to unity that guarantees the numerical stability.Based on Eqs. (3) and (7), we can express the total electromagnetic field in region 0 in terms of the coefficients of the equivalent electric surface currentJ as n
( )
(
)
( ) ( )( )
(
)
( ) ( ) 0 0 1 0 0 1 0 0 ˆ , , ˆ , ˆ , y x y x L N ik y y i k G x inc n n n inc n L n L N ik y y i k G x inc n n n inc n L n E E x y t Q J Q H x y e e H zH x y z Q J Q H x y e e η − − − + =− = − − − + =− = ⎧ = + ⎡ − ⎤ ⎣ ⎦ ⎪⎪ ⎨ ⎪ = + ⎡⎣ − ⎤⎦ ⎪⎩∑ ∑
∑ ∑
. (8)Note that Eq. (8) not only describes the field property in the exterior region of reflection (region 0), but also provides a connection to the fields in the interior region (region 1). The field expressions in the exterior region of transmission (region 4) can also be derived in a similar way without the contribution of the incident field.
When expanding the equivalent surface current on boundary C using the linear bases, 1
care must be taken. For the bases fn plotted in Fig. 2, f0 contains only f0,2, while fN contains
only fN,1. Note that 0 x
ik a N
J =J e based on the Bloch condition.
In order to describe the formulation and to construct the final matrix form in a compact and systematic way, we execute the Galerkin testing procedure [16] first to the field expression on each boundary to get its matrix representation. In other words, on each division segment, a testing local linear basis (5) is multiplied to the electric or magnetic field and integrated over that segment. Then the contributions from various equivalent sources on different boundaries are drawn as basic blocks for constructing the final matrix equation of the whole problem.
In the following, we use three indicators (i, j, p) for specific discussions, where i denotes the boundary on which the testing operators ( E
i
T , H i
T ) are executed, j denotes the boundary on which the equivalent surface currents are located, and p denotes the region in which we apply the field expressions. Note that the equivalent surface currents are related to two indicators (j,
p) and expanded as shown below:
2 , , , , 1 2 , , , 1 ˆ ˆ j j j j j j p p n n s n s n s n s j j j j j p p n n s n s n s J J t f M M z f β α β α = = ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
∑ ∑
∑
∑
, (9)where
(
Jnj,Mnj)
are unknown expansion coefficients and j pβ is an indicator with two possible values (+1, -1) indicating the sign of the current. Another indicator j,
n s
α , with two possible values (1, 0), indicates whether each half basis ,j
n s
f is used or not for a linear basis j
n
f .
For all the boundary fields belonging to regions 1, 2 and 3, we apply Eq. (1) and execute the Galerkin testing procedure to obtain the corresponding matrix representations as follows:
{ }
(
)
( ) ( ){ }
, , , enclosing region enclosing region , , 2 2 2 2 j E p E p i i inc p j p p j p p j p p c j p p p p j p j i ij ij j p H p H p i i inc p j p p j p i T E T E i J M J dl T H T H i M J ωμ ϕ ϕ ϕ ωε ωε ϕ ⎧ ⎡ ⎤ ⎫ ⎪ ′ ′ ′ ′⎪ = ⎨ + ⎢ − × ∇ − ∇ ⋅ ∇ ⎥ ⎬ ⎢ ⎥ ⎪ ⎣ ⎦ ⎪ ⎩ ⎭ ⎡ ⎤ = + ⎣ + ⎦ ′ = + + × ∇∑
∫
∑
b A J B M ( ) ( ) enclosing region enclosing region j p c j p p p j p j i ij ij j p dl ϕ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ ⎫ ⎪ ⎪ ⎡ ⎤ ′⎪ ⎨ ⎬ ⎪ ⎪ ⎣ ⎦ ⎪ ⎩ ⎭ ⎪ ⎪ = + ⎡ + ⎤ ⎣ ⎦ ⎪ ⎩∑
∫
∑
c C J D M , (10) where the matrix elements are given by(
)
{
}
2 2 , , , , , , , , , 1 1 2 2 , , , , , , 1 1 , , ˆ ˆ ˆ 2 ˆ ˆ 2 2 i j i j p i j i i j j j j j j ij mn m r n s c m r m r c p p n s n s p p n s n s p r s p p i j i i j j ij mn m r n s c m r m r c p n s p r s p ij mn m i A t f i t f t f dl dl B t f z f dl dl C α α ωμ β ϕ β ϕ ωε α α β ϕ α = = = = ⎧ ⎡ ⎤ ⎫ ⎪ ′ ′ ′⎪ = ⋅⎨ ⎢ − ∇ ⋅ ∇ ⎥ ⎬ ⎢ ⎥ ⎪ ⎣ ⎦ ⎪ ⎩ ⎭ ′ ′ ⎡ ⎤ = − ⋅ ⎣ × ∇ ⎦ =∑∑
∫
∫
∑∑
∫
∫
{
}
{
}
2 2 , , , , 1 1 2 2 , , , , , 1 1 2 , , , , 1 2 , , , 1 ˆ ˆ 2 ˆ 2 ˆ 2 i j i j i i i j i j j j r n s c m r c n s p n s p r s p i j i j j ij mn m r n s c m r c p p n s p r s p i i i p i m m r m r m r inc c r p i i p i m m r c m r inc r zf t f dl dl D f i f dl dl b t f E dl c zf H dl α β ϕ α α ωε β ϕ α α = = = = = = ⎧ ′ ′ ⎡ ⎤ ⋅ ⎣ × ∇ ⎦ ⎨ ′ ⎡ ⎤ = ⎣ ⎦ = ⋅ = ⋅∑∑
∫
∫
∑∑
∫
∫
∑ ∫
∑ ∫
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ . (11) Similarly for the exterior region, such as region 0, we have the following matrix representations through field expression (8) under the Galerkin testing procedure.{ }
(
)
(
)
( ){ }
(
)
(
)
1 1 1 0 0 1 0 1 1 0 0 0 0 0 1 1 1, 1 1 0 0 1 0 1 1 0 0 0 ˆ , , ˆ , ˆ , x N L i k G x E p E p j p i i inc l n n n inc n L n p p j i i j N i k H p H p j p i i inc n n n inc n n T E T E x y t Q J Q H x y e T H T zH x y z Q J Q H x y e η − + = = = = = = =− = = = = = = = − = = = = = = = ⎧ ⎡ ⎤ ⎫ = ⎨ + ⎣ − ⎦ ⎬ ⎩ ⎭ = + ⎡ ⎤ = + ⎣ − ⎦∑ ∑
∑
b A J ( ) 0 0 1 1 1, 1 x L G x L p p j i i j + =− = = = = = = ⎧ ⎪ ⎪ ⎪ ⎪⎪ ⎨ ⎪ ⎧ ⎫ ⎪ ⎨ ⎬ ⎪ ⎩ ⎭ ⎪ = + ⎪⎩∑
c C J , (12)( ) ( )
(
)
2 0 1 1 1 1, 1, , , , 1 2 0 1 1 1 1, 1, , , , 1 0 1 1 1 0 1, , , , 0 ˆ ˆ ˆ ˆ ˆ , ˆ x x L i k G x p i i i i j mn m r m r m r n r L L i k G x p i i i i j mn m r m r m r n r L p i i i p i m m r m r m r inc n in A t f t Q e dx C t f zQ e dx b t f E x y t Q H α η α α η + = = = = = = = =− + = = = = = = = =− = = = = = = ⎡ ⎤ = ⋅⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = ⋅⎢ ⎥ ⎣ ⎦ = ⋅ −∑
∫
∑
∑
∫
∑
(
)
( )(
)
(
)
( ) 1 1 1 2 0 0 1 0 1 2 0 1 1 1 0 0 1, , , , 0 0 1 0 , ˆ ˆ , ˆ , x x N L i k G x p c n r L n N L i k G x p i i i p p i m m r m r m r inc n inc n r L n x y e dx c α t f zH x y zQ H x y e dx − + = = =− = − + = = = = = = = = =− = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎡ ⎤ ⎪ ⎢ ⎥ ⎪ ⎣ ⎦ ⎪ ⎡ ⎤ ⎪ = ⋅ − ⎢ ⎥ ⎪ ⎣ ⎦ ⎩∑
∫
∑ ∑
∑
∫
∑ ∑
(13)In Eq. (13), the integration of the term 1 ( ) , ( )
x
i k G x i
m r
f = x e + can be analytically calculated. Moreover, it is only related to the summation L
( )
L =−
⋅
∑
. This remarkably saves the computation time comparing with other methods used for considering the periodicity in boundary integral equation formulation.As to matching the boundary conditions, there are three types of matching. (1) The matching for boundaries inside the interior region. (2) The periodic boundary matching. (3) The matching for boundaries between the interior region and the exterior region. With the above matrix representations (Eqs. (10) and (12)) as basic blocks in boundary matching, we can systematically construct the entire matrix equation for the whole problem when programming. The entire matrix equation can finally be expressed in the form
⎡⎢ ⎤ ⎡ ⎤ ⎡ ⎤⎥ ⎢ ⎥ ⎢ ⎥=
⎣ ⎦ ⎣ ⎦ ⎣ ⎦
A B J b
C D M c
, (14)
where matrix A is composed of the sub-matrices A (combinations of ij A with different ijp
p ), J is a column vector consisting of sub-column vectors j
J , and b is also a column vector
with their elements b being combinations of i b with different p . Similar meanings hold for ip
B , C , D , M and c .
Actually, in the final form, the matrix equation has reshaped after eliminating some unnecessary columns or rows of dependent terms due to the nature of periodicity. Note also that once the equivalent surface currents are obtained through matrix inversion, the reflected and transmitted fields of any diffraction order in the exterior region can be readily calculated with Eq. (3) without numerical calculation of the radiation integral. This is one of the advantages of our method.
3. Simulation results
3. A. Numerical verifications
To verify the numerical accuracy of the proposed method, we first consider a reflection-type highly-conducting metallic-lamellar grating with the dielectric constant εγ = -100 (n = i10), as depicted in Fig. 3.
Fig. 3. A reflection-type metallic-lamellar grating with a grating period a, a groove depth d, and a groove width g. The incident wave is p-polarized with an incident angle θi.
The numerical convergence in the calculation of diffraction efficiency of a highly-conducting lamellar grating has been a widely studied issue associated with the CWM, which is based on the concept of plane-wave expansion. In Fig. 4(a), we show the reflectance of the -1 order versus the groove width based on the CWM with the same parameters used in [9]. Both the period a and the groove depth d are set to be 500 nm. The groove width g is taken to be a varying parameter from 5 to 460 nm. The incident field is TE- or p-polarized with the wavelength λ0 = 632.8 nm and the incidence angle θi = 30
o
. Quite many sharp reflectance peaks and dips can be seen in Fig. 4(a), which have been proved spurious and are not the real physical features, even though we have used Li’s Fourier factorization in our calculation [17]. For comparison, Fig. 4(b) shows the result calculated with the proposed method. Here, one can see that the major reflectance variation curve is almost the same as that in Fig. 4(a). The major difference is that those spurious features in Fig. 4(a) disappear. Our result also matches well with that calculated with the rigorous modal method [9]. The superiority of the proposed method to the CWM in this application is due to the fact that we use the plane-wave expansion technique only in the homogeneous regions exterior to the grating region where the rapid field variation is carefully treated with the BIEM. Note that the non-uniform mesh is used in the proposed method for all the boundaries of the problem, except the boundaries between the interior and exterior regions (C1 and C2 in Fig. 1), on which 11 equally sampled
points are taken for describing the equivalent electric surface current, with the same number of the plane waves used for expanding the outgoing waves in the exterior region. For this calculation, about 150 local linear bases in total are used. The non-uniform mesh is especially suitable for modeling the problem with fine geometry. In this case, the small groove width can be described well by using the non-uniform mesh, with the total number of the unknowns not increased much.
(a) (b)
Fig. 4. Reflectance of the -1 order versus the groove width, calculated with (a) the coupled-wave method (CWM) and with (b) the proposed method.
Another accuracy verification is shown in Fig. 5 with the same problem geometry shown in Fig. 3.
(a) (b)
Fig. 5. Zero-order reflectance versus wavelength. Solid lines denote the results of our method and empty squares represent those calculated with the CWM. (a) θi = 0
o
. (b) θi = 80
o
.
The groove width is fixed at 250 nm. The metal is changed into Ag, whose dielectric constant can be described by the Drude model 1 2/[ ( )]
r p i p
ε = −ω ω ω γ+ , with the plasma frequency ωp
= 1.36x1016 s-1 (corresponding to the energy of 9 eV) and the damping frequency γp = ωp/90.
The result of the zero-order reflectance versus the wavelength of a normally incident wave is shown in Fig. 5(a). The solid line indicates the result calculated with the proposed method and the empty squares denote that with the CWM. They match almost perfectly with each other. We further take a large incident angle of 80o for numerical verification with the result shown in Fig. 5(b). In this case, the proposed method also leads to almost exactly the same result as that calculated by the CWM. In this verification, by comparing the results with those of the CWM, we verified the accuracy of the proposed method in calculating the diffraction efficiency for a real metal grating with a large angle of incidence. The CWM also performs well for this case. However, it suffers from the convergence difficulty when dealing with the arbitrarily shaped metallic gratings, which can be easily treated by our method. Other unstructured-mesh numerical techniques, such as the FEM, are suitable for modeling the arbitrarily shaped geometry. However, they might not perform well dealing with grating diffraction at a large angle of incidence if the perfectly matched layer (PML) is used to truncate the computation domain in the perpendicular direction, for the PML does not perform well for an outgoing wave with a large incident angle. In summary, our method is capable of dealing with the arbitrarily shaped metallic gratings, even at a large angle of incidence.
3. B. Simulation results of a wavy structure
After verifying the accuracy of the proposed method, we simulate a wavy layered structure with a substrate-metal-cover-air architecture. The metallic wavy layered structure has recently been used in organic light-emitting devices (OLEDs) for enhancing the emission rate due to the high intensity of the SPP mode and extracting the light via the SPP coupling [18-20]. The coupling mechanisms have been theoretically and experimentally investigated [21, 22]. For an asymmetric substrate-metal-air layered structure, the cross coupling is considered as an important coupling mechanism, in which the two SPP modes at metal/substrate and metal/air interfaces are coupled via the Bragg scattering due to the structure periodicity [19]. However, the cross-coupled SPPs have been shown to occur only over a narrow range of emission wavelengths and angles. Thus an additional dielectric cover layer is introduced for utilizing
the wavy layered structure (substrate-metal-cover-air) to study the enhanced optical transmission by the coupling mechanism of the coupled SPPs. Note that the cover layer is used to symmetrize the layered structure, making the field intensity of each SPP mode appropriately distributed on both interfaces of the metal film. This may result in an efficient coupling effect.
For numerical implementations, the substrate is assumed to be GaN (with the refractive index of 2.5), covered by Ag. The material of the dielectric cover has a refractive index nc,
which is left as a varying parameter. The geometry of the structure is shown in Fig. 6, where a is the grating period, and A is the amplitude of the sinusoidal variation. The thicknesses of the silver layer and the dielectric cover layer are denoted by tm and tc, respectively. A wave is
normally incident toward the +y direction from the substrate.
Fig. 6. Schematic illustration of a wavy layered structure. The grating period is a, and the amplitude of the sinusoidal variation is A. The cover layer, with its thickness denoted by tc, has
a refractive index nc. The thickness of the Ag film is tm.
With this problem geometry, we will examine the zero-order transmittance and power loss versus the wavelength of the incident light for different structure parameters. Here, the power loss is calculated by subtracting all the diffracted power from the incident power. First, a set of standard parameters is chosen as a = 218 nm, A = 20 nm, tm = 44 nm, tc = 100 nm, and nc =
2.5. We first assume that the cover has a refractive index the same as that of the GaN substrate and is thicker than the decay length of the SPP mode. Therefore, the problem geometry becomes symmetric. The assumption of 218 nm for the period a means to excite an SPP mode at λ0 = 650 nm according to the diffraction relation,
kx = n (2π/a) + kSPP. (15)
Here, kx is the in-plane wavevector of the incident wave, kSPP is the wavevector of the SPP
mode [1] at the planar Ag/GaN interface, and n is an appropriate integer. The result from Eq. (15) is an estimate based on the assumptions that the corrugation just causes a small perturbation and the Ag film is thick enough such that the optical activities of its upper and lower interfaces can be regarded independent.
By taking the silver thickness tm as a varying parameter, changing from 44 to 28 nm, the
transmittance and power loss versus wavelength are shown in Fig. 7. In Fig. 7(a), the transmittance peak location blue shifts as tm decreases. However, a shoulder is gradually
formed on the long-wavelength side. The maximum peak reaches 0.78 at tm = 28 nm. At tm =
44 nm, the transmittance peak is located at 656 nm, which is close to the designed wavelength from Eq. (15). In Fig. 7(b), one can see a major and a minor peak in all the power loss spectra. The major loss peak presents a red-shift trend while the minor one shows an opposite trend in decreasing tm. Comparing Figs. 7(a) and (b), one can see that the transmittance peak is always
close to the minor loss peak. Therefore, one can separate the transmittance peak and the major loss peak by narrowing down the Ag layer. In this case, because of the nearly symmetric configuration with the normal incidence, the coupled SPPs dominate the transmission behavior. Each coupled SPP mode, with its field distributed around the two interfaces of the
Ag film, is Bragg scattered into the radiation modes, coupling the light through the Ag film. The dispersion curves of two major coupled SPP modes (even and odd) split further from each other with decreasing metal thickness, accounting for the broadening of the transmittance spectrum and the blue-shifted transmittance peak. Actually, the gradually formed shoulder of in transmittance spectrum is due to another transmittance peak corresponding to one of the two major coupled SPP modes.
(a) (b) (c) (d)
Fig. 7. (a). Zero-order transmittance and (b) power loss versus wavelength. The results are calculated for different Ag-film thicknesses. (c) Magnitude distributions of the magnetic field at the transmittance peak of tm = 44 nm (labeled by (1) in (a)) and the two peaks of tm = 28 nm
(labeled by (2) and (3) in (a)). (d) Distributions of the time-average Poynting vector corresponding to the three plots of (c).
For the transmittance peak of tm = 44 nm (labeled by (1)) and the two peaks of tm = 28 nm
(labeled by (2) and (3)), we illustrate their magnitude distributions of the magnetic field in Fig. 7(c), with the amplitude of the incident magnetic field set at unity. The corresponding
the incoming and the outgoing waves are coupled with the SPP modes via mainly the high-intensity regions around both sides of the Ag film. Such a coupling implies that for the layered wavy structure, the SPP mode indeed plays an important role transferring the energy through the metal film. This is different from the situation of the slit grating structure, in which transmittance is mainly enhanced by the cavity-mode component in the slits and inhibited by the horizontal SPP mode [23].
Then, we take the grating amplitude A as a varying parameter, changing from 24 to 8 nm, to give Fig. 8, in which a slight blue shift of the transmittance peak in decreasing A can be seen. The transmittance peak level reduces considerably when A becomes smaller. The power loss variation is shown in Fig. 8(b).In this case, varying the grating amplitude does not break the symmetry of the structure. Hence, the coupled SPP modes remain symmetric across the Ag film. However, the coupling effect between the coupled SPP modes and the radiation modes will be reduced when the grating amplitude further decreases. As a result, either the incident wave is reflected or its energy is dissipated due to the metallic loss.
(a) (b)
Fig. 8. (a). Zero-order transmittance and (b) power loss versus wavelength. The results are calculated for different amplitudes of the sinusoidal shape.
Next, we change the refractive index of the cover layer nc from 2.5 into 2.7, with other
parameters remaining the same with the aforementioned standard ones. Now because of the problem geometry becomes asymmetric, the coupled SPP modes tend to asymmetrically distribute across the Ag film, resulting in a reduction of the coupling efficiency compared with the result of the symmetric structure. Besides, the increased refractive index of the cover layer may lead to a red shift of the transmittance peak. A possible way to recover the transmittance and blue shift back its peak is reducing the thickness of the cover. This is in the sense of averaging the index of the cover layer and that of the air to an effective value. For this purpose, we vary the cover layer thickness tc from 100 to 50 nm to see the effects on the
transmittance and power loss. The numerical results are shown in Fig. 9. The transmittance peak in Fig. 9(a) is blue-shifted, as expected, when tc is decreased. The peak level reaches a
maximum at about 0.7 when tc = 70 nm and then decreases as tc becomes thinner. In Fig. 9(b),
the major power loss peak is also blue-shifted with decreasing tc. The magnitude distributions
of the magnetic field at the peaks of tc = 100 nm, tc = 70 nm, and tc = 50 nm are shown in Fig.
9(c). The mode distribution is rather symmetric for the case of tc = 70 nm, consisting with the
(a) (b)
(c)
Fig. 9. (a). Zero-order transmittance and (b) power loss versus wavelength. The results are calculated for different thicknesses of the cover layer. (c) Magnitude distributions of the magnetic field at the peaks of tc = 100 nm, tc = 70 nm, and tc = 50 nm.
With the asymmetric structure of nc = 2.7 and a = 218 nm, we change other parameters
into A = 20 nm, tm = 28 nm, and tc = 70 nm, based on the results in Figs. (7)-(9) to optimize
the transmittance level. The results are shown in Fig. 10.
(a) (b) (c)
Fig. 10. (a). Zero-order transmittance and power loss versus wavelength. The structure parameters used are a = 218 nm, A = 20 nm, tm = 28 nm, tc = 70 nm, and nc = 2.7. (b)
Magnitude distribution of the magnetic field at λ0 = 630 nm. (c) Distribution of the
time-average Poynting vector at λ0 = 630 nm.
By adjusting the parameters, the transmittance shown in Fig. 10(a) becomes larger than 0.7 over a wide wavelength range from 610 to 700 nm. The transmittance peak value reaches 0.83 at 630 nm. Figure 10(b) shows the magnitude distribution of the magnetic field at the
illustrated in Fig. 10(c). The incident wave is coupled through the metal film and diffracted into the free-space radiation mode in a similar way mentioned in Fig. 7(c) and (d), even when the cover layer is with a refractive index larger than that of the substrate. Actually, in the sense of averaging the refractive index of the cover layer and that of the air, the mechanism of the coupled SPPs is still effective as long as the refractive index of cover layer is not much larger than that of the substrate. Though we use the sinusoidal shape as the grating geometry to demonstrate the capability of the proposed numerical method, similar results can also be obtained in other sinusoid-like shapes.
4. Conclusions
We have developed a novel hybrid method, based on the BIEM and assisted with the plane-wave expansion technique. This method is flexible and applicable to various 2D grating structures, particularly those with metal/dielectric interfaces. Numerical accuracy has been verified with highly conducting and practical metals. This method can provide quite accurate results in the cases of large incidence angles. The transmittance of a metal/dielectric wavy structure has been investigated with various structure parameters. High transmittance up to 0.83 was obtained. The numerical results also showed that the transmittance enhancement was due to the mechanismof the coupled SPP modes, which was different from the situation of a slit grating structure.
Acknowledgments
This research was supported by the National Science Council, the Republic of China, under the grants NSC 95-2221-E-002-286 and NSC 95-2120-M-002-012, and by US Air Force Scientific Research Office under the contract AOARD-06-4052.
