Open access peer-reviewed chapter

Green Function

Written By

Jing Huang

Submitted: 18 October 2016 Reviewed: 21 February 2017 Published: 14 June 2017

DOI: 10.5772/68028

From the Edited Volume

Recent Studies in Perturbation Theory

Edited by Dimo I. Uzunov

Chapter metrics overview

2,035 Chapter Downloads

View Full Metrics


Both the scalar Green function and the dyadic Green function of an electromagnetic field and the transform from the scalar to dyadic Green function are introduced. The Green function of a transmission line and the propagators are also presented in this chapter.


  • Green function
  • boundary condition
  • scatter
  • propagator
  • convergence

1. Introduction

In 1828, Green introduced a function, which he called a potential, for calculating the distribution of a charge on a surface bounding a region in Rn in the presence of external electromagnetic forces. The Green function has been an interesting topic in modern physics and engineering, especially for the electromagnetic theory in various source distributions (charge, current, and magnetic current), various construct conductors, and dielectric. Even though most problems can be solved without the use of Green functions, the symbolic simplicity with which they could be used to express relationships makes the formulations of many problems simpler and more compact. Moreover, it is easier to conceptualize many problems; especially the dyadic Green function is generalized to layered media of planar, cylindrical, and spherical configurations.


2. Definition of Green function

2.1. Mathematics definition

For the linear operator, there are: L^x=f(t), t > 0;


Rewriting Eq. (1) as:


Defining the Green function as:


So, the solution of Eq. (1) is:


We give several types of Green functions [1]



3. The scalar Green function

3.1. The scalar Green function of an electromagnetic field

The Green function of a wave equation is the solution of the wave equation for a point source [2]. And when the solution to the wave equation due to a point source is known, the solution due to a general source can be obtained by the principle of linear superposition (see Figure 1).

Figure 1.

The radiation of a source s(r) in a volume V.

This is merely a result of the linearity of the wave equation, and that a general source is just a linear superposition of point sources. For example, to obtain the solution to the scalar wave equation in V in Figure 1


we first seek the Green function in the same V, which is the solution to the following equation:


Given g (r, r′), φ(r) can be found easily from the principle of linear superposition, since g (r, r′) is the solution to Eq. (5) with a point source on the right-hand side. To see this more clearly, note that an arbitrary source s(r) is just


which is actually a linear superposition of point sources in mathematical terms. Consequently, the solution to Eq. (5) is just


which is an integral linear superposition of the solution of Eq. (6). Moreover, it can be seen that g(r, r′) ≡ g(r′, r,) from reciprocity irrespective of the shape of V.

To find the solution of Eq. (6) for an unbounded, homogeneous medium, one solves it in spherical coordinates with the origin at r'. By so doing, Eq. (6) becomes


But due to the spherical symmetry of a point source, g(r) must also be spherically symmetric. Then, for r ≠ 0, adopt the proper coordinate origin (the vector r is replaced by the scalar r), the homogeneous, spherically symmetric solution to Eq. (9) is given by


Since sources are absent at infinity, physical grounds then imply that only an outgoing solution can exist; hence,


The constant c is found by matching the singularities at the origin on both sides of Eq. (9). To do this, we substitute Eq. (11) into Eq. (9) and integrate Eq. (9) over a small volume about the origin to yield


Note that the second integral vanishes when ∆V → 0 because dV = 4πr2dr. Moreover, the first integral in Eq. (12) can be converted into a surface integral using Gauss theorem to obtain


or c = 1/(4π).

The solution to Eq. (6) must depend only on rr′. Therefore, in general,


implying that g(r, r') is translationally invariant for unbounded, homogeneous media. Consequently, the solution to Eq. (5), from Eq. (9), is then


Once ϕ(r) and n^ϕ(r) are known on S, then ϕ(r′) away from S could be found


3.2. The scalar Green functions of one-dimensional transmission lines

We consider a transmission line excited by a distributed current source, K(x), as sketched in Figure 2. The line may be finite or infinite, and it may be terminated at either end with impedance or by another line [3]. For a harmonically oscillating current source K(x), the voltage and the current on the line satisfy the following pair of equations:

Figure 2.

Transmission line excited by a distributed current source, K(x).


L and C denote, respectively, the distributed inductance and capacitance of the line.

By eliminating I(x) between Eq. (17) and Eq. (18), there is


where k=ωLC denotes the propagation constant of the line. Eq. (19) has been designated as an inhomogeneous one-dimensional scalar wave equation.

The Green function pertaining to a one-dimensional scalar wave equation of the form of Eq. (19), denoted by g(x, x′), is a solution of the Eq. (9). The solution for g(x, x′) is not completely determined unless there are two boundary conditions which the function must satisfy at the extremities of the spatial domain in which the function is defined. The boundary conditions which must be satisfied by g(x, x′) are the same as those dictated by the original function which we intend to determine, namely, V(x) in the present case. For this reason, the Green functions are classified according to the boundary conditions, which they must obey. Some of the typical ones (for the transmission line) are illustrated in Figure 3.

Figure 3.

Classification of Green functions according to the boundary conditions.

In general, the subscript 0 designates infinite domain so that we have outgoing waves at x±, often called the radiation condition. Subscript 1 means that one of the boundary conditions satisfies the so-called Dirichlet condition, while the other satisfies the radiation condition. When one of the boundary conditions satisfies the so-called Neumann condition, we use subscript 2. Subscript 3 is reserved for the mixed type. Actually, we should have used a double subscript for two distinct boundary conditions. For example, case (b) of Figure 3 should be denoted by g01, indicating that one radiation condition and one Dirichlet condition are involved. With such an understanding, the simplified notation should be acceptable.

In case (d), a superscript becomes necessary because we have two sets of line voltage and current (V1, I1) and (V2, I2) in this problem, and the Green function also has different forms in the two regions. The first superscript denotes the region where this function is defined, and the second superscript denotes the region where the source is located.

Let the domain of x corresponds to (x1, x2). The function g(x, x′) in Eq. (9) can represent any of the three types, g0, g1, and g2, illustrated in Figures 3a–c, respectively. The treatment of case (d) is slightly different, and it will be formulated later.

(a) By multiplying Eq. (19) by g(x, x′) and Eq. (9) by V(x) and taking the difference of the two resultant equations, we obtain


The first term at the right-hand side of the above equation is simply V(xl), and the term at the left-hand side can be simplified by integration by parts, which gives


If we use the unprimed variable x to denote the position of a field point, as usually is the case, Eq. (21) can be changed to [4]


The last identity is due to the symmetrical property of the Green function. The shifting of the primed and unprimed variables is often practiced in our work. For this reason, it is important to point out that g(x′, x), by definition, satisfies the Eq. (9).

The general solutions for Eq. (9) in the two regions (see Figure 3a) are


The choice of the above functions is done with the proper satisfaction of boundary conditions at infinity. At x = x', the function must be continuous, and its derivative is discontinuous.

They are: [g0(x,x)]x0x+0=0, and [dg0(x,x)dx]x0x+0=-1

The physical interpretation of these two conditions is that the voltage at x' is continuous, but the difference of the line currents at x' must be equal to the source current.

(b) The choice of this type of function is done with the proper satisfaction of boundary conditions. At x = x', the function must be continuous, its derivative is discontinuous, and a Dirichlet condition is satisfied at x = 0.


In view of Eq. (24), it can be interpreted as consisting of an incident and a scattered wave; that is


where g1s(x,x)=i2keik(x+x).

Such a notion is not only physically useful, but mathematically it offers a shortcut to finding a composite Green function. It is called as the shortcut method or the method of scattering superposition.

(c) Similarly, the method of scattering superposition suggests that we can start with


To satisfy the Neumann condition at x = 0, we require




(d) In this case, we have two differential equations to start with


It is assumed that the current source is located in region 1 (see Figure 3d). We introduce two Green functions of the third kind, denoted by g(11) (x, x') and g(21) (x, x'). g(21), the first number of the superscript corresponds to the region where the function is defined. The second number corresponds to the region where the source is located; then


At the junction corresponding to x = 0, g(11) and g(21) satisfy the boundary condition that


The last condition corresponds to the physical requirement that the current at the junction must be continuous. Again, by means of the method of scattering superposition, there are


The characteristic impedance of the lines, respectively, is


By the boundary condition, there are


Example: Green function solution of nonlinear Schrodinger equation in the time domain [5].

The nonlinear Schrodinger equation including nonresonant and resonant nonlinear items is:


Where A is the field, β2 and β3 are the second and third order dispersion, respectively. A(z) is the fiber absorption profile. k0=ω0/c, ω0 is the center frequency. Aeff is the effective core area. n is the refractive index.


where g(ω1 + ω2 + ω3) is the Raman gain and f(ω1 + ω2 + ω3) is the Raman nongain coefficient. Г is the attenuation coefficient.

The original nonlinear part is divided into the nonresonant and resonant susceptibility items χNR(3) and χR(3). The solution has the form:


Then, there is:




and taking the operator V^(t) as a perturbation item, the eigenequation n=2kinn!βnnφTn=Eφ is


Assuming E = 1, we get the corresponding characteristic equation:


Its characteristic roots are r1,r2,r3. The solution can be represented as:


where ϕm=exp(irmt),m=1,2,3, and c1,c2,c3 are determined by the initial pulse. The Green function of Eq. (47) is:


Constructing the Green function as:

G0(t,t)={a1φ1+a2φ2+a3φ3, t>tb1φ1+b2φ2+b3φ3, t<tE51

At the point t = t′, there are:


It is reasonable to let b1 = b2 = b3 = 0, then:

a1=φ2φ˙3φ˙2φ3W(t), a2=φ3φ˙1φ˙3φ1W(t), a3=φ1φ˙2φ˙1φ2W(t)E55

Finally, the solution of Eq. (44) can be written with the eigenfunction and Green function:

φ(t)=φ(t)+G0(t,t)V(t)φ(t)dt=ϕ(t)+G0(t,t,E)V(t)ϕ(t)dt+dtG0(t,t,E)V(t)G0(t,t,E)V(t)φ(t)dt=ϕ(t)+G0(t,t,E)V(t)ϕ(t)dt+dtG0(t,t,E)V(t)G0(t,t,E)V(t)ϕ(t)dt++dtG0(t,t)V(t)G0(t,t)V(t)dttimes lG0(tl,tl+1)V(tl+1)φ(tl+1)dtl+1E57

The accuracy can be estimated by the last term of Eq. (57).


4. The dyadic Green function

4.1. The dyadic Green function for the electromagnetic field in a homogeneous isotropic medium

The Green function for the scalar wave equation could be used to find the dyadic Green function for the vector wave equation in a homogeneous, isotropic medium [3]. First, notice that the vector wave equation in a homogeneous, isotropic medium is


Then, by using the fact that ××E(r)=-2E+E and that E=ρ/ε=J/iωε, which follows from the continuity equation, we can rewrite Eq. (58) as


where I^ is an identity operator. In Cartesian coordinates, there are actually three scalar wave equations embedded in the above vector equation, each of which can be solved easily in the manner of Eq. (4). Consequently,


where g(r′r)is the unbounded medium scalar Green function. Moreover, by using the vector identities gf=fg+gf and gF=gF+(g)F, it can be shown that




Hence, Eq. (60) can be rewritten as


It can also be derived using scalar and vector potentials.

Alternatively, Eq. (63) can be written as




is a dyad known as the dyadic Green function for the electric field in an unbounded, homogeneous medium. (A dyad is a 3 × 3 matrix that transforms a vector to a vector. It is also a second rank tensor). Even though Eq. (64) is established for an unbounded, homogeneous medium, such a general relationship also exists in a bounded, homogeneous medium. It could easily be shown from reciprocity that




is the relation between Ji and the electric field produced by Jj. Notice that the above equation implies [6]


Then, by taking transpose of Eq. (66b), Eq. (64) becomes


Alternatively, the dyadic Green function for an unbounded, homogeneous medium can also be written as


By substituting Eq. (67) back into Eq. (58) and writing


we can show quite easily that


Equation (64) or (67), due to the ∇∇ operator inside the integration operating on g(r′r), has a singularity of 1/|r′r|3 when r′r. Consequently, it has to be redefined in this case for it does not converge uniformly, specifically, when r is also in the source region occupied by J(r). Hence, at this point, the evaluation of Eq. (67) in a source region is undefined.

And as the vector analog of Eq. (16)


4.2. The boundary condition

The dyadic Green function is introduced mainly to formulate various canonical electromagnetic problems in a systematic manner to avoid treatments of many special cases which can be treated as one general problem [3, 7, 8]. Some typical problems are illustrated in Figure 4 where (a) shows a current source in the presence of a conducting sphere located in air, (b) shows a conducting cylinder with an aperture which is excited by some source inside the cylinder, (c) shows a rectangular waveguide with a current source placed inside the guide, and (d) shows two semi-infinite isotropic media in contact, such as air and “flat” earth with a current source placed in one of the regions.

Figure 4.

Some typical boundary value problems.

Unless specified otherwise, we assume that for problems involving only one medium such as (a), (b), and (c) the medium is air, then the wave number k is equal to ω(μ0ε0)1/2=2π/λ. The electromagnetic fields in these cases are solutions of the wave Eq. (62) and


The fields must satisfy the boundary conditions required by these problems.

In general, using the notations G^e and G^m to denote, respectively, the electric and the magnetic dyadic Green functions; they are solutions of the dyadic differential equations


is the same as Eq. (70), and there is


(a) and (b): Electric dyadic Green function (the first kind, using the subscript 1 denotes G^e1,G^m1, and the subscript “0” represents the free-space condition that the environment does not have any scattering object) is required to satisfy the dyadic Dirichlet condition on Sd, namely,


So, for (a)


and for (b)


(c) the electric dyadic Green function is required to satisfy the dyadic boundary condition on Sd, namely,


(d) For problems involving two isotropic media such as the configuration shown in Figure 4d, there are two sets of fields [9]. The wave numbers in these two regions are denoted by k1=ω(μ1ε1)1/2 and k2=ω(μ2ε2)1/2. There are four functions for the dyadic Green function of the electric type and another four functions for the magnetic type, denoted, respectively, by G^e11G^e12G^e21 and G^e22, and G^m11G^m12G^m21 and G^m22. The superscript notation in G^e11 means that both the field point and the source point are located in region 1. For G^e21, it means that the field point is located in region 1 and the source point is located in region 2. A current source is located in region 1 only, and the two sets of wave equations are




There are


At the interface, the electromagnetic field and the corresponding dyadic Green function satisfy the following boundary conditions


The electric fields are


5. Vector wave functions, L, M, and N

The vector wave functions are the building blocks of the eigenfunction expansions of various kinds of dyadic Green functions. These functions were first introduced by Hansen [1012] in formulating certain electromagnetic problems.Three kinds of vector wave functions, denoted by L, M, and N, are solutions of the homogeneous vector Helmholtz equation. To derive the eigenfunction expansion of the magnetic dyadic Green functions that are solenoidal and satisfy with the vector wave equation, the L functions are not needed. If we try to find eigenfunction expansion of the electric dyadic Green functions then the L functions are also needed.

A vector wave function, by definition, is an eigenfunction or a characteristic function, which is a solution of the homogeneous vector wave equation ××Fκ2F=0.

There are two independent sets of vector wave functions, which can be constructed using the characteristic function pertaining to a scalar wave equation as the generating function. One kind of vector wave function, called the Cartesian or rectilinear vector wave function, is formed if we let


where ψ1 denotes a characteristic function, which satisfies the scalar wave equation


And c denotes a constant vector, such as x, y, or z. For convenience, we shall designate c as the piloting vector and Ψ as the generating function. Another kind, designated as the spherical vector wavefunction, will be introduced later, whereby the piloting vector is identified as the spherical radial vector R.

Actually, substituting Eq. (91) into Eq. (92), it is


The set of functions so obtained


Ψ2,Ψ3 denote the characteristic functions which also satisfy (92) but may be different from the function used to define M1.

In the following, the expressions for the dyadic Green functions of a rectangular waveguide will be derived asserting to the vector wave functions. The method and the general procedure would apply equally well to other bodies (cylindrical waveguide, circular cylinder in free space, and inhomogeneous media and moving medium).

Figure 5 shows the orientation of the guide with respect to the rectangular coordinate system, and we will choose the unit vector z to represent the piloting vector c.

Figure 5.

A rectangular waveguide.

The scalar wave function


where kx2+ky2+h2=κ2.

the constants kx and ky should have the following characteristic values


The complete expression and the notation for the set of functions M, which satisfy the vector Dirichlet condition are


where Sx=sinkxx,Cx=coskxx, Sy=sinkyy,Cy=coskyy. The subscript “e” attached to Memn is an abbreviation for the word “even,” and “o” for “odd.”

In a similar manner


It is obvious that Memn represents the electric field of the TEmn mode, while Nomn represents that of the TMmn mode.

In summary, the vector wave functions, which can be used to represent the electromagnetic field inside a rectangular waveguide, are of the form




where a(h)=b(h)=1κ2k2, h=±(k2kx2ky2)1/2 and δ0={1,m=0orn=00,m0,n0.

M', N', m', n', h' denote another set of values, which may be distinct or the same as M, N, m, n, h.


6. Retarded and advanced Green functions

Green function is also utilized to solve the Schrödinger equation in quantum mechanics. Being completely equivalent to the Landauer scattering approach, the GF technique has the advantage that it calculates relevant transport quantities (e.g., transmission function) using effective numerical techniques. Besides, the Green function formalism is well adopted for atomic and molecular discrete-level systems and can be easily extended to include inelastic and many-body effects [13, 14].

(A) The definitions of propagators

The time-dependent Schrödinger equation is:


The solution of this equation at time t can be written in terms of the solution at time t′:


where U^(t,t) is called the time-evolution operator.

For the case of a time-independent Hermitian Hamiltonian H^, so that the eigenstates |Ψn(t)=eiEnt/ħ|Ψn with energies En are found from the stationary Schrödinger equation


The eigenfunctions |Ψn are orthogonal and normalized, for discrete energy levels 1:


and form a complete set of states (I^ is the unity operator)


The time-evolution operator for a time-independent Hamiltonian can be written as


This formal solution is difficult to use directly in most cases, but one can obtain the useful eigenstate representation from it. From the identity U^=U^I^ and (107), (109), (110) it follows that


which demonstrates the superposition principle. The wave function at time t is


where Ψn|Ψ(t) are the coefficients of the expansion of the initial function |Ψ(t) on the basis of eigenstates.

It is equivalent and more convenient to introduce two Green operators, also called propagators, retarded G^R(t,t) and advanced G^A(t,t):


so that at t > t′ one has


while at t < t′ it follows


The operators G^R(t,t) at t > t′and G^A(t,t) at t < t′ are the solutions of the equation


with the boundary conditions G^R(t,t)=0 at t < t′, G^A(t,t)=0 at t > t′. Indeed, at t > tEq. (118) satisfies the Schrödinger equation Eq. (105) due to Eq. (117). And integrating Eq. (117) from t=tη to t=t+η where η is an infinitesimally small positive number η=0+, one gets


giving correct boundary condition at t = t′. Thus, if the retarded Green operator G^R(t,t) is known, the time-dependent wave function at any initial condition is found (and makes many other useful things, as we will see below).

For a time-independent Hamiltonian, the Green function is a function of the time difference τ=t-t, and one can consider the Fourier transform


This transform, however, can not be performed in all cases, because G^R(A)(E) includes oscillating terms eiEτ/ħ. To avoid this problem we define the retarded Fourier transform


and the advanced one


where the limit η → 0 is assumed in the end of calculation. With this addition, the integrals are convergent. This definition is equivalent to the definition of a retarded (advanced) function as a function of complex energy variable at the upper (lower) part of the complex plain.

Applying this transform to Eq. (117), the retarded Green operator is


The advanced operator G^A(E) is related to the retarded one through


Using the completeness propertyn|ΨnΨn|=1, there is




Apply the ordinary inverse Fourier transform to G^R(E), the retarded function becomes


Indeed, a simple pole in the complex E plain is at E=Eniη, the residue in this point determines the integral at τ > 0 when the integration contour is closed through the lower half-plane, while at τ < 0 the integration should be closed through the upper half-plane and the integral is zero.

The formalism of retarded Green functions is quite general and can be applied to quantum systems in an arbitrary representation. For example, in the coordinate system Eq. (124) is


(B) Path integral representation of the propagator

In the path integral representation, each path is assigned an amplitude eidtL, L is the Lagrangian function. The propagator is the sum of all the amplitudes associated with the paths connecting xa and xb (Figure 6). Such a summation is an infinite-dimensional integral.

Figure 6.

The total amplitude is the sum of all amplitudes associated with thee paths connecting xa and xb.

The propagator satisfies


Let us divide the time interval [ta, tb] into N equal segments, each of length Δt=(tbta)/N.


where ln[iG(xj,tj,xj1,tj1)]=iΔtL(tj,xj+xj12,xjxj12).

Example: LC circuit-based metamaterials

In this section, we will use the relationship of current and voltage in the LC circuit to build the propagator of the LC circuit field coupled to an atom.

Figure 7 shows the LC-circuit.The following are valid:

Figure 7.

The coupled system, including an LC field and a bipole.




where x=LI, I is the current, V is the voltage, q is the charge quantity, L and C are the inductance and capacitance, respectively. Eq. (132) is equal to a harmonic, and the Lagrangian operator is:


The Lagrangian operator describing the bipole is:


where x is the coordinate of the bipole, ε is the LC field, m is the mass of an electron, and e is the unit of charge. g=1c, and ΩLC=1LC. Defining their action items as:




Taking the coupling effect (exε) into account, the Green function of the coupled system is:


Where x represents the series coordinates x1,x2,…,and so on and ɛ represents ɛ1,ɛ2,…., and so on.


7. The recent applications of the Green function method

7.1. Convergence

In the Green function, the high oscillation of Bessel/Hankel functions in the integrands results in quite time-consuming integrations along the Sommerfeld integration paths (SIP) which ensures that the integrands can satisfy the radiation condition in the direction normal to the interface of a medium. To facilitate the evaluation, the method of moments (MoM) [15], the steepest descent path (SDP) method, and the discrete complex image method (DCIM) [16, 17] are very important methods.

The technique for locating the modes is quite necessary for accurately calculating the spatial Green functions of a layered medium. The path tracking algorithm can obtain all the modes for the configuration shown in Figure 8, even when region 2 is very thick [18]. Like the method in Ref. [19], it does not involve a contour integration and could be extended to more complicated configurations.

Figure 8.

A general configuration with a three-layered medium: region 1 is free space, region 2 is a substrate with thickness h and relative permittivity ɛr1, and region 3 is a half space with relative permittivity ɛr2.

The discrete complex image method (DCIM) has been shown to deteriorate sharply for distances between source and observation points larger than a few wavelengths [20]. So, the total least squares algorithm (TLSA) is applied to the determination of the proper and improper poles of spectral domain multilayered Green’s functions that are closer to the branch point and to the determination of the residues at these poles [21].

The complex-plane for the determination of proper and improper poles is shown in Figure 9. Since half the ellipse is in the proper sheet of the -plane and half the ellipse is in the improper sheet, the poles will not only correctly capture the information of the proper poles but will also capture the information of those improper poles that are closer to the branch point kρ = k0.

Figure 9.

Elliptic path chosen in the complex -plane when applying the total least squares algorithm. The upper half ellipse (solid line) is located in the proper Riemman sheet, and the lower half ellipse (dashed line) is located in the improper sheet.

For the 2-D dielectric photonic crystals as shown in Figure 10, the integral equation is written in terms of the unknown equivalent current sources flowing on the surfaces of the periodic 2-D cylinders. The method of moments is then employed to solve for the unknown current distributions. The required Green function of the problem is represented in terms of a finite summation of complex images. It is shown that when the field-point is far from the periodic sources, it is just sufficient to consider the contribution of the propagating poles in the structure [22]. This will result in a summation of plane waves that has an even smaller size compared with the conventional complex images Green function. This provides an analyzed method for the dielectric periodic structures.

Figure 10.

Typical (a) waveguide and (b) directional coupler in a rectangular lattice.

Others, since the Gaussian function is an eigenfunction of the Hankel transform operator, for the microstrip structures, the spectral Green’s function can be expanded into a Gaussian series [23]. By introducing the mixed-form thin-stratified medium fast-multiple algorithm (MF-TSM-FMA), which includes the multipole expansion and the plane wave expansion in one multilevel tree, the different scales of interaction can be separated by the multilevel nature of the the fast multipole algorithm [24].

The vector wave functions, L, M, and N, are the solutions of the homogeneous vector Helmholtz equation. They can also be used for the analyses of the radiation in multilayer and this method avoids the finite integration in some cases.

7.2. Multilayer structure

The volume integral equation (VIE) can analyze electromagnetic radiation and scattering problems in inhomogeneous objects. By introducing an “impulse response” Green function, and invoking Green theorem, the Helmholtz equation can be cast into an equivalent volume integral equation including the source current or charges distribution. But the number of unknowns is typically large and the equation should be reformulated if there are in contrast both permittivity and permeability. At present, it is utilized to analyse the general scatterers in layered medium [25, 26].

When the inhomogeneity is one dimension, the Green function can be determined analytically in the spectral (Fourier) domain, and the spatial domain counterpart can be obtained by simply inverse Fourier transforming it.

Surface integral equation (SIE) method is another powerful method to handle electromagnetic problems. Similarly, by introducing the Green function, the Helmholtz equation can be cast into an equivalent surface integral equation, where the unknowns are pushed to the boundary of the scatterers [27].

Despite the convergence problem, the locations of the source and observation point may cause the change of Green function form, for example, for a source location either inside or outside the medium, the algebraic form of the Green functions changes as the receiver moves vertically in the direction of stratification from one layer to another [28].

First, we introduce the full-wave computational model [29]. A multilayer structure involving infinitely 1-D periodic chains of parallel circular cylinders in any given layer can be constructed as shown in Figure 11. Each layer consists of a homogeneous slab within which the circular cylinders are embedded. This is the typical aeronautic situation with fiber-reinforced four-layer pile (with fibers orientated at 0°, 45°, −45°, and 90°), but any other arrangement is manageable likewise.

Figure 11.

(a) Sketch of a standard (0, 45, −45, 90) degree, four-layer fiber-reinforced composite laminate as in aeronautics. (b) General two-layer pile of interest exhibiting two different cylinder orientations and associated coordinate systems with geometrical parameters as indicated. (c) Cell defined in the lth layer of multilayered photonic crystals.

In the multilayered photonic crystals, the Rayleig’s method and mode-matching are combined to produce scattering matrices. An S-matrix-based recursive matrix is developed for modeling electromagnetic scattering. Field expansions and the relationship between expansion coefficients are given.

There is a mix treatment for the inhomogeneous and homogeneous multilayered structure [30]. As shown in Figure 12, a substrate is divided into two regions. The top region is laterally inhomogeneous and for the finite-difference method (FDM) or the finite element method (FEM), the volume integral equation, is used. The bottom region is layerwise homogeneous, and the boundary-element methods (BEM) are used. The two regions are connected such as a BEM panel is associated with an FEM node on the interface.

Figure 12.

Substrate is divided into homogeneous and inhomogeneous regions in combined BEM/FEM and BEM/FDM methods.

A Green function was derived for a layerwise uniform substrate and was then used in a layerwise nonuniform substrate with additional boundary conditions applied to the interface. Given that the lateral inhomogeneity is local, volume meshing is used only for the local inhomogeneous regions, BEM meshing is applied to the surfaces of these local regions.

For a field (observation) point in the jth layer and a source point in the kth layer, the Green function has the form:


where the superscripts u and l indicate the upper and lower solutions, respectively, depending on whether the field point (or observation point) is above or below the source point. a and b are the substrate dimensions in the x- and y- directions, respectively, and more details can be found in Refs. [31, 32].

The electromagnetic field in a multilayer structure can be efficiently simplified by the assumption that the multilayer is grounded by a perfect electric conducto (PEC) plane [33, 34]. When the source and the field points are assumed to be inside the dielectric slab, in a layered medium as shown in Figure 13, by applying the boundary conditions, the 1-D Green functions is

Figure 13.

(a) Geometry of an infinite dielectric slab of thickness d grounded by a PEC plane at x = d. (b) Geometry of a finite dielectric slab of thickness 2d and height 2L surrounded by regions □ and □.


where PMC represents the perfect magnetic conductor. The simplified Green function form can be deduced to the cae of (b).

The three-dimensional (3-D) Green function for a continuous, linearly stratified planar media, backed by a PEC ground plane, can also be expressed in terms of a single contour integral involving one-dimensional (1-D) green function. The constructure is shown in Figure 14.

Figure 14.

Representation of the continuous, linearly stratified media by discrete slabs of finite thickness and constant permittivity, ɛp and permeability μp for the pth layer of thickness hp. The thicknesses, permittivities and permeabilities are different for each layer.

The general formulation for a single electric current element has been worked out in detail in Ref. [35] which is based on the appropriate information from Ref. [36].


  1. 1. Abramowitz M, Stegun I. Handbook of Mathematical Functions. New York: Dover Publications; 1965
  2. 2. Chew WC. Waves and Fields in Inhomogenous Media. New York: Wiley-IEEE Press; January 1999
  3. 3. Tai C-T. Dyadic Green Functions in Electromagnetic Theory. 2nd ed. New York: IEEE press; 1994
  4. 4. Nussbaumer HJ. Fast Fourier Transform and Convolution Algorithms. New York: Springer Verlag; 1982
  5. 5. Huang J, Yao J, Xu D. Green function method for the time domain simulation of pulse propagation. Applied Optics. 2014;53:3533–3539
  6. 6. Tai CT. Complementary reciprocity theorems in electromagnetic theory. IEEE Transactions on Antennas and Propagation. 1992;40:675481
  7. 7. Sommerfeld A. Partial Differential Equation. New York: Academic Press; 1949
  8. 8. Kong JA. Electromagnetic Wave Theory. New York: John Wiley & Sons; 1986
  9. 9. Tai CT. Dyadic Green functions for a rectangular waveguide filled with two dielectrics. Journal of Electromagnetic Waves and Applied. 1988;2:245–253
  10. 10. Hansen WW. A new me of expansion in radiation problems. Physical Review. 1935;47:139–143
  11. 11. Hansen WW. Directional characteristics of any antenna over a plane earth. Journal of Applied Physics. 1936;7:460465
  12. 12. Hansen WW. Transformations useful in certain antenna calculations. Journal of Applied Physics. 1937;8:282–286
  13. 13. Di Ventra M. Electrical Transport in Nanoscale Systems. Cambridge: Cambridge University Press; 2008
  14. 14. Ferry DK, Goodnick SM. Transport in Nanostructures. Cambridge: Cambridge University Press; 1997
  15. 15. Michalski KA, Mosig JR. Multilayered media Green’s functions in integral functions. IEEE Transactions on Antennas and Propagation. 1997;45:508–519
  16. 16. Chow YL, Yang JJ, Fang DG, Howard GE. A closed-form spatial Green’s function for the thick microstrip substrate. IEEE Transactions on Microwave Theory and Techniques. 1991;39:588–592
  17. 17. Aksun MI. A robust approach for the derivation of closed-form Green’s functions. IEEE Transactions on Microwave Theory and Techniques. 1996;44:651–658
  18. 18. Hu J, Zhou HX, Song Z, Hong W. Locating all the modes of Green’s function for a three-layered medium based on the path tracking algorithm. IEEE Transactions on Antennas and Propagation. 2009,57:2315–2322
  19. 19. Tsang L, Wu BP. Electromagnetic fields of Hertzian dipoles in layered media of moderate thickness including the effects of all modes. IEEE Antennas and Wireless Propagation Letters. 2007;6:316–319
  20. 20. Shuley NV, Boix RR, Medina F, Horno M. On the fast approximation of Green’s functions in MPIE formulations for planar layered media. IEEE Transactions on Microwave Theory and Techniques. 2002;50:2185–2192
  21. 21. Fructos A L, Boix RR, Raúl Rodríguez-Berral, Mesa F. Efficient determination of the poles and residues of spectral domain multilayered Green’s functions that are relevant in far-field calculations. IEEE Transactions on Antennas and Propagation. 2010;58:218–222
  22. 22. Ameri H, Faraji-Dana R. Green’s function analysis of electromagnetic wave propagation in photonic crystal devices using complex images technique. IEEE Journal of Lightwave Technolgy. 2011;29:298–304
  23. 23. Tajdini MM, Shishegar AA. A novel analysis of microstrip structures using the Gaussian Green’s function method. IEEE Transactions on Antennas and Propagation. 2010;58:88–94
  24. 24. Chen YP, Xiong JL, Chew WC. A mixed-form thin-stratified medium fast-multipole algorithm for both low and mid-frequency problems, IEEE Transactions on Antennas and Propagation. 2011,59:2341–2349
  25. 25. Cui TJ, Chew WC. Fast algorithm for electromagnetic scattering by buried 3-D dielectric objects of large size. IEEE Transactions on Geoscience and Remote Sensing. 1999;37:2597–2608
  26. 26. Millard X, Liu QH. Simulation of near-surface detection of objects in layeredmedia by the BCGS-FFTmethod, IEEE Transactions on Geoscience and Remote Sensing. 2004;42:327–334
  27. 27. Chew WC. Waves and Fields in Inhomogeneous Media. Berlin, Germany: Van Nostrand Reinhold; 1990, Reprinted by New York: IEEE Press; 1995
  28. 28. Chatterjee D. Phase-integral formulation of the single-layer microstrip Green’s function. 2010 URSI International Symposium on Electromagnetic Theory. Berlin: 2010. pp. 970–973
  29. 29. Li C, Lesselier D, Zhong Y. Full-wave computational model of electromagnetic scattering by arbitrarily rotated 1-D periodic multilayer structure. IEEE Transactions on Antennas and Propagation. 2016;64:1047–1060
  30. 30. Xu C, Gharpurey R, Fiez TS, Mayaram K. Extraction of parasitics in inhomogeneous substrates with a new Green function-based method. IEEE Transactions on Computer-Aided Design of Integrated Circuits Systems. 2008;27:1595–1606
  31. 31. Niknejad AM, Gharpurey R, Meyer RG. Numerically stable Green function for modeling and analysis of substrate coupling in integrated circuits. IEEE Transactions on Computer-Aided Design of Integrated Circuits Systems.1998;17:305–315
  32. 32. Xu C, Fiez T, Mayaram K. On the numerical stability of Green’s function for substrate coupling in integrated circuits. IEEE Transactions on Computer-Aided Design of Integrated Circuits Systems. 2005;24:653–658
  33. 33. Parsa A, Paknys R. Interior Green’s function solution for a thick and finite dielectric slab. IEEE Transactions on Antennas and Propagations. 2007;55:3504–3514
  34. 34. Rogier H. New series expansions for the 3-D Green’s function of multilayered media with 1-D periodicity based on perfectly matched layers. IEEE Transactions on Microwave Theory and Techniques. 2007;55:1730–1738
  35. 35. Chatterjee D, Walker SD. Study of Sommerfeld and phase integral approaches for Green’s functions for PEC-terminated inhomogeneous media, technical memorandum NRL/MR/5310-10-9240, January 2010, Naval Research Laboratory, Washington, DC, USA.
  36. 36. Felsen LB, Marcuvuitz N. Radiation and Scattering of Waves, IEEE Classic Reissue. New York, USA: IEEE-Wiley Press; 1994

Written By

Jing Huang

Submitted: 18 October 2016 Reviewed: 21 February 2017 Published: 14 June 2017