Table of Contents Author Guidelines Submit a Manuscript
International Journal of Optics
Volume 2012 (2012), Article ID 532316, 19 pages
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.


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: 𝜀𝜀=1000𝜀2000𝜀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.

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𝐸𝑅,𝐵=𝜅,𝑧exp𝑖𝜅𝑅,𝑧𝑑𝜅(2𝜋)2𝐵𝑅,𝜅,𝑧exp𝑖𝜅(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 𝐵II=𝐵𝜔,𝜅1,𝑔,+exp𝑖𝑘𝑧,1𝑧0+𝐵2,𝑔,+exp𝑖𝑘𝑧,2𝑧0×Θ𝑧𝑧0+𝐵1,𝑏,++𝐵2,𝑏,++𝐵1,𝑔,exp𝑖𝑘𝑧,1𝑧0+𝐵2,𝑔,exp𝑖𝑘𝑧,2𝑧0𝑧×Θ0+𝐵𝑧1,𝑏,+𝐵2,𝑏,+𝐵𝛿𝑧𝑧0𝑅,𝐷exp𝑖𝜅II=𝐷𝜔,𝜅1,𝑔,+exp𝑖𝑘𝑧,1𝑧0+𝐷2,𝑔,+exp𝑖𝑘𝑧,2𝑧0×Θ𝑧𝑧0+𝐷1,𝑏,++𝐷2,𝑏,++𝐷1,𝑔,exp𝑖𝑘𝑧,1𝑧0+𝐷2,𝑔,exp𝑖𝑘𝑧,2𝑧0𝑧×Θ0+𝐷𝑧1,𝑏,+𝐷2,𝑏,+𝐷𝛿𝑧𝑧0𝑅,𝐸exp𝑖𝜅II=𝐸𝜔,𝜅1,𝑔,+exp𝑖𝑘𝑧,1𝑧0+𝐸2,𝑔,+exp𝑖𝑘𝑧,2𝑧0×Θ𝑧𝑧0+𝐸1,𝑏,++𝐸2,𝑏,++𝐸1,𝑔,exp𝑖𝑘𝑧,1𝑧0+𝐸2,𝑔,exp𝑖𝑘𝑧,2𝑧0𝑧×Θ0+𝐸𝑧1,𝑏,+𝐸2,𝑏,+𝐸𝛿𝑧𝑧0𝑅,exp𝑖𝜅(9) where 𝐵𝑖,𝑔,±=𝐵𝑖,𝑔,±𝑅exp𝑖𝜅exp±𝑖𝑘𝑧,𝑖𝑧𝑖,±,𝐵𝑖,𝑏,±=𝐵𝑖,𝑏,±𝑅exp𝑖𝜅exp±𝑖𝑘𝑧,𝑖𝑧𝑖,±,𝐷𝑖,𝑔,±=1𝐵𝑖𝜔𝜇×𝑖,𝑔,±,𝐸𝑖,𝑔,±=𝜀1𝐷𝑖,𝑔,±𝐷𝑖,𝑏,±=1𝐵𝑖𝜔𝜇×𝑖,𝑏,±,𝐸𝑖,𝑏,±=𝜀1𝐷𝑖,𝑏,±.𝑖=1,2,(10)𝐵(𝐷, 𝐸) 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𝑅×exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧,𝐷I=𝐸𝜔,𝜅I=𝐵𝜔,𝜅𝑠,̂𝑒𝑠,I+𝐵𝑝,̂𝑒𝑝,I𝑅×exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧,𝐵III=𝐵𝜔,𝜅𝑠,+𝑠,+III+𝐵𝑝,+𝑝,+III𝑅×exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧,𝐷III=𝐸𝜔,𝜅III=𝐵𝜔,𝜅𝑠,+̂𝑒𝑠,+III+𝐵𝑝,+̂𝑒𝑝,+III𝑅×exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧,(11) where the unit vectors are defined as ̂𝑒𝛼𝑠,±=𝑘𝑠,±×̂𝑧||𝑘𝑠,±||=×̂𝑧𝜅cos𝜃𝜅sin𝜃±𝑘𝑧×001=0,sin𝜃cos𝜃𝛼𝑠,±=𝑘𝑠,±×̂𝑒𝑠,±||𝑘𝑠,±×̂𝑒𝑠,±||=1𝜔𝜀𝑟𝜇𝑟±𝑘𝑧cos𝜃±𝑘𝑧,sin𝜃𝜅𝛼𝑝,±=𝑘𝑝,±×̂𝑧||𝑘𝑝,±||=×̂𝑧𝜅cos𝜃𝜅sin𝜃±𝑘𝑧×001=0,sin𝜃cos𝜃̂𝑒𝛼𝑝,±=𝑘𝑝,±×𝑝,±|||𝑘𝑝,±×𝑝,±|||=1𝜔𝜀𝑟𝜇𝑟𝑘𝑧cos𝜃𝑘𝑧𝜅,sin𝜃(12) where 𝛼=I,II,III denotes the region where the unit vectors belong.

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,𝑔,+exp𝑖𝑘𝑧,1𝑧0+𝐵2,𝑔,+exp𝑖𝑘𝑧,2𝑧0𝐵1,𝑔,exp𝑖𝑘𝑧,1𝑧0𝐵2,𝑔,exp𝑖𝑘𝑧,2𝑧0+̂𝑧×𝐵𝛿𝑧𝑧0𝑅×𝑅𝛿exp𝑖𝜅+𝑖𝜅𝐵exp𝑖𝜅𝑧𝑧0+𝑖𝜔𝜇𝜀𝐸𝛿𝑧𝑧0𝑅exp𝑖𝜅=4𝜋𝑖𝜔𝜇𝑃𝛿𝑧𝑧0𝑅×𝑅𝛿exp𝑖𝜅+4𝜋𝜇𝑖𝜅𝑀exp𝑖𝜅𝑧𝑧0+4𝜋𝜇𝛿𝑧𝑧0𝑅,𝛿̂𝑧×𝑀exp𝑖𝜅(13)𝑧𝑧0×𝐸̂𝑧1,𝑔,+exp𝑖𝑘𝑧,1𝑧0𝐸1,𝑔,exp𝑖𝑘𝑧,1𝑧0+𝐸2,𝑔,+exp𝑖𝑘𝑧,2𝑧0𝐸2,𝑔,exp𝑖𝑘𝑧,2𝑧0+̂𝑧×𝐸𝛿𝑧𝑧0𝑅×𝑅𝛿exp𝑖𝜅+𝑖𝜅𝐸exp𝑖𝜅𝑧𝑧0=𝑖𝜔𝐵𝛿𝑧𝑧0𝑅,exp𝑖𝜅(14) 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 𝑧=): 𝐸II=𝐵𝜔,𝜅;𝑧=0(̂𝑥̂𝑥+̂𝑦̂𝑦)𝑠,̂𝑒𝑠,+𝐵𝑝,̂𝑒𝑝,𝑅,𝐵(̂𝑥̂𝑥+̂𝑦̂𝑦)exp𝑖𝜅II=𝐵𝜔,𝜅;𝑧=0(̂𝑥̂𝑥+̂𝑦̂𝑦)𝑠,𝑠,+𝐵𝑝,𝑝,𝑅,(̂𝑥̂𝑥+̂𝑦̂𝑦)exp𝑖𝜅(15) for 𝑧=0 boundary, where ̂𝑥 and ̂𝑦 are the unit vectors of 𝑥 and 𝑦 directions, 𝐸II=𝐵𝜔,𝜅;𝑧=(̂𝑥̂𝑥+̂𝑦̂𝑦)𝑠,+̂𝑒𝑠,++𝐵𝑝,+̂𝑒𝑝,+𝑅𝐵(̂𝑥̂𝑥+̂𝑦̂𝑦)exp𝑖𝜅exp(𝑖𝑘),II=𝐵𝜔,𝜅;𝑧=(̂𝑥̂𝑥+̂𝑦̂𝑦)𝑠,+𝑠,++𝐵𝑝,+𝑝,+𝑅(̂𝑥̂𝑥+̂𝑦̂𝑦)exp𝑖𝜅exp(𝑖𝑘),(16) 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.(17) 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𝜋𝜔𝜇𝑃̂𝑧𝐵𝜇𝜀𝑝,𝑔,++𝐵𝑝,𝑔,+𝜅𝑘II𝑘𝑧,II=𝑖𝜔𝑘II𝑘𝑧,II𝐵̂𝑒𝑠,±,1𝐵𝜇𝜀𝑠,𝑔,+𝐵𝑠,𝑔,,𝐵=𝑖𝜔𝐵𝜅̂𝑧=0.(18) The fields launching from the source can be gotten by solving (18) and the results are collected as follows: 𝐵𝑠,𝑔,+=2𝜋𝑖𝜔2𝜇𝑘𝑧,II𝜇𝜀̂𝑒𝑠,±II𝑃+2𝜋𝜇𝑖𝜔2𝜇𝜀𝑘𝑧,II𝑠,+II𝑀,𝐵𝑠,𝑔,=2𝜋𝑖𝜔2𝜇𝑘𝑧,II𝜇𝜀̂𝑒𝑠,±II𝑃+2𝜋𝜇𝑖𝜔2𝜇𝜀𝑘𝑧,II𝑠,II𝑀,𝐵𝑝,𝑔,+=2𝜋𝑖𝜇𝜔𝑘II𝑘𝑧,II𝜀𝜇̂𝑒𝑠,±II𝑀2𝜋𝑖𝜔𝑘II𝜇𝑘𝑧,II𝑠,+II𝑃,𝐵𝑝,𝑔,=2𝜋𝑖𝜇𝜔𝑘II𝑘𝑧,II𝜀𝜇̂𝑒𝑠,±II𝑀2𝜋𝑖𝜔𝑘II𝜇𝑘𝑧,II𝑠,II𝑃,(19) 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𝐵𝜇𝜀𝑠,𝑏,+̂𝑒𝑠,+𝑇,II+𝐵𝑠,𝑔,exp𝑖𝑘𝑧,II𝑧0̂𝑒𝑠,𝑇,II+𝐵𝑠,𝑏,̂𝑒𝑠,𝑇,II,𝐵𝑝,̂𝑒𝑝,𝑇,I=1𝐵𝜇𝜀𝑝,𝑏,+̂𝑒𝑝,+𝑇,II+𝐵𝑝,𝑔,exp𝑖𝑘𝑧,II𝑧0̂𝑒𝑝,𝑇,II+𝐵𝑝,𝑏,̂𝑒𝑝,𝑇,II,𝐵𝑝,𝑝,𝑇,I=𝐵𝑝,𝑏,+𝑝,+𝑇,II+𝐵𝑝,𝑔,exp𝑖𝑘𝑧,II𝑧0𝑝,𝑇,II+𝐵𝑝,𝑏,𝑝,𝑇,II,𝐵𝑠,𝑠,𝑇,I=𝐵𝑠,𝑏,+𝑠,+𝑇,II+𝐵𝑠,𝑔,exp𝑖𝑘𝑧,II𝑧0𝑠,𝑇,II+𝐵𝑠,𝑏,𝑠,𝑇,II.(20) Similarly, the boundary equation (16) for 𝑧= gives 𝐵𝑠,+exp𝑖𝑘𝑧,Î𝑒𝑠,+𝑇,III=1𝐵𝜇𝜀𝑠,𝑔,+exp𝑖𝑘𝑧,II𝑧0̂𝑒𝑠,+𝑇,II+𝐵𝑠,𝑏,+exp𝑖𝑘𝑧,IÎ𝑒𝑠,+𝑇,II+𝐵𝑠,𝑏,exp𝑖𝑘𝑧,IÎ𝑒𝑠,𝑇,II,𝐵𝑝,+exp𝑖𝑘𝑧,Î𝑒𝑝,+𝑇,III=1𝐵𝜇𝜀𝑝,𝑔,+exp𝑖𝑘𝑧,II𝑧0̂𝑒𝑝,+𝑇,II+𝐵𝑝,𝑏,+exp𝑖𝑘𝑧,IÎ𝑒𝑝,+𝑇,II+𝐵𝑝,𝑏,exp𝑖𝑘𝑧,IÎ𝑒𝑝,𝑇,II,𝐵𝑝,+exp𝑖𝑘𝑧,𝐼𝑝,+𝑇,III=𝐵𝑝,𝑔,+exp𝑖𝑘𝑧,II𝑧0𝑝,+𝑇,II+𝐵𝑝,𝑏,+exp𝑖𝑘𝑧,II𝑝,+𝑇,II+𝐵𝑝,𝑏,exp𝑖𝑘𝑧,II𝑝,𝑇,II,𝐵𝑠,+exp𝑖𝑘𝑧,I𝑠,+𝑇,III=𝐵𝑠,𝑔,+exp𝑖𝑘𝑧,II𝑧0𝑠,+𝑇,II+𝐵𝑠,𝑏,+exp𝑖𝑘𝑧,II𝑠,+𝑇,II+𝐵𝑠,𝑏,exp𝑖𝑘𝑧,II𝑠,𝑇,II,(21) where ̂𝑒𝑠,±𝑇,I(III),̂𝑒𝑇𝑠,±,II,̂𝑒𝑝,±𝑇,I(III),  ̂𝑒𝑝,±𝑇,II, 𝑝,±𝑇,I(III),𝑝,±𝑇,II,𝑠,±𝑇,I(III),𝑠,±𝑇,II are the tangential parts of these unit vectors in different regions and 𝑘𝑧,II=𝑘𝑧,𝑠=𝑘𝑧,𝑝. 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𝜔,𝜅𝑟𝜇𝜀𝑠II,I𝑒𝑖𝑘𝑧,II𝑧01𝑡𝑠II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+𝑟𝑠II,I𝑡𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑡𝑠II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+×𝑠,𝑅exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧+1𝑟𝜇𝜀𝑝II,I𝑒𝑖𝑘𝑧,II𝑧01𝑡𝑝II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+𝑟𝑝II,I𝑡𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑡𝑝II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+×𝑝,𝑅exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧,𝐸I=1𝜔,𝜅𝑟𝜇𝜀𝑠II,I𝑒𝑖𝑘𝑧,II𝑧01𝑡𝑠II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+𝑟𝑠II,I𝑡𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑡𝑠II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+×̂𝑒𝑠,𝑅exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧+1𝑟𝜇𝜀𝑝II,I𝑒𝑖𝑘𝑧,II𝑧01𝑡𝑝II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+𝑟𝑝II,I𝑡𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑡𝑝II,I2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+×̂𝑒𝑝,𝑅exp𝑖𝜅exp𝑖𝑘𝑧,I𝑧.(22) The fields generated in region III are gotten as 𝐵III=1𝜔,𝜅𝑟𝜀𝜇𝑠II,III𝑡𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)1𝑡𝑠II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0)1𝑡𝑠II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+×𝑠,+𝑒𝑖𝑘𝑧,I𝑅exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧+1𝑟𝜀𝜇𝑝II,III𝑡𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)1𝑡𝑝II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0)1𝑡𝑝II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+×𝑝,+𝑒𝑖𝑘𝑧,I𝑅exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧,𝐸III=1𝜔,𝜅𝑟𝜀𝜇𝑠II,III𝑡𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)1𝑡𝑠II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0)1𝑡𝑠II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑠,𝑔,+×̂𝑒𝑠,+𝑒𝑖𝑘𝑧,I𝑅exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧+1𝑟𝜀𝜇𝑝II,III𝑡𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)1𝑡𝑝II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0)1𝑡𝑝II,III2𝑒𝑖2𝑘𝑧,II𝐵𝑝,𝑔,+×̂𝑒𝑝,+𝑒𝑖𝑘𝑧,I𝑅exp𝑖𝜅exp𝑖𝑘𝑧,III𝑧,(23) 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𝜋𝜇𝑖𝜔𝑘II𝑘𝑧,II𝑡𝑠II,I𝑟𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)𝑠,I𝑠,+II+𝑒𝑖𝑘𝑧,II𝑧0𝑠,I𝑠,II1𝑟𝑠II,I2𝑒𝑖2𝑘𝑧,II+𝑡𝑝II,I𝑒𝑖𝑘𝑧,II𝑧0+𝑟𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑟𝑝II,I2𝑒𝑖2𝑘𝑧,II𝑝,Î𝑒𝑠,±II×𝑒𝑅𝑖𝜅𝑒𝑖𝑘𝑧,I𝑧,𝐺I𝐵𝑃𝜅;𝑧𝑧0=2𝜋𝑖𝜔2𝜇𝑘𝑧,II×𝑡𝑠II,I𝑒𝑖𝑘𝑧,II𝑧0+𝑟𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑟𝑠II,I2𝑒𝑖2𝑘𝑧,II𝑠,Î𝑒𝑠,±II𝑡𝑝II,I𝑒𝑖𝑘𝑧,II𝑧0𝑝,I𝑠,II+𝑟𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)𝑝,I𝑠,+II1𝑟𝑝II,I2𝑒𝑖2𝑘𝑧,II×𝑒𝑅𝑖𝜅𝑒𝑖𝑘𝑧,I𝑧,𝐺I𝐸𝑀𝜅;𝑧𝑧0=2𝜋𝜇𝑖𝜔𝑘II𝑘𝑧,II×𝑟𝑠II,I𝑡𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)̂𝑒𝑠,I𝑠,+II+𝑒𝑖𝑘𝑧,II𝑧0̂𝑒𝑠,I𝑠,II1𝑡𝑠II,I2𝑒𝑖2𝑘𝑧,II+𝑟𝑝II,I𝑒𝑖𝑘𝑧,II𝑧0+𝑡𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑡𝑝II,I2𝑒𝑖2𝑘𝑧,IÎ𝑒𝑝,Î𝑒𝑠,±II×𝑒𝑅𝑖𝜅𝑒𝑖𝑘𝑧,I𝑧,𝐺I𝐸𝑀𝜅;𝑧𝑧0=2𝜋𝑖𝜔2𝜇𝑘𝑧,II𝑡𝑠II,I𝑒𝑖𝑘𝑧,II𝑧0+𝑟𝑠II,I𝑒𝑖𝑘𝑧,II(2𝑧0)1𝑟𝑠II,I2𝑒𝑖2𝑘𝑧,IÎ𝑒𝑠,Î𝑒𝑠,±II𝑡𝑝II,I𝑒𝑖𝑘𝑧,II𝑧0̂𝑒𝑝,I𝑠,II+𝑟𝑝II,I𝑒𝑖𝑘𝑧,II(2𝑧0)̂𝑒𝑝,I𝑠,+II1𝑟𝑝II,I2𝑒𝑖2𝑘𝑧,II×𝑒𝑅𝑖𝜅𝑒𝑖𝑘𝑧,I𝑧,𝐺III𝐵𝑀𝜅;𝑧𝑧0=2𝜋𝜇𝑖𝜔𝑘II𝑘𝑧,II×𝑡𝑠II,III𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)𝑠,+III𝑠,II+𝑒𝑖𝑘𝑧,II(𝑧0)𝑠,+III𝑠,+II1𝑟𝑠II,III2𝑒𝑖2𝑘𝑧,II+𝑡𝑝II,III𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)+𝑒𝑖𝑘𝑧,II(𝑧0)1𝑟𝑝II,III2𝑒𝑖2𝑘𝑧,II𝑝,+IIÎ𝑒𝑠,±II×𝑒𝑖𝜅𝑧,III(𝑧)𝑒𝑖𝑅𝑘,𝐺III𝐵𝑃𝜅;𝑧𝑧0=2𝜋𝑖𝜔2𝜇𝑘𝑧,II×𝑡𝑠II,III𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)+𝑒𝑖𝑘𝑧,II(𝑧0)1𝑟𝑠II,III2𝑒𝑖2𝑘𝑧,II𝑠,+IIÎ𝑒𝑠,±II𝑡𝑝II,III𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)𝑝,+III𝑠,II+𝑒𝑖𝑘𝑧,II(𝑧0)𝑝,+III𝑠,+II1𝑟𝑝II,III2𝑒𝑖2𝑘𝑧,II×𝑒𝑖𝜅𝑧,III(𝑧)𝑒𝑖𝑅𝑘𝐺III𝐸𝑀𝜅𝑧𝑧0=2𝜋𝜇𝑖𝜔𝑘II𝑘𝑧,II×𝑡𝑠II,III𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)̂𝑒𝑠,+III𝑠,II+𝑒𝑖𝑘𝑧,II(𝑧0)̂𝑒𝑠,+III𝑠,+II1𝑟𝑠II,III2𝑒𝑖2𝑘𝑧,II+𝑡𝑝II,III𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)+𝑒𝑖𝑘𝑧,II(𝑧0)1𝑟𝑝II,III2𝑒𝑖2𝑘𝑧,IÎ𝑒𝑝,+IIÎ𝑒𝑠,±II×𝑒𝑖𝜅𝑧,III(𝑧)𝑒𝑖𝑅𝑘,𝐺III𝐸𝑃𝜅𝑧𝑧0=2𝜋𝑖𝜔2𝜇𝑘𝑧,II×𝑡𝑠II,III𝑟𝑠II,III𝑒𝑖𝑘𝑧,II(𝑧0+)+𝑒𝑖𝑘𝑧,II(𝑧0)1𝑟𝑠II,III2𝑒𝑖2𝑘𝑧,IÎ𝑒𝑠,+IIÎ𝑒𝑠,±II𝑡𝑝II,III𝑟𝑝II,III𝑒𝑖𝑘𝑧,II(𝑧0+)̂𝑒𝑝,+III𝑠,II+𝑒𝑖𝑘𝑧,II(𝑧0)̂𝑒𝑝,+III𝑠,+II1𝑟𝑝II,III2𝑒𝑖2𝑘𝑧,II×𝑒𝑖𝜅𝑧,III(𝑧)𝑒𝑖𝑅𝑘.(24) 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 diag(𝜀𝑜,𝜀𝑒,𝜀𝑜), 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.(25) 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𝜅2cos2𝜃+𝑘2𝑧,𝑜±𝑘𝑧,𝑜0,𝜅cos𝜃̂𝑒𝑒,±=1𝜅2cos2𝜃+𝑘2𝑧,𝑒𝜀2𝑒+𝜅2sin2𝜃𝜀2𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑒×𝜀𝑜1𝜅2𝜀sin𝜃cos𝜃𝑒1𝑘2𝑧,𝑒𝜅2cos2𝜃±𝜀𝑜1𝑘𝑧,𝑒.𝜅sin𝜃(26) Thus, we can have the unit vectors for magnetic field and electric displacement vector as 𝑜±=1𝑘𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑜𝜅2𝑘sin𝜃cos𝜃2𝑧,𝑜+𝜅2cos2𝜃𝑘𝑧,𝑜,𝜅sin𝜃𝑒±=1𝜅2cos2𝜃+𝑘2𝑧,𝑒±𝑘𝑧,𝑒0,𝑑𝜅cos𝜃𝑜,±=1𝜅2cos2𝜃+𝑘2𝑧,𝑜±𝑘𝑧,𝑜0,𝑑𝜅cos𝜃𝑒±=1𝑘𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒𝜅2sin𝜃cos𝜃𝑘2𝑧,𝑒𝜅2cos2𝜃±𝑘𝑧,𝑒,𝜅sin𝜃(27) 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: 𝜅2cos2𝜃+𝑘2𝑧,𝑜𝑘𝑜𝐵𝑜,𝑔,𝐵𝑜,𝑔,+=𝑖4𝜋𝜔𝜇𝑃𝑥+𝑖4𝜋𝜇𝜅sin𝜃𝑀𝑧,𝜅2sin𝜃cos𝜃𝑘𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑜𝐵𝑜,𝑔,𝐵𝑜,𝑔,++𝑘𝑧,𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒𝐵𝑒,𝑔,++𝐵𝑒,𝑔,=𝑖4𝜋𝜔𝜇𝑃𝑦𝑖4𝜋𝜇𝜅cos𝜃𝑀𝑧,𝐸𝑧=4𝜋𝜀𝑜𝑃𝑧,𝜅2cos2𝜃+𝑘2𝑧,𝑒𝜔𝜇𝜀𝑒𝐵𝑒,𝑔,+𝐵𝑒,𝑔,+𝑖𝐸𝑧𝜅sin𝜃𝑖𝜔𝐵𝑥𝑘=0,𝑜𝜔𝜇𝜀𝑜𝑘𝑧,𝑜𝑘2𝑧,𝑜+𝜅2cos2𝜃𝐵𝑜,𝑔,++𝐵𝑜,𝑔,+1𝜔𝜇𝜀𝑜𝜅2sin𝜃cos𝜃𝜅2cos2𝜃+𝑘2𝑧,𝑒𝐵𝑒,𝑔,+𝐵𝑒,𝑔,𝑖𝐸𝑧𝜅cos𝜃𝑖𝜔𝐵𝑦𝐵=0,𝑧=0.(28) 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𝑀𝑧,(29) where the coefficients in front of the source are 𝑂1=𝑖2𝜋𝜅𝜅2cos2𝜃+𝑘2𝑧,𝑒𝜔𝜇𝑅cos𝜃𝑘𝑜𝑘𝑧,𝑜+𝑖2𝜋𝜅3𝜔𝜇𝜀𝑒𝑅sin2𝜃cos𝜃𝑘𝑜𝑘𝑧,𝑜𝜀𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑒,𝑂2=𝑖2𝜋𝜔𝑘𝑜𝜇𝜅2cos2𝜃+𝑘2𝑧,𝑜,𝑂3=𝑖2𝜋𝜔2𝜅2cos2𝜃+𝑘2𝑧,𝑒𝜇2𝜀𝑜𝑅𝑘𝑜𝑘𝑧,𝑜,𝑂4=𝑖2𝜋𝜔2𝜅2𝜇2𝜀𝑒𝑅sin𝜃cos𝜃𝑘𝑜𝑘𝑧,𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑒,𝑂5=𝑖2𝜋𝜇𝜅𝑘𝑜sin𝜃𝜅2cos2𝜃+𝑘2𝑧,𝑜,𝑂1=𝑖2𝜋𝜔𝜅2cos2𝜃+𝑘2𝑧,𝑜𝜇𝑘𝑧,𝑒𝑅,𝑂2=𝑖2𝜋𝜔𝜅2𝜇sin𝜃cos𝜃𝑘𝑧,𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑜𝑅,𝑂3=𝑖2𝜋𝜔𝜅𝜇𝜀𝑒sin𝜃𝜀𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑒,𝑂4=𝑖2𝜋𝜔2𝜇2𝜀𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒,𝑂5=𝑖2𝜋𝜅2cos2𝜃+𝑘2𝑧,𝑜𝜅𝜇cos𝜃𝑘𝑧,𝑒𝑅+𝑖2𝜋𝜅3𝜇sin2𝜃cos𝜃𝑘𝑧,𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑜𝑅,(30) where the coefficient 𝑅 in (30) is given as 𝑅=𝑘2𝑧,𝑜+𝜅2cos2𝜃/𝑘2𝑧,𝑒+𝜅2cos2𝜃. 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 𝜅2sin𝜃cos𝜃𝑘𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑜𝐵𝑜,𝑏,++𝐵𝑜,𝑏,+𝐵𝑜,𝑔,𝑒𝑖𝑘𝑧,𝑜𝑧0+𝑘𝑧,𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒𝐵𝑒,𝑏,+𝐵𝑒,𝑏,𝐵𝑒,𝑔,𝑒𝑖𝑘𝑧,𝑒𝑧0=𝐵𝑠,𝑘𝑧,I𝜔cos𝜃+𝐵𝑝,sin𝜃,𝜅2cos2𝜃+𝑘2𝑧,𝑜𝑘𝑜𝐵𝑜,𝑏,++𝐵𝑜,𝑏,+𝐵𝑜,𝑔,𝑒𝑖𝑘𝑧,𝑜𝑧0=𝐵𝑠,𝑘𝑧,I𝜔sin𝜃𝐵𝑝,𝐵cos𝜃,𝑜,𝑏,+𝐵𝑜,𝑏,𝐵𝑜,𝑔,𝑒𝑖𝑘𝑧,𝑜𝑧0𝑘𝑜𝑘𝑧,𝑜𝜔𝜇𝜀𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑜+𝐵𝑒,𝑏,++𝐵𝑒,𝑏,+𝐵𝑒,𝑔,𝑒𝑖𝑘𝑧,𝑒𝑧0𝜅2sin𝜃cos𝜃𝜔𝜇𝜀𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑒=𝐵𝑠,sin𝜃+𝐵𝑝,𝑘𝑧,I𝐵𝜔cos𝜃,𝑒,𝑏,++𝐵𝑒,𝑏,+𝐵𝑒,𝑔,𝑒𝑖𝑘𝑧,𝑒𝑧0𝑘2𝑧,𝑒𝜅2cos2𝜃𝜔𝜇𝜀𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒=𝐵𝑠,cos𝜃+𝐵𝑝,𝑘𝑧,I𝐵𝜔sin𝜃,𝑜,𝑏,+𝑒𝑖𝑘𝑧,𝑜+𝐵𝑜,𝑔,+𝑒𝑖𝑘𝑧,𝑜(𝑧0)+𝐵𝑜,𝑏,𝑒𝑖𝑘𝑧,𝑜×𝜅2sin𝜃cos𝜃𝑘𝑜𝜅2cos2𝜃+𝑘2𝑧,𝑜+𝐵𝑒,𝑏,+𝑒𝑖𝑘𝑧,𝑒+𝐵𝑒,𝑔,+𝑒𝑖𝑘𝑧,𝑒(𝑧0)𝐵𝑒,𝑏,𝑒𝑖𝑘𝑧,𝑒×𝑘𝑧,𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒=𝑘𝑧,III𝐵𝑠,+𝑒𝑖𝑘𝑧,III𝜔cos𝜃+𝐵𝑝,+𝑒𝑖𝑘𝑧,III𝐵sin𝜃𝑜,𝑏,+𝑒𝑖𝑘𝑧,𝑜+𝐵𝑜,𝑔,+𝑒𝑖𝑘𝑧,𝑜(𝑧0)+𝐵𝑜,𝑏,𝑒𝑖𝑘𝑧,𝑜×𝜅2cos2𝜃+𝑘2𝑧,𝑜𝑘𝑜=𝑘𝑧,III𝐵𝑠,+𝑒𝑖𝑘𝑧,III𝜔sin𝜃𝐵𝑝,+𝑒𝑖𝑘𝑧,III𝐵cos𝜃𝑜,𝑏,+𝑒𝑖𝑘𝑧,𝑜+𝐵𝑜,𝑔,+𝑒𝑖𝑘𝑧,𝑜(𝑧0)𝐵𝑜,𝑏,𝑒𝑖𝑘𝑧,𝑜×𝑘𝑜𝑘𝑧,𝑜𝜀𝑜𝜔𝜇𝜅2cos2𝜃+𝑘2𝑧,𝑜+𝐵𝑒,𝑏,+𝑒𝑖𝑘𝑧,𝑒+𝐵𝑒,𝑔,+𝑒𝑖𝑘𝑧,𝑒(𝑧0)+𝐵𝑒,𝑏,𝑒𝑖𝑘𝑧,𝑒×𝜅2sin𝜃cos𝜃𝜀𝑜𝜔𝜇𝜅2cos2𝜃+𝑘2𝑧,𝑒=𝐵𝑠,+𝑒𝑖𝑘𝑧,IIIsin𝜃𝐵𝑝,+𝑒𝑖𝑘𝑧,III𝑘𝑧,III𝐵𝜔cos𝜃𝑒,𝑏,+𝑒𝑖𝑘𝑧,𝑒+𝐵𝑒,𝑔,+𝑒𝑖𝑘𝑧,𝑒(𝑧0)+𝐵𝑒,𝑏,𝑒𝑖𝑘𝑧,𝑒×𝑘2𝑧,𝑒𝜅2cos2𝜃𝜔𝜇𝜀𝑒𝜅2cos2𝜃+𝑘2𝑧,𝑒=𝐵𝑠,+𝑒𝑖𝑘𝑧,IIIcos𝜃𝐵𝑝,+𝑒𝑖𝑘𝑧,III𝑘𝑧,III𝜔sin𝜃.(31) From these equations, the origin of mode conversion is obviously accounted by the equations with both 𝑜- and 𝑒-mode. After solving the fields in regions I and III and linking to the sources with the help of (29), we can get Green functions for these regions:𝐺I𝐵,𝑃=𝑊2𝑂2𝑒𝑖𝑘𝑧,𝑜𝑧0𝑊1𝑂2𝑒𝑖𝑘𝑧,𝑜𝑧0+𝑊5𝑂2𝑒𝑖𝑘𝑧,𝑒𝑧0+𝑊4𝑂2𝑒𝑖𝑘𝑧,𝑒𝑧0𝑠,I𝑊̂𝑥5𝑂1𝑒𝑖𝑘𝑧,𝑒𝑧0+𝑊4𝑂1𝑒𝑖𝑘𝑧,𝑒𝑧0𝑠,I+̂𝑦𝑊2𝑂1𝑒𝑖𝑘𝑧,𝑜𝑧0+𝑊1𝑂1𝑒𝑖𝑘𝑧,𝑜𝑧0𝑊5𝑂3𝑒𝑖𝑘𝑧,𝑒𝑧0+𝑊4𝑂3𝑒𝑖𝑘𝑧,𝑒𝑧0𝑠,I+𝑊̂𝑧3𝑂2𝑒𝑖𝑘𝑧,𝑜𝑧0𝑊1𝑂2𝑒𝑖𝑘𝑧,𝑜𝑧0𝑊7𝑂2𝑒𝑖𝑘𝑧,𝑒𝑧0𝑊6𝑂2𝑒𝑖𝑘𝑧,𝑒𝑧0𝑝,I+𝑊̂𝑥7𝑂1𝑒𝑖𝑘𝑧,𝑒𝑧0+𝑊6𝑂1𝑒𝑖𝑘𝑧,𝑒𝑧0𝑝,I+𝑊̂𝑦3𝑂1𝑒