About this Journal Submit a Manuscript Table of Contents
International Journal of Optics
Volume 2012 (2012), Article ID 532316, 19 pages
doi:10.1155/2012/532316
Research Article

Green Function Formulism for Electromagnetic Wave Generated in Nanostructured Metamaterial of Finite Thickness: Isotropy and Anisotropy

1Department of Physics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
2School of Physical Science and Technology, Suzhou University, No. 1 Shizi Street, Jiangsu, Suzhou 215006, China

Received 31 May 2012; Accepted 2 July 2012

Academic Editor: Zhaolin Lu

Copyright © 2012 Shunbo Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A Green function formulism is developed to calculate the electromagnetic fields generated by sources embedded in nanostructured medium which could be represented by an effective electric permittivity tensor with finite thicknesses. The method begins with the decomposition of the generated mode into the eigenmodes in the medium, which have definite dispersions. To account the interface effect at boundaries, especially the mode conversion at the interface between anisotropic media, mode expansion method is combined into the theoretical framework. Thus, the electromagnetic wave in any given position can be gotten clearly. The formulism can provide conveniences of studying the novel properties of nanostructured metamaterials.

1. Introduction

Artificially constructed metamaterials [1, 2] have unique electromagnetic properties which cannot be obtained in natural materials, such as artificial magnetism, negative refraction [3], and near-field focusing [4]. Recently, active research areas related to metamaterials attracted much attention [5, 6], for instance, high/low epsilon metamaterials, designer dispersion, transformation optics metamaterials, switchable metamaterials, amplifying metamaterials, sensor metamaterials, nonlinear metamaterials, and quantum metamaterials [7]. The origin of novel property of structured metamaterial comes from the modification of effective permittivity and permeability. There are many developed methods to determine the effective permittivity and permeability [8, 9]. Nanostructured metamaterials have wide applications in real world, such as manipulation of thermal radiation [10], plasmonic biosensor [11], artificial magnetic material [12], and superlens imaging [13]. Patterned metal at the interface between two media will generate surface-plasmon polaritons which alter surface waves and subsequently change radiative heat transfer, coherence properties, and Casimir forces [14].

To obtain the electromagnetic property of a nanostructured metamaterial, one would have to calculate Green function of the given system [15, 16]. The conventional formulism [17] can handle homogeneous isotropic media very well, and the interface effects can be trivially accounted by introducing Fresnel coefficients. However, a lot of metamaterials exhibit anisotropy properties due to the special arrangement of structure [18]. Thus, the story becomes quite complex, when anisotropy of media is introduced. One distinguishing phenomenon of anisotropic media is the mode conversion [19, 20] happening at the interface between anisotropic media. Under the situation, the previous formulism is no longer suitable and convenient, since the mode conversion leads to an infinite multiplication for the Fresnel transmission and reflection coefficients. Thus, to study the coherence property of thermally radiative field from anisotropic media consisted by subwavelength structures with finite thickness, we have firstly to build up a general Green function theory working well in both isotropic media and anisotropic media. Exactly the motivation drives us to the present work. Since the development of metamaterial and nanooptics paves the way to design innovative structures, we herein proposed a Green function method to calculate electromagnetic wave generated by a structured metamaterials which have definite permittivity and permeability. In the following sections, we will firstly present our formulism for an arbitrary anisotropic medium with finite thickness in Section 2. To prove the correction of the method, we apply it to an isotropic medium case to recover the results in [17] in the same section. In Section 3, the Green functions linking spatial fields and sources in uniaxial and biaxial medium of finite thickness are calculated by following the general formulism shown in Section 2. A brief conclusion remark is made in the last section.

2. Theoretical Formulism and Its Correction

To begin with, we consider a piece of anisotropic medium with a finite thickness placed in vacuum (Figure 1), so that random sources locate only in the anisotropic medium. In the present work, it is assumed the relative permittivity tensor is anisotropic but its relative permeability is isotropic (the current formulism can be applied to the case that both permittivity and permeability are anisotropic, but we stay at the current situation to achieve a better demonstration of the formulism), and they can be written as: 𝜀 𝜀 = 1 0 0 0 𝜀 2 0 0 0 𝜀 3 , 𝜇 = 𝜇 . ( 1 ) The relative permittivity and relative permeability of vacuum are equal to 1. Once the information is given, for a given frequency the electromagnetic property is totally governed by Maxwell equations [17]: 𝜀 𝑟 𝐸 𝑃 , 𝑟 = 4 𝜋 𝑟 × 𝐵 𝑟 + 𝑖 𝜔 𝜇 𝑟 𝜀 𝑟 𝐸 𝑟 = 𝑖 4 𝜋 𝜔 𝜇 𝑟 𝑃 𝑟 + 4 𝜋 𝜇 𝑟 𝑀 𝐵 𝐸 𝐵 × 𝑟 𝑟 = 0 × 𝑟 𝑖 𝜔 𝑟 = 0 , ( 2 ) where 𝜀 𝑟 represents the relative permittivity tensor in different regions shown in Figure 1, 𝜇 𝑟 means the permeability in different regions, 𝜔 represents 𝜔 / 𝑐 (angular frequency over velocity of light in vacuum), 𝑃 ( 𝑟 ) denotes the random electric dipole source distribution in space, and 𝑀 ( 𝑟 ) is for random magnetic dipole source distribution.

532316.fig.001
Figure 1: Schematic illustration of the geometry of the system under considering shows the thickness ( ) of the anisotropic (isotropic) layer and the coordinate system.
2.1. Propagation Property in Homogeneous Medium

One can extract the information about the eigen-plane-wave modes and their corresponding dispersion relations by doing operations on the homogeneous Maxwell equations (dropping the source terms about 𝑃 and 𝑀 ), and one has 𝐸 × × 𝑟 𝜔 2 𝜇 𝑟 𝜀 𝑟 𝐸 𝑟 = 0 . ( 3 ) Taking into account the plane wave ansatz, we can have the matrix equation determining the eigenmodes for given regions: 𝜔 2 𝜀 𝑥 𝜇 𝑟 𝑘 2 𝑦 𝑘 2 𝑧 𝑘 𝑥 𝑘 𝑦 𝑘 𝑥 𝑘 𝑧 𝑘 𝑥 𝑘 𝑦 𝜔 2 𝜀 𝑦 𝜇 𝑟 𝑘 2 𝑥 𝑘 2 𝑧 𝑘 𝑦 𝑘 𝑧 𝑘 𝑥 𝑘 𝑧 𝑘 𝑦 𝑘 𝑧 𝜔 2 𝜀 𝑧 𝜇 𝑟 𝑘 2 𝑦 𝑘 2 𝑥 × 𝐸 𝑥 𝐸 𝑦 𝐸 𝑧 = 0 , ( 4 ) where 𝜀 𝑥 = 𝜀 1 , 𝜀 𝑦 = 𝜀 2 , 𝜀 𝑧 = 𝜀 3 in region II, 𝜀 𝑥 = 𝜀 𝑦 = 𝜀 𝑧 = 1 in regions I and III, 𝜇 𝑟 = 𝜇 in region II, and 𝜇 𝑟 = 1 in regions I and III. For a nontrivial mode in the medium, the determinant of the coefficient matrix in (4) should vanish. Taking the homogeneous form of the first Maxwell equation in (2) as a constraint, the equation determining eigenmodes is arrived at 𝜔 4 𝜔 2 𝑘 2 𝑥 + 𝑘 2 𝑦 𝜀 𝑧 𝜇 𝑟 + 𝑘 2 𝑥 + 𝑘 2 𝑧 𝜀 𝑦 𝜇 𝑟 + 𝑘 2 𝑧 + 𝑘 2 𝑦 𝜀 𝑥 𝜇 𝑟 + 𝑘 2 𝑥 𝜀 𝑦 𝜀 𝑧 𝜇 𝑟 2 + 𝑘 2 𝑦 𝜀 𝑥 𝜀 𝑧 𝜇 2 𝑟 + 𝑘 2 𝑧 𝜀 𝑥 𝜀 𝑦 𝜇 2 𝑟 × 𝑘 2 𝑥 + 𝑘 2 𝑧 + 𝑘 2 𝑦 = 0 . ( 5 ) From (5) 𝑘 𝑧 ( 𝜔 , 𝑘 𝑥 , 𝑘 𝑦 ) can be solved out as the function of 𝜔 , 𝑘 𝑥 and 𝑘 𝑦 , and there are 4 solutions for 𝑘 𝑧 , namely, { 𝑘 𝑧 , 1 , 𝑘 𝑧 , 1 ; 𝑘 𝑧 , 2 , 𝑘 𝑧 , 2 } ( 𝑘 𝑧 , 1 and 𝑘 𝑧 , 2 are positive real or positive imaginary numbers to denote the wave propagating or decaying along the positive 𝑧 direction.). To get the corresponding polarizations of the eigenmodes, we can solve (4) with the help of the homogeneous form of the first Maxwell equation in (2) and the corresponding dispersion relation of the given eigenmode 𝑘 2 𝑧 , 𝑖 [ 𝑘 𝑧 , 𝑖 ( 𝜔 , 𝑘 𝑥 , 𝑘 𝑦 ) ] 2 = 0 , where the index i is 1 (2), when we solve the corresponding polarization of the first (second) eigenmode. For the convenience, we denote the unit vectors of them by ̂ 𝑒 1 , ± and ̂ 𝑒 2 , ± , and conveniently the corresponding unit vectors of electric displacement vectors and magnetic field vectors are gotten by 𝑑 𝑖 , ± = 𝜀 𝑟 ̂ 𝑒 𝑖 , ± / | 𝜀 𝑟 ̂ 𝑒 𝑖 , ± | and 𝑖 , ± = 𝑘 ± × ̂ 𝑒 𝑖 , ± 𝑘 / | ± × ̂ 𝑒 𝑖 , ± | . Then physics of propagating wave in the medium is quite clear: for propagating wave in the medium with arbitrary polarization its electric displacement vector is always perpendicular to wave vector; from the previous analysis, the electric displacement vector can be always decomposed in the plane perpendicular to the wave vector 𝑘 into two components along 𝑑 1 and 𝑑 2 , respectively. These two components obey different but definite dispersions, once the three parameters ( 𝜔 , 𝑘 𝑥 , and 𝑘 𝑦 ) are known.

2.2. Treatment of Sources

We now restore the source terms in Maxwell equations so that sources enter the picture. The analysis in the last section illuminates that the property of propagating modes in one medium is totally clear, when any three parameters of 𝜔 , 𝑘 𝑥 , 𝑘 𝑦 , and 𝑘 𝑧 are determined. On the other hand, when the spatial distribution and time evolution of the random sources are known, Fourier transformation can always help us to count the modes in a desired way (since we have limited in a given frequency case, 𝜔 is known now): 𝑃 = 𝑟 𝑖 = 1 , 2 𝑑 𝑧 0 𝑒 𝑖 𝑘 𝑧 , 𝑖 𝑧 0 𝑑 𝑘 ( 2 𝜋 ) 3 𝑃 𝑘 𝑖 𝑒 𝑖 𝑘 𝑖 𝑟 = 𝑑 𝑧 0 𝑑 𝜅 ( 2 𝜋 ) 2 𝛿 𝑧 𝑧 0 𝑃 𝜅 𝑒 𝑅 𝑖 𝜅 , 𝑀 = 𝑟 𝑖 = 1 , 2 𝑑 𝑧 0 𝑒 ± 𝑖 𝑘 𝑧 , 𝑖 𝑧 0 𝑑 𝑘 ( 2 𝜋 ) 3 𝑀 𝑘 𝑖 𝑒 𝑖 𝑘 𝑖 𝑟 = 𝑑 𝑧 0 𝑑 𝜅 ( 2 𝜋 ) 2 𝛿 𝑧 𝑧 0 𝑀 𝜅 𝑒 𝑅 𝑖 𝜅 , ( 6 ) where 𝑧 0 [ 0 , ] , since there should be no random sources in vacuum. Here and below it is convenient to put 𝑅 = ( 𝑥 , 𝑦 ) . Such a transformation has clear physical meanings: the modes generated from sources are counted according to their tangential wave vectors, and to make the classification possible the price of cutting the medium into planes along 𝑧 direction is paid. The advantage of the classification is obvious: for each mode generated we can decompose it into two eigenmodes with definite dispersion relations and the same tangential wave vector.

Similarly, the electromagnetic field of a given position in space can be expressed as 𝐸 = 𝑅 , 𝑧 𝑑 𝜅 ( 2 𝜋 ) 2 𝐸 𝑅 , 𝐵 = 𝜅 , 𝑧 e x p 𝑖 𝜅 𝑅 , 𝑧 𝑑 𝜅 ( 2 𝜋 ) 2 𝐵 𝑅 , 𝜅 , 𝑧 e x p 𝑖 𝜅 ( 7 ) where 𝐸 = 𝜅 , 𝑧 𝑑 𝑧 0 𝐺 𝐸 𝑃 𝜅 ; 𝑧 𝑧 0 𝑃 𝜅 , 𝑧 0 + 𝑑 𝑧 0 𝐺 𝐸 𝑀 𝜅 ; 𝑧 𝑧 0 𝑀 𝜅 , 𝑧 0 , 𝐵 = 𝜅 , 𝑧 𝑑 𝑧 0 𝐺 𝐵 𝑃 𝜅 ; 𝑧 𝑧 0 𝑃 𝜅 , 𝑧 0 + 𝑑 𝑧 0 𝐺 𝐵 𝑀 𝜅 ; 𝑧 𝑧 0 𝑀 𝜅 , 𝑧 0 . ( 8 ) In (8), 𝐺 𝐸 𝑃 , 𝐺 𝐵 𝑃 , 𝐺 𝐵 𝑀 , and 𝐺 𝐵 𝑀 are Green functions linking the field in space with the source for a given tangential wave vector 𝜅 . Once getting the relation between fields in space and source, the forms of Green functions are determined spontaneously. Because of the linearity of Maxwell equation, we can focus to one mode with tangential wave vector 𝜅 generated from a source plane labeled as 𝑧 0 , and to account all modes one needs only to integrate all 𝜅 and 𝑧 dimensions of the anisotropic medium.

2.3. Mode Expansion in Each Region

To determine the field generated from source plate locating at 𝑧 0 with tangential wave vector 𝜅 , the mode expansion method would be found useful [21, 22]. Since we have assumed the anisotropy of permeability and isotropy of permittivity, mode expansion according to magnetic field may be more convenient, and they are demonstrated in Figure 2. For the given frequency 𝜔 and tangential wave vector 𝜅 , the fields in region II can be expanded as 𝐵 I I = 𝐵 𝜔 , 𝜅 1 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐵 2 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 × Θ 𝑧 𝑧 0 + 𝐵 1 , 𝑏 , + + 𝐵 2 , 𝑏 , + + 𝐵 1 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐵 2 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 𝑧 × Θ 0 + 𝐵 𝑧 1 , 𝑏 , + 𝐵 2 , 𝑏 , + 𝐵 𝛿 𝑧 𝑧 0 𝑅 , 𝐷 e x p 𝑖 𝜅 I I = 𝐷 𝜔 , 𝜅 1 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐷 2 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 × Θ 𝑧 𝑧 0 + 𝐷 1 , 𝑏 , + + 𝐷 2 , 𝑏 , + + 𝐷 1 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐷 2 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 𝑧 × Θ 0 + 𝐷 𝑧 1 , 𝑏 , + 𝐷 2 , 𝑏 , + 𝐷 𝛿 𝑧 𝑧 0 𝑅 , 𝐸 e x p 𝑖 𝜅 I I = 𝐸 𝜔 , 𝜅 1 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐸 2 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 × Θ 𝑧 𝑧 0 + 𝐸 1 , 𝑏 , + + 𝐸 2 , 𝑏 , + + 𝐸 1 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐸 2 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 𝑧 × Θ 0 + 𝐸 𝑧 1 , 𝑏 , + 𝐸 2 , 𝑏 , + 𝐸 𝛿 𝑧 𝑧 0 𝑅 , e x p 𝑖 𝜅 ( 9 ) where 𝐵 𝑖 , 𝑔 , ± = 𝐵 𝑖 , 𝑔 , ± 𝑅 e x p 𝑖 𝜅 e x p ± 𝑖 𝑘 𝑧 , 𝑖 𝑧 𝑖 , ± , 𝐵 𝑖 , 𝑏 , ± = 𝐵 𝑖 , 𝑏 , ± 𝑅 e x p 𝑖 𝜅 e x p ± 𝑖 𝑘 𝑧 , 𝑖 𝑧 𝑖 , ± , 𝐷 𝑖 , 𝑔 , ± = 1 𝐵 𝑖 𝜔 𝜇 × 𝑖 , 𝑔 , ± , 𝐸 𝑖 , 𝑔 , ± = 𝜀 1 𝐷 𝑖 , 𝑔 , ± 𝐷 𝑖 , 𝑏 , ± = 1 𝐵 𝑖 𝜔 𝜇 × 𝑖 , 𝑏 , ± , 𝐸 𝑖 , 𝑏 , ± = 𝜀 1 𝐷 𝑖 , 𝑏 , ± . 𝑖 = 1 , 2 , ( 1 0 ) 𝐵 ( 𝐷 , 𝐸 ) is to account the local (nonpropagating) field locating on the source plane, and Θ ( 𝑧 𝑧 0 ) is step function to denote the directions of modes launching from the source plane. In regions I and III, the eigenmodes are s- and p-polarized modes, and field can be expanded as 𝐵 I = 𝐵 𝜔 , 𝜅 𝑠 , 𝑠 , I + 𝐵 𝑝 , 𝑝 , I 𝑅 × e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 , 𝐷 I = 𝐸 𝜔 , 𝜅 I = 𝐵 𝜔 , 𝜅 𝑠 , ̂ 𝑒 𝑠 , I + 𝐵 𝑝 , ̂ 𝑒 𝑝 , I 𝑅 × e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 , 𝐵 I I I = 𝐵 𝜔 , 𝜅 𝑠 , + 𝑠 , + I I I + 𝐵 𝑝 , + 𝑝 , + I I I 𝑅 × e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 , 𝐷 I I I = 𝐸 𝜔 , 𝜅 I I I = 𝐵 𝜔 , 𝜅 𝑠 , + ̂ 𝑒 𝑠 , + I I I + 𝐵 𝑝 , + ̂ 𝑒 𝑝 , + I I I 𝑅 × e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 , ( 1 1 ) where the unit vectors are defined as ̂ 𝑒 𝛼 𝑠 , ± = 𝑘 𝑠 , ± × ̂ 𝑧 | | 𝑘 𝑠 , ± | | = × ̂ 𝑧 𝜅 c o s 𝜃 𝜅 s i n 𝜃 ± 𝑘 𝑧 × 0 0 1 = 0 , s i n 𝜃 c o s 𝜃 𝛼 𝑠 , ± = 𝑘 𝑠 , ± × ̂ 𝑒 𝑠 , ± | | 𝑘 𝑠 , ± × ̂ 𝑒 𝑠 , ± | | = 1 𝜔 𝜀 𝑟 𝜇 𝑟 ± 𝑘 𝑧 c o s 𝜃 ± 𝑘 𝑧 , s i n 𝜃 𝜅 𝛼 𝑝 , ± = 𝑘 𝑝 , ± × ̂ 𝑧 | | 𝑘 𝑝 , ± | | = × ̂ 𝑧 𝜅 c o s 𝜃 𝜅 s i n 𝜃 ± 𝑘 𝑧 × 0 0 1 = 0 , s i n 𝜃 c o s 𝜃 ̂ 𝑒 𝛼 𝑝 , ± = 𝑘 𝑝 , ± × 𝑝 , ± | | | 𝑘 𝑝 , ± × 𝑝 , ± | | | = 1 𝜔 𝜀 𝑟 𝜇 𝑟 𝑘 𝑧 c o s 𝜃 𝑘 𝑧 𝜅 , s i n 𝜃 ( 1 2 ) where 𝛼 = I , I I , I I I denotes the region where the unit vectors belong.

532316.fig.002
Figure 2: Schematic illustration of eigenmode expansion in three regions: 𝐵 𝑠 , and 𝐵 𝑝 , denote the eigen-plane-wave modes in region I, and the minus sign is to denote the direction of propagating is along 𝑧 − direction; 𝐵 1 , 𝑔 , ± in region II denotes the first eigenmode generated from the source plane, 𝐵 2 , 𝑔 , ± is to denote the second eigenmode generated from the source plane, 𝐵 1 , 𝑏 , ± in region II expresses the first eigenmode bouncing back at the upper and lower boundaries, and 𝐵 2 , 𝑏 , ± gives the second eigenmode bouncing back at the upper and lower boundaries (the plus and minus signs denotes the propagating along 𝑧 + and 𝑧 − directions); 𝐵 𝑠 , + and 𝐵 𝑝 , + denote the eigen-plane wave modes in region III, and the minus sign is to denote the direction of propagating is along 𝑧 + direction.
2.4. Determination of Generated Field in Each Region and Green Functions

The amplitudes of the four eigenmodes launching from source plane for the given frequency 𝜔 and tangential wave vector 𝜅 can be determined by putting the fields in region II into the second and the last Maxwell equations in (2). The propagating modes in the mode expansion fulfill homogeneous Maxwell equation, and thus the electric fields and the magnetic fields of these propagating modes cancel. The second and the last Maxwell equations in (2), respectively, read 𝛿 𝑧 𝑧 0 × 𝐵 ̂ 𝑧 1 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐵 2 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 𝐵 1 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 𝐵 2 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 + ̂ 𝑧 × 𝐵 𝛿 𝑧 𝑧 0 𝑅 × 𝑅 𝛿 e x p 𝑖 𝜅 + 𝑖 𝜅 𝐵 e x p 𝑖 𝜅 𝑧 𝑧 0 + 𝑖 𝜔 𝜇 𝜀 𝐸 𝛿 𝑧 𝑧 0 𝑅 e x p 𝑖 𝜅 = 4 𝜋 𝑖 𝜔 𝜇 𝑃 𝛿 𝑧 𝑧 0 𝑅 × 𝑅 𝛿 e x p 𝑖 𝜅 + 4 𝜋 𝜇 𝑖 𝜅 𝑀 e x p 𝑖 𝜅 𝑧 𝑧 0 + 4 𝜋 𝜇 𝛿 𝑧 𝑧 0 𝑅 , 𝛿 ̂ 𝑧 × 𝑀 e x p 𝑖 𝜅 ( 1 3 ) 𝑧 𝑧 0 × 𝐸 ̂ 𝑧 1 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 𝐸 1 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 1 𝑧 0 + 𝐸 2 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 𝐸 2 , 𝑔 , e x p 𝑖 𝑘 𝑧 , 2 𝑧 0 + ̂ 𝑧 × 𝐸 𝛿 𝑧 𝑧 0 𝑅 × 𝑅 𝛿 e x p 𝑖 𝜅 + 𝑖 𝜅 𝐸 e x p 𝑖 𝜅 𝑧 𝑧 0 = 𝑖 𝜔 𝐵 𝛿 𝑧 𝑧 0 𝑅 , e x p 𝑖 𝜅 ( 1 4 ) where we have used Θ ( 𝑧 𝑧 0 ) = ̂ 𝑧 𝛿 ( 𝑧 𝑧 0 ) and 𝛿 ( 𝑧 𝑧 0 ) = ̂ 𝑧 𝛿 ( 𝑧 𝑧 0 ) . By requiring the coefficients in front of 𝛿 ( 𝑧 𝑧 0 ) and 𝛿 ( 𝑧 𝑧 0 ) equal for the two sides, the amplitudes 𝐵 𝑖 , 𝑔 , ± can be solved out. The fields generated in the regions outside of the medium can be determined by imposing the condition of continuity of tangential field at two boundaries ( 𝑧 = 0 and 𝑧 = ): 𝐸 I I = 𝐵 𝜔 , 𝜅 ; 𝑧 = 0 ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) 𝑠 , ̂ 𝑒 𝑠 , + 𝐵 𝑝 , ̂ 𝑒 𝑝 , 𝑅 , 𝐵 ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) e x p 𝑖 𝜅 I I = 𝐵 𝜔 , 𝜅 ; 𝑧 = 0 ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) 𝑠 , 𝑠 , + 𝐵 𝑝 , 𝑝 , 𝑅 , ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) e x p 𝑖 𝜅 ( 1 5 ) for 𝑧 = 0 boundary, where ̂ 𝑥 and ̂ 𝑦 are the unit vectors of 𝑥 and 𝑦 directions, 𝐸 I I = 𝐵 𝜔 , 𝜅 ; 𝑧 = ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) 𝑠 , + ̂ 𝑒 𝑠 , + + 𝐵 𝑝 , + ̂ 𝑒 𝑝 , + 𝑅 𝐵 ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) e x p 𝑖 𝜅 e x p ( 𝑖 𝑘 ) , I I = 𝐵 𝜔 , 𝜅 ; 𝑧 = ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) 𝑠 , + 𝑠 , + + 𝐵 𝑝 , + 𝑝 , + 𝑅 ( ̂ 𝑥 ̂ 𝑥 + ̂ 𝑦 ̂ 𝑦 ) e x p 𝑖 𝜅 e x p ( 𝑖 𝑘 ) , ( 1 6 ) for 𝑧 = boundary.

2.5. Correction of the Formulism: Isotropic Case

To prove the correction of the formulism, in this subsection we apply it to isotropic medium and compare the results gotten with well-known results in [17]. Equations from (13) to (16) are solved under the current case to get the Green functions. Since now the relative permittivity in region II then reduces to a scalar (we denote it by 𝜀 ), the eigenmodes in the region are also 𝑠 - and 𝑝 -polarized modes. The dispersions for 𝑠 - and 𝑝 -polarized modes are the same and given as 𝜔 2 𝜀 𝜇 = 𝜅 2 + 𝑘 2 𝑧 in region II and 𝜔 2 = 𝜅 2 + 𝑘 2 𝑧 in region I and III. The unit vectors of electric field and magnetic field are given in (12). The further convenience is that ̂ 𝑒 𝑠 ± , 𝜅 and ̂ 𝑧 build up a right-hand coordinate, which can simplify the solution of (13) and (14) by decomposing the fields according to the three directions. After equating the coefficients in front of 𝛿 ( 𝑧 𝑧 0 ) in (13) and (14), the following relations are found: ̂ 𝑧 × 𝐵 = 4 𝜋 𝜇 ̂ 𝑧 × 𝑀 , ̂ 𝑧 × 𝐸 = 0 . ( 1 7 ) The local magnetic source relates with local magnetic field as 𝐵 ̂ 𝑒 𝑠 ± = 4 𝜋 𝜇 𝑀 ̂ 𝑒 𝑠 ± and 𝐵 𝜅 = 4 𝜋 𝜇 𝑀 𝜅 and ̂ 𝑒 𝑠 , ± and 𝜅 components of local electric field vanish ( 𝐸 ̂ 𝑒 𝑠 ± = 𝐸 𝜅 = 0 ). From the equality of the coefficients in front of 𝛿 ( 𝑧 𝑧 0 ) , we can get the other equations, and we write them according to the components along ̂ 𝑒 𝑠 ± , 𝜅 , and ̂ 𝑧 directions for (13), and (14) respectively, 1 𝐵 𝜀 𝜇 𝑠 , 𝑔 , + + 𝐵 𝑠 , 𝑔 , = 4 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 𝑃 ̂ 𝑒 𝑠 , ± 4 𝜋 𝜇 𝑖 𝜅 𝜔 𝑘 𝑧 𝑀 , 𝐵 ̂ 𝑧 𝑝 , 𝑔 , + 𝐵 𝑝 , 𝑔 , , , 1 = 4 𝜋 𝑖 𝜔 𝜇 𝑃 𝜅 𝑖 𝜔 𝜀 𝜇 𝐸 ̂ 𝑧 = 𝑖 4 𝜋 𝜔 𝜇 𝑃 ̂ 𝑧 𝐵 𝜇 𝜀 𝑝 , 𝑔 , + + 𝐵 𝑝 , 𝑔 , + 𝜅 𝑘 I I 𝑘 𝑧 , I I = 𝑖 𝜔 𝑘 I I 𝑘 𝑧 , I I 𝐵 ̂ 𝑒 𝑠 , ± , 1 𝐵 𝜇 𝜀 𝑠 , 𝑔 , + 𝐵 𝑠 , 𝑔 , , 𝐵 = 𝑖 𝜔 𝐵 𝜅 ̂ 𝑧 = 0 . ( 1 8 ) The fields launching from the source can be gotten by solving (18) and the results are collected as follows: 𝐵 𝑠 , 𝑔 , + = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I 𝜇 𝜀 ̂ 𝑒 𝑠 , ± I I 𝑃 + 2 𝜋 𝜇 𝑖 𝜔 2 𝜇 𝜀 𝑘 𝑧 , I I 𝑠 , + I I 𝑀 , 𝐵 𝑠 , 𝑔 , = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I 𝜇 𝜀 ̂ 𝑒 𝑠 , ± I I 𝑃 + 2 𝜋 𝜇 𝑖 𝜔 2 𝜇 𝜀 𝑘 𝑧 , I I 𝑠 , I I 𝑀 , 𝐵 𝑝 , 𝑔 , + = 2 𝜋 𝑖 𝜇 𝜔 𝑘 I I 𝑘 𝑧 , I I 𝜀 𝜇 ̂ 𝑒 𝑠 , ± I I 𝑀 2 𝜋 𝑖 𝜔 𝑘 I I 𝜇 𝑘 𝑧 , I I 𝑠 , + I I 𝑃 , 𝐵 𝑝 , 𝑔 , = 2 𝜋 𝑖 𝜇 𝜔 𝑘 I I 𝑘 𝑧 , I I 𝜀 𝜇 ̂ 𝑒 𝑠 , ± I I 𝑀 2 𝜋 𝑖 𝜔 𝑘 I I 𝜇 𝑘 𝑧 , I I 𝑠 , I I 𝑃 , ( 1 9 ) where the results are written in a compact form with the help of the unit vectors defined and the same with (2.37) in [17], and thus Green functions in infinite bulk isotropic medium gotten by present method are correct.

To determine the fields in regions I and III, the boundary equations (15) and (16) should be solved. Isotropy of the medium in region II provides great simplification: since the orthogonality between tangential wave vectors ̂ 𝑒 𝑇 𝑠 , ± = 𝑇 𝑝 , ± and 𝑇 𝑠 , ± = ̂ 𝑒 𝑇 𝑝 , ± always holds, it is not necessary to arrange the equations in 𝑥 - 𝑦 - 𝑧 coordinate. Taking the advantage, the boundary equation (15) for 𝑧 = 0 becomes 𝐵 𝑠 , ̂ 𝑒 𝑠 , 𝑇 , I = 1 𝐵 𝜇 𝜀 𝑠 , 𝑏 , + ̂ 𝑒 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑔 , e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑠 , 𝑇 , I I + 𝐵 𝑠 , 𝑏 , ̂ 𝑒 𝑠 , 𝑇 , I I , 𝐵 𝑝 , ̂ 𝑒 𝑝 , 𝑇 , I = 1 𝐵 𝜇 𝜀 𝑝 , 𝑏 , + ̂ 𝑒 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑔 , e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑝 , 𝑇 , I I + 𝐵 𝑝 , 𝑏 , ̂ 𝑒 𝑝 , 𝑇 , I I , 𝐵 𝑝 , 𝑝 , 𝑇 , I = 𝐵 𝑝 , 𝑏 , + 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑔 , e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑝 , 𝑇 , I I + 𝐵 𝑝 , 𝑏 , 𝑝 , 𝑇 , I I , 𝐵 𝑠 , 𝑠 , 𝑇 , I = 𝐵 𝑠 , 𝑏 , + 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑔 , e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑠 , 𝑇 , I I + 𝐵 𝑠 , 𝑏 , 𝑠 , 𝑇 , I I . ( 2 0 ) Similarly, the boundary equation (16) for 𝑧 = gives 𝐵 𝑠 , + e x p 𝑖 𝑘 𝑧 , I ̂ 𝑒 𝑠 , + 𝑇 , I I I = 1 𝐵 𝜇 𝜀 𝑠 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑏 , + e x p 𝑖 𝑘 𝑧 , I I ̂ 𝑒 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑏 , e x p 𝑖 𝑘 𝑧 , I I ̂ 𝑒 𝑠 , 𝑇 , I I , 𝐵 𝑝 , + e x p 𝑖 𝑘 𝑧 , I ̂ 𝑒 𝑝 , + 𝑇 , I I I = 1 𝐵 𝜇 𝜀 𝑝 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑏 , + e x p 𝑖 𝑘 𝑧 , I I ̂ 𝑒 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑏 , e x p 𝑖 𝑘 𝑧 , I I ̂ 𝑒 𝑝 , 𝑇 , I I , 𝐵 𝑝 , + e x p 𝑖 𝑘 𝑧 , 𝐼 𝑝 , + 𝑇 , I I I = 𝐵 𝑝 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑏 , + e x p 𝑖 𝑘 𝑧 , I I 𝑝 , + 𝑇 , I I + 𝐵 𝑝 , 𝑏 , e x p 𝑖 𝑘 𝑧 , I I 𝑝 , 𝑇 , I I , 𝐵 𝑠 , + e x p 𝑖 𝑘 𝑧 , I 𝑠 , + 𝑇 , I I I = 𝐵 𝑠 , 𝑔 , + e x p 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑏 , + e x p 𝑖 𝑘 𝑧 , I I 𝑠 , + 𝑇 , I I + 𝐵 𝑠 , 𝑏 , e x p 𝑖 𝑘 𝑧 , I I 𝑠 , 𝑇 , I I , ( 2 1 ) where ̂ 𝑒 𝑠 , ± 𝑇 , I ( I I I ) , ̂ 𝑒 𝑇 𝑠 , ± , I I , ̂ 𝑒 𝑝 , ± 𝑇 , I ( I I I ) ,   ̂ 𝑒 𝑝 , ± 𝑇 , I I , 𝑝 , ± 𝑇 , I ( I I I ) , 𝑝 , ± 𝑇 , I I , 𝑠 , ± 𝑇 , I ( I I I ) , 𝑠 , ± 𝑇 , I I are the tangential parts of these unit vectors in different regions and 𝑘 𝑧 , I I = 𝑘 𝑧 , 𝑠 = 𝑘 𝑧 , 𝑝 . Once 𝐵 𝑠 , and 𝐵 𝑝 , ( 𝐵 𝑠 , + and 𝐵 𝑝 , + ) are solved from (20) and (21), the generated fields in region I (III) from the source plane are determined. The fields in region I are given as follows: 𝐵 I = 1 𝜔 , 𝜅 𝑟 𝜇 𝜀 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 1 𝑡 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + 𝑟 𝑠 I I , I 𝑡 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑡 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + × 𝑠 , 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 + 1 𝑟 𝜇 𝜀 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 1 𝑡 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + 𝑟 𝑝 I I , I 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑡 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + × 𝑝 , 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 , 𝐸 I = 1 𝜔 , 𝜅 𝑟 𝜇 𝜀 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 1 𝑡 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + 𝑟 𝑠 I I , I 𝑡 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑡 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + × ̂ 𝑒 𝑠 , 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 + 1 𝑟 𝜇 𝜀 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 1 𝑡 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + 𝑟 𝑝 I I , I 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑡 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + × ̂ 𝑒 𝑝 , 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I 𝑧 . ( 2 2 ) The fields generated in region III are gotten as 𝐵 I I I = 1 𝜔 , 𝜅 𝑟 𝜀 𝜇 𝑠 I I , I I I 𝑡 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 1 𝑡 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑡 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + × 𝑠 , + 𝑒 𝑖 𝑘 𝑧 , I 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 + 1 𝑟 𝜀 𝜇 𝑝 I I , I I I 𝑡 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 1 𝑡 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑡 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + × 𝑝 , + 𝑒 𝑖 𝑘 𝑧 , I 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 , 𝐸 I I I = 1 𝜔 , 𝜅 𝑟 𝜀 𝜇 𝑠 I I , I I I 𝑡 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 1 𝑡 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑡 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑠 , 𝑔 , + × ̂ 𝑒 𝑠 , + 𝑒 𝑖 𝑘 𝑧 , I 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 + 1 𝑟 𝜀 𝜇 𝑝 I I , I I I 𝑡 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 1 𝑡 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑡 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝐵 𝑝 , 𝑔 , + × ̂ 𝑒 𝑝 , + 𝑒 𝑖 𝑘 𝑧 , I 𝑅 e x p 𝑖 𝜅 e x p 𝑖 𝑘 𝑧 , I I I 𝑧 , ( 2 3 ) where 𝑡 𝑠 𝑖 , 𝑗 = ( 𝑘 𝑧 , 𝑖 𝑘 𝑧 , 𝑗 ) / ( 𝑘 𝑧 , 𝑖 + 𝑘 𝑧 , 𝑗 ) , 𝑡 𝑠 𝑖 , 𝑗 = 2 𝑘 𝑧 , 𝑖 / ( 𝑘 𝑧 , 𝑖 + 𝑘 𝑧 , 𝑗 ) , 𝑡 𝑝 𝑖 , 𝑗 = ( 𝑘 𝑧 , 𝑖 𝜀 𝑗 𝜇 𝑗 𝑘 𝑧 , 𝑗 𝜀 𝑖 𝜇 𝑖 ) / ( 𝑘 𝑧 , 𝑖 𝜀 𝑗 𝜇 𝑗 + 𝑘 𝑧 , 𝑗 𝜀 𝑖 𝜇 𝑖 )   , 𝑟 𝑝 𝑖 , 𝑗 = 2 𝑘 𝑧 , 𝑖 𝜀 𝑗 𝜇 𝑗 𝜀 𝑖 𝜇 𝑖 / ( 𝑘 𝑧 , 𝑖 𝜀 𝑗 𝜇 𝑗 + 𝑘 𝑧 , 𝑗 𝜀 𝑖 𝜇 𝑖 ) are Fresnel coefficients. Since 𝐵 𝑠 , 𝑔 , ± and 𝐵 𝑝 , 𝑔 , ± have been related with the sources, the fields in regions I and III in fact can spontaneously relate with the sources. Thus, Green functions linking the fields in these regions and source plane are ready to be gotten by following definition (8): 𝐺 I 𝐵 𝑀 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝜇 𝑖 𝜔 𝑘 I I 𝑘 𝑧 , I I 𝑡 𝑠 I I , I 𝑟 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 𝑠 , I 𝑠 , + I I + 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑠 , I 𝑠 , I I 1 𝑟 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I + 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 + 𝑟 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑟 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝑝 , I ̂ 𝑒 𝑠 , ± I I × 𝑒 𝑅 𝑖 𝜅 𝑒 𝑖 𝑘 𝑧 , I 𝑧 , 𝐺 I 𝐵 𝑃 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I × 𝑡 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 + 𝑟 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑟 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝑠 , I ̂ 𝑒 𝑠 , ± I I 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 𝑝 , I 𝑠 , I I + 𝑟 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 𝑝 , I 𝑠 , + I I 1 𝑟 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I × 𝑒 𝑅 𝑖 𝜅 𝑒 𝑖 𝑘 𝑧 , I 𝑧 , 𝐺 I 𝐸 𝑀 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝜇 𝑖 𝜔 𝑘 I I 𝑘 𝑧 , I I × 𝑟 𝑠 I I , I 𝑡 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) ̂ 𝑒 𝑠 , I 𝑠 , + I I + 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑠 , I 𝑠 , I I 1 𝑡 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I + 𝑟 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 + 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑡 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I ̂ 𝑒 𝑝 , I ̂ 𝑒 𝑠 , ± I I × 𝑒 𝑅 𝑖 𝜅 𝑒 𝑖 𝑘 𝑧 , I 𝑧 , 𝐺 I 𝐸 𝑀 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I 𝑡 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 + 𝑟 𝑠 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) 1 𝑟 𝑠 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I ̂ 𝑒 𝑠 , I ̂ 𝑒 𝑠 , ± I I 𝑡 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I 𝑧 0 ̂ 𝑒 𝑝 , I 𝑠 , I I + 𝑟 𝑝 I I , I 𝑒 𝑖 𝑘 𝑧 , I I ( 2 𝑧 0 ) ̂ 𝑒 𝑝 , I 𝑠 , + I I 1 𝑟 𝑝 I I , I 2 𝑒 𝑖 2 𝑘 𝑧 , I I × 𝑒 𝑅 𝑖 𝜅 𝑒 𝑖 𝑘 𝑧 , I 𝑧 , 𝐺 I I I 𝐵 𝑀 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝜇 𝑖 𝜔 𝑘 I I 𝑘 𝑧 , I I × 𝑡 𝑠 I I , I I I 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 𝑠 , + I I I 𝑠 , I I + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 𝑠 , + I I I 𝑠 , + I I 1 𝑟 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I + 𝑡 𝑝 I I , I I I 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑟 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝑝 , + I I I ̂ 𝑒 𝑠 , ± I I × 𝑒 𝑖 𝜅 𝑧 , I I I ( 𝑧 ) 𝑒 𝑖 𝑅 𝑘 , 𝐺 I I I 𝐵 𝑃 𝜅 ; 𝑧 𝑧 0 = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I × 𝑡 𝑠 I I , I I I 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑟 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I 𝑠 , + I I I ̂ 𝑒 𝑠 , ± I I 𝑡 𝑝 I I , I I I 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) 𝑝 , + I I I 𝑠 , I I + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 𝑝 , + I I I 𝑠 , + I I 1 𝑟 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I × 𝑒 𝑖 𝜅 𝑧 , I I I ( 𝑧 ) 𝑒 𝑖 𝑅 𝑘 𝐺 I I I 𝐸 𝑀 𝜅 𝑧 𝑧 0 = 2 𝜋 𝜇 𝑖 𝜔 𝑘 I I 𝑘 𝑧 , I I × 𝑡 𝑠 I I , I I I 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) ̂ 𝑒 𝑠 , + I I I 𝑠 , I I + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) ̂ 𝑒 𝑠 , + I I I 𝑠 , + I I 1 𝑟 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I + 𝑡 𝑝 I I , I I I 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑟 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I ̂ 𝑒 𝑝 , + I I I ̂ 𝑒 𝑠 , ± I I × 𝑒 𝑖 𝜅 𝑧 , I I I ( 𝑧 ) 𝑒 𝑖 𝑅 𝑘 , 𝐺 I I I 𝐸 𝑃 𝜅 𝑧 𝑧 0 = 2 𝜋 𝑖 𝜔 2 𝜇 𝑘 𝑧 , I I × 𝑡 𝑠 I I , I I I 𝑟 𝑠 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) 1 𝑟 𝑠 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I ̂ 𝑒 𝑠 , + I I I ̂ 𝑒 𝑠 , ± I I 𝑡 𝑝 I I , I I I 𝑟 𝑝 I I , I I I 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 + ) ̂ 𝑒 𝑝 , + I I I 𝑠 , I I + 𝑒 𝑖 𝑘 𝑧 , I I ( 𝑧 0 ) ̂ 𝑒 𝑝 , + I I I 𝑠 , + I I 1 𝑟 𝑝 I I , I I I 2 𝑒 𝑖 2 𝑘 𝑧 , I I × 𝑒 𝑖 𝜅 𝑧 , I I I ( 𝑧 ) 𝑒 𝑖 𝑅 𝑘 . ( 2 4 ) These results exactly recover (4.18) in [17], and thus the correction of current formulism is proved.

3. Green Functions for Sources in Anisotropic Medium

The formulism developed in the last section is then applied to obtain Green functions linking sources in anisotropic medium and fields in space. The advantage of mode expansion method would be obviously shown in the following study.

3.1. Sources in Uniaxial Crystal

Under the situation, the permittivity tensor in region II would be read as d i a g ( 𝜀 𝑜 , 𝜀 𝑒 , 𝜀 𝑜 ) , where without losing generality we have assumed the anisotropy is along 𝑦 direction (so that mode conversion still happens at the interface). Then the equation determining eigenmodes (5) becomes 𝜔 2 𝑘 2 𝑥 𝜀 𝑜 𝜇 𝑘 2 𝑦 𝜀 𝑜 𝜇 𝑘 2 𝑧 𝜀 𝑜 𝜇 𝜔 2 𝑘 2 𝑥 𝜀 𝑒 𝜇 𝑘 2 𝑦 𝜀 𝑜 𝜇 𝑘 2 𝑧 𝜀 𝑒 𝜇 = 0 . ( 2 5 ) Thus, there are two dispersions gotten: the mode with dispersion 𝜔 2 𝑘 2 𝑥 / 𝜀 𝑜 𝜇 𝑘 2 𝑦 / 𝜀 𝑜 𝜇 𝑘 2 𝑧 / 𝜀 𝑜 𝜇 = 0 is known as 𝑜 -light, and the mode with dispersion 𝜔 2 𝑘 2 𝑥 / 𝜀 𝑒 𝜇 𝑘 2 𝑦 / 𝜀 𝑜 𝜇 𝑘 2 𝑧 / 𝜀 𝑒 𝜇 = 0 is known as 𝑒 -light. The corresponding unit vector for electric field can be solved from (4) and the constraint (the homogeneous form of the first Maxwell equation in (2)) and given as ̂ 𝑒 𝑜 , ± = 1 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 ± 𝑘 𝑧 , 𝑜 0 , 𝜅 c o s 𝜃 ̂ 𝑒 𝑒 , ± = 1 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝜀 2 𝑒 + 𝜅 2 s i n 2 𝜃 𝜀 2 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 × 𝜀 𝑜 1 𝜅 2 𝜀 s i n 𝜃 c o s 𝜃 𝑒 1 𝑘 2 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 ± 𝜀 𝑜 1 𝑘 𝑧 , 𝑒 . 𝜅 s i n 𝜃 ( 2 6 ) Thus, we can have the unit vectors for magnetic field and electric displacement vector as 𝑜 ± = 1 𝑘 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝜅 2 𝑘 s i n 𝜃 c o s 𝜃 2 𝑧 , 𝑜 + 𝜅 2 c o s 2 𝜃 𝑘 𝑧 , 𝑜 , 𝜅 s i n 𝜃 𝑒 ± = 1 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 ± 𝑘 𝑧 , 𝑒 0 , 𝑑 𝜅 c o s 𝜃 𝑜 , ± = 1 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 ± 𝑘 𝑧 , 𝑜 0 , 𝑑 𝜅 c o s 𝜃 𝑒 ± = 1 𝑘 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝜅 2 s i n 𝜃 c o s 𝜃 𝑘 2 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 ± 𝑘 𝑧 , 𝑒 , 𝜅 s i n 𝜃 ( 2 7 ) where 𝑘 𝑜 = 𝜅 2 + 𝑘 2 𝑧 , 𝑜 and 𝑘 𝑒 = 𝜅 2 + 𝑘 2 𝑧 , 𝑒 . Following the mode expansion denoted in Section 2.3 the equations determining the fields launching from source plane can be gotten by applying (13) and (14) under current situation; the equality of coefficients in front of 𝛿 ( 𝑧 𝑧 0 ) and 𝛿 ( 𝑧 𝑧 0 ) provides the equations determining these fields and constraints about the local field ( 𝐵 𝑥 , 𝐵 𝑦 ) = 4 𝜋 𝜇 ( 𝑀 𝑥 , 𝑀 𝑦 ) and ( 𝐸 𝑥 , 𝐸 𝑦 ) = 0 . We summarize them from (13) and (14), respectively, for 𝑥 , 𝑦 , and 𝑧 directions: 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝑘 𝑜 𝐵 𝑜 , 𝑔 , 𝐵 𝑜 , 𝑔 , + = 𝑖 4 𝜋 𝜔 𝜇 𝑃 𝑥 + 𝑖 4 𝜋 𝜇 𝜅 s i n 𝜃 𝑀 𝑧 , 𝜅 2 s i n 𝜃 c o s 𝜃 𝑘 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝐵 𝑜 , 𝑔 , 𝐵 𝑜 , 𝑔 , + + 𝑘 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝐵 𝑒 , 𝑔 , + + 𝐵 𝑒 , 𝑔 , = 𝑖 4 𝜋 𝜔 𝜇 𝑃 𝑦 𝑖 4 𝜋 𝜇 𝜅 c o s 𝜃 𝑀 𝑧 , 𝐸 𝑧 = 4 𝜋 𝜀 𝑜 𝑃 𝑧 , 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝜔 𝜇 𝜀 𝑒 𝐵 𝑒 , 𝑔 , + 𝐵 𝑒 , 𝑔 , + 𝑖 𝐸 𝑧 𝜅 s i n 𝜃 𝑖 𝜔 𝐵 𝑥 𝑘 = 0 , 𝑜 𝜔 𝜇 𝜀 𝑜 𝑘 𝑧 , 𝑜 𝑘 2 𝑧 , 𝑜 + 𝜅 2 c o s 2 𝜃 𝐵 𝑜 , 𝑔 , + + 𝐵 𝑜 , 𝑔 , + 1 𝜔 𝜇 𝜀 𝑜 𝜅 2 s i n 𝜃 c o s 𝜃 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝐵 𝑒 , 𝑔 , + 𝐵 𝑒 , 𝑔 , 𝑖 𝐸 𝑧 𝜅 c o s 𝜃 𝑖 𝜔 𝐵 𝑦 𝐵 = 0 , 𝑧 = 0 . ( 2 8 ) After solving (28), the fields launching from the source plane are determined: 𝐵 𝑜 , 𝑔 , + = 𝑂 1 𝑃 𝑧 + 𝑂 2 𝑃 𝑥 + 𝑂 3 𝑀 𝑦 𝑂 4 𝑀 𝑥 𝑂 5 𝑀 𝑧 , 𝐵 𝑜 , 𝑔 , = 𝑂 1 𝑃 𝑧 𝑂 2 𝑃 𝑥 + 𝑂 3 𝑀 𝑦 𝑂 4 𝑀 𝑥 + 𝑂 5 𝑀 𝑧 , 𝐵 𝑒 , 𝑔 , + = 𝑂 1 𝑃 𝑦 + 𝑂 2 𝑃 𝑥 + 𝑂 3 𝑃 𝑧 + 𝑂 4 𝑀 𝑥 𝑂 5 𝑀 𝑧 , 𝐵 𝑒 , 𝑔 , = 𝑂 1 𝑃 𝑦 + 𝑂 2 𝑃 𝑥 𝑂 3 𝑃 𝑧 𝑂 4 𝑀 𝑥 𝑂 5 𝑀 𝑧 , ( 2 9 ) where the coefficients in front of the source are 𝑂 1 = 𝑖 2 𝜋 𝜅 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝜔 𝜇 𝑅 c o s 𝜃 𝑘 𝑜 𝑘 𝑧 , 𝑜 + 𝑖 2 𝜋 𝜅 3 𝜔 𝜇 𝜀 𝑒 𝑅 s i n 2 𝜃 c o s 𝜃 𝑘 𝑜 𝑘 𝑧 , 𝑜 𝜀 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 , 𝑂 2 = 𝑖 2 𝜋 𝜔 𝑘 𝑜 𝜇 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 , 𝑂 3 = 𝑖 2 𝜋 𝜔 2 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝜇 2 𝜀 𝑜 𝑅 𝑘 𝑜 𝑘 𝑧 , 𝑜 , 𝑂 4 = 𝑖 2 𝜋 𝜔 2 𝜅 2 𝜇 2 𝜀 𝑒 𝑅 s i n 𝜃 c o s 𝜃 𝑘 𝑜 𝑘 𝑧 , 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 , 𝑂 5 = 𝑖 2 𝜋 𝜇 𝜅 𝑘 𝑜 s i n 𝜃 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 , 𝑂 1 = 𝑖 2 𝜋 𝜔 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝜇 𝑘 𝑧 , 𝑒 𝑅 , 𝑂 2 = 𝑖 2 𝜋 𝜔 𝜅 2 𝜇 s i n 𝜃 c o s 𝜃 𝑘 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝑅 , 𝑂 3 = 𝑖 2 𝜋 𝜔 𝜅 𝜇 𝜀 𝑒 s i n 𝜃 𝜀 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 , 𝑂 4 = 𝑖 2 𝜋 𝜔 2 𝜇 2 𝜀 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 , 𝑂 5 = 𝑖 2 𝜋 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝜅 𝜇 c o s 𝜃 𝑘 𝑧 , 𝑒 𝑅 + 𝑖 2 𝜋 𝜅 3 𝜇 s i n 2 𝜃 c o s 𝜃 𝑘 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝑅 , ( 3 0 ) where the coefficient 𝑅 in (30) is given as 𝑅 = 𝑘 2 𝑧 , 𝑜 + 𝜅 2 c o s 2 𝜃 / 𝑘 2 𝑧 , 𝑒 + 𝜅 2 c o s 2 𝜃 . The fields in the regions outside of the anisotropic medium can be determined by solving boundary equations, and they can be written for 𝑧 = 0 and 𝑧 = , respectively, as 𝜅 2 s i n 𝜃 c o s 𝜃 𝑘 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝐵 𝑜 , 𝑏 , + + 𝐵 𝑜 , 𝑏 , + 𝐵 𝑜 , 𝑔 , 𝑒 𝑖 𝑘 𝑧 , 𝑜 𝑧 0 + 𝑘 𝑧 , 𝑒 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 𝐵 𝑒 , 𝑏 , + 𝐵 𝑒 , 𝑏 , 𝐵 𝑒 , 𝑔 , 𝑒 𝑖 𝑘 𝑧 , 𝑒 𝑧 0 = 𝐵 𝑠 , 𝑘 𝑧 , I 𝜔 c o s 𝜃 + 𝐵 𝑝 , s i n 𝜃 , 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 𝑘 𝑜 𝐵 𝑜 , 𝑏 , + + 𝐵 𝑜 , 𝑏 , + 𝐵 𝑜 , 𝑔 , 𝑒 𝑖 𝑘 𝑧 , 𝑜 𝑧 0 = 𝐵 𝑠 , 𝑘 𝑧 , I 𝜔 s i n 𝜃 𝐵 𝑝 , 𝐵 c o s 𝜃 , 𝑜 , 𝑏 , + 𝐵 𝑜 , 𝑏 , 𝐵 𝑜 , 𝑔 , 𝑒 𝑖 𝑘 𝑧 , 𝑜 𝑧 0 𝑘 𝑜 𝑘 𝑧 , 𝑜 𝜔 𝜇 𝜀 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑜 + 𝐵 𝑒 , 𝑏 , + + 𝐵 𝑒 , 𝑏 , + 𝐵 𝑒 , 𝑔 , 𝑒 𝑖 𝑘 𝑧 , 𝑒 𝑧 0 𝜅 2 s i n 𝜃 c o s 𝜃 𝜔 𝜇 𝜀 𝑜 𝜅 2 c o s 2 𝜃 + 𝑘 2 𝑧 , 𝑒 = 𝐵 𝑠 ,