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

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: 𝜀⃡=⎛⎜⎜⎜⎜⎝𝜀¹0𝜀0₂000𝜀0₃⎞⎟⎟⎟⎟⎠,𝜇=𝜇.(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𝜋∇⋅𝑃⃗𝑟⃗,𝜀⃡_𝑟⋅𝐸⃗𝑟⃗=+−𝑖44𝜋𝜋𝜇𝜔_𝑟𝜇_𝑟𝑃⃗𝑟⃗∇×𝐸∇⃗⋅𝑟𝐵⃗⃗−𝑟𝑖⃗𝜔=𝐵⃗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 ∇×∇×𝐸⃗𝑟⃗−𝜔²𝜇_𝑟𝜀⃡_𝑟⋅𝐸⃗𝑟⃗=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𝑥×⎛⎜⎜⎜⎜⎝𝐸_𝑥⎞⎟⎟⎟⎟⎠𝐸_𝑦𝐸_𝑧⎞⎟⎟⎟⎟⎠=0,(4) where 𝜀_𝑥=𝜀₁,𝜀_𝑦=𝜀₂,𝜀_𝑧=𝜀₃ 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 𝜔⁴−𝜔²𝑘^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𝑧,𝑖−[𝑘_𝑧,𝑖(𝜔,𝑘_𝑥,𝑘_𝑦)]²=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 𝑑₁ and 𝑑₂, 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}(2𝑑𝜋𝑘)⃗³𝑃_⃗𝑘𝑖𝑒^{𝑖⃗𝑘𝑖⋅⃗𝑟}=𝑑𝑧₀(2𝑑𝜋𝜅)⃗²𝑀𝑟⃗=_𝑖=1,2𝛿𝑧−𝑧₀𝑃_⃗𝜅𝑒^{⃗𝑅𝑖⃗𝜅⋅},𝑑𝑧⁰𝑒^{±𝑖𝑘𝑧,𝑖𝑧0}(2𝑑𝜋𝑘)⃗³𝑀_⃗𝑘𝑖𝑒^{𝑖⃗𝑘𝑖⋅⃗𝑟}=𝑑𝑧₀(2𝑑𝜋𝜅)⃗²𝛿𝑧−𝑧₀𝑀_⃗𝜅𝑒^{⃗𝑅𝑖⃗𝜅⋅},(6) where 𝑧₀∈[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𝜋)²𝐵⃗𝑅⃗,𝑧=𝑑𝜅⃗𝐸⃗𝜅⃗,𝑧exp𝑖𝜅⃗⋅𝑅⃗,(2𝜋)²𝐵⃗𝜅⃗,𝑧exp𝑖𝜅⃗⋅𝑅⃗,(7) where 𝐸⃗𝜅⃗,𝑧=+𝑑𝑧𝑑₀𝐺𝑧₀⃡_𝐸𝑃𝜅⃗;𝑧−𝑧₀𝑃⃗𝜅⃗,𝑧₀𝐵⃗𝜅⃗,𝑧=𝑑𝑧₀𝐺⃡_𝐸𝑀𝜅⃗;𝑧−𝑧₀𝑀𝜅⃗,𝑧₀,+𝑑𝐺𝑧₀⃡_𝐵𝑃𝜅⃗;𝑧−𝑧₀𝑃⃗𝜅⃗,𝑧₀𝐺⃡_𝐵𝑀𝜅⃗;𝑧−𝑧₀𝑀𝜅⃗,𝑧₀.(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 𝑧₀, 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 𝑧₀ 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,𝑔,+𝑧e−xp𝑧₀−𝑖𝑘_𝑧,1𝑧₀+𝐵_2,𝑔,+exp−𝑖𝑘_𝑧,2𝑧₀+𝐵_1,𝑏,++𝐵_2,𝑏,+×+Θ𝐵𝑧_1,𝑔,−0exp𝑖𝑘_𝑧,1𝑧₀+𝐵_2,𝑔,−exp𝑖𝑘_𝑧,2𝑧₀−𝑧+𝐵_1,𝑏,−+𝐵_2,𝑏,−𝐷+𝐵⃗𝛿𝑧−𝑧₀exp𝑖𝜅⃗⋅_II𝑅⃗,𝜔,𝜅⃗=×𝐷Θ𝑧_1,𝑔,+e−x𝑧p₀−𝑖𝑘_𝑧,1𝑧₀+𝐷_2,𝑔,+exp−𝑖𝑘_𝑧,2𝑧₀+𝐷_1,𝑏,++𝐷_2,𝑏,++×Θ𝐷𝑧_1,𝑔,−0exp𝑖𝑘_𝑧,1𝑧₀+𝐷_2,𝑔,−exp𝑖𝑘_𝑧,2𝑧₀−𝑧+𝐷_1,𝑏,−+𝐷_2,𝑏,−𝐸+𝐷⃗𝛿𝑧−𝑧₀exp𝑖𝜅⃗⋅_II𝑅⃗,𝜔,𝜅⃗=×𝐸Θ_1,𝑔,+𝑧e−xp𝑧₀−𝑖𝑘_𝑧,1𝑧₀+𝐸_2,𝑔,+exp−𝑖𝑘_𝑧,2𝑧₀+𝐸_1,𝑏,++𝐸_2,𝑏,++×Θ𝐸𝑧_1,𝑔,−0exp𝑖𝑘_𝑧,1𝑧₀+𝐸_2,𝑔,−exp𝑖𝑘_𝑧,2𝑧₀−𝑧+𝐸_1,𝑏,−+𝐸_2,𝑏,−+𝐸⃗𝛿𝑧−𝑧₀exp𝑖𝜅⃗⋅𝑅⃗,(9) where 𝐵𝐵_{𝑖,𝑔,±𝑖,𝑏,±}=𝐵_𝑖,𝑔,±exp𝑖𝜅⃗⋅𝑅⃗exp±𝑖𝑘_𝑧,𝑖𝑧ℎ_𝑖,±,=𝐵𝐷_𝑖,𝑏,±ex_𝑖,𝑔,±p𝑖𝜅⃗⋅𝑅⃗exp±𝑖𝑘_𝑧,𝑖𝑧ℎ_𝑖,±,=−1𝑖𝜔𝜇∇×𝐵^𝑖,𝑔,±𝐸_𝑖,𝑔,±,=𝜀⃡⁻¹⋅𝐷_𝑖,𝑔,±𝐷_𝑖,𝑏,±=−1𝑖𝜔𝜇∇×𝐵^𝑖,𝑏,±𝐸_𝑖,𝑏,±,=𝜀⃡⁻¹⋅𝐷_𝑖,𝑏,±.𝑖=1,2,(10)𝐵⃗(𝐷⃗, 𝐸⃗) is to account the local (nonpropagating) field locating on the source plane, and Θ(𝑧−𝑧₀) 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𝜔,𝜅_𝑠,−⃗×exp𝑒𝑖̂𝜅⃗_𝑠,−^I⋅+𝑅⃗𝐵ex_𝑝,−p𝑒̂−_𝑝,−^I𝑖𝑘_𝑧,I𝐵_III𝑧,×ex𝜔p,𝜅⃗𝑖𝜅=⃗⋅𝑅𝐵⃗e^𝑠,+xℎp𝑖𝑘^𝑠,+III+_𝑧,III𝐵_𝑝,+ℎ_𝑝,+^III𝐷_III𝑧,=𝐵𝜔,𝜅⃗=𝐸_III𝜔,𝜅_𝑠,+⃗×exp𝑒𝑖̂𝜅⃗⋅_𝑠,+^III𝑅+⃗𝐵exp_𝑝,+𝑒̂𝑖𝑘_𝑝,+^III_𝑧,III𝑧,(11) where the unit vectors are defined as 𝑒̂_𝛼𝑠,±=𝑘^𝑠,±×𝑧̂||𝑘_𝑠,±×𝑧||̂=⎛⎜⎜⎜⎜⎝𝜅𝜅±csoi𝑘ns_𝑧𝜃𝜃ℎ⎞⎟⎟⎟⎟⎠×⎛⎜⎜⎜⎜⎝001⎞⎟⎟⎟⎟⎠=⎛⎜⎜⎜⎜⎝−scin0os𝜃𝜃_𝛼𝑠,±⎞⎟⎟⎟⎟⎠,=𝑘^𝑠,±×𝑒̂_𝑠,±||𝑘_𝑠,±×𝑒̂_𝑠,±||=𝜔√1𝜀_𝑟𝜇_𝑟⎛⎜⎜⎜⎜⎝±±𝑘𝑘_𝑧𝑧cos𝜃ℎ−s𝜅_𝛼𝑝,±in𝜃⎞⎟⎟⎟⎟⎠,=𝑘^𝑝,±×𝑧̂||𝑘_𝑝,±×𝑧||̂=⎛⎜⎜⎜⎜⎝𝜅𝜅±csoi𝑘ns_𝑧𝜃𝜃𝑒̂_𝛼𝑝,±⎞⎟⎟⎟⎟⎠×⎛⎜⎜⎜⎜⎝001⎞⎟⎟⎟⎟⎠=⎛⎜⎜⎜⎜⎝−scin0os𝜃𝜃⎞⎟⎟⎟⎟⎠,=−𝑘^𝑝,±×ℎ_𝑝,±|||𝑘_𝑝,±×ℎ_𝑝,±|||=𝜔√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 𝛿𝑧×−𝑧₀𝐵𝑧̂_1,𝑔,+−𝐵e_1,𝑔,−xp−𝑖𝑘_𝑧,1𝑧₀+𝐵_2,𝑔,+exp−𝑖𝑘_𝑧,2𝑧₀+𝑧̂×𝐵⃗𝛿_exp𝑖𝑘_𝑧,1𝑧₀−𝐵_2,𝑔,−exp𝑖𝑘_𝑧,2𝑧₀×𝐵⃗exp𝑖𝜅⃗𝑧⋅𝑅−⃗𝑧₀𝛿e𝑧x−p𝑧₀𝑖𝜅⃗⋅𝑅⃗+𝑖𝜅⃗+𝑖𝜔𝜇𝜀⃡⋅𝐸⃗𝛿𝑧−𝑧₀=−4𝜋𝑖𝜔𝜇𝑃⃗𝛿𝑧−𝑧₀exp𝑖𝜅⃗⋅𝑅⃗×𝑀exp𝑖𝜅⃗⋅𝑅⃗𝛿𝑧ex−p𝑧₀𝑖𝜅⃗⋅𝑅⃗+4𝜋𝜇𝑖𝜅⃗+4𝜋𝜇𝛿_𝛿𝑧−𝑧₀𝑧−𝑧₀𝑧̂×𝑀exp𝑖𝜅⃗⋅𝑅⃗,(13)×𝐸𝑧̂_1,𝑔,++𝐸e_2,𝑔,+xp−𝑖𝑘_𝑧,1𝑧₀−𝐸_1,𝑔,−exp𝑖𝑘_𝑧,1𝑧₀+𝑧̂×𝐸⃗𝛿_exp−𝑖𝑘_𝑧,2𝑧₀−𝐸_2,𝑔,−exp𝑖𝑘_𝑧,2𝑧₀×𝐸⃗exp𝑖𝜅𝑧⃗−⋅𝑅𝑧⃗₀𝛿ex𝑧p−𝑧𝑖𝜅₀⃗⋅𝑅⃗+𝑖𝜅⃗=𝑖𝜔𝐵⃗𝛿𝑧−𝑧₀exp𝑖𝜅⃗⋅𝑅⃗,(14) where we have used ∇Θ(𝑧−𝑧₀)=𝑧𝛿̂(𝑧−𝑧₀) and ∇𝛿(𝑧−𝑧₀)=𝑧𝛿̂′(𝑧−𝑧₀). By requiring the coefficients in front of 𝛿(𝑧−𝑧₀) and 𝛿′(𝑧−𝑧₀) 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𝑧 in region II and 𝜔²=𝜅²+𝑘^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 𝛿′(𝑧−𝑧₀) 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 𝛿(𝑧−𝑧₀), we can get the other equations, and we write them according to the components along 𝑒̂_𝑠±, 𝜅, and 𝑧̂ directions for (13), and (14) respectively, 1√𝜀𝜇𝐵^𝑠,𝑔,++𝐵_{𝑠,𝑔,−}=4𝜋𝑖𝜔²𝜇𝑘_𝑧𝑃𝑒̂_𝑠,±−4𝜋𝜇𝑖𝜅𝜔𝑘_𝑧𝐵𝑀𝑧_𝑝,𝑔,+̂,1𝑖𝜔𝜀−𝜇𝐸𝐵𝑧̂=_{𝑝,𝑔,−}−=𝑖4−𝜋4𝜔𝜋𝜇𝑖𝜔𝑃𝜇𝑧𝑃̂,𝜅,√𝜇𝜀𝐵^𝑝,𝑔,++𝐵_{𝑝,𝑔,−}+𝜅𝑘^II𝑘_𝑧,II=𝑖𝜔𝑘^II𝑘_𝑧,II1𝐵𝑒̂_𝑠,±,√𝜇𝜀𝐵^𝑠,𝑔,+𝐵−𝑧̂𝐵=0._{𝑠,𝑔,−}=𝑖𝜔𝐵𝜅,(18) The fields launching from the source can be gotten by solving (18) and the results are collected as follows: 𝐵_𝑠,𝑔,+=2𝜋𝑖𝜔²𝜇𝑘_𝑧,II√𝜇𝜀𝑒̂_𝑠,±^II⋅𝑃⃗+2𝜋𝜇𝑖𝜔²𝜇𝜀𝑘_𝑧,IIℎ^𝑠,+II𝐵_{𝑠,𝑔,−}⋅𝑀,=2𝜋𝑖𝜔²𝜇𝑘_𝑧,II√𝜇𝜀𝑒̂_𝑠,±^II⋅𝑃⃗+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𝑧₀𝑒̂_{𝑠,−𝑇,}^II+𝐵_{𝑠,𝑏,−}𝑒̂_{𝑠,−𝑇,}^II,𝑒̂_{𝑝,−𝑇,}^I=1√𝜇𝜀𝐵^𝑝,𝑏,+𝐵_𝑝,−𝑒̂_𝑝,+𝑇,^II+𝐵_{𝑝,𝑔,−}exp𝑖𝑘_𝑧,II𝑧₀𝑒̂_{𝑝,−𝑇,}^II+𝐵_{𝑝,𝑏,−}𝑒̂_{𝑝,−𝑇,}^II,ℎ_{𝑝,−𝑇,}^I𝐵=_𝑠,−𝐵_𝑝,𝑏,+ℎ_𝑝,+𝑇,^II+𝐵_{𝑝,𝑔,−}exp𝑖𝑘_𝑧,II𝑧₀ℎ_{𝑝,−𝑇,}^II+𝐵_{𝑝,𝑏,−}ℎ_{𝑝,−𝑇,}^II,ℎ_{𝑠,−𝑇,}^I=𝐵_𝑠,𝑏,+ℎ_𝑠,+𝑇,^II+𝐵_{𝑠,𝑔,−}exp𝑖𝑘_𝑧,II𝑧₀ℎ_{𝑠,−𝑇,}^II+𝐵_{𝑠,𝑏,−}ℎ_{𝑠,−𝑇,}^II.(20) Similarly, the boundary equation (16) for 𝑧=ℎ gives 𝐵_𝑠,+exp𝑖𝑘_𝑧,Iℎ𝑒̂_𝑠,+𝑇,^III=1√𝜇𝜀𝐵^𝑠,𝑔,+exp𝑖𝑘_𝑧,IIℎ−𝑧₀𝑒̂_𝑠,+𝑇,^II+𝐵_𝑠,𝑏,+exp𝑖𝑘_𝑧,IIℎ𝑒̂_𝑠,+𝑇,^II𝐵_𝑝,++𝐵_{𝑠,𝑏,−}exp−𝑖𝑘_𝑧,IIℎ𝑒̂_{𝑠,−𝑇,}^II,exp𝑖𝑘_𝑧,Iℎ𝑒̂_𝑝,+𝑇,^III=1√𝜇𝜀𝐵^𝑝,𝑔,+exp𝑖𝑘_𝑧,IIℎ−𝑧₀𝑒̂_𝑝,+𝑇,^II+𝐵_𝑝,𝑏,+exp𝑖𝑘_𝑧,IIℎ𝑒̂_𝑝,+𝑇,^II𝐵_𝑝,++𝐵_{𝑝,𝑏,−}exp−𝑖𝑘_𝑧,IIℎ𝑒̂_{𝑝,−𝑇,}^II,exp𝑖𝑘_𝑧,𝐼ℎℎ_𝑝,+𝑇,^III=𝐵_𝑝,𝑔,+exp𝑖𝑘_𝑧,IIℎ−𝑧₀ℎ_𝑝,+𝑇,^II𝐵_𝑠,++𝐵_𝑝,𝑏,+exp𝑖𝑘_𝑧,IIℎℎ_𝑝,+𝑇,^II+𝐵_{𝑝,𝑏,−}exp−𝑖𝑘_𝑧,IIℎℎ_{𝑝,−𝑇,}^II,exp𝑖𝑘_𝑧,Iℎℎ_𝑠,+𝑇,^III=𝐵_𝑠,𝑔,+exp𝑖𝑘_𝑧,IIℎ−𝑧₀ℎ_𝑠,+𝑇,^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𝑧0}1−𝑡^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑠,𝑔,−}+𝑟^𝑠II,I𝑡^𝑠_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑡^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑠,𝑔,+×ℎ_𝑠,−+1exp𝑖𝜅⃗⋅𝑅⃗exp−𝑖𝑘_𝑧,I𝑧√𝜇𝜀𝑟^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}1−𝑡^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑝,𝑔,−}+𝑟^𝑝II,I𝑡^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑡^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑝,𝑔,+×ℎ_𝑝,−𝐸_Iexp𝑖𝜅⃗⋅𝑅⃗exp−𝑖𝑘_𝑧,I𝑧,=𝜔,𝜅1⃗√𝜇𝜀𝑟^𝑠II,I𝑒^{𝑖𝑘𝑧,II𝑧0}1−𝑡^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑠,𝑔,−}+𝑟^𝑠II,I𝑡^𝑠_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑡^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑠,𝑔,+×𝑒̂_𝑠,−+1exp𝑖𝜅⃗⋅𝑅⃗exp−𝑖𝑘_𝑧,I𝑧√𝜇𝜀𝑟^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}1−𝑡^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑝,𝑔,−}+𝑟^𝑝II,I𝑡^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑡^𝑝_II,I²𝑒^{𝑖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,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑠,𝑔,−}+𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑡^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑠,𝑔,+×ℎ_𝑠,++1𝑒^{−𝑖𝑘𝑧,Iℎ}exp𝑖𝜅⃗⋅𝑅⃗exp𝑖𝑘_𝑧,III𝑧√𝜀𝜇𝑟^𝑝II,III𝑡^𝑝_II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}1−𝑡^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑝,𝑔,−}+𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑡^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑝,𝑔,+×ℎ_𝑝,+𝐸_III𝑒^{−𝑖𝑘𝑧,Iℎ}exp𝑖𝜅⃗⋅𝑅⃗exp𝑖𝑘_𝑧,III𝑧,=𝜔,1𝜅⃗√𝜀𝜇𝑟^𝑠II,III𝑡^𝑠_II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}1−𝑡^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑠,𝑔,−}+𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑡^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_𝑠,𝑔,+×𝑒̂_𝑠,++1𝑒^{−𝑖𝑘𝑧,Iℎ}exp𝑖𝜅⃗⋅𝑅⃗exp𝑖𝑘_𝑧,III𝑧√𝜀𝜇𝑟^𝑝II,III𝑡^𝑝_II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}1−𝑡^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝐵_{𝑝,𝑔,−}+𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑡^𝑝_II,III²𝑒^{𝑖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_𝐵𝑀2𝜋𝜇𝜅𝑖⃗;𝜔𝑧𝑘−_II𝑧₀𝑘_𝑧,II⎛⎜⎜⎝𝑡^𝑠_II,I𝑟^𝑠II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}ℎ_𝑠,−^Iℎ_𝑠,+^II+𝑒^{𝑖𝑘𝑧,II𝑧0}ℎ_𝑠,−^Iℎ_𝑠,−^II1−𝑟^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}+𝑡^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}+𝑟^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑟^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}ℎ_𝑝,−^I𝑒̂_𝑠,±^II×𝑒_{⃗𝑅𝑖⃗𝜅⋅}⎞⎟⎟⎠𝐺⃡_I𝐵𝑃𝑒^{−𝑖𝑘𝑧,I𝑧},=2𝜋𝜅𝑖𝜔⃗;₂𝑧−𝑧₀𝜇𝑘_𝑧,II×⎛⎜⎜⎝𝑡^𝑠II,I𝑒^{𝑖𝑘𝑧,II𝑧0}+𝑟^𝑠_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑟^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}ℎ_𝑠,−^I𝑒̂_𝑠,±^II−𝑡^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}ℎ_𝑝,−^Iℎ_𝑠,−^II+𝑟^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}ℎ_𝑝,−^Iℎ_𝑠,+^II1−𝑟^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}×𝑒_{⃗𝑅𝑖⃗𝜅⋅}⎞⎟⎟⎠𝐺⃡_I𝐸𝑀𝑒^{−𝑖𝑘𝑧,I𝑧},=2𝜋𝜇𝜅𝑖⃗;𝜔𝑧𝑘−_II𝑧₀𝑘_𝑧,II×⎡⎢⎢⎣𝑟^𝑠II,I𝑡^𝑠II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}𝑒̂_𝑠,−^Iℎ_𝑠,+^II+𝑒^{𝑖𝑘𝑧,II𝑧0}𝑒̂_𝑠,−^Iℎ_𝑠,−^II1−𝑡^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}+𝑟^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}+𝑡^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑡^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝑒̂_𝑝,−^I𝑒̂_𝑠,±^II×𝑒_{⃗𝑅𝑖⃗𝜅⋅}⎤⎥⎥⎦𝐺⃡_I𝐸𝑀𝑒^{−𝑖𝑘𝑧,I𝑧},=2𝜋𝑖𝜅𝜔⃗;₂𝑧−𝑧₀𝜇𝑘_𝑧,II⎡⎢⎢⎣𝑡^𝑠_II,I𝑒^{𝑖𝑘𝑧,II𝑧0}+𝑟^𝑠_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}1−𝑟^𝑠_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}𝑒̂_𝑠,−^I𝑒̂_𝑠,±^II−𝑡^𝑝II,I𝑒^{𝑖𝑘𝑧,II𝑧0}𝑒̂_𝑝,−^Iℎ_𝑠,−^II+𝑟^𝑝_II,I𝑒^{𝑖𝑘𝑧,II(2ℎ−𝑧0)}𝑒̂_𝑝,−^Iℎ_𝑠,+^II1−𝑟^𝑝_II,I²𝑒^{𝑖2𝑘𝑧,IIℎ}×𝑒_{⃗𝑅𝑖⃗𝜅⋅}⎤⎥⎥⎦𝐺⃡_III𝐵𝑀𝑒^{−𝑖𝑘𝑧,I𝑧},=2𝜋𝜇𝜅𝑖⃗;𝜔𝑧𝑘−_II𝑧₀𝑘_𝑧,II×⎡⎢⎢⎣𝑡^𝑠II,III𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}ℎ_𝑠,+^IIIℎ_𝑠,−^II+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}ℎ_𝑠,+^IIIℎ_𝑠,+^II1−𝑟^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}+𝑡^𝑝II,III𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑟^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}ℎ_𝑝,+^III𝑒̂_𝑠,±^II×𝑒_{𝑖𝜅𝑧,III(𝑧−ℎ)}⎤⎥⎥⎦𝐺⃡_III𝐵𝑃𝑒^{𝑖⃗⃗𝑅𝑘⋅},=2𝜋𝜅𝑖𝜔⃗;₂𝑧−𝑧₀𝜇𝑘_𝑧,II×⎡⎢⎢⎣𝑡^𝑠II,III𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑟^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}ℎ_𝑠,+^III𝑒̂_𝑠,±^II−𝑡^𝑝II,III𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}ℎ_𝑝,+^IIIℎ_𝑠,−^II+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}ℎ_𝑝,+^IIIℎ_𝑠,+^II1−𝑟^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}×𝑒_{𝑖𝜅𝑧,III(𝑧−ℎ)}⎤⎥⎥⎦𝑒^{𝑖⃗⃗𝑅𝑘⋅}𝐺=⃡^III_𝐸𝑀2𝜋𝜇𝜅𝑖⃗𝑧𝜔−𝑘_II𝑧₀𝑘_𝑧,II×⎡⎢⎢⎣𝑡^𝑠II,III𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}𝑒̂_𝑠,+^IIIℎ_𝑠,−^II+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}𝑒̂_𝑠,+^IIIℎ_𝑠,+^II1−𝑟^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}+𝑡^𝑝II,III𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑟^𝑝_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝑒̂_𝑝,+^III𝑒̂_𝑠,±^II×𝑒_{𝑖𝜅𝑧,III(𝑧−ℎ)}⎤⎥⎥⎦𝐺⃡_III𝐸𝑃𝑒^{𝑖⃗⃗𝑅𝑘⋅},=2𝜋𝜅𝑖𝜔⃗𝑧₂−𝑧₀𝜇𝑘_𝑧,II×⎡⎢⎢⎣𝑡^𝑠II,III𝑟^𝑠II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}1−𝑟^𝑠_II,III²𝑒^{𝑖2𝑘𝑧,IIℎ}𝑒̂_𝑠,+^III𝑒̂_𝑠,±^II−𝑡^𝑝II,III𝑟^𝑝II,III𝑒^{𝑖𝑘𝑧,II(𝑧0+ℎ)}𝑒̂_𝑝,+^IIIℎ_𝑠,−^II+𝑒^{𝑖𝑘𝑧,II(ℎ−𝑧0)}𝑒̂_𝑝,+^IIIℎ_𝑠,+^II1−𝑟^𝑝_II,III²𝑒^{𝑖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𝑧_𝑒𝜇=0.(25) Thus, there are two dispersions gotten: the mode with dispersion 𝜔²−𝑘^2𝑥/𝜀_𝑜𝜇−𝑘^2𝑦/𝜀_𝑜𝜇−𝑘^2𝑧/𝜀_𝑜𝜇=0 is known as 𝑜-light, and the mode with dispersion 𝜔²−𝑘^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𝜅²cos²𝜃+𝑘^2𝑧,𝑜⎛⎜⎜⎜⎜⎝±𝑘_𝑧,𝑜𝑒̂_𝑒,±−𝜅c0os𝜃⎞⎟⎟⎟⎟⎠,=1𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝜀^2𝑒+𝜅²sin²𝜃𝜀^2𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑒×⎛⎜⎜⎜⎜⎝𝜀^𝑜−1𝜀𝜅²sin𝜃cos𝜃_𝑒−1±−𝜀𝑘^2𝑧,𝑒_𝑜−1−𝜅²cos²𝜃𝑘_𝑧,𝑒𝜅sin𝜃⎞⎟⎟⎟⎟⎠.(26) Thus, we can have the unit vectors for magnetic field and electric displacement vector as ℎ_𝑜±=1𝑘_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜⎛⎜⎜⎜⎜⎝𝑘−𝜅²sin𝜃cos𝜃_2𝑧,𝑜∓𝑘+_𝑧,𝑜𝜅²cos²𝜃ℎ𝜅sin𝜃_𝑒±⎞⎟⎟⎟⎟⎠,=1𝜅²cos²𝜃+𝑘^2𝑧,𝑒⎛⎜⎜⎜⎜⎝±𝑘_𝑧,𝑒𝑑−𝜅c0os𝜃_𝑜,±⎞⎟⎟⎟⎟⎠,=1𝜅²cos²𝜃+𝑘^2𝑧,𝑜⎛⎜⎜⎜⎜⎝±𝑘_𝑧,𝑜𝑑−𝜅c0os𝜃_𝑒±⎞⎟⎟⎟⎟⎠,=1𝑘_𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑒⎛⎜⎜⎜⎜⎝−𝑘𝜅²_2𝑧,𝑒sin𝜃cos𝜃±𝑘−_𝑧,𝑒𝜅²cos²𝜃𝜅sin𝜃⎞⎟⎟⎟⎟⎠,(27) where 𝑘_𝑜=𝜅²+𝑘^2𝑧,𝑜 and 𝑘_𝑒=𝜅²+𝑘^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 𝛿(𝑧−𝑧₀) and 𝛿′(𝑧−𝑧₀) 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: 𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑘_𝑜𝐵_{𝑜,𝑔,−}=−𝑖4𝜋𝜔𝜇𝑃_𝑥−𝐵_𝑜,𝑔,+𝜅₂+𝑖4𝜋𝜇𝜅sin𝜃𝑀_𝑧,sin𝜃cos𝜃𝑘_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜+𝑘_𝑧,𝑒𝐵_{𝑜,𝑔,−}−𝐵_𝑜,𝑔,+𝜅²cos²𝜃+𝑘^2𝑧,𝑒=−𝑖4𝜋𝜔𝜇𝑃_𝑦𝐵_𝑒,𝑔,++𝐵_{𝑒,𝑔,−}𝐸_𝑧−𝑖4𝜋𝜇𝜅cos𝜃𝑀_𝑧,=−4𝜋𝜀_𝑜𝑃_𝑧,𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝜔𝜇𝜀_𝑒+𝑖𝐸_𝑧𝐵_𝑒,𝑔,+−𝐵_{𝑒,𝑔,−}𝑘𝜅sin𝜃−𝑖𝜔𝐵_𝑥=0,_𝑜𝜔𝜇𝜀_𝑜𝑘_𝑧,𝑜𝑘^2𝑧,𝑜+𝜅²cos²𝜃𝐵^𝑜,𝑔,++1+𝐵_{𝑜,𝑔,−}𝜔𝜇𝜀_𝑜𝜅²sin𝜃cos𝜃𝜅²cos²𝜃+𝑘^2𝑧,𝑒−𝑖𝐸_𝑧𝐵_𝑒,𝑔,+−𝐵_{𝑒,𝑔,−}𝜅cos𝜃−𝐵𝑖𝜔𝐵_𝑦=0,_𝑧=0.(28) After solving (28), the fields launching from the source plane are determined: 𝐵𝐵_{𝑜,𝑔,+𝑜,𝑔,−}=−𝑂₁𝑃_𝑧+𝑂₂𝑃_𝑥+𝑂₃𝑀_𝑦−𝑂₄𝑀_𝑥−𝑂₅𝑀_𝑧,𝐵_𝑒,𝑔,+=−𝑂₁𝑃_𝑧−𝑂₂𝑃_𝑥+𝑂₃𝑀_𝑦−𝑂₄𝑀_𝑥+𝑂₅𝑀_𝑧,𝐵_{𝑒,𝑔,−}=−𝑂₁𝑃_𝑦+𝑂₂𝑃_𝑥+𝑂₃𝑃_𝑧+𝑂₄𝑀_𝑥−𝑂₅𝑀_𝑧,=−𝑂₁𝑃_𝑦+𝑂₂𝑃_𝑥−𝑂₃𝑃_𝑧−𝑂₄𝑀_𝑥−𝑂₅𝑀_𝑧,(29) where the coefficients in front of the source are 𝑂₁=⎛⎜⎜⎜⎝𝑖2𝜋𝜅𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝜔𝜇𝑅cos𝜃𝑘_𝑜𝑘_𝑧,𝑜+𝑖2𝜋𝜅³𝜔𝜇𝜀_𝑒𝑅sin²𝜃cos𝜃𝑘_𝑜𝑘_𝑧,𝑜𝜀_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝑂₂⎞⎟⎟⎟⎠,=𝑖2𝜋𝜔𝑘^𝑜𝜇𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑂₃,=𝑖2𝜋𝜔²𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝜇²𝜀_𝑜𝑅𝑘_𝑜𝑘_𝑧,𝑜𝑂₄,=𝑖2𝜋𝜔²𝜅²𝜇²𝜀_𝑒𝑅sin𝜃cos𝜃𝑘_𝑜𝑘_𝑧,𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝑂₅,=𝑖2𝜋𝜇𝜅𝑘^𝑜sin𝜃𝜅²cos²𝜃+𝑘^2𝑧,𝑜,𝑂₁=𝑖2𝜋𝜔𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝜇𝑘_𝑧,𝑒𝑅,𝑂₂=𝑖2𝜋𝜔𝜅²𝜇sin𝜃cos𝜃𝑘_𝑧,𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑅,𝑂₃=𝑖2𝜋𝜔𝜅𝜇𝜀^𝑒sin𝜃𝜀_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑒,𝑂₄=𝑖2𝜋𝜔²𝜇²𝜀_𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑒,𝑂₅=⎡⎢⎢⎢⎣𝑖2𝜋𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝜅𝜇cos𝜃+𝑖2𝜋𝜅₃𝑘_𝑧,𝑒𝑅𝜇sin²𝜃cos𝜃𝑘_𝑧,𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑅⎤⎥⎥⎥⎦,(30) where the coefficient 𝑅 in (30) is given as 𝑅=𝑘^2𝑧,𝑜+𝜅²cos²𝜃/𝑘^2𝑧,𝑒+𝜅²cos²𝜃. 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 −𝜅²sin𝜃cos𝜃𝑘_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜+𝑘_𝑧,𝑒𝐵_𝑜,𝑏,++𝐵_{𝑜,𝑏,−}+𝐵_{𝑜,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑜𝑧0}𝜅²cos²𝜃+𝑘^2𝑧,𝑒=−𝐵_𝑠,−𝐵_𝑒,𝑏,+−𝐵_{𝑒,𝑏,−}−𝐵_{𝑒,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑒𝑧0}𝑘_𝑧,I𝜔cos𝜃+𝐵^𝑝,−sin𝜃,𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑘_𝑜𝐵_𝑜,𝑏,+=−𝐵_𝑠,−+𝐵_{𝑜,𝑏,−}+𝐵_{𝑜,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑜𝑧0}𝑘_𝑧,I𝜔sin𝜃−𝐵^𝑝,−𝐵cos𝜃,_𝑜,𝑏,+−𝐵_{𝑜,𝑏,−}−𝐵_{𝑜,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑜𝑧0}𝑘^𝑜𝑘_𝑧,𝑜𝜔𝜇𝜀_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜+𝐵_𝑒,𝑏,++𝐵_{𝑒,𝑏,−}+𝐵_{𝑒,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑒𝑧0}𝜅²sin𝜃cos𝜃𝜔𝜇𝜀_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑒=𝐵_𝑠,−sin𝜃+𝐵_𝑝,−𝑘_𝑧,I𝐵𝜔cos𝜃,_𝑒,𝑏,++𝐵_{𝑒,𝑏,−}+𝐵_{𝑒,𝑔,−}𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑘^2𝑧,𝑒−𝜅²cos²𝜃𝜔𝜇𝜀_𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑒=−𝐵_𝑠,−cos𝜃+𝐵_𝑝,−𝑘_𝑧,I𝐵𝜔sin𝜃,_𝑜,𝑏,+×𝑒−𝜅^{𝑖𝑘𝑧,𝑜ℎ}₂+𝐵_𝑜,𝑔,+𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝐵_{𝑜,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑜ℎ}sin𝜃cos𝜃𝑘_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜×+𝐵_𝑒,𝑏,+𝑒𝑘^{𝑖𝑘𝑧,𝑒ℎ}_𝑧,𝑒+𝐵_𝑒,𝑔,+𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}−𝐵_{𝑒,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑒ℎ}𝜅²cos²𝜃+𝑘^2𝑧,𝑒=𝑘^𝑧,III𝐵_𝑠,+𝑒^{𝑖𝑘𝑧,IIIℎ}𝜔cos𝜃+𝐵^𝑝,+𝐵𝑒^{𝑖𝑘𝑧,IIIℎ}sin𝜃_𝑜,𝑏,+×𝑒^{𝑖𝑘𝑧,𝑜ℎ}+𝐵_𝑜,𝑔,+𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝐵_{𝑜,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑜ℎ}𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑘_𝑜=𝑘^𝑧,III𝐵_𝑠,+𝑒^{𝑖𝑘𝑧,IIIℎ}𝜔sin𝜃−𝐵^𝑝,+𝐵𝑒^{𝑖𝑘𝑧,IIIℎ}cos𝜃_𝑜,𝑏,+×𝑒^{𝑖𝑘𝑧,𝑜ℎ}+𝐵𝑘_{𝑜,𝑔,+𝑜}𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝐵_{𝑜,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑜ℎ}𝑘_𝑧,𝑜𝜀_𝑜𝜔𝜇𝜅²cos²𝜃+𝑘^2𝑧,𝑜×+𝐵𝜅_𝑒,𝑏,+2𝑒^{𝑖𝑘𝑧,𝑒ℎ}+𝐵_𝑒,𝑔,+𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝐵_{𝑒,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑒ℎ}sin𝜃cos𝜃𝜀_𝑜𝜔𝜇𝜅²cos²𝜃+𝑘^2𝑧,𝑒=𝐵_𝑠,+𝑒^{𝑖𝑘𝑧,IIIℎ}sin𝜃−𝐵_𝑝,+𝑒^{𝑖𝑘𝑧,IIIℎ}𝑘_𝑧,III𝐵𝜔cos𝜃_𝑒,𝑏,+×𝑒−^{𝑖𝑘𝑧,𝑒ℎ}𝑘+_2𝑧,𝑒𝐵_𝑒,𝑔,+𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝐵_{𝑒,𝑏,−}𝑒^{−𝑖𝑘𝑧,𝑒ℎ}−𝜅²cos²𝜃𝜔𝜇𝜀_𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑒=−𝐵_𝑠,+𝑒^{𝑖𝑘𝑧,IIIℎ}cos𝜃−𝐵_𝑝,+𝑒^{𝑖𝑘𝑧,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_𝐵,𝑃=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎢⎣−𝑊²−𝑂𝑊₂𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₂𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₅𝑂₂𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑠,−I𝑥₅̂+−𝑊₂𝑂₁𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑠,−I𝑦̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₁𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₃𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑠,−I𝑧₃⎤⎥⎥⎥⎥⎥⎦̂+𝑂𝑊₂𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₂𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₇𝑂₂𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑝,−I𝑥₇̂+𝑊𝑂₁𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₆𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑝,−I𝑦₃̂×exp𝑖𝜅⃗⋅𝑅𝑂⃗¹e𝑒xp^{𝑖𝑘𝑧,𝑜𝑧0}−+𝑖𝑘𝑊_1𝑧,I𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₃𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}⎤⎥⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐸,𝑃𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎢⎣−𝑊²−𝑂𝑊₂𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₂𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₅𝑂₂𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑠,−I𝑥₅̂+−𝑊₂𝑂₁𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑠,−I𝑦̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₁𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₃𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑠,−I𝑧₃⎤⎥⎥⎥⎥⎥⎦̂+𝑂₂𝑊𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₂𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₇𝑂₂𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑝,−I𝑥₇̂+𝑊𝑂₁𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₆𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑝,−I𝑦₃̂×exp𝑖𝜅⃗⋅𝑅𝑂⃗₁𝑒exp^{𝑖𝑘𝑧,𝑜𝑧0}−+𝑖𝑊𝑘_1𝑧,I𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₃𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑝,−I𝑧⎤⎥⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐵,𝑀𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣−𝑊²+𝑂𝑊₄𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₄𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₄𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑠,−I𝑥₂̂+𝑊𝑂₃𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}ℎ_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎣𝑊𝑂₅𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₅𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₄𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑠,−I𝑧₃⎤⎥⎥⎥⎥⎦̂−𝑂𝑊₄𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₄𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₄𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑝,−I𝑥₃̂+−𝑊₃𝑂₃𝑒^{𝑖𝑘𝑧,𝑜𝑧0}𝑀_𝑦+𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}ℎ_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑂⃗₅𝑒exp^{𝑖𝑘𝑧,𝑜𝑧0}−+𝑖𝑊𝑘_1𝑧,I𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₅𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₆𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}ℎ^𝑝,−I𝑧⎤⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐸,𝑀𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣−𝑊²+𝑂𝑊₄𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₄𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₄𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₄𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑠,−I𝑥₂̂+𝑊𝑂₃𝑒^{𝑖𝑘𝑧,𝑜𝑧0}−𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}𝑒̂_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₅𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}−𝑊₅𝑂₅𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₄𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑠,−I𝑧₃̂⎤⎥⎥⎥⎥⎦−𝑂₄𝑊𝑒^{𝑖𝑘𝑧,𝑜𝑧0}+𝑊₁𝑂₄𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₄𝑒^{𝑖𝑘𝑧,𝑒𝑧0}−𝑊₆𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑝,−I𝑥₃̂+−𝑊₃𝑂₃𝑒^{𝑖𝑘𝑧,𝑜𝑧0}𝑀_𝑦+𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}𝑒̂_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅⃗𝑂e₅x𝑒p−^{𝑖𝑘𝑧,𝑜𝑧0}𝑖+𝑘𝑊_1𝑧,I𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜𝑧0}+𝑊₇𝑂₅𝑒^{𝑖𝑘𝑧,𝑒𝑧0}+𝑊₆𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒𝑧0}𝑒̂^𝑝,−I𝑧⎤⎥⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_III𝐵,𝑃𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎢⎣−𝑊¹−𝑂𝑊₂𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₂𝑂₂𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₄𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₂𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑠,+III𝑥₄̂+−𝑊₁𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₁𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑠,+III𝑦̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₂𝑂₁𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₃𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑠,+III𝑧₁̂⎤⎥⎥⎥⎥⎥⎦−𝑂₂𝑊𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₃𝑂₂𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₆𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₂𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑝,+III𝑥₆̂+𝑊𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₁𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑝,+III𝑦₁̂×exp𝑖𝜅⃗⋅𝑅𝑂⃗₁e𝑒xp𝑖𝑘^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}_𝑧,III+𝑊₃𝑂₁𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₃𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑝,+III𝑧̂⎤⎥⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐸,𝑃,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎢⎣−𝑊¹−𝑂𝑊₂𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₂𝑂₂𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₄𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₂𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑠,+III𝑥₄̂+−𝑊₁𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₁𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑠,+III𝑦̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₁𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₂𝑂₁𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₃𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑠,+III𝑧₁⎤⎥⎥⎥⎥⎥⎦̂−𝑂₂𝑊𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₃𝑂₂𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₆𝑂₂𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₂𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑝,+III𝑥₆̂+𝑊𝑂₁𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₁𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑝,+III𝑦₁̂×exp𝑖𝜅⃗⋅𝑂𝑅₁⃗𝑒exp𝑖𝑘^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+_𝑧,III𝑊₃𝑂₁𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₃𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₃𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑝,+III𝑧⎤⎥⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐵,𝑀,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣−𝑊¹+𝑂𝑊₄𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₂𝑂₄𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₄𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₄𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑠,+III𝑥₁̂+−𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₂𝑂₃𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}ℎ_𝑠,+^III𝑦̂+⎡⎢⎢⎢⎢⎣𝑊𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₂𝑂₅𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}𝑒^{−𝑖𝑘𝑧ℎ}+𝑊₁𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₅𝑂₅𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑠,+III𝑧₁⎤⎥⎥⎥⎥⎦̂−𝑂₄𝑊𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₃𝑂₄𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₄𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑝,+III𝑥₁̂+−𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₃𝑂₃𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}ℎ_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑂⃗₅e𝑒xp𝑖𝑘^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}_𝑧,III+𝑊₃𝑂₅𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}−𝑊₇𝑂₅𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}ℎ^𝑝,+III𝑧⎤⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐸,𝑀,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣−𝑊¹+𝑂₄𝑊𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₂𝑂₄𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₄𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₅𝑂₄𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑠,+III𝑥₁̂+−𝑊₄𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₂𝑂₃𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑦̂+⎡⎢⎢⎢⎢⎢⎣𝑊𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₂𝑂₅𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}𝑒^{−𝑖𝑘𝑧ℎ}+𝑊₁𝑂₅𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₅𝑂₅𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑠,+III𝑧₁̂⎤⎥⎥⎥⎥⎦−𝑂₄𝑊𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₃𝑂₄𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₄𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}+𝑊₇𝑂₄𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑝,+III𝑥₁̂+−𝑊₁𝑂₃𝑒^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}+𝑊₃𝑂₃𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅⃗exp𝑂₅𝑒𝑖𝑘^{−𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}_𝑧,III+𝑊₃𝑂₅𝑒^{𝑖𝑘𝑧,𝑜(ℎ−𝑧0)}−𝑊₆𝑂₅𝑒^{−𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}−𝑊₇𝑂₅𝑒^{𝑖𝑘𝑧,𝑒(ℎ−𝑧0)}𝑒̂^𝑝,+III𝑧⎤⎥⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠(𝑧−ℎ).(32)

The green functions in (32) describe the relation between the fields generated outside of the uniaxial crystal and the source embedded in uniaxial crystal. (The parameters (𝑊_𝑖,𝑖=1,2,…,7) in (32) will be given in Appendix A.)

3.2. Sources in Biaxial Crystal

In the subsection, the Green functions linking the fields in region I and III and sources embedded in biaxial crystal in region II are given. The relative permittivity tensor now has the form diag(𝜀₁,𝜀₂,𝜀₃) with 𝜀₁≠𝜀₂≠𝜀₃. Thus, the equation determining the dispersions of eigenmodes now becomes 𝜔⁴−𝜔²𝑘^2𝑥+𝑘^2𝑦𝜀₃𝜇+𝑘^2𝑥+𝑘^2𝑧𝜀₂𝜇+𝑘^2𝑧+𝑘^2𝑦+𝑘_2𝑥𝜀₁𝜇𝜀₂𝜀₃𝜇²+𝑘^2𝑦𝜀₁𝜀₃𝜇²+𝑘^2𝑧𝜀₁𝜀₃𝜇²𝑘^2𝑥+𝑘^2𝑧+𝑘^2𝑦=0.(33) The dispersion can be gotten as followes: 𝑘^2𝑧,1=(−𝑏+√𝑏²−4𝑎𝑐)/2𝑎 and 𝑘^2𝑧,2=(−𝑏−√𝑏²−4𝑎𝑐)/2𝑎 (each of them denotes two modes propagating (or decaying) along positive and negative 𝑧 direction), where 𝑎, 𝑏, and 𝑐 are, respectively, the coefficients in front of the fourth-order, second-order and zeroth-order terms of 𝑘_𝑧, and they have the following concrete forms: 𝑎=1𝑐=𝜀𝑘₁𝜀_32𝑥𝜇,𝑏=𝜀₂𝑘𝜀^2𝑥₃𝜇+𝜀₁𝑘𝜀^2𝑦₃𝜇−𝜀𝜔₂𝜇²−𝜀𝜔₁𝜇²,+𝑘^2𝑦𝑘^2𝑥+𝑘_2𝑥𝜀₂𝜀₃𝜇+𝜀₁𝑘𝜀^2𝑦₃𝜇𝜀₃1𝜇+𝜀𝜔₂𝜇²+𝑘^2𝑦𝜀₃1𝜇+𝜀𝜔₁𝜇²+𝜔⁴.(34) The corresponding unit vectors of electric fields of eigenmodes can be solved as [23] 𝑒̂_𝑖,±=𝑋^2𝑖+1+𝑍^2𝑖^−1/2𝑋_𝑖,1,±𝑍_𝑖,(35) where 𝑋_𝑖=𝜅²cos𝜃sin𝜃𝜅²+𝑘^2𝑧,𝑖−𝜔²𝜀₃𝜇𝜔²𝜀₁𝜇−𝑘^2𝑧,𝑖−𝜅²sin²𝜃𝜔²𝜀₃𝜇−𝜅²−𝜅²cos²𝜃𝑘^2𝑧,𝑖𝑍_𝑖=𝑘^𝑧,𝑖𝜅sin𝜃𝜅²+𝑘^2𝑧,𝑖−𝜔²𝜀₁𝜇𝜔²𝜀₁𝜇−𝜅²sin²𝜃−𝑘^2𝑧,𝑖𝜔²𝜀₃𝜇−𝜅²−𝜅²cos²𝜃𝑘^2𝑧,𝑖.(36) In (36) 𝑘_𝑧,𝑖 is a positive square root from 𝑘^2𝑧,𝑖 so that it denotes the propagating or decaying wave along positive 𝑧 direction. Thus, the unit vectors of electric displacement vector and magnetic field can be gotten asℎ𝑑_{𝑖,±𝑖,±}=𝜀²¹𝑋^2𝑖+𝜀²²+𝜀²³𝑍^2𝑖^−1/2𝜀₁𝑋_𝑖,𝜀₂,±𝜀₃𝑍_𝑖,=1𝑁_𝑖±𝜅sin𝜃𝑍_𝑖where𝑁_𝑖−𝑘_𝑧,𝑖,±𝑘_𝑧,𝑖𝑋_𝑖−𝑍_𝑖𝜅cos𝜃,𝜅cos𝜃−𝜅sin𝜃𝑋_𝑖,=𝜅sin𝜃𝑍_𝑖−𝑘_𝑧,𝑖²+𝑘_𝑧,𝑖𝑋_𝑖−𝑍_𝑖𝜅cos𝜃²+𝜅cos𝜃−𝜅sin𝜃𝑋_𝑖².(37)

Then we can determine the fields launching from the source plane with the help of these unit vectors defined. The general equations (13) and (14) can be mapped to concrete forms with each equation gotten giving the field equality in different directions. The relations, (𝐵_𝑥,𝐵_𝑦)=4𝜋𝜇(𝑀_𝑥,𝑀_𝑦) and (𝐸_𝑥,𝐸_𝑦)=0, are found again by following the fact of vanishing of the coefficients in front of 𝛿′ function. The equations determining the fields launching from the source plane are: 𝑘₁𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹^1/2+𝑘₂𝐵_1,𝑔,−−𝐵_1,𝑔,+𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²²^1/2+𝑖𝐸_𝑧𝐵_2,𝑔,−−𝐵_2,𝑔,+𝜅𝑋s₁in𝜃=𝑖𝜔𝐵_𝑥,𝑘₁𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹^1/2+𝑋₂𝐵_1,𝑔,−−𝐵_1,𝑔,+𝑘₂𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²²^1/2−𝑖𝐸_𝑧𝐵_2,𝑔,−−𝐵_2,𝑔,+𝜅cos𝜃=𝑖𝐵𝜔_𝑧𝐵_𝑦,−𝑘=0_𝑧,1𝑋₁−𝑍₁𝜅cos𝜃𝑁₁𝐵_1,𝑔,+−𝑘_𝑧,2+𝐵_1,𝑔,−𝑋₂−𝑍₂𝜅cos𝜃𝑁₂𝐵_2,𝑔,+=−𝑖4𝜋𝜔𝜇𝑃_𝑥+𝐵_2,𝑔,−𝜅sin𝜃𝑍₁+𝑖4𝜋𝜇𝜅sin𝜃𝑀_𝑧,−𝑘_𝑧,1𝑁₁𝐵_1,𝑔,++𝜅sin𝜃𝑍₂+𝐵_1,𝑔,−−𝑘_𝑧,2𝑁₂𝐵_2,𝑔,+=−𝑖4𝜋𝜔𝜇𝑃_𝑦+𝐵_2,𝑔,−𝐸_𝑧−𝑖4𝜋𝜇𝜅cos𝜃𝑀_𝑧,=−4𝜋𝜀₃𝑃_𝑧,(38) where 𝑘₁=𝜅²+𝑘^2𝑧,1 and 𝑘₂=𝜅²+𝑘^2𝑧,2. The solutions of (35) and (36) give the relations between the fields launching from the source plane and the sources as 𝐵_1,𝑔,−=12𝑉¹⁰𝐵_1,𝑔,+𝑃_𝑥+𝑉₁₁𝑃_𝑦+𝑉₁₂𝑀_𝑧+𝑉₁𝑀_𝑥+𝑉₂𝑀_𝑦+𝑉₃𝑃_𝑧=12𝑉¹⁰𝐵_2,𝑔,−𝑃_𝑥+𝑉₁₁𝑃_𝑦+𝑉₁₂𝑀_𝑧−𝑉₁𝑀_𝑥−𝑉₂𝑀_𝑦−𝑉₃𝑃_𝑧=12𝑉⁷𝐵_2,𝑔,+𝑃_𝑥+𝑉₈𝑃_𝑦+𝑉₉𝑀_𝑧+𝑉₄𝑀_𝑥+𝑉₅𝑀_𝑦+𝑉₆𝑃_𝑧=12𝑉⁷𝑃_𝑥+𝑉₈𝑃_𝑦+𝑉₉𝑀_𝑧−𝑉₄𝑀_𝑥−𝑉₅𝑀_𝑦−𝑉₆𝑃_𝑧.(39) The details of the definitions of parameters 𝑉_𝑖(𝑖=1,2,…,12) in (39) are given in Appendix B. To determine the fields generated in regions I and III, the boundary equations should be solved, and we firstly write down their concrete forms from the continuity of the tangential fields for the interfaces at 𝑧=0 and 𝑧=ℎ, respectively (we have used the fact that 𝑘_𝑧=𝑘_𝑧,I=𝑘_𝑧,III in the following equations, when 𝜅⃗ is the same), 𝑆₁𝐵+𝑆_1,𝑏,+2−𝐵_1,𝑏,−−𝐵_1,𝑔,−𝑒^{𝑖𝑘𝑧,1𝑧0}=−𝐵_𝑠,−𝐵_2,𝑏,+−𝐵_2,𝑏,−−𝐵_2,𝑔,−𝑒^{𝑖𝑘𝑧,2𝑧0}𝑘_𝑧𝜔cos𝜃+𝐵^𝑝,−𝑆sin𝜃,₃𝐵+𝑆_1,𝑏,+4−𝐵_1,𝑏,−−𝐵_1,𝑔,−𝑒^{𝑖𝑘𝑧,1𝑧0}=−𝐵_𝑠,−𝐵_2,𝑏,+−𝐵_2,𝑏,−−𝐵_2,𝑔,−𝑒^{𝑖𝑘𝑧,2𝑧0}𝑘_𝑧𝜔sin𝜃−𝐵^𝑝,−𝑆cos𝜃,₅𝑋+₁𝑆𝐵_1,𝑏,+6+𝐵_1,𝑏,−+𝐵_1,𝑔,−𝑒^{𝑖𝑘𝑧,1𝑧0}=𝐵_𝑠,−𝑋₂𝐵_2,𝑏,++𝐵_2,𝑏,−+𝐵_2,𝑔,−𝑒^{𝑖𝑘𝑧,2𝑧0}sin𝜃+𝐵_𝑝,−𝑘_𝑧𝑆𝜔cos𝜃,₅𝐵+𝑆_1,𝑏,+6+𝐵_1,𝑏,−+𝐵_1,𝑔,−𝑒^{𝑖𝑘𝑧,1𝑧0}=−𝐵_𝑠,−𝐵_2,𝑏,++𝐵_2,𝑏,−+𝐵_2,𝑔,−𝑒^{𝑖𝑘𝑧,2𝑧0}cos𝜃+𝐵_𝑝,−𝑘_𝑧𝑆𝜔sin𝜃,₁𝐵+𝑆_1,𝑏,+2𝑒^{𝑖𝑘𝑧,1ℎ}−𝐵_1,𝑏,−𝑒^{−𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑔,+𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}=𝑘_𝑧𝐵_2,𝑏,+𝑒^{𝑖𝑘𝑧,2ℎ}−𝐵_2,𝑏,−𝑒^{−𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑔,+𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}cos𝜃𝜔𝐵^𝑠,+𝑆₃𝑒^{𝑖𝑘𝑧ℎ}+sin𝜃𝐵_𝑝,+𝑒^{𝑖𝑘𝑧ℎ},𝐵+𝑆_1,𝑏,+4𝑒^{𝑖𝑘𝑧,1ℎ}−𝐵_1,𝑏,−𝑒^{−𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑔,+𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}=𝑘_𝑧𝐵_2,𝑏,+𝑒^{𝑖𝑘𝑧,2ℎ}−𝐵_2,𝑏,−𝑒^{−𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑔,+𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}sin𝜃𝜔𝐵^𝑠,+𝑆₅𝑒^{𝑖𝑘𝑧ℎ}−cos𝜃𝐵_𝑝,+𝑒^{𝑖𝑘𝑧ℎ},𝑋+₁𝑆𝐵_1,𝑏,+6𝑒^{𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑏,−𝑒^{−𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑔,+𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}=sin𝜃𝑋𝐵₂_𝑠,+𝐵_2,𝑏,+𝑒^{𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑏,−𝑒^{−𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑔,+𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒^{𝑖𝑘𝑧ℎ}−𝑘^𝑧cos𝜃𝜔𝐵^𝑝,+𝑒^{𝑖𝑘𝑧ℎ}𝑆₅𝐵+𝑆_1,𝑏,+6𝑒^{𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑏,−𝑒^{−𝑖𝑘𝑧,1ℎ}+𝐵_1,𝑔,+𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}=−cos𝐵𝜃𝐵_{2,𝑏,+𝑠,+}𝑒^{𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑏,−𝑒^{−𝑖𝑘𝑧,2ℎ}+𝐵_2,𝑔,+𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒^{𝑖𝑘𝑧ℎ}−𝑘^𝑧sin𝜃𝜔𝐵^𝑝,+𝑒^{𝑖𝑘𝑧ℎ},(40) where the parameters are given as 𝑆₁=𝜅sin𝜃𝑍¹−𝑘_𝑧,1𝑁₁,𝑆₂=𝜅sin𝜃𝑍²−𝑘_𝑧,2𝑁₂𝑆₃,=𝑘^𝑧,1𝑋₁−𝑍₁𝜅cos𝜃𝑁₁,𝑆₄=𝑘^𝑧,2𝑋₂−𝑍₂𝜅cos𝜃𝑁₂𝑆₅,=𝑘¹𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹𝑆₆,=𝑘²𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²².(41) The fields in regions I and III can be solved out from (40) and (41). Consequently, the Green functions relating them to the sources in biaxial crystal can be gotten with the help of (39) and the definitions given in (8):𝐺^I_𝐵,𝑃=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈²+𝑈𝑉₁₀𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₀𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₇𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₇𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑠,−^I𝑥₂̂+𝑈𝑉₁₁𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₈𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₈𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₃₃𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₃𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₆𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₆𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ^𝑠,−I𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₀₃𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₀𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₇𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₇𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑥̂+−𝑈₃𝑉₁₁𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₈𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₈𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅⃗𝑉e₃x𝑒p−^{𝑖𝑘𝑧,1𝑧0}𝑖+𝑘𝑈_1𝑧,I𝑉₃𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₆𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₆𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐸,𝑃𝑧,=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈²+𝑈𝑉₁₀𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₀𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₇𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₇𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑥₂̂+𝑈𝑉₁₁𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₈𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₈𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₃₃𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₃𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₆𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₆𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₀₃𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₀𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₇𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₇𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑥̂+−𝑈₃𝑉₁₁𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₈𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₈𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅⃗e𝑉x₃p𝑒−𝑖^{𝑖𝑘𝑧,1𝑧0}𝑘+_𝑧,I𝑈₁𝑉₃𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₆𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₆𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐵,𝑀𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈²+𝑈𝑉₁𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₁𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₄𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₄𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑠,−^I𝑥₂̂+𝑈𝑉₂𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₅𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₅𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₁₂₃𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₂𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₉𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₉𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑠,−^I𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₃𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₁𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₄𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₄𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑥̂+−𝑈₃𝑉₂𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₅𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₅𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₁₂e𝑒xp^{𝑖𝑘𝑧,1𝑧0}−−𝑖𝑘𝑈_1𝑧,I𝑉₁₂𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₉𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₉𝑒^{−𝑖𝑘𝑧,2𝑧0}ℎ_𝑝,−^I𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_I𝐸,𝑀𝑧,=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈²+𝑈𝑉₁𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₁𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₄𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₄𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑥₂̂+𝑈𝑉₂𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₅𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₄𝑉₅𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑦₂̂+⎡⎢⎢⎢⎢⎣−𝑉𝑈₁₂₃𝑒^{𝑖𝑘𝑧,1𝑧0}−𝑈₁𝑉₁₂𝑒^{−𝑖𝑘𝑧,1𝑧0}+𝑈₅𝑉₉𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₄𝑉₉𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑠,−^I𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₃𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₁𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₄𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₄𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑥̂+−𝑈₃𝑉₂𝑒^{𝑖𝑘𝑧,1𝑧0}+𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₅𝑒^{𝑖𝑘𝑧,2𝑧0}+𝑈₆𝑉₅𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₁₂e𝑒xp^{𝑖𝑘𝑧,1𝑧0}−−𝑖𝑘𝑈_1𝑧,I𝑉₁₂𝑒^{−𝑖𝑘𝑧,1𝑧0}−𝑈₇𝑉₉𝑒^{𝑖𝑘𝑧,2𝑧0}−𝑈₆𝑉₉𝑒^{−𝑖𝑘𝑧,2𝑧0}𝑒̂_𝑝,−^I𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺_III𝐵,𝑃𝑧,=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈¹+𝑈𝑉₁₀𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₀𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₇𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₇𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑥₁̂+𝑈𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₈𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₈𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑦₁̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₃₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₃𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₆𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₆𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₀₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₃𝑉₁₀𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₇𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₇𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑥̂+−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₃𝑉₁₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₈𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₈𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₃e𝑒xp𝑖𝑘^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}_𝑧,III+𝑈₃𝑉₃𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₆𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₆𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐸,𝑃,=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈¹+𝑈𝑉₁₀𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₀𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₇𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₇𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑥₁̂+𝑈𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₈𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₈𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑦₁̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₃₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₃𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₆𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₆𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑧̂⎤⎥⎥⎥⎥⎦+−𝑈𝑉₁₀₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₃𝑉₁₀𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₇𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₇𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑥̂+−𝑈₁𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₃𝑉₁₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₈𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₈𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₃𝑒exp𝑖𝑘^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}_𝑧,III+𝑈₃𝑉₃𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₆𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₆𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑧̂⎤⎥⎥⎥⎥⎦⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐵,𝑀,=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈¹+𝑈𝑉₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₄𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₄𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑥₁̂+𝑈𝑉₂𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₅𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₅𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑦₁̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₁₂₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₉𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₉𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑠,+^III𝑧⎤⎥⎥⎥⎥⎦̂+−𝑈𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₃𝑉₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₄𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₄𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑥̂+−𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₃𝑉₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₅𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₅𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₁₂e𝑒xp𝑖𝑘^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}_𝑧,III−𝑈₃𝑉₁₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₉𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₉𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}ℎ_𝑝,+^III𝑧⎤⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠𝐺(𝑧−ℎ)_III𝐸,𝑀,=12⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝⎡⎢⎢⎢⎢⎣𝑈¹+𝑈𝑉₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₄𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₄𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑥₁̂+𝑈𝑉₂𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₂𝑉₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₅𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₅𝑉₅𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑦₁̂+⎡⎢⎢⎢⎢⎣−𝑈𝑉₁₂₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}−𝑈₂𝑉₁₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₄𝑉₉𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₅𝑉₉𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑠,+^III𝑧⎤⎥⎥⎥⎥⎦̂+−𝑈𝑉₁₁𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₃𝑉₁𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₄𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₄𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑥̂+−𝑈₁𝑉₂𝑒^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₃𝑉₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₅𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}−𝑈₇𝑉₅𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑦̂×exp𝑖𝜅⃗⋅𝑅𝑉⃗₁₂e𝑒xp𝑖𝑘^{−𝑖𝑘𝑧,1(ℎ−𝑧0)}_𝑧,III−𝑈₃𝑉₁₂𝑒^{𝑖𝑘𝑧,1(ℎ−𝑧0)}+𝑈₆𝑉₉𝑒^{−𝑖𝑘𝑧,2(ℎ−𝑧0)}+𝑈₇𝑉₉𝑒^{𝑖𝑘𝑧,2(ℎ−𝑧0)}𝑒̂_𝑝,+^III𝑧⎤⎥⎥⎥⎥⎦̂⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠(𝑧−ℎ).(42)

Green functions gotten in (42) give the relations between the fields generated in region I (or region III) and the sources embedded in biaxial crystal (region II). The definitions of parameters (𝑈_𝑖,𝑖=1,2,…,7) in (42) will be given in Appendix C.

4. Conclusion Remark

In summary, a new Green function formulism is developed by combining eigenmode expansion method and conventional formulism. The correction of the formulism is supported by recovering well-known results gotten by conventional formulism. The powerful eigenmode expansion method removes the difficulties brought by the interface between anisotropic media. Given to the advantage, the Green functions relating the fields generated outside of anisotropic medium and the sources inside were calculated and shown. They may have potential applications in studying the electromagnetic response and its spatial distribution of anisotropic metamaterial medium, which has well-defined effective medium tensor.

Appendices

A. Definitions of Parameters 𝑊_𝑖 in (32)

The parameters 𝑊_𝑖 come out, when one solves the equation groups (31). Thus, it can be related to the coefficients in front of the field amplitudes. In the solution processes, they firstly relate to parameters 𝐼_𝑖,𝑖=1,2,…,6, and 𝐹_𝑖,𝑖=1,2,…,9, as 𝑊₁=𝐼³𝐼₅−𝐼₂𝐼₆𝐼₁𝐼₅−𝐼₂𝐼₄𝑒^{𝑖𝑘𝑧,𝑜ℎ}𝑊₂,=𝐹¹𝐼₅−𝐹₄𝐼₂𝑊₃𝐹₇𝐼₁𝐼₅−𝐼₂𝐼₄𝑒^{−𝑖𝑘𝑧,𝑜ℎ},=𝐹¹𝐼₄−𝐹₄𝐼₁𝑊₄𝐹₇𝐼₂𝐼₄−𝐼₁𝐼₅𝑒^{−𝑖𝑘𝑧,𝑜ℎ},=𝐼²𝐼₁𝐼₅−𝐼₂𝐼₄𝑊₅𝑒^{𝑖𝑘𝑧,𝑜ℎ},=𝐼⁵𝐼₁𝐼₅−𝐼₂𝐼₄𝑊₆𝑒^{−𝑖𝑘𝑧,𝑜ℎ},=𝐼¹𝐼₂𝐼₄−𝐼₁𝐼₅𝑊₇𝑒^{𝑖𝑘𝑧,𝑜ℎ},=𝐼⁴𝐼₂𝐼₄−𝐼₁𝐼₅𝑒^{−𝑖𝑘𝑧,𝑜ℎ}.(A.1) The definitions of parameters 𝐼_𝑖,𝑖=1,2,…,6, are the following: 𝐼₁=𝐹₂−𝐹¹𝐹₈𝐹₇𝐼₂,=𝐹₃−𝐹¹𝐹₉𝐹₇𝐼₃,=𝐹¹𝐷₆𝐹₇𝐷₅+𝐷²𝐷₅𝐼₄,=𝐹₅−𝐹⁴𝐹₈𝐹₇𝐼₅,=𝐹₆−𝐹⁴𝐹₉𝐹₇,𝐼₆=𝐹⁴𝐷₆𝐹₇𝐷₅−𝐷¹𝐷₅.(A.2) The definitions of parameters 𝐹_𝑖,𝑖=1,2,…,9, are given as 𝐹₁=𝐷²𝐷₇𝐷₅𝐹₂+𝐷₄,=𝐷₁𝑒^{−𝑖𝑘𝑧,𝑒ℎ}+𝐷²𝐷₆𝐷₅𝐹₃𝑒^{𝑖𝑘𝑧,𝑜ℎ},=𝐷₃𝑒^{−𝑖𝑘𝑧,𝑒ℎ}−𝐷²𝐷₈𝐷₅𝐹₄𝑒^{𝑖𝑘𝑧,𝑜ℎ},=−𝐷¹𝐷₇𝐷₅𝐹₅−𝐷₃,=−𝐷₂𝑒^{𝑖𝑘𝑧,𝑒ℎ}−𝐷¹𝐷₆𝐷₅𝐹₆𝑒^{𝑖𝑘𝑧,𝑜ℎ},=−𝐷₄𝑒^{𝑖𝑘𝑧,𝑒ℎ}+𝐷¹𝐷₈𝐷₅𝐹₇𝑒^{𝑖𝑘𝑧,𝑜ℎ},=𝐷₈−𝐷⁶𝐷₇𝐷₅𝐹₈,=𝐷₅𝑒^{−𝑖𝑘𝑧,𝑜ℎ}−𝐷²⁶𝐷₅𝐹₉𝑒^{𝑖𝑘𝑧,𝑜ℎ},=−𝐷₇𝑒^{−𝑖𝑘𝑧,𝑜ℎ}+𝐷⁶𝐷₈𝐷₅𝑒^{𝑖𝑘𝑧,𝑜ℎ}.(A.3) The parameters 𝐷_𝑖,𝑖=1,2,…,8, relate to the coefficients in (31) directly (since the media in region I and region III are the same, we used the fact that 𝑘_𝑧=𝑘_𝑧,I=𝑘_𝑧,III, when 𝜅⃗ is the same): 𝐷₁=𝑘^𝑧2𝐴₂𝐴₃𝜔𝐴³cos𝜃+𝐴₁sin𝜃+𝜔𝜇𝜀^𝑒𝑅2𝑘_𝑜𝐴₃𝐷cos𝜃,₂=𝑘^𝑧2𝐴₂𝐴₃𝜔𝐴³cos𝜃+𝐴₁sin𝜃−𝜔𝜇𝜀^𝑒𝑅2𝑘_𝑜𝐴₃𝐷cos𝜃,₃=12𝐴₂𝐴₃𝐴₃sin𝜃−𝐴₁cos𝜃+𝑘^𝑧𝜇𝜀_𝑒𝑅2𝑘_𝑜𝐴₃𝐷sin𝜃,₄=12𝐴₂𝐴₃𝐴₃sin𝜃−𝐴₁cos𝜃−𝑘^𝑧𝜇𝜀_𝑒𝑅2𝑘_𝑜𝐴₃𝐷sin𝜃,₅=𝑘^𝑧,𝑒𝜔𝜇𝜀_𝑒√𝜇𝜀_𝑜𝑅²2𝑘_𝑧,𝑜𝑘_𝑜𝐴₂𝐴₃𝑘^𝑜sin𝜃𝜔𝜇𝜀_𝑒𝑅𝐴³−𝑅cos𝜃√𝜇𝜀_𝑜𝐴₁+𝑘^𝑧sin𝜃𝐷2𝐴₃𝜔₆,=𝑘^𝑧,𝑒𝜔𝜇𝜀_𝑒√𝜇𝜀_𝑜𝑅²2𝑘_𝑧,𝑜𝑘_𝑜𝐴₂𝐴₃𝑘^𝑜sin𝜃𝜔𝜇𝜀_𝑒𝑅𝐴³−𝑅cos𝜃√𝜇𝜀_𝑜𝐴₁−𝑘^𝑧sin𝜃𝐷2𝐴₃𝜔₇,=𝑘^𝑧,𝑒𝜔𝜇𝜀_𝑒√𝜇𝜀_𝑜𝑅²2𝑘_𝑧,𝑜𝑘_𝑜𝐴₂𝐴₃𝑘^𝑜𝑘_𝑧cos𝜃𝜔²𝜇𝜀_𝑒𝑅𝐴³+𝑅𝑘^𝑧sin𝜃𝜔√𝜇𝜀_𝑜𝐴₁+c2o𝐴s₃𝜃𝐷₈,=𝑘^𝑧,𝑒𝜔𝜇𝜀_𝑒√𝜇𝜀_𝑜𝑅²2𝑘_𝑧,𝑜𝑘_𝑜𝐴₂𝐴₃𝑘^𝑜𝑘_𝑧cos𝜃𝜔²𝜇𝜀_𝑒𝑅𝐴³+𝑅𝑘^𝑧sin𝜃𝜔√𝜇𝜀_𝑜𝐴₁−c2o𝐴s₃𝜃,(A.4) where the parameters 𝐴_𝑖(𝑖=1,2,3) are given as 𝐴₁=𝜅²sin𝜃cos𝜃𝑘_𝑜𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝐴₂,=𝑘^𝑧,𝑒𝜅²cos²𝜃+𝑘^2𝑧,𝑒𝐴₃,=𝜅²cos²𝜃+𝑘^2𝑧,𝑜𝑘_𝑜.(A.5)

B. Definitions of Parameters 𝑉_𝑖 in (39)

The parameters come from the processes of solving equation groups (38), and they can relate to the coefficients of them as 𝑉₁=−𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹^1/2𝑖4𝜋𝜔𝜇𝑋₂𝑉₂𝑘₁𝑋₁−𝑋₂,=𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹^1/2𝑖4𝜋𝜔𝜇𝑉₃𝑘₁𝑋₁−𝑋₂,=−𝑖4𝜋𝜅𝜔𝜇𝜀²¹𝑋²¹+𝜀²²+𝜀²³𝑍²¹^1/2cos𝜃+𝑋₁sin𝜃𝑉₄𝑘₁𝜀₃𝑋₁−𝑋₂,=−𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²²^1/2𝑖4𝜋𝜇𝜔𝑋₁𝑉₅𝑘₂𝑋₂−𝑋₁,=𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²²^1/2𝑖4𝜋𝜇𝜔𝑉₆𝑘₂𝑋₂−𝑋₁,=−𝑖4𝜋𝜅𝜔𝜇𝜀²¹𝑋²²+𝜀²²+𝜀²³𝑍²²^1/2cos𝜃+𝑋₁sin𝜃𝑉₇𝜀₃𝑘₂𝑋₂−𝑋₁,=−𝑖4𝜋𝜔𝜇𝑁²𝜅sin𝜃𝑍₁−𝑘_𝑧,1𝑉₈𝒜,=−𝑖4𝜋𝜔𝜇𝑁²𝑘_𝑧,1𝑋₁−𝑍₁𝜅cos𝜃𝑉₉𝒜,=𝑁²𝑖4𝜋𝜇𝜅sin𝜃𝜅sin𝜃𝑍₁−𝑘_𝑧,1−ℬ𝑉₁₀𝒜,=−𝑖4𝜋𝜔𝜇𝑁¹𝜅sin𝜃𝑍₂−𝑘_𝑧,2𝑉₁₁𝒞,=−𝑖4𝜋𝜔𝜇𝑁¹𝑘_𝑧,2𝑋₂−𝑍₂𝜅cos𝜃𝑉₁₂𝒞,=𝑁¹𝑖4𝜋𝜇𝜅sin𝜃𝜅sin𝜃𝑍₂−𝑘_𝑧,2−𝒟𝒞,(B.1) where 𝒜 denotes (𝜅sin𝜃𝑍₂−𝑘_𝑧,2)(𝑘_𝑧,1𝑋₁−𝑍₁𝜅cos𝜃)−(𝑘_𝑧,2𝑋₂−𝑍₂𝜅cos𝜃)(𝜅sin𝜃𝑍₁−𝑘_𝑧,1), ℬ denotes 𝑖4𝜋𝜇𝜅cos𝜃(𝑘_𝑧,1𝑋₁−𝑍₁𝜅cos𝜃), 𝒞 denotes (𝜅sin𝜃𝑍₁−𝑘_𝑧,1)(𝑘_𝑧,2𝑋₂−𝑍₂𝜅cos𝜃)−(𝑘_𝑧,1𝑋₁−𝑍₁𝜅cos𝜃)(𝜅sin𝜃𝑍₂−𝑘_𝑧,2), and 𝒟 denotes 𝑖4𝜋𝜇𝜅cos𝜃(𝑘_𝑧,2𝑋₂−𝑍₂𝜅cos𝜃).

C. Definitions of Parameters 𝑈_𝑖 in (42)

The parameters 𝑈_𝑖 are gotten in the processes of solving equation groups (40) and (41), and they can be related to the coefficients in front of the field amplitudes. They firstly relate to parameters 𝐼^𝑖,𝑖=1,2,…,6, and 𝐹^𝑖,𝑖=1,2,…,9, which will defined later: 𝑈₁=𝐼^3𝐼^5−𝐼^2𝐼^6𝐼^1𝐼^5−𝐼^2𝐼^4𝑈₂𝑒^{𝑖𝑘𝑧,1ℎ},=𝐹^1𝐼^5−𝐹^4𝐼^2𝐹₇𝐼^1𝐼^5−𝐼^2𝐼^4𝑒^{−𝑖𝑘𝑧,1ℎ}𝑈₃,=𝐹^1𝐼^4−𝐹^4𝐼^1𝐹^7𝐼^2𝐼^4−𝐼^1𝐼^5𝑒^{−𝑖𝑘𝑧,1ℎ}𝑈₄,=𝐼^2𝐼^1𝐼^5−𝐼^2𝐼^4𝑈₅,=𝐼^5𝐼^1𝐼^5−𝐼^2𝐼^4𝑈₆𝑒^{−𝑖𝑘𝑧,1ℎ},=𝐼^1𝐼^2𝐼^4−𝐼^1𝐼^5𝑈₇𝑒^{𝑖𝑘𝑧,1ℎ},=𝐼^4𝐼^2𝐼^4−𝐼^1𝐼^5𝑒^{−𝑖𝑘𝑧,1ℎ}.(C.1) The parameters 𝐼^𝑖,𝑖=1,2,…,6, are given as functions of 𝐹^𝑖,𝑖=1,2,…,9, and 𝑄_𝑖,𝑖=1,2,…,8𝐼^1=𝐹^2−𝐹^1𝐹^8𝐹^7𝐼_2,=𝐹^3−𝐹^1𝐹^9𝐹^7𝐼_3,=𝐹^1𝑄₁𝐹^7𝑄₃+𝑄⁵𝑄₃𝐼_4,=𝐹^5−𝐹^4𝐹^8𝐹^7𝐼_5,=𝐹^6−𝐹^4𝐹^9𝐹^7𝐼_6,=𝐹^4𝑄₁𝐹^7𝑄₃−𝑄⁷𝑄₃.(C.2) The parameters 𝐹^𝑖,𝑖=1,2,…,9, are defined as the functions of 𝑄_𝑖,𝑖=1,2,…,8, 𝐹^1=𝑄⁴𝑄₅𝑄₃𝐹_2+𝑄₆,=𝑄₇𝑒^{−𝑖𝑘𝑧,2ℎ}+𝑄⁵𝑄₁𝑄₃𝐹_3𝑒^{𝑖𝑘𝑧,1ℎ},=𝑄₈𝑒^{−𝑖𝑘𝑧,2ℎ}−𝑄⁵𝑄₂𝑄₃𝐹_4𝑒^{𝑖𝑘𝑧,1ℎ},=−𝑄⁷𝑄₄𝑄₃𝐹_5−𝑄₈,=−𝑄₅𝑒^{𝑖𝑘𝑧,2ℎ}−𝑄⁷𝑄₁𝑄₃𝐹_6𝑒^{𝑖𝑘𝑧,1ℎ},=−𝑄₆𝑒^{𝑖𝑘𝑧,2ℎ}+𝑄⁷𝑄₂𝑄₃𝐹_7𝑒^{𝑖𝑘𝑧,1ℎ},=𝑄₂−𝑄¹𝑄₄𝑄₃𝐹_8,=𝑄₃𝑒^{−𝑖𝑘𝑧,1ℎ}−𝑄²¹𝑄₃𝐹_9𝑒^{𝑖𝑘𝑧,1ℎ},=−𝑄₄𝑒^{−𝑖𝑘𝑧,1ℎ}+𝑄¹𝑄₂𝑄₃𝑒^{𝑖𝑘𝑧,1ℎ}.(C.3) The parameters 𝑄_𝑖,𝑖=1,2,…,8, relate with the coefficients in (41) and (42) directly: 𝑄₁=−𝑋¹cos𝜃+sin𝜃𝑋₂+𝑋₂−𝑋₁sin𝜃2𝑆₅𝑋₂−𝑋₁𝑋₁+𝑆⁴cos𝜃−𝑆₂sin𝜃𝑘_𝑧𝑄2𝑆₂𝑆₃−𝑆₁𝑆₄𝜔₂,=−𝑆⁴sin𝜃+𝑆₂cos𝜃+𝑋2₁𝑆₂𝑆₃−𝑆₁𝑆₄sin𝜃−cos𝜃𝑋₂+𝑋₂−𝑋₁cos𝜃2𝑆₅𝑋₂−𝑋₁𝑋₁𝑘_𝑧𝑄𝜔₃,=𝑋²−𝑋₁sin𝜃−𝑋₂sin𝜃+cos𝜃𝑋₁−𝑆₄2𝑆₅𝑋₁𝑋₂−𝑋₁cos𝜃−𝑆₂sin𝜃𝑄2𝑆₂𝑆₃−𝑆₁𝑆₄𝑘𝜔₄^𝑧,=cos𝜃𝑋²−𝑋₁−𝑋₂cos𝜃−𝑋₁sin𝜃+𝑆2𝑆₅𝑋₁𝑋₂−𝑋₁𝑘𝜔₄^𝑧sin𝜃+𝑆₂cos𝜃𝑄₅2𝑆₂𝑆₃−𝑆₁𝑆₄=−−𝑆²𝑆₃−𝑆₁𝑆₄cos𝜃−𝑆₁𝑆₄cos𝜃−𝑆₂sin𝜃2𝑆₂𝑆₃−𝑆₁𝑆₄𝑆₂𝑘_𝑧+𝑋𝜔₁cos𝜃+sin𝜃𝑄₆2𝑆₆𝑋₂−𝑋₁,=cos𝜃−𝑋¹sin𝜃+2𝑆𝑆₆₂𝑋₂−𝑋₁𝑘𝜔^𝑧𝑆₃−𝑆₁𝑆₄sin𝜃+𝑆₁𝑆₄sin𝜃+𝑆₂cos𝜃𝑄₇2𝑆₂𝑆₃−𝑆₁𝑆₄𝑆₂,=cos𝜃𝑆²𝑆₃−𝑆₁𝑆₄+𝑆₁𝑆₄cos𝜃−𝑆₂sin𝜃+sin𝜃+cos𝜃𝑋2𝑆₂₁𝑆₂𝑆₃−𝑆₁𝑆₄𝑘𝜔^𝑧𝑄₈2𝑆₆𝑋₂−𝑋₁,=sin𝜃𝑆²𝑆₃−𝑆₁𝑆₄+𝑆₁𝑆₄sin𝜃+𝑆₂cos𝜃−cos𝜃−𝑋₁2𝑆₂𝑆₂𝑆₃−𝑆₁𝑆₄sin𝜃2𝑆₆𝑋₂−𝑋₁𝑘𝜔^𝑧.(C.4)

Acknowledgments

W. Wen would like to thank the support for this work by the Hong Kong RGC Grant HKUST2/CRF/11G. B. Hou would like to acknowledge the support from the National Natural Science Foundation of China (Grant no. 11104198) and from a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).