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

Authors: Shunbo Li; Xiao Xiao; Bo Hou; Weijia Wen
Journal: International Journal of Optics
 (2012)

---

## 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.


---

## Body

## 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 thicknessh 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)ε⃡=(ε1000ε2000ε3),μ=μ.
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]:
(2)∇·(ε⃡r·E⃑(r⃑))=-4π∇·P⃑(r⃑),∇×B⃑(r⃑)+iω~μrε⃡r·E⃑(r⃑)=-i4πω~μrP⃑(r⃑)+4πμr∇×M⃑(r⃑)∇·B⃑(r⃑)=0∇×E⃑(r⃑)-iω~B⃑(r⃑)=0,
where ε⃡r represents the relative permittivity tensor in different regions shown in Figure 1, μr means the permeability in different regions, ω~ represents ω/c (angular frequency over velocity of light in vacuum), P⃑(r⃑) denotes the random electric dipole source distribution in space, and M⃑(r⃑) is for random magnetic dipole source distribution.Figure 1
Schematic illustration of the geometry of the system under considering shows the thickness (h) 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 aboutP⃑ and M⃑), and one has
(3)∇×(∇×E⃑(r⃑))-ω~2μrε⃡r·E⃑(r⃑)=0.
Taking into account the plane wave ansatz, we can have the matrix equation determining the eigenmodes for given regions:
(4)(ω~2εxμr-ky2-kz2kxkykxkzkxkyω~2εyμr-kx2-kz2kykzkxkzkykzω~2εzμr-ky2-kx2)×(ExEyEz)=0,
where εx=ε1,εy=ε2,εz=ε3 in region II, εx=εy=εz=1 in regions I and III, μr=μ in region II, and μr=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
(5)ω~4-ω~2(kx2+ky2εzμr+kx2+kz2εyμr+kz2+ky2εxμr)+(kx2εyεzμr2+ky2εxεzμr2+kz2εxεyμr2)×(kx2+kz2+ky2)=0.
From (5) kz(ω~,kx,ky) can be solved out as the function of ω~, kx and ky, and there are 4 solutions for kz, namely, {kz,1,-kz,1;kz,2,-kz,2} (kz,1 and kz,2 are positive real or positive imaginary numbers to denote the wave propagating or decaying along the positive z 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 kz,i2-[kz,i(ω~,kx,ky)]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 e^1,± and e^2,±, and conveniently the corresponding unit vectors of electric displacement vectors and magnetic field vectors are gotten by d^i,±=ε⃡r·e^i,±/|ε⃡r·e^i,±| and h^i,±=k⃑±×e^i,±/|k⃑±×e^i,±|. 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 k⃑ into two components along d^1 and d^2, respectively. These two components obey different but definite dispersions, once the three parameters (ω~, kx, and ky) 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ω~, kx, ky, and kz 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):
(6)P⃑(r⃑)=∑i=1,2∫dz0e∓ikz,iz0∫dk⃑(2π)3P⃑k⃑ieik⃑i·r⃑=∫dz0∫dκ⃑(2π)2δ(z-z0)P⃑κ⃑eiκ⃑·R⃑,M⃑(r⃑)=∑i=1,2∫dz0e±ikz,iz0∫dk⃑(2π)3M⃑k⃑ieik⃑i·r⃑=∫dz0∫dκ⃑(2π)2δ(z-z0)M⃑κ⃑eiκ⃑·R⃑,
where z0∈[0,h], since there should be no random sources in vacuum. Here and below it is convenient to put R⃑=(x,y). 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 z 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(7)E⃑(R⃑,z)=∫dκ⃑(2π)2E⃑(κ⃑,z)exp(iκ⃑·R⃑),B⃑(R⃑,z)=∫dκ⃑(2π)2B⃑(κ⃑,z)exp(iκ⃑·R⃑),
where
(8)E⃑(κ⃑,z)=∫dz0G⃡EP(κ⃑;z-z0)P⃑(κ⃑,z0)+∫dz0G⃡EM(κ⃑;z-z0)M⃑(κ⃑,z0),B⃑(κ⃑,z)=∫dz0G⃡BP(κ⃑;z-z0)P⃑(κ⃑,z0)+∫dz0G⃡BM(κ⃑;z-z0)M⃑(κ⃑,z0).
In (8), G⃡EP, G⃡BP, G⃡BM, and G⃡BM 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 z0, and to account all modes one needs only to integrate all κ⃑ and z dimensions of the anisotropic medium.


### 2.3. Mode Expansion in Each Region

To determine the field generated from source plate locating atz0 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
(9)B⃑II(ω~,κ⃑)=[B⃑1,g,+exp(-ikz,1z0)+B⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+B⃑1,b,++B⃑2,b,++[B⃑1,g,-exp(ikz,1z0)+B⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+B⃑1,b,-+B⃑2,b,-+B⃑δ(z-z0)exp(iκ⃑·R⃑),D⃑II(ω~,κ⃑)=[D⃑1,g,+exp(-ikz,1z0)+D⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+D⃑1,b,++D⃑2,b,++[D⃑1,g,-exp(ikz,1z0)+D⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+D⃑1,b,-+D⃑2,b,-+D⃑δ(z-z0)exp(iκ⃑·R⃑),E⃑II(ω~,κ⃑)=[E⃑1,g,+exp(-ikz,1z0)+E⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+E⃑1,b,++E⃑2,b,++[E⃑1,g,-exp(ikz,1z0)+E⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+E⃑1,b,-+E⃑2,b,-+E⃑δ(z-z0)exp(iκ⃑·R⃑),
where
(10)B⃑i,g,±=Bi,g,±exp(iκ⃑·R⃑)exp(±ikz,iz)h^i,±,B⃑i,b,±=Bi,b,±exp(iκ⃑·R⃑)exp(±ikz,iz)h^i,±,D⃑i,g,±=-1iω~μ∇×B⃑i,g,±,E⃑i,g,±=ε⃡-1·D⃑i,g,±D⃑i,b,±=-1iω~μ∇×B⃑i,b,±,E⃑i,b,±=ε⃡-1·D⃑i,b,±.i=1,2,B⃑(D⃑, E⃑) is to account the local (nonpropagating) field locating on the source plane, and Θ(z-z0) 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
(11)B⃑I(ω~,κ⃑)=(Bs,-h^s,-I+Bp,-h^p,-I)B⃑I(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(-ikz,Iz),D⃑I(ω~,κ⃑)=E⃑I(ω~,κ⃑)D⃑I(ω~,κ⃑)=(Bs,-e^s,-I+Bp,-e^p,-I)D⃑I(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(-ikz,Iz),B⃑III(ω~,κ⃑)=(Bs,+h^s,+III+Bp,+h^p,+III)B⃑III(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(ikz,IIIz),D⃑III(ω~,κ⃑)=E⃑III(ω~,κ⃑)D⃑III(ω~,κ⃑)=(Bs,+e^s,+III+Bp,+e^p,+III)D⃑III(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(ikz,IIIz),
where the unit vectors are defined as
(12)e^s,±α=k⃑s,±×z^|k⃑s,±×z^|=(κcosθκsinθ±kz)×(001)=(sinθ-cosθ0),h^s,±α=k⃑s,±×e^s,±|k⃑s,±×e^s,±|=1ω~εrμr(±kzcosθ±kzsinθ-κ),h^p,±α=k⃑p,±×z^|k⃑p,±×z^|=(κcosθκsinθ±kz)×(001)=(sinθ-cosθ0),e^p,±α=-k⃑p,±×h^p,±|k⃑p,±×h^p,±|=1ω~εrμr(∓kzcosθ∓kzsinθκ),
where α=I,  II,  III denotes the region where the unit vectors belong.Figure 2
Schematic illustration of eigenmode expansion in three regions:B⃑s,- and B⃑p,- denote the eigen-plane-wave modes in region I, and the minus sign is to denote the direction of propagating is along z− direction; B⃑1,g,± in region II denotes the first eigenmode generated from the source plane, B⃑2,g,± is to denote the second eigenmode generated from the source plane, B⃑1,b,± in region II expresses the first eigenmode bouncing back at the upper and lower boundaries, and B⃑2,b,± gives the second eigenmode bouncing back at the upper and lower boundaries (the plus and minus signs denotes the propagating along z+ and z− directions); B⃑s,+ and B⃑p,+ denote the eigen-plane wave modes in region III, and the minus sign is to denote the direction of propagating is along z+ 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
(13)δ(z-z0)z^×[B⃑1,g,+exp(-ikz,1z0)+B⃑2,g,+exp(-ikz,2z0)-B⃑1,g,-exp(ikz,1z0)-B⃑2,g,-exp(ikz,2z0)]+z^×B⃑δ′(z-z0)exp(iκ⃑·R⃑)+iκ⃑×B⃑exp(iκ⃑·R⃑)δ(z-z0)+iω~με⃡·E⃑δ(z-z0)exp(iκ⃑·R⃑)=-4πiω~μP⃑δ(z-z0)exp(iκ⃑·R⃑)+4πμiκ⃑×M⃑exp(iκ⃑·R⃑)δ(z-z0)+4πμδ′(z-z0)z^×M⃑exp(iκ⃑·R⃑),(14)δ(z-z0)z^×[E⃑1,g,+exp(-ikz,1z0)-E⃑1,g,-exp(ikz,1z0)+E⃑2,g,+exp(-ikz,2z0)-E⃑2,g,-exp(ikz,2z0)]+z^×E⃑δ′(z-z0)exp(iκ⃑·R⃑)+iκ⃑×E⃑exp(iκ⃑·R⃑)δ(z-z0)=iω~B⃑δ(z-z0)exp(iκ⃑·R⃑),
where we have used ∇Θ(z-z0)=z^δ(z-z0) and ∇δ(z-z0)=z^δ'(z-z0). By requiring the coefficients in front of δ(z-z0) and δ'(z-z0) equal for the two sides, the amplitudes Bi,g,± 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 (z=0 and z=h):
(15)E⃑II(ω~,κ⃑;z=0)·(x^x^+y^y^)=(Bs,-e^s,-+Bp,-e^p,-)·(x^x^+y^y^)exp(iκ⃑·R⃑),B⃑II(ω~,κ⃑;z=0)·(x^x^+y^y^)=(Bs,-h^s,-+Bp,-h^p,-)·(x^x^+y^y^)exp(iκ⃑·R⃑),
for z=0 boundary, where x^ and y^ are the unit vectors of x and y directions,
(16)E⃑II(ω~,κ⃑;z=h)·(x^x^+y^y^)=(Bs,+e^s,++Bp,+e^p,+)·(x^x^+y^y^)exp(iκ⃑·R⃑)exp(ikh),B⃑II(ω~,κ⃑;z=h)·(x^x^+y^y^)=(Bs,+h^s,++Bp,+h^p,+)·(x^x^+y^y^)exp(iκ⃑·R⃑)exp(ikh),
for z=h 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 s- and p-polarized modes. The dispersions for s- and p-polarized modes are the same and given as ω~2εμ=κ2+kz2 in region II and ω~2=κ2+kz2 in region I and III. The unit vectors of electric field and magnetic field are given in (12). The further convenience is that e^s±, κ^ and z^ 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 δ'(z-z0) in (13) and (14), the following relations are found:
(17)z^×B⃑=4πμz^×M⃑,z^×E⃑=0.
The local magnetic source relates with local magnetic field as Be^s±=4πμMe^s± and Bκ^=4πμMκ^ and e^s,± and κ^ components of local electric field vanish (Ee^s±=Eκ^=0). From the equality of the coefficients in front of δ(z-z0), we can get the other equations, and we write them according to the components along e^s±, κ^, and z^ directions for (13), and (14) respectively,
(18)1εμ(Bs,g,++Bs,g,-)=4πiω~2μkzPe^s,±-4πμiκω~kzMz^,Bp,g,+-Bp,g,-=-4πiω~μPκ^,iωεμEz^=-i4πω~μPz^,1με(Bp,g,++Bp,g,-)+κkIIkz,II=iω~kIIkz,IIBe^s,±,1με(Bs,g,+-Bs,g,-)=iω~Bκ^,Bz^=0.
The fields launching from the source can be gotten by solving (18) and the results are collected as follows:
(19)Bs,g,+=(2πiω~2μkz,IIμε(e^s,±II·P⃑)+2πμiω~2μεkz,II(h^s,+II·M⃑)),Bs,g,-=(2πiω~2μkz,IIμε(e^s,±II·P⃑)+2πμiω~2μεkz,II(h^s,-II·M⃑)),Bp,g,+=(2πiμω~kIIkz,IIεμ(e^s,±II·M⃑)-2πiω~kIIμkz,II(h^s,+II·P⃑)),Bp,g,-=(2πiμω~kIIkz,IIεμ(e^s,±II·M⃑)-2πiω~kIIμkz,II(h^s,-II·P⃑)),
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 e^s,±T=h^p,±T and h^s,±T=-e^p,±T always holds, it is not necessary to arrange the equations in x-y-z coordinate. Taking the advantage, the boundary equation (15) for z=0 becomes
(20)Bs,-e^s,-T,I=1με(Bs,b,+e^s,+T,II+Bs,g,-exp(ikz,IIz0)e^s,-T,II+Bs,b,-e^s,-T,II),Bp,-e^p,-T,I=1με(Bp,b,+e^p,+T,II+Bp,g,-exp(ikz,IIz0)e^p,-T,II+Bp,b,-e^p,-T,II),Bp,-h^p,-T,I=Bp,b,+h^p,+T,II+Bp,g,-exp(ikz,IIz0)h^p,-T,II+Bp,b,-h^p,-T,II,Bs,-h^s,-T,I=Bs,b,+h^s,+T,II+Bs,g,-exp(ikz,IIz0)h^s,-T,II+Bs,b,-h^s,-T,II.
Similarly, the boundary equation (16) for z=h gives
(21)Bs,+exp(ikz,Ih)e^s,+T,III=1με(Bs,g,+exp[ikz,II(h-z0)]e^s,+T,II+Bs,b,+exp(ikz,IIh)e^s,+T,II+Bs,b,-exp(-ikz,IIh)e^s,-T,II),Bp,+exp(ikz,Ih)e^p,+T,III=1με(Bp,g,+exp[ikz,II(h-z0)]e^p,+T,II+Bp,b,+exp(ikz,IIh)e^p,+T,II+Bp,b,-exp(-ikz,IIh)e^p,-T,II),Bp,+exp(ikz,Ih)h^p,+T,III=Bp,g,+exp[ikz,II(h-z0)]h^p,+T,II+Bp,b,+exp(ikz,IIh)h^p,+T,II+Bp,b,-exp(-ikz,IIh)h^p,-T,II,Bs,+exp(ikz,Ih)h^s,+T,III=Bs,g,+exp[ikz,II(h-z0)]h^s,+T,II+Bs,b,+exp(ikz,IIh)h^s,+T,II+Bs,b,-exp(-ikz,IIh)h^s,-T,II,
where e^s,±T,I(III),e^s,±T,II,e^p,±T,I(III),  e^p,±T,II, h^p,±T,I(III),h^p,±T,II,h^s,±T,I(III),h^s,±T,II are the tangential parts of these unit vectors in different regions and kz,II=kz,s=kz,p. Once Bs,- and Bp,-(Bs,+ and Bp,+) 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:
(22)B⃑I(ω~,κ⃑)=1με[rII,Iseikz,IIz01-tII,Is2ei2kz,IIhBs,g,-+rII,IstII,Iseikz,II(2h-z0)1-tII,Is2ei2kz,IIhBs,g,+]×h^s,-exp(iκ⃑·R⃑)exp(-ikz,Iz)+1με[rII,Ipeikz,IIz01-tII,Ip2ei2kz,IIhBp,g,-+rII,IptII,Ipeikz,II(2h-z0)1-tII,Ip2ei2kz,IIhBp,g,+]×h^p,-exp(iκ⃑·R⃑)exp(-ikz,Iz),E⃑I(ω~,κ⃑)=1με[rII,Iseikz,IIz01-tII,Is2ei2kz,IIhBs,g,-+rII,IstII,Iseikz,II(2h-z0)1-tII,Is2ei2kz,IIhBs,g,+]×e^s,-exp(iκ⃑·R⃑)exp(-ikz,Iz)+1με[rII,Ipeikz,IIz01-tII,Ip2ei2kz,IIhBp,g,-+rII,IptII,Ipeikz,II(2h-z0)1-tII,Ip2ei2kz,IIhBp,g,+]×e^p,-exp(iκ⃑·R⃑)exp(-ikz,Iz).
The fields generated in region III are gotten as
(23)B⃑III(ω~,κ⃑)=1εμ[rII,IIIstII,IIIseikz,II(z0+h)1-tII,IIIs2ei2kz,IIhBs,g,-+rII,IIIseikz,II(h-z0)1-tII,IIIs2ei2kz,IIhBs,g,+]×h^s,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz)+1εμ[rII,IIIptII,IIIpeikz,II(z0+h)1-tII,IIIp2ei2kz,IIhBp,g,-+rII,IIIpeikz,II(h-z0)1-tII,IIIp2ei2kz,IIhBp,g,+]×h^p,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz),E⃑III(ω~,κ⃑)=1εμ[rII,IIIstII,IIIseikz,II(z0+h)1-tII,IIIs2ei2kz,IIhBs,g,-+rII,IIIseikz,II(h-z0)1-tII,IIIs2ei2kz,IIhBs,g,+]×e^s,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz)+1εμ[rII,IIIptII,IIIpeikz,II(z0+h)1-tII,IIIp2ei2kz,IIhBp,g,-+rII,IIIpeikz,II(h-z0)1-tII,IIIp2ei2kz,IIhBp,g,+]×e^p,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz),
where ti,js=(kz,i-kz,j)/(kz,i+kz,j), ti,js=2kz,i/(kz,i+kz,j), ti,jp=(kz,iεjμj-kz,jεiμi)/(kz,iεjμj+kz,jεiμi)  ,ri,jp=2kz,iεjμjεiμi/(kz,iεjμj+kz,jεiμi) are Fresnel coefficients. Since Bs,g,± and Bp,g,± 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):
(24)G⃡BMI(κ⃑;z-z0)=2πμiω~kIIkz,II(tII,Is(rII,Iseikz,II(2h-z0)h^s,-Ih^s,+II+eikz,IIz0h^s,-Ih^s,-II)1-rII,Is2ei2kz,IIh+tII,Ip(eikz,IIz0+rII,Ipeikz,II(2h-z0))1-rII,Ip2ei2kz,IIhh^p,-Ie^s,±II)×eiκ⃑·R⃑e-ikz,Iz,G⃡BPI(κ⃑;z-z0)=2πiω~2μkz,II×(tII,Is(eikz,IIz0+rII,Iseikz,II(2h-z0))1-rII,Is2ei2kz,IIhh^s,-Ie^s,±II-tII,Ip(eikz,IIz0h^p,-Ih^s,-II+rII,Ipeikz,II(2h-z0)h^p,-Ih^s,+II)1-rII,Ip2ei2kz,IIh)×eiκ⃑·R⃑e-ikz,Iz,G⃡EMI(κ⃑;z-z0)=2πμiω~kIIkz,II×[rII,Is(tII,Iseikz,II(2h-z0)e^s,-Ih^s,+II+eikz,IIz0e^s,-Ih^s,-II)1-tII,Is2ei2kz,IIh+rII,Ip(eikz,IIz0+tII,Ipeikz,II(2h-z0))1-tII,Ip2ei2kz,IIhe^p,-Ie^s,±II]×eiκ⃑·R⃑e-ikz,Iz,G⃡EMI(κ⃑;z-z0)=2πiω~2μkz,II[tII,Is(eikz,IIz0+rII,Iseikz,II(2h-z0))1-rII,Is2ei2kz,IIhe^s,-Ie^s,±II-tII,Ip(eikz,IIz0e^p,-Ih^s,-II+rII,Ipeikz,II(2h-z0)e^p,-Ih^s,+II)1-rII,Ip2ei2kz,IIh]×eiκ⃑·R⃑e-ikz,Iz,G⃡BMIII(κ⃑;z-z0)=2πμiω~kIIkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)h^s,+IIIh^s,-II+eikz,II(h-z0)h^s,+IIIh^s,+II)1-rII,IIIs2ei2kz,IIh+tII,IIIp(rII,IIIpeikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIp2ei2kz,IIhh^p,+IIIe^s,±II]×eiκz,III(z-h)eik⃑·R⃑,G⃡BPIII(κ⃑;z-z0)=2πiω~2μkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIs2ei2kz,IIhh^s,+IIIe^s,±II-tII,IIIp(rII,IIIpeikz,II(z0+h)h^p,+IIIh^s,-II+eikz,II(h-z0)h^p,+IIIh^s,+II)1-rII,IIIp2ei2kz,IIh]×eiκz,III(z-h)eik⃑·R⃑G⃡EMIII(κ⃑z-z0)=2πμiω~kIIkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)e^s,+IIIh^s,-II+eikz,II(h-z0)e^s,+IIIh^s,+II)1-rII,IIIs2ei2kz,IIh+tII,IIIp(rII,IIIpeikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIp2ei2kz,IIhe^p,+IIIe^s,±II]×eiκz,III(z-h)eik⃑·R⃑,G⃡EPIII(κ⃑z-z0)=2πiω~2μkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIs2ei2kz,IIhe^s,+IIIe^s,±II-tII,IIIp(rII,IIIpeikz,II(z0+h)e^p,+IIIh^s,-II+eikz,II(h-z0)e^p,+IIIh^s,+II)1-rII,IIIp2ei2kz,IIh]×eiκz,III(z-h)eik⃑·R⃑.
These results exactly recover (4.18) in [17], and thus the correction of current formulism is proved.


## 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 aboutP⃑ and M⃑), and one has
(3)∇×(∇×E⃑(r⃑))-ω~2μrε⃡r·E⃑(r⃑)=0.
Taking into account the plane wave ansatz, we can have the matrix equation determining the eigenmodes for given regions:
(4)(ω~2εxμr-ky2-kz2kxkykxkzkxkyω~2εyμr-kx2-kz2kykzkxkzkykzω~2εzμr-ky2-kx2)×(ExEyEz)=0,
where εx=ε1,εy=ε2,εz=ε3 in region II, εx=εy=εz=1 in regions I and III, μr=μ in region II, and μr=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
(5)ω~4-ω~2(kx2+ky2εzμr+kx2+kz2εyμr+kz2+ky2εxμr)+(kx2εyεzμr2+ky2εxεzμr2+kz2εxεyμr2)×(kx2+kz2+ky2)=0.
From (5) kz(ω~,kx,ky) can be solved out as the function of ω~, kx and ky, and there are 4 solutions for kz, namely, {kz,1,-kz,1;kz,2,-kz,2} (kz,1 and kz,2 are positive real or positive imaginary numbers to denote the wave propagating or decaying along the positive z 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 kz,i2-[kz,i(ω~,kx,ky)]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 e^1,± and e^2,±, and conveniently the corresponding unit vectors of electric displacement vectors and magnetic field vectors are gotten by d^i,±=ε⃡r·e^i,±/|ε⃡r·e^i,±| and h^i,±=k⃑±×e^i,±/|k⃑±×e^i,±|. 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 k⃑ into two components along d^1 and d^2, respectively. These two components obey different but definite dispersions, once the three parameters (ω~, kx, and ky) 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ω~, kx, ky, and kz 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):
(6)P⃑(r⃑)=∑i=1,2∫dz0e∓ikz,iz0∫dk⃑(2π)3P⃑k⃑ieik⃑i·r⃑=∫dz0∫dκ⃑(2π)2δ(z-z0)P⃑κ⃑eiκ⃑·R⃑,M⃑(r⃑)=∑i=1,2∫dz0e±ikz,iz0∫dk⃑(2π)3M⃑k⃑ieik⃑i·r⃑=∫dz0∫dκ⃑(2π)2δ(z-z0)M⃑κ⃑eiκ⃑·R⃑,
where z0∈[0,h], since there should be no random sources in vacuum. Here and below it is convenient to put R⃑=(x,y). 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 z 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(7)E⃑(R⃑,z)=∫dκ⃑(2π)2E⃑(κ⃑,z)exp(iκ⃑·R⃑),B⃑(R⃑,z)=∫dκ⃑(2π)2B⃑(κ⃑,z)exp(iκ⃑·R⃑),
where
(8)E⃑(κ⃑,z)=∫dz0G⃡EP(κ⃑;z-z0)P⃑(κ⃑,z0)+∫dz0G⃡EM(κ⃑;z-z0)M⃑(κ⃑,z0),B⃑(κ⃑,z)=∫dz0G⃡BP(κ⃑;z-z0)P⃑(κ⃑,z0)+∫dz0G⃡BM(κ⃑;z-z0)M⃑(κ⃑,z0).
In (8), G⃡EP, G⃡BP, G⃡BM, and G⃡BM 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 z0, and to account all modes one needs only to integrate all κ⃑ and z dimensions of the anisotropic medium.


## 2.3. Mode Expansion in Each Region

To determine the field generated from source plate locating atz0 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
(9)B⃑II(ω~,κ⃑)=[B⃑1,g,+exp(-ikz,1z0)+B⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+B⃑1,b,++B⃑2,b,++[B⃑1,g,-exp(ikz,1z0)+B⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+B⃑1,b,-+B⃑2,b,-+B⃑δ(z-z0)exp(iκ⃑·R⃑),D⃑II(ω~,κ⃑)=[D⃑1,g,+exp(-ikz,1z0)+D⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+D⃑1,b,++D⃑2,b,++[D⃑1,g,-exp(ikz,1z0)+D⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+D⃑1,b,-+D⃑2,b,-+D⃑δ(z-z0)exp(iκ⃑·R⃑),E⃑II(ω~,κ⃑)=[E⃑1,g,+exp(-ikz,1z0)+E⃑2,g,+exp(-ikz,2z0)]×Θ(z-z0)+E⃑1,b,++E⃑2,b,++[E⃑1,g,-exp(ikz,1z0)+E⃑2,g,-exp(ikz,2z0)]×Θ(z0-z)+E⃑1,b,-+E⃑2,b,-+E⃑δ(z-z0)exp(iκ⃑·R⃑),
where
(10)B⃑i,g,±=Bi,g,±exp(iκ⃑·R⃑)exp(±ikz,iz)h^i,±,B⃑i,b,±=Bi,b,±exp(iκ⃑·R⃑)exp(±ikz,iz)h^i,±,D⃑i,g,±=-1iω~μ∇×B⃑i,g,±,E⃑i,g,±=ε⃡-1·D⃑i,g,±D⃑i,b,±=-1iω~μ∇×B⃑i,b,±,E⃑i,b,±=ε⃡-1·D⃑i,b,±.i=1,2,B⃑(D⃑, E⃑) is to account the local (nonpropagating) field locating on the source plane, and Θ(z-z0) 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
(11)B⃑I(ω~,κ⃑)=(Bs,-h^s,-I+Bp,-h^p,-I)B⃑I(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(-ikz,Iz),D⃑I(ω~,κ⃑)=E⃑I(ω~,κ⃑)D⃑I(ω~,κ⃑)=(Bs,-e^s,-I+Bp,-e^p,-I)D⃑I(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(-ikz,Iz),B⃑III(ω~,κ⃑)=(Bs,+h^s,+III+Bp,+h^p,+III)B⃑III(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(ikz,IIIz),D⃑III(ω~,κ⃑)=E⃑III(ω~,κ⃑)D⃑III(ω~,κ⃑)=(Bs,+e^s,+III+Bp,+e^p,+III)D⃑III(ω~,κ⃑)=×exp(iκ⃑·R⃑)exp(ikz,IIIz),
where the unit vectors are defined as
(12)e^s,±α=k⃑s,±×z^|k⃑s,±×z^|=(κcosθκsinθ±kz)×(001)=(sinθ-cosθ0),h^s,±α=k⃑s,±×e^s,±|k⃑s,±×e^s,±|=1ω~εrμr(±kzcosθ±kzsinθ-κ),h^p,±α=k⃑p,±×z^|k⃑p,±×z^|=(κcosθκsinθ±kz)×(001)=(sinθ-cosθ0),e^p,±α=-k⃑p,±×h^p,±|k⃑p,±×h^p,±|=1ω~εrμr(∓kzcosθ∓kzsinθκ),
where α=I,  II,  III denotes the region where the unit vectors belong.Figure 2
Schematic illustration of eigenmode expansion in three regions:B⃑s,- and B⃑p,- denote the eigen-plane-wave modes in region I, and the minus sign is to denote the direction of propagating is along z− direction; B⃑1,g,± in region II denotes the first eigenmode generated from the source plane, B⃑2,g,± is to denote the second eigenmode generated from the source plane, B⃑1,b,± in region II expresses the first eigenmode bouncing back at the upper and lower boundaries, and B⃑2,b,± gives the second eigenmode bouncing back at the upper and lower boundaries (the plus and minus signs denotes the propagating along z+ and z− directions); B⃑s,+ and B⃑p,+ denote the eigen-plane wave modes in region III, and the minus sign is to denote the direction of propagating is along z+ 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
(13)δ(z-z0)z^×[B⃑1,g,+exp(-ikz,1z0)+B⃑2,g,+exp(-ikz,2z0)-B⃑1,g,-exp(ikz,1z0)-B⃑2,g,-exp(ikz,2z0)]+z^×B⃑δ′(z-z0)exp(iκ⃑·R⃑)+iκ⃑×B⃑exp(iκ⃑·R⃑)δ(z-z0)+iω~με⃡·E⃑δ(z-z0)exp(iκ⃑·R⃑)=-4πiω~μP⃑δ(z-z0)exp(iκ⃑·R⃑)+4πμiκ⃑×M⃑exp(iκ⃑·R⃑)δ(z-z0)+4πμδ′(z-z0)z^×M⃑exp(iκ⃑·R⃑),(14)δ(z-z0)z^×[E⃑1,g,+exp(-ikz,1z0)-E⃑1,g,-exp(ikz,1z0)+E⃑2,g,+exp(-ikz,2z0)-E⃑2,g,-exp(ikz,2z0)]+z^×E⃑δ′(z-z0)exp(iκ⃑·R⃑)+iκ⃑×E⃑exp(iκ⃑·R⃑)δ(z-z0)=iω~B⃑δ(z-z0)exp(iκ⃑·R⃑),
where we have used ∇Θ(z-z0)=z^δ(z-z0) and ∇δ(z-z0)=z^δ'(z-z0). By requiring the coefficients in front of δ(z-z0) and δ'(z-z0) equal for the two sides, the amplitudes Bi,g,± 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 (z=0 and z=h):
(15)E⃑II(ω~,κ⃑;z=0)·(x^x^+y^y^)=(Bs,-e^s,-+Bp,-e^p,-)·(x^x^+y^y^)exp(iκ⃑·R⃑),B⃑II(ω~,κ⃑;z=0)·(x^x^+y^y^)=(Bs,-h^s,-+Bp,-h^p,-)·(x^x^+y^y^)exp(iκ⃑·R⃑),
for z=0 boundary, where x^ and y^ are the unit vectors of x and y directions,
(16)E⃑II(ω~,κ⃑;z=h)·(x^x^+y^y^)=(Bs,+e^s,++Bp,+e^p,+)·(x^x^+y^y^)exp(iκ⃑·R⃑)exp(ikh),B⃑II(ω~,κ⃑;z=h)·(x^x^+y^y^)=(Bs,+h^s,++Bp,+h^p,+)·(x^x^+y^y^)exp(iκ⃑·R⃑)exp(ikh),
for z=h 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 s- and p-polarized modes. The dispersions for s- and p-polarized modes are the same and given as ω~2εμ=κ2+kz2 in region II and ω~2=κ2+kz2 in region I and III. The unit vectors of electric field and magnetic field are given in (12). The further convenience is that e^s±, κ^ and z^ 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 δ'(z-z0) in (13) and (14), the following relations are found:
(17)z^×B⃑=4πμz^×M⃑,z^×E⃑=0.
The local magnetic source relates with local magnetic field as Be^s±=4πμMe^s± and Bκ^=4πμMκ^ and e^s,± and κ^ components of local electric field vanish (Ee^s±=Eκ^=0). From the equality of the coefficients in front of δ(z-z0), we can get the other equations, and we write them according to the components along e^s±, κ^, and z^ directions for (13), and (14) respectively,
(18)1εμ(Bs,g,++Bs,g,-)=4πiω~2μkzPe^s,±-4πμiκω~kzMz^,Bp,g,+-Bp,g,-=-4πiω~μPκ^,iωεμEz^=-i4πω~μPz^,1με(Bp,g,++Bp,g,-)+κkIIkz,II=iω~kIIkz,IIBe^s,±,1με(Bs,g,+-Bs,g,-)=iω~Bκ^,Bz^=0.
The fields launching from the source can be gotten by solving (18) and the results are collected as follows:
(19)Bs,g,+=(2πiω~2μkz,IIμε(e^s,±II·P⃑)+2πμiω~2μεkz,II(h^s,+II·M⃑)),Bs,g,-=(2πiω~2μkz,IIμε(e^s,±II·P⃑)+2πμiω~2μεkz,II(h^s,-II·M⃑)),Bp,g,+=(2πiμω~kIIkz,IIεμ(e^s,±II·M⃑)-2πiω~kIIμkz,II(h^s,+II·P⃑)),Bp,g,-=(2πiμω~kIIkz,IIεμ(e^s,±II·M⃑)-2πiω~kIIμkz,II(h^s,-II·P⃑)),
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 e^s,±T=h^p,±T and h^s,±T=-e^p,±T always holds, it is not necessary to arrange the equations in x-y-z coordinate. Taking the advantage, the boundary equation (15) for z=0 becomes
(20)Bs,-e^s,-T,I=1με(Bs,b,+e^s,+T,II+Bs,g,-exp(ikz,IIz0)e^s,-T,II+Bs,b,-e^s,-T,II),Bp,-e^p,-T,I=1με(Bp,b,+e^p,+T,II+Bp,g,-exp(ikz,IIz0)e^p,-T,II+Bp,b,-e^p,-T,II),Bp,-h^p,-T,I=Bp,b,+h^p,+T,II+Bp,g,-exp(ikz,IIz0)h^p,-T,II+Bp,b,-h^p,-T,II,Bs,-h^s,-T,I=Bs,b,+h^s,+T,II+Bs,g,-exp(ikz,IIz0)h^s,-T,II+Bs,b,-h^s,-T,II.
Similarly, the boundary equation (16) for z=h gives
(21)Bs,+exp(ikz,Ih)e^s,+T,III=1με(Bs,g,+exp[ikz,II(h-z0)]e^s,+T,II+Bs,b,+exp(ikz,IIh)e^s,+T,II+Bs,b,-exp(-ikz,IIh)e^s,-T,II),Bp,+exp(ikz,Ih)e^p,+T,III=1με(Bp,g,+exp[ikz,II(h-z0)]e^p,+T,II+Bp,b,+exp(ikz,IIh)e^p,+T,II+Bp,b,-exp(-ikz,IIh)e^p,-T,II),Bp,+exp(ikz,Ih)h^p,+T,III=Bp,g,+exp[ikz,II(h-z0)]h^p,+T,II+Bp,b,+exp(ikz,IIh)h^p,+T,II+Bp,b,-exp(-ikz,IIh)h^p,-T,II,Bs,+exp(ikz,Ih)h^s,+T,III=Bs,g,+exp[ikz,II(h-z0)]h^s,+T,II+Bs,b,+exp(ikz,IIh)h^s,+T,II+Bs,b,-exp(-ikz,IIh)h^s,-T,II,
where e^s,±T,I(III),e^s,±T,II,e^p,±T,I(III),  e^p,±T,II, h^p,±T,I(III),h^p,±T,II,h^s,±T,I(III),h^s,±T,II are the tangential parts of these unit vectors in different regions and kz,II=kz,s=kz,p. Once Bs,- and Bp,-(Bs,+ and Bp,+) 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:
(22)B⃑I(ω~,κ⃑)=1με[rII,Iseikz,IIz01-tII,Is2ei2kz,IIhBs,g,-+rII,IstII,Iseikz,II(2h-z0)1-tII,Is2ei2kz,IIhBs,g,+]×h^s,-exp(iκ⃑·R⃑)exp(-ikz,Iz)+1με[rII,Ipeikz,IIz01-tII,Ip2ei2kz,IIhBp,g,-+rII,IptII,Ipeikz,II(2h-z0)1-tII,Ip2ei2kz,IIhBp,g,+]×h^p,-exp(iκ⃑·R⃑)exp(-ikz,Iz),E⃑I(ω~,κ⃑)=1με[rII,Iseikz,IIz01-tII,Is2ei2kz,IIhBs,g,-+rII,IstII,Iseikz,II(2h-z0)1-tII,Is2ei2kz,IIhBs,g,+]×e^s,-exp(iκ⃑·R⃑)exp(-ikz,Iz)+1με[rII,Ipeikz,IIz01-tII,Ip2ei2kz,IIhBp,g,-+rII,IptII,Ipeikz,II(2h-z0)1-tII,Ip2ei2kz,IIhBp,g,+]×e^p,-exp(iκ⃑·R⃑)exp(-ikz,Iz).
The fields generated in region III are gotten as
(23)B⃑III(ω~,κ⃑)=1εμ[rII,IIIstII,IIIseikz,II(z0+h)1-tII,IIIs2ei2kz,IIhBs,g,-+rII,IIIseikz,II(h-z0)1-tII,IIIs2ei2kz,IIhBs,g,+]×h^s,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz)+1εμ[rII,IIIptII,IIIpeikz,II(z0+h)1-tII,IIIp2ei2kz,IIhBp,g,-+rII,IIIpeikz,II(h-z0)1-tII,IIIp2ei2kz,IIhBp,g,+]×h^p,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz),E⃑III(ω~,κ⃑)=1εμ[rII,IIIstII,IIIseikz,II(z0+h)1-tII,IIIs2ei2kz,IIhBs,g,-+rII,IIIseikz,II(h-z0)1-tII,IIIs2ei2kz,IIhBs,g,+]×e^s,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz)+1εμ[rII,IIIptII,IIIpeikz,II(z0+h)1-tII,IIIp2ei2kz,IIhBp,g,-+rII,IIIpeikz,II(h-z0)1-tII,IIIp2ei2kz,IIhBp,g,+]×e^p,+e-ikz,Ihexp(iκ⃑·R⃑)exp(ikz,IIIz),
where ti,js=(kz,i-kz,j)/(kz,i+kz,j), ti,js=2kz,i/(kz,i+kz,j), ti,jp=(kz,iεjμj-kz,jεiμi)/(kz,iεjμj+kz,jεiμi)  ,ri,jp=2kz,iεjμjεiμi/(kz,iεjμj+kz,jεiμi) are Fresnel coefficients. Since Bs,g,± and Bp,g,± 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):
(24)G⃡BMI(κ⃑;z-z0)=2πμiω~kIIkz,II(tII,Is(rII,Iseikz,II(2h-z0)h^s,-Ih^s,+II+eikz,IIz0h^s,-Ih^s,-II)1-rII,Is2ei2kz,IIh+tII,Ip(eikz,IIz0+rII,Ipeikz,II(2h-z0))1-rII,Ip2ei2kz,IIhh^p,-Ie^s,±II)×eiκ⃑·R⃑e-ikz,Iz,G⃡BPI(κ⃑;z-z0)=2πiω~2μkz,II×(tII,Is(eikz,IIz0+rII,Iseikz,II(2h-z0))1-rII,Is2ei2kz,IIhh^s,-Ie^s,±II-tII,Ip(eikz,IIz0h^p,-Ih^s,-II+rII,Ipeikz,II(2h-z0)h^p,-Ih^s,+II)1-rII,Ip2ei2kz,IIh)×eiκ⃑·R⃑e-ikz,Iz,G⃡EMI(κ⃑;z-z0)=2πμiω~kIIkz,II×[rII,Is(tII,Iseikz,II(2h-z0)e^s,-Ih^s,+II+eikz,IIz0e^s,-Ih^s,-II)1-tII,Is2ei2kz,IIh+rII,Ip(eikz,IIz0+tII,Ipeikz,II(2h-z0))1-tII,Ip2ei2kz,IIhe^p,-Ie^s,±II]×eiκ⃑·R⃑e-ikz,Iz,G⃡EMI(κ⃑;z-z0)=2πiω~2μkz,II[tII,Is(eikz,IIz0+rII,Iseikz,II(2h-z0))1-rII,Is2ei2kz,IIhe^s,-Ie^s,±II-tII,Ip(eikz,IIz0e^p,-Ih^s,-II+rII,Ipeikz,II(2h-z0)e^p,-Ih^s,+II)1-rII,Ip2ei2kz,IIh]×eiκ⃑·R⃑e-ikz,Iz,G⃡BMIII(κ⃑;z-z0)=2πμiω~kIIkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)h^s,+IIIh^s,-II+eikz,II(h-z0)h^s,+IIIh^s,+II)1-rII,IIIs2ei2kz,IIh+tII,IIIp(rII,IIIpeikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIp2ei2kz,IIhh^p,+IIIe^s,±II]×eiκz,III(z-h)eik⃑·R⃑,G⃡BPIII(κ⃑;z-z0)=2πiω~2μkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIs2ei2kz,IIhh^s,+IIIe^s,±II-tII,IIIp(rII,IIIpeikz,II(z0+h)h^p,+IIIh^s,-II+eikz,II(h-z0)h^p,+IIIh^s,+II)1-rII,IIIp2ei2kz,IIh]×eiκz,III(z-h)eik⃑·R⃑G⃡EMIII(κ⃑z-z0)=2πμiω~kIIkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)e^s,+IIIh^s,-II+eikz,II(h-z0)e^s,+IIIh^s,+II)1-rII,IIIs2ei2kz,IIh+tII,IIIp(rII,IIIpeikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIp2ei2kz,IIhe^p,+IIIe^s,±II]×eiκz,III(z-h)eik⃑·R⃑,G⃡EPIII(κ⃑z-z0)=2πiω~2μkz,II×[tII,IIIs(rII,IIIseikz,II(z0+h)+eikz,II(h-z0))1-rII,IIIs2ei2kz,IIhe^s,+IIIe^s,±II-tII,IIIp(rII,IIIpeikz,II(z0+h)e^p,+IIIh^s,-II+eikz,II(h-z0)e^p,+IIIh^s,+II)1-rII,IIIp2ei2kz,IIh]×eiκz,III(z-h)eik⃑·R⃑.
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 asdiag(εo,εe,εo), where without losing generality we have assumed the anisotropy is along y direction (so that mode conversion still happens at the interface). Then the equation determining eigenmodes (5) becomes
(25)(ω~2-kx2εoμ-ky2εoμ-kz2εoμ)(ω~2-kx2εeμ-ky2εoμ-kz2εeμ)=0.
Thus, there are two dispersions gotten: the mode with dispersion ω~2-kx2/εoμ-ky2/εoμ-kz2/εoμ=0 is known as o-light, and the mode with dispersion ω~2-kx2/εeμ-ky2/εoμ-kz2/εeμ=0 is known as e-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
(26)e^o,±=1κ2cos2θ+kz,o2(±kz,o0-κcosθ),e^e,±=1(κ2cos2θ+kz,e2)εe2+κ2sin2θεo2κ2cos2θ+kz,e2×(εo-1κ2sinθcosθεe-1(-kz,e2-κ2cos2θ)±εo-1kz,eκsinθ).
Thus, we can have the unit vectors for magnetic field and electric displacement vector as
(27)h^o±=1koκ2cos2θ+kz,o2(-κ2sinθcosθkz,o2+κ2cos2θ∓kz,oκsinθ),h^e±=1κ2cos2θ+kz,e2(±kz,e0-κcosθ),d^o,±=1κ2cos2θ+kz,o2(±kz,o0-κcosθ),d^e±=1keκ2cos2θ+kz,e2(κ2sinθcosθ-kz,e2-κ2cos2θ±kz,eκsinθ),
where ko=κ2+kz,o2 and ke=κ2+kz,e2. 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 δ(z-z0) and δ'(z-z0) provides the equations determining these fields and constraints about the local field (Bx,By)=4πμ(Mx,My) and (Ex,Ey)=0. We summarize them from (13) and (14), respectively, for x, y, and z directions:
(28)κ2cos2θ+kz,o2ko(Bo,g,--Bo,g,+)=-i4πω~μPx+i4πμκsinθMz,κ2sinθcosθkoκ2cos2θ+kz,o2(Bo,g,--Bo,g,+)+kz,eκ2cos2θ+kz,e2(Be,g,++Be,g,-)=-i4πω~μPy-i4πμκcosθMz,Ez=-4πεoPz,κ2cos2θ+kz,e2ω~μεe(Be,g,+-Be,g,-)+iEzκsinθ-iω~Bx=0,koω~μεokz,okz,o2+κ2cos2θ(Bo,g,++Bo,g,-)+1ω~μεoκ2sinθcosθκ2cos2θ+kz,e2(Be,g,+-Be,g,-)-iEzκcosθ-iω~By=0,Bz=0.
After solving (28), the fields launching from the source plane are determined:
(29)Bo,g,+=-O1Pz+O2Px+O3My-O4Mx-O5Mz,Bo,g,-=-O1Pz-O2Px+O3My-O4Mx+O5Mz,Be,g,+=-O-1Py+O-2Px+O-3Pz+O-4Mx-O-5Mz,Be,g,-=-O-1Py+O-2Px-O-3Pz-O-4Mx-O-5Mz,
where the coefficients in front of the source are
(30)O1=(i2πκκ2cos2θ+kz,e2ω~μRcosθkokz,o+i2πκ3ω~μεeRsin2θcosθkokz,oεoκ2cos2θ+kz,e2),O2=i2πω~koμκ2cos2θ+kz,o2,O3=i2πω~2κ2cos2θ+kz,e2μ2εoRkokz,o,O4=i2πω~2κ2μ2εeRsinθcosθkokz,oκ2cos2θ+kz,e2,O5=i2πμκkosinθκ2cos2θ+kz,o2,O-1=i2πω~κ2cos2θ+kz,o2μkz,eR,O-2=i2πω~κ2μsinθcosθkz,eκ2cos2θ+kz,o2R,O-3=i2πω~κμεesinθεoκ2cos2θ+kz,e2,O-4=i2πω~2μ2εeκ2cos2θ+kz,e2,O-5=[i2πκ2cos2θ+kz,o2κμcosθkz,eR+i2πκ3μsin2θcosθkz,eκ2cos2θ+kz,o2R],
where the coefficient R in (30) is given as R=kz,o2+κ2cos2θ/kz,e2+κ2cos2θ. The fields in the regions outside of the anisotropic medium can be determined by solving boundary equations, and they can be written for z=0 and z=h, respectively, as
(31)-κ2sinθcosθkoκ2cos2θ+kz,o2(Bo,b,++Bo,b,-+Bo,g,-eikz,oz0)+kz,eκ2cos2θ+kz,e2(Be,b,+-Be,b,--Be,g,-eikz,ez0)=-Bs,-kz,Iω~cosθ+Bp,-sinθ,κ2cos2θ+kz,o2ko(Bo,b,++Bo,b,-+Bo,g,-eikz,oz0)=-Bs,-kz,Iω~sinθ-Bp,-cosθ,(Bo,b,+-Bo,b,--Bo,g,-eikz,oz0)kokz,oω~μεoκ2cos2θ+kz,o2+(Be,b,++Be,b,-+Be,g,-eikz,ez0)κ2sinθcosθω~μεoκ2cos2θ+kz,e2=Bs,-sinθ+Bp,-kz,Iω~cosθ,(Be,b,++Be,b,-+Be,g,-eikz,ez0)(-kz,e2-κ2cos2θ)ω~μεeκ2cos2θ+kz,e2=-Bs,-cosθ+Bp,-kz,Iω~sinθ,(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)+Bo,b,-e-ikz,oh)×-κ2sinθcosθkoκ2cos2θ+kz,o2+(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)-Be,b,-e-ikz,eh)×kz,eκ2cos2θ+kz,e2=kz,IIIBs,+eikz,IIIhω~cosθ+Bp,+eikz,IIIhsinθ(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)+Bo,b,-e-ikz,oh)×κ2cos2θ+kz,o2ko=kz,IIIBs,+eikz,IIIhω~sinθ-Bp,+eikz,IIIhcosθ(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)-Bo,b,-e-ikz,oh)×kokz,oεoω~μκ2cos2θ+kz,o2+(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)+Be,b,-e-ikz,eh)×κ2sinθcosθεoω~μκ2cos2θ+kz,e2=Bs,+eikz,IIIhsinθ-Bp,+eikz,IIIhkz,IIIω~cosθ(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)+Be,b,-e-ikz,eh)×-kz,e2-κ2cos2θω~μεeκ2cos2θ+kz,e2=-Bs,+eikz,IIIhcosθ-Bp,+eikz,IIIhkz,IIIω~sinθ.
From these equations, the origin of mode conversion is obviously accounted by the equations with both o- and e-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:(32)GB,PI=([(-W2O2eikz,oz0-W1O2e-ikz,oz0+W5O-2eikz,ez0+W4O-2e-ikz,ez0)h^s,-Ix^-(W5O-1eikz,ez0+W4O-1e-ikz,ez0)h^s,-Iy^+(-W2O1eikz,oz0+W1O1e-ikz,oz0-W5O-3eikz,ez0+W4O-3e-ikz,ez0)h^s,-Iz^]+[(W3O2eikz,oz0-W1O2e-ikz,oz0-W7O-2eikz,ez0-W6O-2e-ikz,ez0)h^p,-Ix^+(W7O-1eikz,ez0+W6O-1e-ikz,ez0)h^p,-Iy^+(W3O1eikz,oz0+W1O1e-ikz,oz0+W7O-3eikz,ez0-W6O-3e-ikz,ez0)])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,PI=([(-W2O2eikz,oz0-W1O2e-ikz,oz0+W5O-2eikz,ez0+W4O-2e-ikz,ez0)e^s,-Ix^-(W5O-1eikz,ez0+W4O-1e-ikz,ez0)e^s,-Iy^+(-W2O1eikz,oz0+W1O1e-ikz,oz0-W5O-3eikz,ez0+W4O-3e-ikz,ez0)e^s,-Iz^]+[(W3O2eikz,oz0-W1O2e-ikz,oz0-W7O-2eikz,ez0-W6O-2e-ikz,ez0)e^p,-Ix^+(W7O-1eikz,ez0+W6O-1e-ikz,ez0)e^p,-Iy^+(W3O1eikz,oz0+W1O1e-ikz,oz0+W7O-3eikz,ez0-W6O-3e-ikz,ez0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,MI=([(-W2O4eikz,oz0+W1O4e-ikz,oz0-W5O-4eikz,ez0+W4O-4e-ikz,ez0)h^s,-Ix^+(W2O3eikz,oz0-W1O3e-ikz,oz0)h^s,-Iy^+(W2O5eikz,oz0+W1O5e-ikz,oz0-W5O-5eikz,ez0-W4O-5e-ikz,ez0)h^s,-Iz^]+[(W3O4eikz,oz0+W1O4e-ikz,oz0+W7O-4eikz,ez0-W6O-4e-ikz,ez0)h^p,-Ix^-(W3O3eikz,oz0My+W1O3e-ikz,oz0)h^p,-Iy^+(-W3O5eikz,oz0+W1O5e-ikz,oz0+W7O-5eikz,ez0+W6O-5e-ikz,ez0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,MI=([(-W2O4eikz,oz0+W1O4e-ikz,oz0-W5O-4eikz,ez0+W4O-4e-ikz,ez0)e^s,-Ix^+(W2O3eikz,oz0-W1O3e-ikz,oz0)e^s,-Iy^+(W2O5eikz,oz0+W1O5e-ikz,oz0-W5O-5eikz,ez0-W4O-5e-ikz,ez0)e^s,-Iz^]+[(W3O4eikz,oz0+W1O4e-ikz,oz0+W7O-4eikz,ez0-W6O-4e-ikz,ez0)e^p,-Ix^-(W3O3eikz,oz0My+W1O3e-ikz,oz0)e^p,-Iy^+(-W3O5eikz,oz0+W1O5e-ikz,oz0+W7O-5eikz,ez0+W6O-5e-ikz,ez0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,PIII=([(-W1O2e-ikz,o(h-z0)-W2O2eikz,o(h-z0)+W4O-2e-ikz,e(h-z0)+W5O-2eikz,e(h-z0))h^s,+IIIx^-(W4O-1e-ikz,e(h-z0)+W5O-1eikz,e(h-z0))h^s,+IIIy^+(-W1O1e-ikz,o(h-z0)+W2O1eikz,o(h-z0)-W4O-3e-ikz,e(h-z0)+W5O-3eikz,e(h-z0))h^s,+IIIz^]+[(W1O2e-ikz,o(h-z0)-W3O2eikz,o(h-z0)+W6O-2e-ikz,e(h-z0)+W7O-2eikz,e(h-z0))h^p,+IIIx^-(W6O-1e-ikz,e(h-z0)+W7O-1eikz,e(h-z0))h^p,+IIIy^+(W1O1e-ikz,o(h-z0)+W3O1eikz,o(h-z0)-W6O-3e-ikz,e(h-z0)+W7O-3eikz,e(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,PIII=([(-W1O2e-ikz,o(h-z0)-W2O2eikz,o(h-z0)+W4O-2e-ikz,e(h-z0)+W5O-2eikz,e(h-z0))e^s,+IIIx^-(W4O-1e-ikz,e(h-z0)+W5O-1eikz,e(h-z0))e^s,+IIIy^+(-W1O1e-ikz,o(h-z0)+W2O1eikz,o(h-z0)-W4O-3e-ikz,e(h-z0)+W5O-3eikz,e(h-z0))e^s,+IIIz^]+[(W1O2e-ikz,o(h-z0)-W3O2eikz,o(h-z0)+W6O-2e-ikz,e(h-z0)+W7O-2eikz,e(h-z0))e^p,+IIIx^-(W6O-1e-ikz,e(h-z0)+W7O-1eikz,e(h-z0))e^p,+IIIy^+(W1O1e-ikz,o(h-z0)+W3O1eikz,o(h-z0)-W6O-3e-ikz,e(h-z0)+W7O-3eikz,e(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GB,MIII=([(-W1O4e-ikz,o(h-z0)+W2O4eikz,o(h-z0)-W4O-4e-ikz,e(h-z0)+W5O-4eikz,e(h-z0))h^s,+IIIx^+(W1O3e-ikz,o(h-z0)-W2O3eikz,o(h-z0))h^s,+IIIy^+(-W4O-5e-ikz,e(h-z0)+W2O5eikz,o(h-z0)e-ikzh+W1O5e-ikz,o(h-z0)-W5O-5eikz,e(h-z0))h^s,+IIIz^]+[(W1O4e-ikz,o(h-z0)+W3O4eikz,o(h-z0)-W6O-4e-ikz,e(h-z0)+W7O-4eikz,e(h-z0))h^p,+IIIx^-(W1O3e-ikz,o(h-z0)+W3O3eikz,o(h-z0))h^p,+IIIy^+(-W1O5e-ikz,o(h-z0)+W3O5eikz,o(h-z0)-W6O-5e-ikz,e(h-z0)-W7O-5eikz,e(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,MIII=([(-W1O4e-ikz,o(h-z0)+W2O4eikz,o(h-z0)-W4O-4e-ikz,e(h-z0)+W5O-4eikz,e(h-z0))e^s,+IIIx^+(W1O3e-ikz,o(h-z0)-W2O3eikz,o(h-z0))e^s,+IIIy^+(-W4O-5e-ikz,e(h-z0)+W2O5eikz,o(h-z0)e-ikzh+W1O5e-ikz,o(h-z0)-W5O-5eikz,e(h-z0))e^s,+IIIz^]+[(W1O4e-ikz,o(h-z0)+W3O4eikz,o(h-z0)-W6O-4e-ikz,e(h-z0)+W7O-4eikz,e(h-z0))e^p,+IIIx^-(W1O3e-ikz,o(h-z0)+W3O3eikz,o(h-z0))e^p,+IIIy^+(-W1O5e-ikz,o(h-z0)+W3O5eikz,o(h-z0)-W6O-5e-ikz,e(h-z0)-W7O-5eikz,e(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)).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 (Wi,i=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 formdiag(ε1,ε2,ε3) with ε1≠ε2≠ε3. Thus, the equation determining the dispersions of eigenmodes now becomes
(33)ω~4-ω~2(kx2+ky2ε3μ+kx2+kz2ε2μ+kz2+ky2ε1μ)+(kx2ε2ε3μ2+ky2ε1ε3μ2+kz2ε1ε3μ2)(kx2+kz2+ky2)=0.
The dispersion can be gotten as followes: kz,12=(-b+b2-4ac)/2a and kz,22=(-b-b2-4ac)/2a (each of them denotes two modes propagating (or decaying) along positive and negative z direction), where a, b, and c are, respectively, the coefficients in front of the fourth-order, second-order and zeroth-order terms of kz, and they have the following concrete forms:
(34)a=1ε1ε3μ,b=kx2ε2ε3μ+ky2ε1ε3μ-ω~2ε2μ-ω~2ε1μ,c=(kx2+ky2)(kx2ε2ε3μ+ky2ε1ε3μ)+kx2(1ε3μ+ω~2ε2μ)+ky2(1ε3μ+ω~2ε1μ)+ω~4.
The corresponding unit vectors of electric fields of eigenmodes can be solved as [23]
(35)e^i,±=(Xi2+1+Zi2)-1/2(Xi,1,±Zi),
where
(36)Xi=κ2cosθsinθ(κ2+kz,i2-ω~2ε3μ)(ω~2ε1μ-kz,i2-κ2sin2θ)(ω~2ε3μ-κ2)-κ2cos2θkz,i2Zi=kz,iκsinθ(κ2+kz,i2-ω~2ε1μ)(ω~2ε1μ-κ2sin2θ-kz,i2)(ω~2ε3μ-κ2)-κ2cos2θkz,i2.
In (36) kz,i is a positive square root from kz,i2 so that it denotes the propagating or decaying wave along positive z direction. Thus, the unit vectors of electric displacement vector and magnetic field can be gotten as(37)d^i,±=(ε12Xi2+ε22+ε32Zi2)-1/2(ε1Xi,ε2,±ε3Zi),h^i,±=1Ni(±(κsinθZi-kz,i),±(kz,iXi-Ziκcosθ),(κcosθ-κsinθXi)),whereNi=(κsinθZi-kz,i)2+(kz,iXi-Ziκcosθ)2+(κcosθ-κsinθXi)2.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, (Bx,By)=4πμ(Mx,My) and (Ex,Ey)=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:
(38)k1ω~μ(ε12X12+ε22+ε32Z12)1/2[B1,g,--B1,g,+]+k2ω~μ(ε12X22+ε22+ε32Z22)1/2[B2,g,--B2,g,+]+iEzκsinθ=iω~Bx,X1k1ω~μ(ε12X12+ε22+ε32Z12)1/2[B1,g,--B1,g,+]+X2k2ω~μ(ε12X22+ε22+ε32Z22)1/2[B2,g,--B2,g,+]-iEzκcosθ=iω~By,Bz=0-(kz,1X1-Z1κcosθ)N1(B1,g,++B1,g,-)-(kz,2X2-Z2κcosθ)N2(B2,g,++B2,g,-)=-i4πω~μPx+i4πμκsinθMz,(κsinθZ1-kz,1)N1(B1,g,++B1,g,-)+(κsinθZ2-kz,2)N2(B2,g,++B2,g,-)=-i4πω~μPy-i4πμκcosθMz,Ez=-4πε3Pz,
where k1=κ2+kz,12 and k2=κ2+kz,22. The solutions of (35) and (36) give the relations between the fields launching from the source plane and the sources as
(39)B1,g,-=12(V10Px+V11Py+V12Mz+V1Mx+V2My+V3Pz)B1,g,+=12(V10Px+V11Py+V12Mz-V1Mx-V2My-V3Pz)B2,g,-=12(V7Px+V8Py+V9Mz+V4Mx+V5My+V6Pz)B2,g,+=12(V7Px+V8Py+V9Mz-V4Mx-V5My-V6Pz).
The details of the definitions of parameters Vi(i=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 z=0 and z=h, respectively (we have used the fact that kz=kz,I=kz,III in the following equations, when κ⃑ is the same),
(40)S1(B1,b,+-B1,b,--B1,g,-eikz,1z0)+S2(B2,b,+-B2,b,--B2,g,-eikz,2z0)=-Bs,-kzω~cosθ+Bp,-sinθ,S3(B1,b,+-B1,b,--B1,g,-eikz,1z0)+S4(B2,b,+-B2,b,--B2,g,-eikz,2z0)=-Bs,-kzω~sinθ-Bp,-cosθ,S5X1(B1,b,++B1,b,-+B1,g,-eikz,1z0)+S6X2(B2,b,++B2,b,-+B2,g,-eikz,2z0)=Bs,-sinθ+Bp,-kzω~cosθ,S5(B1,b,++B1,b,-+B1,g,-eikz,1z0)+S6(B2,b,++B2,b,-+B2,g,-eikz,2z0)=-Bs,-cosθ+Bp,-kzω~sinθ,S1(B1,b,+eikz,1h-B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S2(B2,b,+eikz,2h-B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=kzcosθω~Bs,+eikzh+sinθBp,+eikzh,S3(B1,b,+eikz,1h-B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S4(B2,b,+eikz,2h-B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=kzsinθω~Bs,+eikzh-cosθBp,+eikzh,S5X1(B1,b,+eikz,1h+B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S6X2(B2,b,+eikz,2h+B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=sinθBs,+eikzh-kzcosθω~Bp,+eikzhS5(B1,b,+eikz,1h+B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S6(B2,b,+eikz,2h+B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=-cosθBs,+eikzh-kzsinθω~Bp,+eikzh,
where the parameters are given as
(41)S1=κsinθZ1-kz,1N1,S2=κsinθZ2-kz,2N2,S3=kz,1X1-Z1κcosθN1,S4=kz,2X2-Z2κcosθN2,S5=k1ω~με12X12+ε22+ε32Z12,S6=k2ω~με12X22+ε22+ε32Z22.
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):(42)GB,PI=12([(U2V10eikz,1z0-U1V10e-ikz,1z0+U5V7eikz,2z0+U4V7e-ikz,2z0)h^s,-Ix^+(U2V11eikz,1z0-U1V11e-ikz,1z0+U5V8eikz,2z0+U4V8e-ikz,2z0)h^s,-Iy^+(U2V3eikz,1z0+U1V3e-ikz,1z0+U5V6eikz,2z0-U4V6e-ikz,2z0)h^s,-Iz^]+[(-U3V10eikz,1z0-U1V10e-ikz,1z0-U7V7eikz,2z0-U6V7e-ikz,2z0)h^p,-Ix^+(-U3V11eikz,1z0-U1V11e-ikz,1z0-U7V8eikz,2z0-U6V8e-ikz,2z0)h^p,-Iy^+(-U3V3eikz,1z0+U1V3e-ikz,1z0-U7V6eikz,2z0+U6V6e-ikz,2z0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,PI=12([(U2V10eikz,1z0-U1V10e-ikz,1z0+U5V7eikz,2z0+U4V7e-ikz,2z0)e^s,-Ix^+(U2V11eikz,1z0-U1V11e-ikz,1z0+U5V8eikz,2z0+U4V8e-ikz,2z0)e^s,-Iy^+(U2V3eikz,1z0+U1V3e-ikz,1z0+U5V6eikz,2z0-U4V6e-ikz,2z0)e^s,-Iz^]+[(-U3V10eikz,1z0-U1V10e-ikz,1z0-U7V7eikz,2z0-U6V7e-ikz,2z0)e^p,-Ix^+(-U3V11eikz,1z0-U1V11e-ikz,1z0-U7V8eikz,2z0-U6V8e-ikz,2z0)e^p,-Iy^+(-U3V3eikz,1z0+U1V3e-ikz,1z0-U7V6eikz,2z0+U6V6e-ikz,2z0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,MI=([(U2V1eikz,1z0+U1V1e-ikz,1z0+U5V4eikz,2z0-U4V4e-ikz,2z0)h^s,-Ix^+(U2V2eikz,1z0+U1V2e-ikz,1z0+U5V5eikz,2z0-U4V5e-ikz,2z0)h^s,-Iy^+(U2V12eikz,1z0-U1V12e-ikz,1z0+U5V9eikz,2z0+U4V9e-ikz,2z0)h^s,-Iz^]+[(-U3V1eikz,1z0+U1V1e-ikz,1z0-U7V4eikz,2z0+U6V4e-ikz,2z0)h^p,-Ix^+(-U3V2eikz,1z0+U1V2e-ikz,1z0-U7V5eikz,2z0+U6V5e-ikz,2z0)h^p,-Iy^+(-U3V12eikz,1z0-U1V12e-ikz,1z0-U7V9eikz,2z0-U6V9e-ikz,2z0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,MI=([(U2V1eikz,1z0+U1V1e-ikz,1z0+U5V4eikz,2z0-U4V4e-ikz,2z0)e^s,-Ix^+(U2V2eikz,1z0+U1V2e-ikz,1z0+U5V5eikz,2z0-U4V5e-ikz,2z0)e^s,-Iy^+(U2V12eikz,1z0-U1V12e-ikz,1z0+U5V9eikz,2z0+U4V9e-ikz,2z0)e^s,-Iz^]+[(-U3V1eikz,1z0+U1V1e-ikz,1z0-U7V4eikz,2z0+U6V4e-ikz,2z0)e^p,-Ix^+(-U3V2eikz,1z0+U1V2e-ikz,1z0-U7V5eikz,2z0+U6V5e-ikz,2z0)e^p,-Iy^+(-U3V12eikz,1z0-U1V12e-ikz,1z0-U7V9eikz,2z0-U6V9e-ikz,2z0e^p,-Iz^)])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,PIII=12([(U1V10e-ikz,1(h-z0)-U2V10eikz,1(h-z0)+U4V7e-ikz,2(h-z0)+U5V7eikz,2(h-z0))h^s,+IIIx^+(U1V11e-ikz,1(h-z0)-U2V11eikz,1(h-z0)+U4V8e-ikz,2(h-z0)+U5V8eikz,2(h-z0))h^s,+IIIy^+(U1V3e-ikz,1(h-z0)+U2V3eikz,1(h-z0)+U4V6e-ikz,2(h-z0)-U5V6eikz,2(h-z0))h^s,+IIIz^]+[(-U1V10e-ikz,1(h-z0)-U3V10eikz,1(h-z0)+U6V7e-ikz,2(h-z0)+U7V7eikz,2(h-z0))h^p,+IIIx^+(-U1V11e-ikz,1(h-z0)-U3V11eikz,1(h-z0)+U6V8e-ikz,2(h-z0)+U7V8eikz,2(h-z0))h^p,+IIIy^+(-U1V3e-ikz,1(h-z0)+U3V3eikz,1(h-z0)+U6V6e-ikz,2(h-z0)-U7V6eikz,2(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,PIII=12([(U1V10e-ikz,1(h-z0)-U2V10eikz,1(h-z0)+U4V7e-ikz,2(h-z0)+U5V7eikz,2(h-z0))e^s,+IIIx^+(U1V11e-ikz,1(h-z0)-U2V11eikz,1(h-z0)+U4V8e-ikz,2(h-z0)+U5V8eikz,2(h-z0))e^s,+IIIy^+(U1V3e-ikz,1(h-z0)+U2V3eikz,1(h-z0)+U4V6e-ikz,2(h-z0)-U5V6eikz,2(h-z0))e^s,+IIIz^]+[(-U1V10e-ikz,1(h-z0)-U3V10eikz,1(h-z0)+U6V7e-ikz,2(h-z0)+U7V7eikz,2(h-z0))e^p,+IIIx^+(-U1V11e-ikz,1(h-z0)-U3V11eikz,1(h-z0)+U6V8e-ikz,2(h-z0)+U7V8eikz,2(h-z0))e^p,+IIIy^+(-U1V3e-ikz,1(h-z0)+U3V3eikz,1(h-z0)+U6V6e-ikz,2(h-z0)-U7V6eikz,2(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GB,MIII=12([(U1V1e-ikz,1(h-z0)+U2V1eikz,1(h-z0)+U4V4e-ikz,2(h-z0)-U5V4eikz,2(h-z0))h^s,+IIIx^+(U1V2e-ikz,1(h-z0)+U2V2eikz,1(h-z0)+U4V5e-ikz,2(h-z0)-U5V5eikz,2(h-z0))h^s,+IIIy^+(U1V12e-ikz,1(h-z0)-U2V12eikz,1(h-z0)+U4V9e-ikz,2(h-z0)+U5V9eikz,2(h-z0))h^s,+IIIz^]+[(-U1V1e-ikz,1(h-z0)+U3V1eikz,1(h-z0)+U6V4e-ikz,2(h-z0)-U7V4eikz,2(h-z0))h^p,+IIIx^+(-U1V2e-ikz,1(h-z0)+U3V2eikz,1(h-z0)+U6V5e-ikz,2(h-z0)-U7V5eikz,2(h-z0))h^p,+IIIy^+(-U1V12e-ikz,1(h-z0)-U3V12eikz,1(h-z0)+U6V9e-ikz,2(h-z0)+U7V9eikz,2(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,MIII=12([(U1V1e-ikz,1(h-z0)+U2V1eikz,1(h-z0)+U4V4e-ikz,2(h-z0)-U5V4eikz,2(h-z0))e^s,+IIIx^+(U1V2e-ikz,1(h-z0)+U2V2eikz,1(h-z0)+U4V5e-ikz,2(h-z0)-U5V5eikz,2(h-z0))e^s,+IIIy^+(U1V12e-ikz,1(h-z0)-U2V12eikz,1(h-z0)+U4V9e-ikz,2(h-z0)+U5V9eikz,2(h-z0))e^s,+IIIz^]+[(-U1V1e-ikz,1(h-z0)+U3V1eikz,1(h-z0)+U6V4e-ikz,2(h-z0)-U7V4eikz,2(h-z0))e^p,+IIIx^+(-U1V2e-ikz,1(h-z0)+U3V2eikz,1(h-z0)+U6V5e-ikz,2(h-z0)-U7V5eikz,2(h-z0))e^p,+IIIy^+(-U1V12e-ikz,1(h-z0)-U3V12eikz,1(h-z0)+U6V9e-ikz,2(h-z0)+U7V9eikz,2(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)).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 (Ui,i=1,2,…,7) in (42) will be given in Appendix C.


## 3.1. Sources in Uniaxial Crystal

Under the situation, the permittivity tensor in region II would be read asdiag(εo,εe,εo), where without losing generality we have assumed the anisotropy is along y direction (so that mode conversion still happens at the interface). Then the equation determining eigenmodes (5) becomes
(25)(ω~2-kx2εoμ-ky2εoμ-kz2εoμ)(ω~2-kx2εeμ-ky2εoμ-kz2εeμ)=0.
Thus, there are two dispersions gotten: the mode with dispersion ω~2-kx2/εoμ-ky2/εoμ-kz2/εoμ=0 is known as o-light, and the mode with dispersion ω~2-kx2/εeμ-ky2/εoμ-kz2/εeμ=0 is known as e-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
(26)e^o,±=1κ2cos2θ+kz,o2(±kz,o0-κcosθ),e^e,±=1(κ2cos2θ+kz,e2)εe2+κ2sin2θεo2κ2cos2θ+kz,e2×(εo-1κ2sinθcosθεe-1(-kz,e2-κ2cos2θ)±εo-1kz,eκsinθ).
Thus, we can have the unit vectors for magnetic field and electric displacement vector as
(27)h^o±=1koκ2cos2θ+kz,o2(-κ2sinθcosθkz,o2+κ2cos2θ∓kz,oκsinθ),h^e±=1κ2cos2θ+kz,e2(±kz,e0-κcosθ),d^o,±=1κ2cos2θ+kz,o2(±kz,o0-κcosθ),d^e±=1keκ2cos2θ+kz,e2(κ2sinθcosθ-kz,e2-κ2cos2θ±kz,eκsinθ),
where ko=κ2+kz,o2 and ke=κ2+kz,e2. 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 δ(z-z0) and δ'(z-z0) provides the equations determining these fields and constraints about the local field (Bx,By)=4πμ(Mx,My) and (Ex,Ey)=0. We summarize them from (13) and (14), respectively, for x, y, and z directions:
(28)κ2cos2θ+kz,o2ko(Bo,g,--Bo,g,+)=-i4πω~μPx+i4πμκsinθMz,κ2sinθcosθkoκ2cos2θ+kz,o2(Bo,g,--Bo,g,+)+kz,eκ2cos2θ+kz,e2(Be,g,++Be,g,-)=-i4πω~μPy-i4πμκcosθMz,Ez=-4πεoPz,κ2cos2θ+kz,e2ω~μεe(Be,g,+-Be,g,-)+iEzκsinθ-iω~Bx=0,koω~μεokz,okz,o2+κ2cos2θ(Bo,g,++Bo,g,-)+1ω~μεoκ2sinθcosθκ2cos2θ+kz,e2(Be,g,+-Be,g,-)-iEzκcosθ-iω~By=0,Bz=0.
After solving (28), the fields launching from the source plane are determined:
(29)Bo,g,+=-O1Pz+O2Px+O3My-O4Mx-O5Mz,Bo,g,-=-O1Pz-O2Px+O3My-O4Mx+O5Mz,Be,g,+=-O-1Py+O-2Px+O-3Pz+O-4Mx-O-5Mz,Be,g,-=-O-1Py+O-2Px-O-3Pz-O-4Mx-O-5Mz,
where the coefficients in front of the source are
(30)O1=(i2πκκ2cos2θ+kz,e2ω~μRcosθkokz,o+i2πκ3ω~μεeRsin2θcosθkokz,oεoκ2cos2θ+kz,e2),O2=i2πω~koμκ2cos2θ+kz,o2,O3=i2πω~2κ2cos2θ+kz,e2μ2εoRkokz,o,O4=i2πω~2κ2μ2εeRsinθcosθkokz,oκ2cos2θ+kz,e2,O5=i2πμκkosinθκ2cos2θ+kz,o2,O-1=i2πω~κ2cos2θ+kz,o2μkz,eR,O-2=i2πω~κ2μsinθcosθkz,eκ2cos2θ+kz,o2R,O-3=i2πω~κμεesinθεoκ2cos2θ+kz,e2,O-4=i2πω~2μ2εeκ2cos2θ+kz,e2,O-5=[i2πκ2cos2θ+kz,o2κμcosθkz,eR+i2πκ3μsin2θcosθkz,eκ2cos2θ+kz,o2R],
where the coefficient R in (30) is given as R=kz,o2+κ2cos2θ/kz,e2+κ2cos2θ. The fields in the regions outside of the anisotropic medium can be determined by solving boundary equations, and they can be written for z=0 and z=h, respectively, as
(31)-κ2sinθcosθkoκ2cos2θ+kz,o2(Bo,b,++Bo,b,-+Bo,g,-eikz,oz0)+kz,eκ2cos2θ+kz,e2(Be,b,+-Be,b,--Be,g,-eikz,ez0)=-Bs,-kz,Iω~cosθ+Bp,-sinθ,κ2cos2θ+kz,o2ko(Bo,b,++Bo,b,-+Bo,g,-eikz,oz0)=-Bs,-kz,Iω~sinθ-Bp,-cosθ,(Bo,b,+-Bo,b,--Bo,g,-eikz,oz0)kokz,oω~μεoκ2cos2θ+kz,o2+(Be,b,++Be,b,-+Be,g,-eikz,ez0)κ2sinθcosθω~μεoκ2cos2θ+kz,e2=Bs,-sinθ+Bp,-kz,Iω~cosθ,(Be,b,++Be,b,-+Be,g,-eikz,ez0)(-kz,e2-κ2cos2θ)ω~μεeκ2cos2θ+kz,e2=-Bs,-cosθ+Bp,-kz,Iω~sinθ,(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)+Bo,b,-e-ikz,oh)×-κ2sinθcosθkoκ2cos2θ+kz,o2+(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)-Be,b,-e-ikz,eh)×kz,eκ2cos2θ+kz,e2=kz,IIIBs,+eikz,IIIhω~cosθ+Bp,+eikz,IIIhsinθ(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)+Bo,b,-e-ikz,oh)×κ2cos2θ+kz,o2ko=kz,IIIBs,+eikz,IIIhω~sinθ-Bp,+eikz,IIIhcosθ(Bo,b,+eikz,oh+Bo,g,+eikz,o(h-z0)-Bo,b,-e-ikz,oh)×kokz,oεoω~μκ2cos2θ+kz,o2+(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)+Be,b,-e-ikz,eh)×κ2sinθcosθεoω~μκ2cos2θ+kz,e2=Bs,+eikz,IIIhsinθ-Bp,+eikz,IIIhkz,IIIω~cosθ(Be,b,+eikz,eh+Be,g,+eikz,e(h-z0)+Be,b,-e-ikz,eh)×-kz,e2-κ2cos2θω~μεeκ2cos2θ+kz,e2=-Bs,+eikz,IIIhcosθ-Bp,+eikz,IIIhkz,IIIω~sinθ.
From these equations, the origin of mode conversion is obviously accounted by the equations with both o- and e-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:(32)GB,PI=([(-W2O2eikz,oz0-W1O2e-ikz,oz0+W5O-2eikz,ez0+W4O-2e-ikz,ez0)h^s,-Ix^-(W5O-1eikz,ez0+W4O-1e-ikz,ez0)h^s,-Iy^+(-W2O1eikz,oz0+W1O1e-ikz,oz0-W5O-3eikz,ez0+W4O-3e-ikz,ez0)h^s,-Iz^]+[(W3O2eikz,oz0-W1O2e-ikz,oz0-W7O-2eikz,ez0-W6O-2e-ikz,ez0)h^p,-Ix^+(W7O-1eikz,ez0+W6O-1e-ikz,ez0)h^p,-Iy^+(W3O1eikz,oz0+W1O1e-ikz,oz0+W7O-3eikz,ez0-W6O-3e-ikz,ez0)])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,PI=([(-W2O2eikz,oz0-W1O2e-ikz,oz0+W5O-2eikz,ez0+W4O-2e-ikz,ez0)e^s,-Ix^-(W5O-1eikz,ez0+W4O-1e-ikz,ez0)e^s,-Iy^+(-W2O1eikz,oz0+W1O1e-ikz,oz0-W5O-3eikz,ez0+W4O-3e-ikz,ez0)e^s,-Iz^]+[(W3O2eikz,oz0-W1O2e-ikz,oz0-W7O-2eikz,ez0-W6O-2e-ikz,ez0)e^p,-Ix^+(W7O-1eikz,ez0+W6O-1e-ikz,ez0)e^p,-Iy^+(W3O1eikz,oz0+W1O1e-ikz,oz0+W7O-3eikz,ez0-W6O-3e-ikz,ez0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,MI=([(-W2O4eikz,oz0+W1O4e-ikz,oz0-W5O-4eikz,ez0+W4O-4e-ikz,ez0)h^s,-Ix^+(W2O3eikz,oz0-W1O3e-ikz,oz0)h^s,-Iy^+(W2O5eikz,oz0+W1O5e-ikz,oz0-W5O-5eikz,ez0-W4O-5e-ikz,ez0)h^s,-Iz^]+[(W3O4eikz,oz0+W1O4e-ikz,oz0+W7O-4eikz,ez0-W6O-4e-ikz,ez0)h^p,-Ix^-(W3O3eikz,oz0My+W1O3e-ikz,oz0)h^p,-Iy^+(-W3O5eikz,oz0+W1O5e-ikz,oz0+W7O-5eikz,ez0+W6O-5e-ikz,ez0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,MI=([(-W2O4eikz,oz0+W1O4e-ikz,oz0-W5O-4eikz,ez0+W4O-4e-ikz,ez0)e^s,-Ix^+(W2O3eikz,oz0-W1O3e-ikz,oz0)e^s,-Iy^+(W2O5eikz,oz0+W1O5e-ikz,oz0-W5O-5eikz,ez0-W4O-5e-ikz,ez0)e^s,-Iz^]+[(W3O4eikz,oz0+W1O4e-ikz,oz0+W7O-4eikz,ez0-W6O-4e-ikz,ez0)e^p,-Ix^-(W3O3eikz,oz0My+W1O3e-ikz,oz0)e^p,-Iy^+(-W3O5eikz,oz0+W1O5e-ikz,oz0+W7O-5eikz,ez0+W6O-5e-ikz,ez0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,PIII=([(-W1O2e-ikz,o(h-z0)-W2O2eikz,o(h-z0)+W4O-2e-ikz,e(h-z0)+W5O-2eikz,e(h-z0))h^s,+IIIx^-(W4O-1e-ikz,e(h-z0)+W5O-1eikz,e(h-z0))h^s,+IIIy^+(-W1O1e-ikz,o(h-z0)+W2O1eikz,o(h-z0)-W4O-3e-ikz,e(h-z0)+W5O-3eikz,e(h-z0))h^s,+IIIz^]+[(W1O2e-ikz,o(h-z0)-W3O2eikz,o(h-z0)+W6O-2e-ikz,e(h-z0)+W7O-2eikz,e(h-z0))h^p,+IIIx^-(W6O-1e-ikz,e(h-z0)+W7O-1eikz,e(h-z0))h^p,+IIIy^+(W1O1e-ikz,o(h-z0)+W3O1eikz,o(h-z0)-W6O-3e-ikz,e(h-z0)+W7O-3eikz,e(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,PIII=([(-W1O2e-ikz,o(h-z0)-W2O2eikz,o(h-z0)+W4O-2e-ikz,e(h-z0)+W5O-2eikz,e(h-z0))e^s,+IIIx^-(W4O-1e-ikz,e(h-z0)+W5O-1eikz,e(h-z0))e^s,+IIIy^+(-W1O1e-ikz,o(h-z0)+W2O1eikz,o(h-z0)-W4O-3e-ikz,e(h-z0)+W5O-3eikz,e(h-z0))e^s,+IIIz^]+[(W1O2e-ikz,o(h-z0)-W3O2eikz,o(h-z0)+W6O-2e-ikz,e(h-z0)+W7O-2eikz,e(h-z0))e^p,+IIIx^-(W6O-1e-ikz,e(h-z0)+W7O-1eikz,e(h-z0))e^p,+IIIy^+(W1O1e-ikz,o(h-z0)+W3O1eikz,o(h-z0)-W6O-3e-ikz,e(h-z0)+W7O-3eikz,e(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GB,MIII=([(-W1O4e-ikz,o(h-z0)+W2O4eikz,o(h-z0)-W4O-4e-ikz,e(h-z0)+W5O-4eikz,e(h-z0))h^s,+IIIx^+(W1O3e-ikz,o(h-z0)-W2O3eikz,o(h-z0))h^s,+IIIy^+(-W4O-5e-ikz,e(h-z0)+W2O5eikz,o(h-z0)e-ikzh+W1O5e-ikz,o(h-z0)-W5O-5eikz,e(h-z0))h^s,+IIIz^]+[(W1O4e-ikz,o(h-z0)+W3O4eikz,o(h-z0)-W6O-4e-ikz,e(h-z0)+W7O-4eikz,e(h-z0))h^p,+IIIx^-(W1O3e-ikz,o(h-z0)+W3O3eikz,o(h-z0))h^p,+IIIy^+(-W1O5e-ikz,o(h-z0)+W3O5eikz,o(h-z0)-W6O-5e-ikz,e(h-z0)-W7O-5eikz,e(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,MIII=([(-W1O4e-ikz,o(h-z0)+W2O4eikz,o(h-z0)-W4O-4e-ikz,e(h-z0)+W5O-4eikz,e(h-z0))e^s,+IIIx^+(W1O3e-ikz,o(h-z0)-W2O3eikz,o(h-z0))e^s,+IIIy^+(-W4O-5e-ikz,e(h-z0)+W2O5eikz,o(h-z0)e-ikzh+W1O5e-ikz,o(h-z0)-W5O-5eikz,e(h-z0))e^s,+IIIz^]+[(W1O4e-ikz,o(h-z0)+W3O4eikz,o(h-z0)-W6O-4e-ikz,e(h-z0)+W7O-4eikz,e(h-z0))e^p,+IIIx^-(W1O3e-ikz,o(h-z0)+W3O3eikz,o(h-z0))e^p,+IIIy^+(-W1O5e-ikz,o(h-z0)+W3O5eikz,o(h-z0)-W6O-5e-ikz,e(h-z0)-W7O-5eikz,e(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)).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 (Wi,i=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 formdiag(ε1,ε2,ε3) with ε1≠ε2≠ε3. Thus, the equation determining the dispersions of eigenmodes now becomes
(33)ω~4-ω~2(kx2+ky2ε3μ+kx2+kz2ε2μ+kz2+ky2ε1μ)+(kx2ε2ε3μ2+ky2ε1ε3μ2+kz2ε1ε3μ2)(kx2+kz2+ky2)=0.
The dispersion can be gotten as followes: kz,12=(-b+b2-4ac)/2a and kz,22=(-b-b2-4ac)/2a (each of them denotes two modes propagating (or decaying) along positive and negative z direction), where a, b, and c are, respectively, the coefficients in front of the fourth-order, second-order and zeroth-order terms of kz, and they have the following concrete forms:
(34)a=1ε1ε3μ,b=kx2ε2ε3μ+ky2ε1ε3μ-ω~2ε2μ-ω~2ε1μ,c=(kx2+ky2)(kx2ε2ε3μ+ky2ε1ε3μ)+kx2(1ε3μ+ω~2ε2μ)+ky2(1ε3μ+ω~2ε1μ)+ω~4.
The corresponding unit vectors of electric fields of eigenmodes can be solved as [23]
(35)e^i,±=(Xi2+1+Zi2)-1/2(Xi,1,±Zi),
where
(36)Xi=κ2cosθsinθ(κ2+kz,i2-ω~2ε3μ)(ω~2ε1μ-kz,i2-κ2sin2θ)(ω~2ε3μ-κ2)-κ2cos2θkz,i2Zi=kz,iκsinθ(κ2+kz,i2-ω~2ε1μ)(ω~2ε1μ-κ2sin2θ-kz,i2)(ω~2ε3μ-κ2)-κ2cos2θkz,i2.
In (36) kz,i is a positive square root from kz,i2 so that it denotes the propagating or decaying wave along positive z direction. Thus, the unit vectors of electric displacement vector and magnetic field can be gotten as(37)d^i,±=(ε12Xi2+ε22+ε32Zi2)-1/2(ε1Xi,ε2,±ε3Zi),h^i,±=1Ni(±(κsinθZi-kz,i),±(kz,iXi-Ziκcosθ),(κcosθ-κsinθXi)),whereNi=(κsinθZi-kz,i)2+(kz,iXi-Ziκcosθ)2+(κcosθ-κsinθXi)2.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, (Bx,By)=4πμ(Mx,My) and (Ex,Ey)=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:
(38)k1ω~μ(ε12X12+ε22+ε32Z12)1/2[B1,g,--B1,g,+]+k2ω~μ(ε12X22+ε22+ε32Z22)1/2[B2,g,--B2,g,+]+iEzκsinθ=iω~Bx,X1k1ω~μ(ε12X12+ε22+ε32Z12)1/2[B1,g,--B1,g,+]+X2k2ω~μ(ε12X22+ε22+ε32Z22)1/2[B2,g,--B2,g,+]-iEzκcosθ=iω~By,Bz=0-(kz,1X1-Z1κcosθ)N1(B1,g,++B1,g,-)-(kz,2X2-Z2κcosθ)N2(B2,g,++B2,g,-)=-i4πω~μPx+i4πμκsinθMz,(κsinθZ1-kz,1)N1(B1,g,++B1,g,-)+(κsinθZ2-kz,2)N2(B2,g,++B2,g,-)=-i4πω~μPy-i4πμκcosθMz,Ez=-4πε3Pz,
where k1=κ2+kz,12 and k2=κ2+kz,22. The solutions of (35) and (36) give the relations between the fields launching from the source plane and the sources as
(39)B1,g,-=12(V10Px+V11Py+V12Mz+V1Mx+V2My+V3Pz)B1,g,+=12(V10Px+V11Py+V12Mz-V1Mx-V2My-V3Pz)B2,g,-=12(V7Px+V8Py+V9Mz+V4Mx+V5My+V6Pz)B2,g,+=12(V7Px+V8Py+V9Mz-V4Mx-V5My-V6Pz).
The details of the definitions of parameters Vi(i=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 z=0 and z=h, respectively (we have used the fact that kz=kz,I=kz,III in the following equations, when κ⃑ is the same),
(40)S1(B1,b,+-B1,b,--B1,g,-eikz,1z0)+S2(B2,b,+-B2,b,--B2,g,-eikz,2z0)=-Bs,-kzω~cosθ+Bp,-sinθ,S3(B1,b,+-B1,b,--B1,g,-eikz,1z0)+S4(B2,b,+-B2,b,--B2,g,-eikz,2z0)=-Bs,-kzω~sinθ-Bp,-cosθ,S5X1(B1,b,++B1,b,-+B1,g,-eikz,1z0)+S6X2(B2,b,++B2,b,-+B2,g,-eikz,2z0)=Bs,-sinθ+Bp,-kzω~cosθ,S5(B1,b,++B1,b,-+B1,g,-eikz,1z0)+S6(B2,b,++B2,b,-+B2,g,-eikz,2z0)=-Bs,-cosθ+Bp,-kzω~sinθ,S1(B1,b,+eikz,1h-B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S2(B2,b,+eikz,2h-B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=kzcosθω~Bs,+eikzh+sinθBp,+eikzh,S3(B1,b,+eikz,1h-B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S4(B2,b,+eikz,2h-B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=kzsinθω~Bs,+eikzh-cosθBp,+eikzh,S5X1(B1,b,+eikz,1h+B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S6X2(B2,b,+eikz,2h+B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=sinθBs,+eikzh-kzcosθω~Bp,+eikzhS5(B1,b,+eikz,1h+B1,b,-e-ikz,1h+B1,g,+eikz,1(h-z0))+S6(B2,b,+eikz,2h+B2,b,-e-ikz,2h+B2,g,+eikz,2(h-z0))=-cosθBs,+eikzh-kzsinθω~Bp,+eikzh,
where the parameters are given as
(41)S1=κsinθZ1-kz,1N1,S2=κsinθZ2-kz,2N2,S3=kz,1X1-Z1κcosθN1,S4=kz,2X2-Z2κcosθN2,S5=k1ω~με12X12+ε22+ε32Z12,S6=k2ω~με12X22+ε22+ε32Z22.
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):(42)GB,PI=12([(U2V10eikz,1z0-U1V10e-ikz,1z0+U5V7eikz,2z0+U4V7e-ikz,2z0)h^s,-Ix^+(U2V11eikz,1z0-U1V11e-ikz,1z0+U5V8eikz,2z0+U4V8e-ikz,2z0)h^s,-Iy^+(U2V3eikz,1z0+U1V3e-ikz,1z0+U5V6eikz,2z0-U4V6e-ikz,2z0)h^s,-Iz^]+[(-U3V10eikz,1z0-U1V10e-ikz,1z0-U7V7eikz,2z0-U6V7e-ikz,2z0)h^p,-Ix^+(-U3V11eikz,1z0-U1V11e-ikz,1z0-U7V8eikz,2z0-U6V8e-ikz,2z0)h^p,-Iy^+(-U3V3eikz,1z0+U1V3e-ikz,1z0-U7V6eikz,2z0+U6V6e-ikz,2z0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,PI=12([(U2V10eikz,1z0-U1V10e-ikz,1z0+U5V7eikz,2z0+U4V7e-ikz,2z0)e^s,-Ix^+(U2V11eikz,1z0-U1V11e-ikz,1z0+U5V8eikz,2z0+U4V8e-ikz,2z0)e^s,-Iy^+(U2V3eikz,1z0+U1V3e-ikz,1z0+U5V6eikz,2z0-U4V6e-ikz,2z0)e^s,-Iz^]+[(-U3V10eikz,1z0-U1V10e-ikz,1z0-U7V7eikz,2z0-U6V7e-ikz,2z0)e^p,-Ix^+(-U3V11eikz,1z0-U1V11e-ikz,1z0-U7V8eikz,2z0-U6V8e-ikz,2z0)e^p,-Iy^+(-U3V3eikz,1z0+U1V3e-ikz,1z0-U7V6eikz,2z0+U6V6e-ikz,2z0)e^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,MI=([(U2V1eikz,1z0+U1V1e-ikz,1z0+U5V4eikz,2z0-U4V4e-ikz,2z0)h^s,-Ix^+(U2V2eikz,1z0+U1V2e-ikz,1z0+U5V5eikz,2z0-U4V5e-ikz,2z0)h^s,-Iy^+(U2V12eikz,1z0-U1V12e-ikz,1z0+U5V9eikz,2z0+U4V9e-ikz,2z0)h^s,-Iz^]+[(-U3V1eikz,1z0+U1V1e-ikz,1z0-U7V4eikz,2z0+U6V4e-ikz,2z0)h^p,-Ix^+(-U3V2eikz,1z0+U1V2e-ikz,1z0-U7V5eikz,2z0+U6V5e-ikz,2z0)h^p,-Iy^+(-U3V12eikz,1z0-U1V12e-ikz,1z0-U7V9eikz,2z0-U6V9e-ikz,2z0)h^p,-Iz^])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GE,MI=([(U2V1eikz,1z0+U1V1e-ikz,1z0+U5V4eikz,2z0-U4V4e-ikz,2z0)e^s,-Ix^+(U2V2eikz,1z0+U1V2e-ikz,1z0+U5V5eikz,2z0-U4V5e-ikz,2z0)e^s,-Iy^+(U2V12eikz,1z0-U1V12e-ikz,1z0+U5V9eikz,2z0+U4V9e-ikz,2z0)e^s,-Iz^]+[(-U3V1eikz,1z0+U1V1e-ikz,1z0-U7V4eikz,2z0+U6V4e-ikz,2z0)e^p,-Ix^+(-U3V2eikz,1z0+U1V2e-ikz,1z0-U7V5eikz,2z0+U6V5e-ikz,2z0)e^p,-Iy^+(-U3V12eikz,1z0-U1V12e-ikz,1z0-U7V9eikz,2z0-U6V9e-ikz,2z0e^p,-Iz^)])×exp(iκ⃑·R⃑)exp(-ikz,Iz),GB,PIII=12([(U1V10e-ikz,1(h-z0)-U2V10eikz,1(h-z0)+U4V7e-ikz,2(h-z0)+U5V7eikz,2(h-z0))h^s,+IIIx^+(U1V11e-ikz,1(h-z0)-U2V11eikz,1(h-z0)+U4V8e-ikz,2(h-z0)+U5V8eikz,2(h-z0))h^s,+IIIy^+(U1V3e-ikz,1(h-z0)+U2V3eikz,1(h-z0)+U4V6e-ikz,2(h-z0)-U5V6eikz,2(h-z0))h^s,+IIIz^]+[(-U1V10e-ikz,1(h-z0)-U3V10eikz,1(h-z0)+U6V7e-ikz,2(h-z0)+U7V7eikz,2(h-z0))h^p,+IIIx^+(-U1V11e-ikz,1(h-z0)-U3V11eikz,1(h-z0)+U6V8e-ikz,2(h-z0)+U7V8eikz,2(h-z0))h^p,+IIIy^+(-U1V3e-ikz,1(h-z0)+U3V3eikz,1(h-z0)+U6V6e-ikz,2(h-z0)-U7V6eikz,2(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,PIII=12([(U1V10e-ikz,1(h-z0)-U2V10eikz,1(h-z0)+U4V7e-ikz,2(h-z0)+U5V7eikz,2(h-z0))e^s,+IIIx^+(U1V11e-ikz,1(h-z0)-U2V11eikz,1(h-z0)+U4V8e-ikz,2(h-z0)+U5V8eikz,2(h-z0))e^s,+IIIy^+(U1V3e-ikz,1(h-z0)+U2V3eikz,1(h-z0)+U4V6e-ikz,2(h-z0)-U5V6eikz,2(h-z0))e^s,+IIIz^]+[(-U1V10e-ikz,1(h-z0)-U3V10eikz,1(h-z0)+U6V7e-ikz,2(h-z0)+U7V7eikz,2(h-z0))e^p,+IIIx^+(-U1V11e-ikz,1(h-z0)-U3V11eikz,1(h-z0)+U6V8e-ikz,2(h-z0)+U7V8eikz,2(h-z0))e^p,+IIIy^+(-U1V3e-ikz,1(h-z0)+U3V3eikz,1(h-z0)+U6V6e-ikz,2(h-z0)-U7V6eikz,2(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GB,MIII=12([(U1V1e-ikz,1(h-z0)+U2V1eikz,1(h-z0)+U4V4e-ikz,2(h-z0)-U5V4eikz,2(h-z0))h^s,+IIIx^+(U1V2e-ikz,1(h-z0)+U2V2eikz,1(h-z0)+U4V5e-ikz,2(h-z0)-U5V5eikz,2(h-z0))h^s,+IIIy^+(U1V12e-ikz,1(h-z0)-U2V12eikz,1(h-z0)+U4V9e-ikz,2(h-z0)+U5V9eikz,2(h-z0))h^s,+IIIz^]+[(-U1V1e-ikz,1(h-z0)+U3V1eikz,1(h-z0)+U6V4e-ikz,2(h-z0)-U7V4eikz,2(h-z0))h^p,+IIIx^+(-U1V2e-ikz,1(h-z0)+U3V2eikz,1(h-z0)+U6V5e-ikz,2(h-z0)-U7V5eikz,2(h-z0))h^p,+IIIy^+(-U1V12e-ikz,1(h-z0)-U3V12eikz,1(h-z0)+U6V9e-ikz,2(h-z0)+U7V9eikz,2(h-z0))h^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)),GE,MIII=12([(U1V1e-ikz,1(h-z0)+U2V1eikz,1(h-z0)+U4V4e-ikz,2(h-z0)-U5V4eikz,2(h-z0))e^s,+IIIx^+(U1V2e-ikz,1(h-z0)+U2V2eikz,1(h-z0)+U4V5e-ikz,2(h-z0)-U5V5eikz,2(h-z0))e^s,+IIIy^+(U1V12e-ikz,1(h-z0)-U2V12eikz,1(h-z0)+U4V9e-ikz,2(h-z0)+U5V9eikz,2(h-z0))e^s,+IIIz^]+[(-U1V1e-ikz,1(h-z0)+U3V1eikz,1(h-z0)+U6V4e-ikz,2(h-z0)-U7V4eikz,2(h-z0))e^p,+IIIx^+(-U1V2e-ikz,1(h-z0)+U3V2eikz,1(h-z0)+U6V5e-ikz,2(h-z0)-U7V5eikz,2(h-z0))e^p,+IIIy^+(-U1V12e-ikz,1(h-z0)-U3V12eikz,1(h-z0)+U6V9e-ikz,2(h-z0)+U7V9eikz,2(h-z0))e^p,+IIIz^])×exp(iκ⃑·R⃑)exp(ikz,III(z-h)).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 (Ui,i=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.



---

*Source: 532316-2012-09-06.xml*