InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Electrical and Electronic Engineering » "Solutions and Applications of Scattering, Propagation, Radiation and Emission of Electromagnetic Waves", book edited by Ahmed Kishk, ISBN 978-953-51-0838-2, Published: November 14, 2012 under CC BY 3.0 license. © The Author(s).

# Analytical Grounds for Modern Theory of Two- Dimensionally Periodic Gratings

By L. G. Velychko, Yu. K. Sirenko and E. D. Vinogradova
DOI: 10.5772/51007

Article top

## Overview

Figure 1. Geometry of a two-dimensionally periodic grating.

Figure 2. A two-dimensionally periodic grating between two dielectric half-spaces as element of a multi-layered structure.

Figure 3. Construction of sets of the values v p ( 1 ) , u p ( j ) , and λ p ( p = 0, ± 1, ± 2,... ) from sets of the values v n m ( 1, E ) , u n m ( j , E ) , v n m ( 1, H ) , u n m ( j , H ) , and λ n m ( m , n = 0, ± 1, ± 2,... ): (a) p = 0,1,2,... ; (b) p = − 1, − 2, − 3,... .

Figure 4. Schematic drawing of a double-layered structure.

Figure 5. On determination of propagation directions for spatial harmonics of the field formed by a two-dimensionally periodic structure.

Figure 6. Natural domain of variation of the spectral parameter k : the first sheet of the surface K .

# Analytical Grounds for Modern Theory of Two-Dimensionally Periodic Gratings

L. G. Velychko1, Yu. K. Sirenko1 and E. D. Vinogradova2

## 1. Introduction

Rigorous models of one-dimensionally periodic diffraction gratings made their appearance in the 1970s, when the corresponding theoretical problems had been considered in the context of classical mathematical disciplines such as mathematical physics, computational mathematics, and the theory of differential and integral equations. Periodic structures are currently the objects of undiminishing attention. They are among the most called-for dispersive elements providing efficient polarization, frequency and spatial signal selection. Fresh insights into the physics of wave processes in diffraction gratings are being implemented into radically new devices operating in gigahertz, terahertz, and optical ranges, into new materials with inclusions ranging in size from micro- to nanometers, and into novel circuits for in-situ man-made and natural material measurements.

However, the potentialities of classical two-dimensional models [1-7] are limited. Both theory and applications invite further investigation of three-dimensional, vector models of periodic structures in increasing frequency. In our opinion these models should be based on time-domain (TD) representations and implemented numerically by the mesh methods [8,9]. It follows from the well-known facts: (i) TD-approaches are free from the idealizations inherent in the frequency domain; (ii) they are universal owing to minimal restrictions imposed on geometrical and material parameters of the objects under study; (iii) they allow explicit computational schemes, which do not require inversion of any operators and call for an adequate run time when implementing on present-day computers; (iv) they result in data easy convertible into a standard set of frequency-domain characteristics. To this must be added that in recent years the local and nonlocal exact absorbing conditions (EAC) have been derived and tested [6,7]. They allow one to replace an open initial boundary value problem occurring in the electrodynamic theory of gratings with a closed problem. In addition, the efficient fast Fourier transform accelerated finite-difference schemes with EAC for characterizing different resonant structures have been constructed and implemented [10].

It is evident that the computational scheme solving a grating problem must be stable and convergent, computational error must be predictable, while the numerical results are bound to be unambiguously treated in physical terms. To comply with these requirements, it is important to carry out theoretical analysis at each stage of the modeling (formulation of boundary value and initial boundary value problems, determination of the correctness classes for them, study of qualitative characteristics of singularities of analytical continuation for solutions of model boundary value problems into a domain of complex-valued frequencies, etc.).

In the present work, we present a series of analytical results providing the necessary theoretical background to the numerical solution of initial boundary value problems as applied to two-dimensionally periodic structures. Section 1 is an Introduction. In Section 2 we give general information required to formulate a model problem in electrodynamic theory of gratings. Sections 3 and 4 are devoted to correct and efficient truncation of the computational space in the problems describing spatial-temporal electromagnetic wave transformation in two-dimensionally periodic structures. Some important characteristics and properties of transient and steady-state fields in regular parts of the rectangular Floquet channel are discussed in Sections 5 and 7. In Section 6, the method of transformation operators (the TD-analog of the generalized scattering matrix method) is described; by applying this method the computational resources can be optimized when calculating a multi-layered periodic structure or a structure on a thick substrate. In Section 8, elements of spectral theory for two-dimensionally periodic gratings are given in view of its importance to physical analysis of resonant scattering of pulsed and monochromatic waves by open periodic resonators.

## 2. Fundamental Equations, Domain of Analysis, Initial and Boundary Conditions

Space-time and space-frequency transformations of electromagnetic waves in diffraction gratings, waveguide units, open resonators, radiators, etc. are described by the solutions of initial boundary value problems and boundary value problems for Maxwell’s equations. In this chapter, we consider the problems of electromagnetic theory of gratings resulting from the following system of Maxwell’s equations for waves propagating in stationary, locally inhomogeneous, isotropic, and frequency dispersive media [9,11]:

 rotH→(g,t)=η0−1∂[E→(g,t)+χε(g,t)∗E→(g,t)]∂t+χσ(g,t)∗E→(g,t)+j→(g,t), (1)
 rotE→(g,t)=−η0∂[H→(g,t)+χμ(g,t)∗H→(g,t)]∂t, (2)

where

g={x,y,z} is the point in a three-dimensional space R3 ;

x , y , and z are the Cartesian coordinates;

E(g,t)={Ex,Ey,Ez} and H(g,t)={Hx,Hy,Hz} are the electric and magnetic field vectors;

η0=(μ0/ε0)1/2 is the intrinsic impedance of free space;

ε0 and μ0 are permittivity and permeability of free space;

j(g,t) is the extraneous current density vector;

χε(g,t) , χμ(g,t) , and χσ(g,t) are the electric, magnetic, and specific conductivity susceptibilities; f1(t)f2(t)=f1(tτ)f2(τ)dτ stands for the convolution operation.

We use the SI system of units. From here on we shall use the term “time” for the parameter t , which is measured in meters, to mean the product of the natural time and the velocity of light in vacuum.

With no frequency dispersion in the domain GR3 , for the points gG we have

χε(g,t)=δ(t)[ε(g)1],χμ(g,t)=δ(t)[μ(g)1],χσ(g,t)=δ(t)σ(g),

where δ(t) is the Dirac delta-function; ε(g) , μ(g) , and σ(g) are the relative permittivity, relative permeability, and specific conductivity of a locally inhomogeneous medium, respectively. Then equations (1) and (2) take the form:

 rotH→(g,t)=η0−1ε(g)∂E→(g,t)∂t+σ(g)E→(g,t)+j→(g,t), (3)
 rotE→(g,t)=−η0μ(g)∂H→(g,t)∂t. (4)

In vacuum, where ε(g)=μ(g)=1 and σ(g)=0 , they can be rewritten in the form of the following vector problems [6]:

or

 [Δ−∂2∂t2]H→(g,t)=F→H(g,t), η0−1∂∂tE→(g,t)=rotH→(g,t)−j→(g,t),F→H(g,t)=−rotj→(g,t). (6)

By Δ we denote the Laplace operator. As shown in [6], the operator graddivE can be omitted in (5) from the following reasons. By denoting the volume density of induced and external electric charge through ρ1(g,t) and ρ2(g,t) , we can write graddivE=ε01grad(ρ1+ρ2) . In homogeneous medium, where ε and σ are positive and non-negative constants, we have ρ1(g,t)=ρ1(g,0)exp(tσ/ε) , and if ρ1(g,0)=0 , then ρ1(g,t)=0 for any t>0 . The remaining term ε01gradρ2 can be moved to the right-hand side of (5) as a part of the function defining current sources of the electric field.

To formulate the initial boundary value problem for hyperbolic equations (1)-(6) [12], initial conditions at t=0 and boundary conditions on the external and internal boundaries of the domain of analysis Q should be added. In 3-D vector or scalar problems, the domain Q is a part of the R3 -space bounded by the surfaces S that are the boundaries of the domains intS , filled with a perfect conductor: Q=R3\intS¯ . In the so-called open problems, the domain of analysis may extend to infinity along one or more spatial coordinates.

The system of boundary conditions for initial boundary value problems is formulated in the following way [11]:

• on the perfectly conducting surface S the tangential component of the electric field vector is zero at all times t

 Εtg(g,t)|g∈S=0 for t≥0; (7)

the normal component of the magnetic field vector on S is equal to zero ( Hnr(g,t)|gS=0 ), and the function Htg(g,t)|gS defines the so-called surface currents generated on S by the external electromagnetic field;

• on the surfaces Sε,μ,σ , where material properties of the medium have discontinuities, as well as all over the domain Q , the tangential components Etg(g,t) and Htg(g,t) of the electric and magnetic field vectors must be continuous;

• in the vicinity of singular points of the boundaries of Q , i.e. the points where the tangents and normals are undetermined, the field energy density must be spatially integrable;

• if the domain Q is unbounded and the field {E(g,t),H(g,t)} is generated by the sources having bounded supports in Q then for any finite time interval (0,T) one can construct a closed virtual boundary MQ sufficiently removed from the sources such that

 {E→(g,t),H→(g,t)}|g∈M,t∈(0,T)=0. (8)

The initial state of the system is determined by the initial conditions at t=0 . The reference states E(g,0) and H(g,0) in the system (1), (2) or the system (3), (4) are the same as E(g,0) and [E(g,t)/t]|t=0 ( H(g,0) and [H(g,t)/t]|t=0 ) in the differential forms of the second order (in the terms of t ), to which (1), (2) or (3), (4) are transformed if the vector H (vector E ) is eliminated (see, for example, system (5), (6)). Thus, (5) should be complemented with the initial conditions

 Ε→(g,0)=φ→(g),∂∂tΕ→(g,t)|t=0=ψ→(g),g∈Q¯. (9)

The functions φ(g) , ψ(g) , and F(g,t) (called the instantaneous and current source functions) usually have limited support in the closure of the domain Q . It is the practice to divide current sources into hard and soft [9]: soft sources do not have material supports and thus they are not able to scatter electromagnetic waves. Instantaneous sources are obtained from the pulsed wave Ui(g,t) exciting an electrodynamic structure: φ(g)=Ui(g,0) and ψ(g)=[Ui(g,t)/t]|t=0 . The pulsed signal Ui(g,t) itself should satisfy the corresponding wave equation and the causality principle. It is also important to demand that the pulsed signal has not yet reached the scattering boundaries by the moment t=0 .

The latter is obviously impossible if infinite structures (for example, gratings) are illuminated by plane pulsed waves that propagate in the direction other than the normal to certain infinite boundary. Such waves are able to run through a part of the scatterer’s surface by any moment of time. As a result a mathematically correct modeling of the process becomes impossible: the input data required for the initial boundary value problem to be set are defined, as a matter of fact, by the solution of this problem.

## 3. Time Domain: Initial Boundary Value Problems

The vector problem describing the transient states of the field nearby the gratings whose geometry is presented in Figure 1 can be written in the form

 {rotH→(g,t)=η0−1∂[E→(g,t)+χε(g,t)∗E→(g,t)]∂t+χσ(g,t)∗E→(g,t)+j→(g,t),rotE→(g,t)=−η0∂[H→(g,t)+χμ(g,t)∗H→(g,t)]∂t, g={x,y,z}∈Q, t>0Ε→(g,0)=φ→E(g), H→(g,0)=φ→H(g), g∈Q¯Εtg(g,t)|g∈S=0, Hnr(g,t)|g∈S=0, t≥0 . (10)

Here, Q¯ is the closure of Q , χε,μ,σ(g,t) are piecewise continuous functions and the surfaces S are assumed to be sufficiently smooth. From this point on it will be also assumed that the continuity conditions for tangential components of the field vectors are satisfied, if required. The domain of analysis Q=R3\intS¯ occupies a great deal of the R3 -space. The problem formulated for that domain can be resolved analytically or numerically only in two following cases.

#### Figure 1.

Geometry of a two-dimensionally periodic grating.

• The problem (10) degenerates into a conventional Cauchy problem ( intS¯= , the medium is homogeneous and nondispersive, while the supports of the functions F(g,t) , φ(g) , and ψ(g) are bounded). With some inessential restrictions for the source functions, the classical and generalized solution of the Cauchy problem does exist; it is unique and is described by the well-known Poisson formula [12].

• The functions F(g,t) , φ(g) , and ψ(g) have the same displacement symmetry as the periodic structure. In this case, the domain of analysis can be reduced to QN={gQ:0<x<lx;0<y<ly} , by adding to problem (10) periodicity conditions [7] on lateral surfaces of the rectangular Floquet channel R={gR3:0<x<lx;0<y<ly} .

The domain of analysis can also be reduced to QN in a more general case. The objects of analysis are in this case not quite physical (complex-valued sources and waves). However, by simple mathematical transformations, all the results can be presented in the customary, physically correct form. There are several reasons (to one of them we have referred at the end of Section 3) why the modeling of physically realizable processes in the electromagnetic theory of gratings should start with the initial boundary value problems for the images fN(g,t,Φx,Φy) of the functions f(g,t) describing the actual sources:

 f(g,t)=∫−∞∞∫−∞∞f˜(z,t,Φx,Φy)exp(2πiΦxxlx)exp(2πiΦyyly)dΦxdΦy=∫−∞∞∫−∞∞fN(g,t,Φx,Φy)dΦxdΦy (11)

From (11) it follows that

fN{fNx}(x+lx,y,z,t,Φx,Φy)=e2πiΦxfN{fNx}(x,y,z,t,Φx,Φy) ,

fN{fNy}(x,y+ly,z,t,Φx,Φy)=e2πiΦyfN{fNy}(x,y,z,t,Φx,Φy)

or, in other symbols,

D[fN](x+lx,y)=e2πiΦxD[fN](x,y) , D[fN](x,y+ly)=e2πiΦyD[fN](x,y) .

The use of the foregoing conditions truncates the domain of analysis to the domain QN , which is a part of the Floquet channel R , and allows us to rewrite problem (10) in the form

 E→(g,t)=∫−∞∞∫−∞∞E→N(g,t,Φx,Φy)dΦxdΦy,H→(g,t)=∫−∞∞∫−∞∞H→N(g,t,Φx,Φy)dΦxdΦy (12)

and

 {rotH→N(g,t)=η0−1∂[E→N(g,t)+χε(g,t)∗E→N(g,t)]∂t+χσ(g,t)∗E→N(g,t)+j→N(g,t),rotE→N(g,t)=−η0∂[H→N(g,t)+χμ(g,t)∗H→N(g,t)]∂t,g∈QN,t>0Ε→N(g,0)=φ→EN(g),H→N(g,0)=φ→HN(g),g∈Q¯ND[E→N(H→N)](lx,y)=e2πiΦxD[E→N(H→N)](0,y),0≤y≤lyD[E→N(H→N)](x,ly)=e2πiΦyD[E→N(H→N)](x,0),0≤x≤lxΕtgN(g,t)|g∈S=0,HnrN(g,t)|g∈S=0,t≥0. (13)

It is known [6-8] that initial boundary value problems for the above discussed equations can be formulated such that they are uniquely solvable in the Sobolev space W21(QT) , where QT=Q×(0,T) and 0tT . On this basis we suppose in the subsequent discussion that the problem (13) for all t[0,T] has also a generalized solution from the space W21(QN,T) and that the uniqueness theorem is true in this space. Here symbol × stands for the operation of direct product of two sets, (0,T) and [0,T] are open and closed intervals, Wmn(G) is the set of all elements f(g) from the space Lm(G) whose generalized derivatives up to the order n inclusive also belong to Lm(G) . Lm(G) is the space of the functions f(g)={fx,fy,fz} (for gG ) such that the functions |f...(g)|m are integrable on the domain G .

## 4. Exact Absorbing Conditions for the Rectangular Floquet Channel

In this section, we present analytical results relative to the truncation of the computational space in open 3-D initial boundary value problems of the electromagnetic theory of gratings. In Section 3, by passing on to some special transforms of the functions describing physically realizable sources, the problem for infinite gratings have been reduced to that formulated in the rectangular Floquet channel R or, in other words, in the rectangular waveguide with quasi-periodic boundary conditions. Now we perform further reduction of the domain QN to the region QLN={gQN:|z|<L} (all the sources and inhomogeneities of the Floquet channel R are supposedly located in this domain). For this purpose the exact absorbing conditions [6,7,10,13,14] for the artificial boundaries L± ( z=±L ) of the domain QLN will be constructed such that their inclusion into (13) does not change the correctness class of the problem and its solution EN(g,t) , HN(g,t) .

From here on we omit the superscripts N in (13). By applying the technique similar to that described in [13,14], represent the solution E(g,t) of (13) in the closure of the domains A={gR:z>L} and B={gR:z<L} in the following form:

 E→(g,t)=∑n,m=−∞∞u→nm±(z,t)μnm(x,y),{x,y}∈R¯z,t≥0, (14)

where the superscript ‘ + ’ corresponds to zL and ‘ ’ to zL and the following notation is used:

Rz=(0<x<lx)×(0<y<ly) ;

{μnm(x,y)} ( n,m=0,±1,±2,... ) is the complete in L2(Rz) orthonormal system of the functions μnm(x,y)=(lxly)1/2exp(iαnx)exp(iβmy) ;

αn=2π(Φx+n)/lx , βm=2π(Φy+m)/ly , and λnm2=αn2+βm2 .

The space-time amplitudes unm±(z,t) satisfy the equations

 {[−∂2∂t2+∂2∂z2−λnm2]u→nm±(z,t)=0,t>0u→nm±(z,0)=0,∂∂tu→nm±(z,t)|t=0=0,{z≥Lz≤−L}. (15)

Equations (14) and (15) are obtained by separating variables in the homogeneous boundary value problems for the equation [Δ2/t2]E(g,t)=0 (see formula (5)) and taking into account that in the domains A and B we have graddivE(g,t)=0 and FE(g,t)=0 . It is also assumed that the field generated by the current and instantaneous sources located in QL has not yet reached the boundaries L± by the moment of time t=0 .

The solutions unm±(z,t) of the vector problems (15), as well as in the case of scalar problems [13,14], can be written as

 u→nm±(±L,t)=∓∫0tJ0[λnm(t−τ)]u→nm±′(±L,τ)dτ,t≥0 . (16)

The above formula represents nonlocal EAC for the space-time amplitudes of the field E(g,t) in the cross-sections z=±L of the Floquet channel R . The exact nonlocal and local absorbing conditions for the field E(g,t) on the artificial boundaries L± follow immediately from (16) and (14):

 E→(x,y,±L,t)=∓∑n,m=−∞∞{∫0tJ0[λnm(t−τ)][∫0lx∫0ly∂E→(x˜,y˜,z,τ)∂z|z=±Lμnm∗(x˜,y˜)dx˜dy˜]dτ}μnm(x,y),{x,y}∈R¯z, t≥0 (17)

and

 E→(x,y,±L,t)=2π∫0π/2∂W→E(x,y,t,φ)∂tdφ,{x,y}∈R¯z,t≥0{[∂2∂t2−sin2φ(∂2∂x2+∂2∂y2)]W→E(x,y,t,φ)=∓∂E→(g,t)∂z|z=±L,{x,y}∈Rz,t>0W→E(x,y,t,φ)|t=0=∂W→E(x,y,t,φ)∂t|t=0=0,{x,y}∈R¯zD[W→E](lx,y)=e2πiΦxD[W→E](0,y),0≤y≤lyD[W→E](x,ly)=e2πiΦyD[W→E](x,0),0≤x≤lx,t≥0. (18)

Here, unm±(±L,τ)=unm±(z,τ)/z|z=±L , J0(t) is the zero-order Bessel function, the superscript ‘ ’ stands for the complex conjugation operation, WE(x,y,t,φ) is some auxiliary function, where the numerical parameter φ lies in the range 0φπ/2 .

It is obvious that the magnetic field vector H(g,t) of the pulsed waves U(g,t)={E(g,t),H(g,t)} outgoing towards the domains A and B satisfies similar boundary conditions on L± . The boundary conditions for E(g,t) and H(g,t) (nonlocal or local) taken together reduce the computational space for the problem (13) to the domain QL (a part of the Floquet channel R ) that contains all the sources and obstacles.

Now suppose that in addition to the sources j(g,t) , φE(g) , and φH(g) , there exist sources jA(g,t) , φEA(g) , and φHA(g) located in A and generating some pulsed wave Ui(g,t)={Ei(g,t),Hi(g,t)} being incident on the boundary L+ at times t>0 . The field Ui(g,t) is assumed to be nonzero only in the domain A . Since the boundary conditions (17), (18) remain valid for any pulsed wave outgoing through L± towards z=± [13,14], then the total field {E(g,t),H(g,t)} is the solution of the initial boundary value problem (13) in the domain QL with the boundary conditions (17) or (18) on L and the following conditions on the artificial boundary L+ :

 E→s(x,y,L,t)=−∑n,m=−∞∞{∫0tJ0[λnm(t−τ)][∫0lx∫0ly∂E→s(x˜,y˜,z,τ)∂z|z=Lμnm∗(x˜,y˜)dx˜dy˜]dτ}μnm(x,y),{x,y}∈R¯z, t≥0 (19)

or

 E→s(x,y,L,t)=2π∫0π/2∂W→E(x,y,t,φ)∂tdφ,{x,y}∈R¯z,t≥0​​{[∂2∂t2−sin2φ(∂2∂x2+∂2∂y2)]W→E(x,y,t,φ)=−∂E→s(g,t)∂z|z=L,{x,y}∈Rz,t>0W→E(x,y,t,φ)|t=0=∂W→E(x,y,t,φ)∂t|t=0=0,{x,y}∈R¯zD[W→E](lx,y)=e2πiΦxD[W→E](0,y),0≤y≤lyD[W→E](x,ly)=e2πiΦyD[W→E](x,0),0≤x≤lx,t≥0. (20)

Here Us(g,t)={Es(g,t),Hs(g,t)}=U(g,t)Ui(g,t) ( gA , t>0 ) is the pulsed wave outgoing towards z=+ . It is generated by the incident wave Ui(g,t) (‘reflection’ from the virtual boundary L+ ) and the sources j(g,t) , φE(g) , and φH(g) .

## 5. Some Important Characteristics of Transient Fields in the Rectangular Floquet Channel

For numerical implementation of the computational schemes involving boundary conditions like (19) or (20), the function Ui(g,t) for t[0,T] and its normal derivative with respect to the boundary L+ are to be known. To obtain the required data for the wave Ui(g,t) generated by a given set of sources jA(g,t) , φEA(g) , and φHA(g) , the following initial boundary value problem for a regular hollow Floquet channel R are to be solved:

The function ρ2A(g,t) here determines the volume density of foreign electric charge.

First we determine the longitudinal components Ezi and Hzi of the field {Ei,Hi} at all points g of the domain R for all times t>0 . Let us consider the scalar initial boundary value problems following from (21):

 {[−∂2∂t2+Δ]{EziHzi}={Fz,EAFz,HA},g∈R,t>0{Ezi(g,0)Hzi(g,0)}={φz,EAφz,HA},{∂Ezi(g,t)/∂t|t=0∂Hzi(g,t)/∂t|t=0}={ψz,EAψz,HA},g∈R¯D[Ezi(Hzi)](lx,y)=e2πiΦxD[Ezi(Hzi)](0,y),0≤y≤lyD[Ezi(Hzi)](x,ly)=e2πiΦyD[Ezi(Hzi)](x,0),0≤x≤lx,t≥0. (22)

By separating of the transverse variables x and y in (22) represent the solution of the problem as

 {Ezi(g,t)Hzi(g,t)}=∑n,m=−∞∞{vnm(z,E)(z,t)vnm(z,H)(z,t)}μnm(x,y) (23)

To determine the scalar functions vnm(z,E)(z,t) and vnm(z,H)(z,t) , we have to invert the following Cauchy problems for the one-dimensional Klein-Gordon equations:

 {[−∂2∂t2+∂2∂z2−λnm]{vnm(z,E)(z,t)vnm(z,H)(z,t)}={Fnm(z,E)AFnm(z,H)A},t>0,−∞

Here Fnm(z,E)A , φnm(z,E)A , ψnm(z,E)A and Fnm(z,H)A , φnm(z,H)A , ψnm(z,H)A are the amplitudes of the Fourier transforms of the functions Fz,EA , φz,EA , ψz,EA and Fz,HA , φz,HA , ψz,HA in the basic set {μnm(x,y)} .

Let us continue analytically the functions vnm(z,E)(z,t) , vnm(z,H)(z,t) and Fnm(z,E)A , Fnm(z,H)A by zero on the semi-axis t0 and pass on to the generalized formulation of the Cauchy problem (24) [12]:

 B(λnm)[vnm(z,E)(z,t)vnm(z,H)(z,t)]≡[−∂2∂t2+∂2∂z2−λnm2]{vnm(z,E)(z,t)vnm(z,H)(z,t)}={Fnm(z,E)AFnm(z,H)A}−δ(1)(t){φnm(z,E)Aφnm(z,H)A}−δ(t){ψnm(z,E)Aψnm(z,H)A}={fnm(z,E)fnm(z,H)},−∞

where δ(t) and δ(m)(t) are the Dirac delta-function and its derivative of the order m . Taking into account the properties of the fundamental solution G(z,t,λ)=(1/2)χ(t|z|)J0(λt2z2) of the operator B(λ) [6,13,14] ( χ(t) is the Heaviside step function), the solutions vnm(z,E)(z,t) and vnm(z,H)(z,t) of equations (25) can be written as

 {vnm(z,E)(z,t)vnm(z,H)(z,t)}=G(z,t,λnm)∗{fnm(z,E)fnm(z,H)}=−12∫−∞t−|z−ω|∫−∞∞J0(λnm(t−τ)2−(z−ω)2)({Fnm(z,E)AFnm(z,H)A}−δ(1)(t){φnm(z,E)Aφnm(z,H)A}−δ(t){ψnm(z,E)Aψnm(z,H)A})dωdτ,−∞

Relations (23) and (26) completely determine the longitudinal components of the field {Ei,Hi} .

Outside the bounded domain enclosing all the sources, in the domain GR , where the waves generated by these sources propagate freely, the following relations [6,14] are valid:

 {E→i=(∂2UE∂x∂z−∂2UH∂y∂t)x→+(∂2UE∂y∂z+∂2UH∂x∂t)y→+(∂2UE∂z2−∂2UE∂t2)z→η0H→i=(∂2UE∂y∂t+∂2UH∂x∂z)x→+(−∂2UE∂x∂t+∂2UH∂y∂z)y→+(∂2UH∂z2−∂2UH∂t2)z→ (27)

in which

 UE,H(g,t)=∑n,m=−∞∞unmE,H(z,t)μnm(x,y) (28)

are the scalar Borgnis functions such that [Δ2/t2][UE,H(g,t)/t]=0 . Equations (23), (26)-(28) determine the field {Ei,Hi} at all points g of the domain G for all times t>0 . Really, since at the time point t=0 the domain G is undisturbed, then we have [Δ2/t2]UE,H=0 ( gG , t>0 ). Hence, in view of (27), (28), it follows:

Ez=2UEz22UEt2=(2UEx2+2UEy2)=n,m=λnm2unmEμnm,

η0Hz=2UHz22UHt2=(2UHx2+2UHy2)=n,m=λnm2unmHμnm

and (see representation (23))

 unmE(z,t)=(λnm)−2vnm(z,E)(z,t),unH(z,t)=η0(λnm)−2vnm(z,H)(z,t). (29)

Hence the functions UE,H(g,t) as well as the transverse components of the field {Ei,Hi} are determined.

The foregoing suggests the following important conclusion: the fields generated in the reflection zone (the domain A ) and transmission zone (the domain B ) of a periodic structure are uniquely determined by their longitudinal (directed along z -axis) components and can be represented in the following form (see also formulas (14) and (23)). For the incident wave we have

 {Ezi(g,t)Hzi(g,t)}=∑n,m=−∞∞{vnm(z,E)(z,t)vnm(z,H)(z,t)}μnm(x,y),g∈A¯,t≥0, (30)

for the reflected wave Us(g,t) (which coincides with the total field U(g,t) if Ui(g,t)0 ) we have

 {Ezs(g,t)orEz(g,t)Hzs(g,t)orHz(g,t)}=∑n,m=−∞∞{unm(z,E)+(z,t)unm(z,H)+(z,t)}μnm(x,y),g∈A¯,t≥0 (31)

and for the transmitted wave (coinciding in the domain B with the total field U(g,t) ) we can write

 {Ez(g,t)Hz(g,t)}=∑n,m=−∞∞{unm(z,E)−(z,t)unm(z,H)−(z,t)}μnm(x,y),g∈B¯,t≥0. (32)

In applied problems, the most widespread are situations where a periodic structure is excited by one of the partial components of TE -wave (with Ezi(g,t)=0 ) or TM -wave (with Hzi(g,t)=0 ) [7]. Consider, for example, a partial wave of order pq . Then we have

Ui(g,t)=Upq(H)i(g,t):Hzi(g,t)=vpq(z,H)(z,t)μpq(x,y)

or

Ui(g,t)=Upq(E)i(g,t):Ezi(g,t)=vpq(z,E)(z,t)μpq(x,y) .

The excitation of this kind is implemented in our models in the following way. The time function vpq(z,H)(L,t) or vpq(z,E)(L,t) is defined on the boundary L+ . This function determines the width of the pulse Ui(g,t) , namely, the frequency range [K1,K2] such that for all frequencies k from this range ( k=2π/λ , λ is the wavelength in free space) the value

γ=|v˜pq(z,HorE)(L,k)|maxk[K1;K2]|v˜pq(z,HorE)(L,k)|

where v˜pq(z,HorE)(L,k) is the spectral amplitude of the pulse vpq(z,HorE)(L,t) , exceeds some given value γ=γ0 . All spectral characteristics f˜(k) are obtainable from the temporal characteristics f(t) by applying the Laplace transform

 f˜(k)=∫0∞f(t)eiktdt↔f(t)=12π∫iα−∞iα+∞f˜(k)e−iktdk,0≤α≤Imk. (33)

For numerical implementation of the boundary conditions (19) and (20) and for calculating space-time amplitudes of the transverse components of the wave Ui(g,t) in the cross-section z=L of the Floquet channel (formulas (27) and (29)), the function (vpq(z,HorE))(L,t) are to be determined. To do this, we apply the following relation [7,14]:

 v→pq(HorE)(L,t)=∫0tJ0[λpq(t−τ)](v→pq(HorE))′(L,t)dτ,t≥0. (34)

which is valid for all the amplitudes of the pulsed wave Ui(g,t) outgoing towards z= and does not violate the causality principle.

## 6. Transformation Operator Method

### 6.1. Evolutionary basis of a signal and transformation operators

Let us place an arbitrary periodic structure of finite thickness between two homogeneous dielectric half-spaces z1=zL>0 (with ε=ε1 ) and z2=zL>0 (with ε=ε2 ). Let also a local coordinate system gj={xj,yj,zj} be associated with each of these half-spaces (Figure 2). Assume that the distant sources located in the domain A of the upper half-space generate a primary wave U1i(g,t)={E1i(g,t),H1i(g,t)} being incident on the artificial boundary L+ (on the plane z1=0 ) as viewed from z1= .

Denote by Ujs(g,t)={Ejs(g,t),Hjs(g,t)} the fields resulting from scattering of the primary wave U1i(g,t) in the domains A (where the total field is U(g,t)={E(g,t),H(g,t)}=U1s(g,t)+U1i(g,t) ) and B (where U(g,t)=U2s(g,t) ). In Section 5, we have shown that the fields under consideration are uniquely determined by their longitudinal components, which can be given, for example, as:

 {Ezi(g,t)Hzi(g,t)}=∑n,m=−∞∞{vnm(1,E)(z1,t)vnm(1,H)(z1,t)}μnm(x,y),z1≥0,t≥0 (35)
 {Ezs(g,t)Hzs(g,t)}=∑n,m=−∞∞{unm(j,E)(zj,t)unm(j,H)(zj,t)}μnm(x,y),zj≥0,t≥0,j=1,2 (36)

(see also formulas (30)-(32)). Here, as before, {μnm(x,y)}n,m= is the complete (in L2(Rz) ) orthonormal system of transverse eigenfunctions of the Floquet channel R (see Section 4), while the space-time amplitudes unm(j,E)(zj,t) and unm(j,H)(zj,t) are determined by the solutions of the following problems (see also problem (15)) for the one-dimensional Klein-Gordon equations:

 {[−εj∂2∂t2+∂2∂zj2−λnm2]unm(j,EorH)(zj,t)=0,t>0unm(j,EorH)(zj,0)=0,∂∂tunm(j,EorH)(zj,t)|t=0=0,j=1,2,n,m=0,±1,±2,.... (37)

Compose from the functions vnm(1,E)(z1,t) , vnm(1,H)(z1,t) , unm(j,E)(zj,t) , unm(j,H)(zj,t) and the eigenvalues λnm ( n,m=0,±1,±2,... ) the sets v(1)(z1,t)={vp(1)(z1,t)}p= , u(j)(zj,t)={up(j)(zj,t)}p= , and {λp}p= such that their members are defined according to the rules depicted in Figure 3. The sets v(1)(z1,t) and u(j)(zj,t) are said to be evolutionary bases of signals U1i(g,t) and Ujs(g,t) . They describe completely and unambiguously transformation of the corresponding nonsine waves in the regular Floquet channels A and B filled with dielectric.

### Figure 2.

A two-dimensionally periodic grating between two dielectric half-spaces as element of a multi-layered structure.

Let us introduce by the relations

 up(j)′(0,t)≡∂∂zjup(j)(zj,t)|zj=0=∫0t∑q=−∞∞[SpqAA(t−τ)δj1+SpqBA(t−τ)δj2]vq(1)(0,τ)dτ,t≥0,p=0,±1,±2,...,j=1,2 (38)
 u(j)′(0,t)={up(j)′(0,t)}p=[SAAδj1+SBAδj2][v(1)(0,τ)],t≥0,j=1,2 (39)

the boundary (on the boundaries zj=0 ) transformation operators SAA and SBA of the evolutionary basis v(1)(z1,t) of the wave U1i(g,t) incoming from the domain A . Here δmn stands for the Kronecker delta, the operators’ elements SnmXY specify the space-time energy transformation from the domain Y into the domain X and from the mode of order m into the mode of order n .

It is evident that the operators SAA and SBA working in the space of evolutionary bases are intrinsic characteristics of the periodic structure placed between two dielectric half-spaces. They totalize an impact of the structure on elementary excitations composing any incident signal U1i(g,t) . Thus for vq(1)(0,t)=δqrδ(tη) , where r is an integer and η>0 , we have up(1)(0,t)=SprAA(tη) and up(2)(0,t)=SprBA(tη) . We use this example with an abstract nonphysical signal by methodological reasons in order to associate the transformation operators’ components SprAA(tτ) and SprBA(tτ) with an ‘elementary excitation’.

### Figure 3.

Construction of sets of the values vp(1) , up(j) , and λp ( p=0,±1,±2,... ) from sets of the values vnm(1,E) , unm(j,E) , vnm(1,H) , unm(j,H) , and λnm ( m,n=0,±1,±2,... ): (a) p=0,1,2,... ; (b) p=1,2,3,... .

The operators SAA and SBA determine all the features of transient states on the upper and bottom boundaries of the layer enclosing the periodic structure. Secondary waves outgoing from these boundaries propagate freely in the regular Floquet channels A and B therewith undergoing deformations (see, for example, [6]). The space-time amplitudes up(j)(zj,t) of the partial components of these waves (the elements of the evolutionary bases of the signals Ujs(g,t) ) vary differently for different values of p and j . These variations on any finite sections of the Floquet channels A и B are described by the diagonal transporting operators Z0z1A and Z0z2B acting according the rule:

 u(j)(zj,t)={up(j)(zj,t)}=[Z0→z1Aδj1+Z0→z2Bδj2][u(j)′(0,τ)],j=1,2. (40)

The structure of the operators given by (40) can be detailed by the formula

 up(j)(zj,t)=−1εj∫0J0[λp(t−τ)2εj−zj2]χ(t−τεj−zj)up(j)′(0,τ)dτ,t≥0,zj≥0,p=0,±1,±2,...,j=1,2, (41)

which reflects general properties of solutions of homogeneous problems (37), i.e. the solutions that satisfy zero initial conditions and are free from the components propagating in the direction of decreasing zj . The derivation technique for (41) is discussed at length in [6,13,14].

### 6.2. Equations of the operator method in the problems for multilayer periodic structures

The operators SAA and SBA completely define properties of the periodic structure excited from the channel A . By analogy with (38) we can determine transformation operators SBB and SAB for evolutionary basis v(2)(z2,t)={vp(2)(z2,t)}p= of the wave U2i(g,t)={E2i(g,t),H2i(g,t)} incident onto the boundary z2=0 from the channel B :

 up(j)′(0,t)=∫0t∑m=−∞∞[SpqAB(t−τ)δj1+SpqBB(t−τ)δj2]vq(2)(0,τ)dτ,t≥0,p=0,±1,±2,...,j=1,2. (42)

Let us construct the algorithm for calculating scattering characteristics of a multilayer structure consisting of two-dimensionally periodic gratings, for which the operators SAA , SBA , SpqAB , and SpqBB are known. Consider a double-layer structure, whose geometry is given in Figure 4. Two semi-transparent periodic gratings I and II are separated by a dielectric layer of finite thickness M (here ε=ε2(I)=ε1(II) ) and placed between the upper and the bottom dielectric half-spaces with the permittivity ε1(I) and ε2(II) , respectively. Let also a pulsed wave like (35) be incident onto the boundary z1(I)=0 from the Floquet channel A .

Retaining previously accepted notation (the evident changes are conditioned by the presence of two different gratings I and II), represent the solution of the corresponding initial boundary value problem in the regular domains A , B , and C in a symbolic form

U(A)=p=[vp(1)(z1(I),t)+up(1)(z1(I),t)]μp(x,y),U(B)=p=[up(2)(z2(I),t)+up(1)(z1(II),t)]μp(x,y),U(C)=p=up(2)(z2(II),t)μp(x,y).

The first terms in the square brackets correspond to the waves propagating towards the domain C , while the second ones correspond to the waves propagating towards the domain A (Figure 4). The set {μp(x,y)}p= is formed from the functions μnm(x,y) , (n,m=0,±1,±2,... ), while the set {λp}p= is composed from the values λnm , (n,m=0,±1,±2,... ) (Figure 3).

### Figure 4.

Schematic drawing of a double-layered structure.

By denoting

u(j)(I)zj(I)u(j)(zj(I),t)|zj(I)=0,u(j)(I)={up(j)(zj(I),t)}|zj(I)=0,

according to formulas (38)-(42), we construct the following system of operator equations:

 {u(1)′(I)=SAA(I)[v(1)(I)]+SAB(I)Zz1(II)=0→MB[u(1)′(II)]u(2)′(I)=SBA(I)[v(1)(I)]+SBB(I)Zz1(II)=0→MB[u(1)′(II)]u(1)′(II)=SBB(II)Zz2(I)=0→MB[u(2)′(I)]u(2)′(II)=SCB(II)Zz2(I)=0→MB[u(2)′(I)]. (43)

Equations (43) clearly represent step-by-step response of the complex structure on the excitation by the signal U1i(g,t) with the evolutionary basis v(1)(z1(I),t)={vp(1)(z1(I),t)}p= (or simply v(1)(I) ). Тhus, for example, the first equation can be interpreted as follows. A signal u(1)(I) (the secondary field in A ) is a sum of two signals, where the first signal is a result of the reflection of the incident signal v(1)(I) by the grating I , while another one is determined by the signal u(1)(II) being deformed during propagation in the channel B and interaction with the grating I .

By method of elimination the system (43) is reduced to the operator equation of the second kind

 u(2)′(I)=SBA(I)[v(1)(I)]+SBB(I)Zz1(II)=0→MBSBB(II)Zz2(I)=0→MB[u(2)′(I)] (44)

and some formulas for calculating the electromagnetic field components in all regions of the two-layered structure. The observation time t for the unknown function u(2)(I) from the left-hand side of equation (44) strictly greater of any moment of time τ for the function u(2)(I) in the right-hand side of the equation (owing to finiteness of wave velocity). Therefore equation (44) can be inverted explicitly in the framework of standard algorithm of step-by-step progression through time layers. Upon realization of this scheme and calculation of the boundary operators by (38), (42), the two-layered structure can be used as ‘elementary’ unit of more complex structures.

Turning back to (38)-(42), we see that the operators entering these equations act differently that their analogues in the frequency domain, where the boundary operators relate a pair ‘field field’. Reasoning from the structure of the transport operators Z0z1A and Z0z2B (formulas (40) and (41)), we relate a pair ‘field directional derivative with respect to the propagation direction’ to increase numerical efficiency of the corresponding computational algorithms.

## 7. Some Important Properties of Steady-State Fields in the Rectangular Floquet Channel

### 7.1. Excitation by a TM -wave

Let a grating (Figure 1) be excited form the domain A by a pulsed TM -wave Ui(g,t)=Upq(E)i(g,t):Ezi(g,t)=vpq(z,E)(z,t)μpq(x,y) and the region QL is free from the sources j(g,t) , φE(g) , and φH(g) . The field generated in the domains A and B is determined completely by their longitudinal components. They can be represented in the form of (31), (32). Define steady-state fields {E˜(g,k),H˜(g,k)} (see formula (33) with Imk=0 ) corresponding to the pulsed fields {Ei,Hi} , {Es,Hs} in A and the pulsed field {E,H} in B , by their z -components:

 {E˜zi(g,k)H˜zi(g,k)}={v˜pq(z,E)(k)0}e−iΓpq(z−L)μpq(x,y),g∈A¯ (45)
 {E˜zs(g,k)H˜zs(g,k)}=∑n,m=−∞∞{u˜nm(z,E)+(k)u˜nm(z,H)+(k)}eiΓnm(z−L)μnm(x,y),g∈A¯ (46)
 {E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{u˜nm(z,E)−(k)u˜nm(z,H)−(k)}e−iΓnm(z+L)μnm(x,y),g∈B¯ (47)

where the following notation is used: v˜pq(z,E)(k)vpq(z,E)(L,t) , u˜nm(z,EorH)±(k)unm(z,EorH)±(±L,t) , Γnm=(k2λnm2)1/2 , ReΓnmRek0 , ImΓnm0 [7].

The amplitudes u˜nm(z,EorH)±(k) form the system of the so-called scattering coefficients of the grating, namely, the reflection coefficients

 Rpq(E)nm(H)=u˜nm(z,H)+(k)v˜pq(z,E)(k), Rpq(E)nm(E)=u˜nm(z,E)+(k)v˜pq(z,E)(k), n,m=0,±1,±2,... (48)

specifying efficiency of transformation of pq -th harmonic of a monochromatic TM -wave into of order nm -th harmonics of the scattered field {E˜s,H˜s} in the reflection zone, and the transmission coefficients

 Tpq(E)nm(H)=u˜nm(z,H)−(k)v˜pq(z,E)(k), Tpq(E)nm(E)=u˜nm(z,E)−(k)v˜pq(z,E)(k), n,m=0,±1,±2,... (49)

determining the efficiency of excitation of the transmitted harmonics in the domain B .

These coefficients are related by the energy balance equations

 ∑n,m=−∞∞1λnm2[(|Rpq(E)nm(E)|2+|Tpq(E)nm(E)|2)±η02(|Rpq(E)nm(H)|2+|Tpq(E)nm(H)|2)]{ReΓnmImΓnm}=1λpq2{ReΓpq+2ImΓpqImRpq(E)pq(E)ImΓpq−2ReΓpqImRpq(E)pq(E)}∓1ε0{W1W2},p,q=0,±1,±2,…, (50)
 W1=ε0η0k∫QLσ(g,k)|E→˜(g,k)|2dg,W2=∫QL[μ0μ(g,k)|H→˜(g,k)|2−ε0ε(g,k)|E→˜(g,k)|2]dg (51)

They follow from the complex power theorem (Poynting theorem) in the integral form [11]

 ∮SL([E→˜×H→˜∗]⋅ds→)=∫QLdiv[E→˜×H→˜∗]dg=ikη0∫QLμ|H→˜|2dg−ikη0∫QLε|E→˜|2dg−∫QLσ|E→˜|2dg (52)

where ε(g,k)1=χ˜ε(g,k)χε(g,t) , μ(g,k)1=χ˜μ(g,k)χμ(g,t) , σ(g,k)=χ˜σ(g,k)χσ(g,t) , ds is the vector element of the surface SL bounding the domain QL . Equations (50)-(52) have been derived starting from the following boundary value problem for a diffraction grating illuminated by a plane TM -wave U˜pq(E)i(g,k):E˜zi(g,k)=exp[iΓpq(zL)]μpq(x,y) :

 {η0rotH→˜(g,k)=−ikε¯(g,k)E→˜(g,k),rotE→˜(g,k)=ikη0μ(g,k)H→˜(g,k),g∈QLD[E→˜(H→˜)](lx,y)=e2πiΦxD[E→˜(H→˜)](0,y),0≤y≤ly,|z|
 {E˜z(g,k)H˜z(g,k)}={10}e−iΓpq(z−L)μpq(x,y)+∑n,m=−∞∞{Rpq(E)nm(E)(k)Rpq(E)nm(H)(k)}eiΓnm(z−L)μnm(x,y),g∈A¯,{E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Tpq(E)nm(E)(k)Tpq(E)nm(H)(k)}e−iΓnm(z+L)μnm(x,y),g∈B¯. (54)

When deriving (50), (51) we have also used the equations relating z -components of the eigenmode of the Floquet channel

 U→˜(g,k):E˜z(g,k)=Ae±iΓzμ(x,y) and H˜z(g,k)=Be±iΓzμ(x,y) (55)

(subscripts nm are omitted) with its longitudinal components:

 E˜x=−βkη0λ2H˜z∓αΓλ2E˜z,E˜y=αkη0λ2H˜z∓βΓλ2E˜z,H˜x=∓αΓλ2H˜z+βkη0λ2E˜z,H˜y=∓βΓλ2H˜z−αkη0λ2E˜z. (56)

Here, ε¯(g,k)=ε(g,k)+iη0σ(g,k)/k , μ(x,y)=(lxly)1/2exp(iαx)exp(iβy) , Γ=k2λ2 , λ2=α2+β2 .

According to the Lorentz lemma [11], the fields {E˜(1),H˜(1)} and {E˜(2),H˜(2)} resulting from the interaction of a grating with two plane TM -waves

U˜pq(E)i(1)(g,k):E˜zi(1)(g,k)=exp[iΓpq(Φx,Φy)(zL)]μpq(x,y,Φx,Φy) and

U˜r,s(E)i(2)(g,k):E˜zi(2)(g,k)=exp[iΓr,s(Φx,Φy)(zL)]μr,s(x,y,Φx,Φy) ,

satisfy the following equation

 ∮SL(([E→˜(1)×H→˜(2)]−[E→˜(2)×H→˜(1)])⋅ds→)=0. (57)

From (57), using (54) and (56), we obtain

 Rpq(E)rs(E)(Φx,Φy)λp,q2(Φx,Φy)Γpq(Φx,Φy)=R−r,−s(E)−p,−q(E)(−Φx,−Φy)λ−r,−s2(−Φx,−Φy)Γ−r,−s(−Φx,−Φy), p,q,r,s=0,±1,±2,... (58)

– the reciprocity relations, which are of considerable importance in the physical analysis of wave scattering by periodic structures as well as when testing numerical algorithms for boundary problems (53), (54).

Assume now that the first wave U˜pq(E)i(1)(g,k): :E˜zi(1)(g,k)=exp[iΓpq(Φx,Φy)(zL)]μpq(x,y,Φx,Φy)=U˜pq(E)i(1)(g,k,A) be incident on the grating from the domain A , as in the case considered above, while another wave U˜r,s(E)i(2)(g,k):E˜zi(2)(g,k,B)=exp[iΓr,s(Φx,Φy)(z+L)]μr,s(x,y,Φx,Φy) is incident from B . Both of these waves satisfy equation (57), whence we have

 Tpq(E)rs(E)(Φx,Φy,A)λp,q2(Φx,Φy)Γpq(Φx,Φy)=T−r,−s(E)−p,−q(E)(−Φx,−Φy,B)λ−r,−s2(−Φx,−Φy)Γ−r,−s(−Φx,−Φy),p,q,r,s=0,±1,±2,... (59)

### 7.2. Excitation by a TE -wave

Let a grating be excited form the domain A by a pulsed TE -wave Ui(g,t)=Upq(H)i(g,t):Hzi(g,t)=vpq(z,H)(z,t)μpq(x,y) and the region QL is free from the sources j(g,t) , φE(g) , and φH(g) . The field generated in the domains A and B is determined completely by their longitudinal components. They can be represented in the form of (31), (32). Define steady-state fields {E˜(g,k),H˜(g,k)} corresponding to the pulsed fields {Ei,Hi} , {Es,Hs} in A and the pulsed field {E,H} in B , by their z -components as was done for the TM -case (see equations (45)-(47)). Introduce the scattering coefficients Rpq(H)nm(E) , Rpq(H)nm(H) , Tpq(H)nm(E) , and Tpq(H)nm(H) by the relations like (48). These coefficients can be determined from the problems

 {η0rotH→˜(g,k)=−ikε¯(g,k)E→˜(g,k),rotE→˜(g,k)=ikη0μ(g,k)H→˜(g,k),g∈QLD[E→˜(H→˜)](lx,y)=e2πiΦxD[E→˜(H→˜)](0,y),0≤y≤ly,|z|
 {E˜z(g,k)H˜z(g,k)}={01}e−iΓpq(z−L)μpq(x,y)+∑n,m=−∞∞{Rpq(H)nm(E)(k)Rpq(H)nm(H)(k)}eiΓnm(z−L)μnm(x,y),g∈A¯,{E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Tpq(H)nm(E)(k)Tpq(H)nm(H)(k)}e−iΓnm(z+L)μnm(x,y),g∈B¯ (61)

and satisfy the following relations, which are corollaries from the Poynting theorem and the Lorentz lemma:

 ∑n,m=−∞∞1λnm2[(|Rpq(H)nm(H)|2+|Tpq(H)nm(H)|2)±1η02(|Rpq(H)nm(E)|2+|Tpq(H)nm(E)|2)]{ReΓnmImΓnm}==1λpq2{ReΓpq+2ImΓpqImRpq(H)pq(H)ImΓpq−2ReΓpqImRpq(H)pq(H)}−1μ0{W1W2},p,q=0,±1,±2,… (62)

and

 Rpq(H)rs(H)(Φx,Φy)λp,q2(Φx,Φy)Γpq(Φx,Φy)=R−r,−s(H)−p,−q(H)(−Φx,−Φy)λ−r,−s2(−Φx,−Φy)Γ−r,−s(−Φx,−Φy), p,q,r,s=0,±1,±2,... (63)
 Tpq(H)rs(H)(Φx,Φy,A)λp,q2(Φx,Φy)Γpq(Φx,Φy)=T−r,−s(H)−p,−q(H)(−Φx,−Φy,B)λ−r,−s2(−Φx,−Φy)Γ−r,−s(−Φx,−Φy),p,q,r,s=0,±1,±2,... (64)

### 7.3. General properties of the grating’s secondary field

Let now k be a real positive frequency parameter, and let an arbitrary semi-transparent grating (Figure 1) be excited from the domain A by a homogeneous TM - or TE -wave

 U→˜pq(EorH)i(g,k):{E˜zi(g,k) or H˜zi(g,k)}=e−iΓpq(z−L)μpq(x,y),p,q:ImΓpq=0. (65)

The terms of infinite series in (54) and (61) are z -components of nm -th harmonics of the scattered field for the domains A and B . The complex amplitudes Rpq(EorH)nm(EorH) and Tpq(EorH)nm(EorH) are the functions of k , Φx , Φy , as well as of the geometry and material parameters of the grating. Every harmonic for which ImΓnm=0 and ReΓnm>0 is a homogeneous plane wave propagating away from the grating along the vector knm : kx=αn , ky=βm , kz=Γnm (in A ; Figure 5) or kz=Γnm (in B ). The frequencies k such that Γnm(k)=0 ( k=knm±=±|λnm| ) are known as threshold frequency or sliding points [1-6]. At those points, a spatial harmonic of order nm with ImΓnm>0 are transformed into a propagating homogeneous pane wave.

It is obvious that the propagation directions knm of homogeneous harmonics of the secondary field depends on their order nm , on the values of k and on the directing vector of the incident wave kpqi : kxi=αp , kyi=βq , kzi=Γpq . According to (50) and (62), we can write the following formulas for the values, which determine the ‘energy content’ of harmonics, or in other words, the relative part of the energy directed by the structure into the relevant spatial radiation channel:

 (WR)pqnm=(|Rpq(E)nm(E)|2+η02|Rpq(E)nm(H)|2)ReΓnmλnm2λpq2Γpq=(WR)pq(E)nm(E)+(WR)pq(E)nm(H),(WT)pqnp=(|Tpq(E)nm(E)|2+η02|Tpq(E)nm(H)|2)ReΓnmλnm2λpq2Γpq=(WT)pq(E)np(E)+(WT)pq(E)np(H) (66)

(for TM -case) and

### Figure 5.

On determination of propagation directions for spatial harmonics of the field formed by a two-dimensionally periodic structure.

 (WR)pqnm=(|Rpq(H)nm(H)|2+1η02|Rpq(H)nm(E)|2)ReΓnmλnm2λpq2Γpq=(WR)pq(H)nm(H)+(WR)pq(H)nm(E),(WT)pqnp=(|Tpq(H)nm(H)|2+1η02|Tpq(H)nm(E)|2)ReΓnmλnm2λpq2Γpq=(WT)pq(H)nm(H)+(WT)pq(H)nm(E) (67)

(for TE -case). The channel corresponding to the nm -th harmonic will be named ‘open’ if ImΓnm=0 . The regime with a single open channel ( nm=pq ) will be called the single-mode regime.

Since |kpqi|=|knm|=k , the nm -th harmonic of the secondary field in the reflection zone propagates in opposition to the incident wave only if αn=αp and βm=βq or, in other notation, if

 n=−2Φx−p and m=−2Φy−q (68)

Generation of the nonspecularly reflected mode of this kind is termed the auto-collimation.

The amplitudes Rpq(EorH)nm(EorH) or Tpq(EorH)nm(EorH) are not all of significance for the physical analysis. In the far-field zone, the secondary field is formed only by the propagating harmonics of the orders nm such that ReΓnm0 . However, the radiation field in the immediate proximity of the grating requires a consideration of the contribution of damped harmonics (nm:ImΓnm>0 ). Moreover, in some situations (resonance mode) this contribution is the dominating one [6].

### 7.4. The simplest corollaries of the reciprocity relations and the energy conservation law

Let us formulate several corollaries of the relations (50), (58), (59), and (62)-(64) basing on the results presented in [3] and [7] for one-dimensionally periodic gratings and assuming that ε(g,k)0 , μ(g,k)0 , and σ(g,k)0 .

• The upper lines in (50) and (62) represent the energy conservation law for propagating waves. If ImΓpq=0 , the energy of the scattered field is clearly related with the energy of the incident wave. The energy of the wave U˜pq(EorH)i(g,k) is partially absorbed by the grating (only if W10 ), and the remaining part is distributed between spatial TM - and TE -harmonics propagating in the domains A and B (the wave is reradiating into the directions z=± ). If a plane inhomogeneous wave be incident on a grating ( ImΓpq>0 ), the total energy is defined by the imaginary part of reflection coefficient Rpq(EorH)pq(EorH) , which in this case is nonnegative.

• The relations in the bottom lines in (50), (62) limit the values of n,m=|Rpq(E)nm(E)|2λnm2ImΓnm , n,m=|Tpq(E)nm(E)|2λnm2ImΓnm , etc. and determine thereby the class of infinite sequences

 l¯2={a={anm}nm=−∞∞:∑nm=−∞∞|anm|2n2+m2∞} (69)

or energetic space, to which amplitudes of the scattered harmonics Rpq(E)nm(E) , Tpq(E)nm(E) , etc. belong.

• It follows from (58), (59), (63), and (64) that for all semi-transparent and reflecting gratings we can write

 (WR)00(EorH)00(EorH)(Φx,Φy)=(WR)00(EorH)00(EorH)(−Φx,−Φy),(WT)00(EorH)00(EorH)(Φx,Φy,A)=(WT)00(EorH)00(EorH)(−Φx,−Φy,B). (70)

The first equation in (70) proves that the efficiency of transformation of the TM - or TE -wave into the specular reflected wave of the same polarization remains unchanged if the grating is rotated in the plane x0y about z -axis through 180° . The efficiency of transformation into the principal transmitted wave of the same polarizations does not also vary with the grating rotation about the axis lying in the plane x0y and being normal to the vector k00 (Figure 5).

• When r=s=p=q=0 we derive from (58), (59), (63), and (64) that

 R00(EorH)00(EorH)(Φx,Φy)=R00(EorH)00(EorH)(−Φx,−Φy),T00(EorH)00(EorH)(Φx,Φy,A)=T00(EorH)00(EorH)(−Φx,−Φy,B). (71)

That means that even if a semi-transparent or reflecting grating is non symmetric with respect to the any planes, the reflection and transmission coefficients entering (71) do not depend on the proper changes in the angles of incidence of the primary wave.

• Relations (50), (58) allow the following regularities to be formulated for ideal (σ(g,k)0 ) asymmetrical reflecting gratings. Let the parameters k , Φx , and Φy be such that ReΓ00(Φx,Φy)>0 and ReΓnm(Φx,Φy)=0 for n,m0 . If the incident wave is an inhomogeneous plane wave U˜±p,±q(E)i(g,k,±Φx,±Φy) , then

 (|R±p,±q(E)00(E)(±Φx,±Φy)|2+η02|R±p,±q(E)00(H)(±Φx,±Φy)|2)ReΓ00(±Φx,±Φy)λ002(±Φx,±Φy)==2ImR±p,±q(E)±p,±q(E)(±Φx,±Φy)ImΓ±p,±q(±Φx,±Φy)λ±p,±q2(±Φx,±Φy). (72)

Since Rpq(E)pq(E)(Φx,Φy)=Rp,q(E)p,q(E)(Φx,Φy) , we derive from (72)

 |Rp,q(E)00(E)(Φx,Φy)|2+η02|Rp,q(E)00(H)(Φx,Φy)|2==|R−p,−q(E)00(E)(−Φx,−Φy)|2+η02|R−p,−q(E)00(H)(−Φx,−Φy)|2. (73)

It is easy to realize a physical meaning of the equation (73) and of similar relation for TE -case, which may be of interest for diffraction electronics. If a grating is excited by a damped harmonic, the efficiency of transformation into the unique propagating harmonic of spatial spectrum is unaffected by the structure rotation in the plane x0y about z -axis through 180° . The above-stated corollaries have considerable utility in testing numerical results and making easier their physical interpretation. The use of these corollaries may considerably reduce amount of calculations.

## 8. Elements of Spectral Theory for Two-Dimensionally Periodic Gratings

The spectral theory of gratings studies singularities of analytical continuation of solutions of boundary value problems formulated in the frequency domain (see, for example, problems (53), (54) and (60), (61)) into the domain of complex-valued (nonphysical) values of real parameters (like frequency, propagation constants, etc.) and the role of these singularities in resonant and anomalous modes in monochromatic and pulsed wave scattering. The fundamental results of this theory for one-dimensionally periodic gratings are presented in [4,6,7]. We present some elements of the spectral theory for two-dimensionally periodic structures, which follow immediately form the results obtained in the previous sections. The frequency k acts as a spectral parameter; a two-dimensionally periodic grating is considered as an open periodic resonator.

### 8.1. Canonical Green function

Let a solution G˜0(g,p,k) of the scalar problem

 {[Δg+k2][G˜0(g,p,k)]=δ(g−p),g={xg,yg,zg}∈R,p={xp,yp,zp}∈QLD[G˜0](lx,yg)=e2πiΦxD[G˜0](0,yg),0≤yg≤ly,|zg|≤LD[G˜0](xg,ly)=e2πiΦyD[G˜0](xg,0),0≤xg≤lx,|zg|≤LG˜0(g,p,k)=∑n,m=−∞∞{Anm(p,k)Bnm(p,k)}e±iΓnm(zg∓L)μnm(xg,yg),g∈{A¯B¯} (74)

is named the canonical Green function for 2-D periodic gratings. In the case of the elementary periodic structure with the absence of any material scatterers, the problems of this kind but with arbitrary right-hand parts of the Helmholtz equation are formulated for the monochromatic waves generated by quasi-periodic current sources located in the region |z|<L .

Let us construct G˜0(g,p,k) as a superposition of free-space Green functions:

 G˜0(g,p,k)=−14π∑n,m=−∞∞exp[ik|g−pnm|]|g−pnm|e2πinΦxe2πimΦy,pnm={xp+nlx,yp+mly,zp}. (75)

By using in (75) the Poisson summation formula [15] and the tabulated integrals [16]

exp(ipx2+a2)x2+a2eibxdx=πiH0(1)(a|p2b2|) and H0(1)(px2+a2)eibxdx=2exp(iap2b2)p2b2 ,

where H0(1)(x) is the Hankel function of the first kind, we obtain

 G˜0(g,p,k)=−i2lxly∑n,m=−∞∞ei[αn(xg−xp)+βm(yg−yp)]exp[i|zg−zp|Γnm]Γnm (76)

The surface K of analytic continuation of the canonical Green function (76) into the domain of complex-valued k is an infinite-sheeted Riemann surface consisting of the complex planes kC with cuts along the lines (Rek)2(Imk)2λnm2=0 , n,m=0,±1,±2,... , Imk0 (Figure 6). The first (physical) sheet Ck of the surface K is uniquely determined by the radiation conditions for G˜0(g,p,k) in the domains A and B , i.e. by the choice of ReΓnmRek0 and ImΓnm0 on the axis Imk=0 . On this sheet, in the domain 0argk<π , we have ImΓnm>0 , while ReΓnm0 for 0argkπ/2 and ReΓnm0 for π/2argkπ . In the domain 3π/2argk<2π for finite number of functions Γnm(k) (with n and m such that (Rek)2(Imk)2λnm2>0 ), the inequalities ImΓnm<0 and ReΓnm>0 hold; for the rest of these functions we have ImΓnm>0 and ReΓnm0 . In the domain π<argk3π/2 , the situation is similar only the signs of ReΓnm are opposite. On the subsequent sheets (each of them with its own pair {k;Γnm(k)} ), the signs (root branches) of Γnm(k) are opposite to those they have on the first sheet for a finite number of n and m . The cuts (solid lines in Fig. 6) originate from the real algebraic branch points knm±=±|λnm| .

### Figure 6.

Natural domain of variation of the spectral parameter k : the first sheet of the surface K .

In the vicinity of some fixed point KK the function G˜0(g,p,k) can be expanded into a Loran series in terms of the local variable [17]

κ={kK;K{knm±}kK;K{knm±} .

Therefore, this function is meromorphic on the surface K . Calculating the residuals Resk=k¯G˜0(g,p,k) at the simple poles k¯{knm±} , we obtain nontrivial solutions of homogeneous ( U˜i(g,k)0 ) canonical ( ε¯(g,k)1 , μ(g,k)1 , intS¯= ) problems (53), (54) and (60), (61):

 E→˜(g,knm±)={E˜x,E˜y,E˜z}; E˜x,yorz=ax,yorzexp[i(αnx+βmy)] and H→˜(g,knm±)=(iknm±η0)−1rotE→˜(g,knm±) (77)

where ax,yorz are the arbitrary constants. These solutions determine free oscillations in the space stratified by the following conditions:

 D[E→˜(H→˜)](x+lx,y)=e2πiΦxD[E→˜(H→˜)](x,y),D[E→˜(H→˜)](x,y+ly)=e2πiΦyD[E→˜(H→˜)](x,y). (78)

### 8.2. Spectrum qualitative characteristics

Let a set Ωk of the points {k¯j}jK such that for all k{k¯j}j the homogeneous (spectral) problem

 {η0rotH→˜(g,k)=−ikε¯(g,k)E→˜(g,k),rotE→˜(g,k)=ikη0μ(g,k)H→˜(g,k),g∈QLD[E→˜(H→˜)](lx,y)=e2πiΦxD[E→˜(H→˜)](0,y),0≤y≤ly,|z|
 {E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Anm(E)(k)Anm(H)(k)}eiΓnm(z−L)μnm(x,y),g∈A¯{E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Bnm(E)(k)Bnm(H)(k)}e−iΓnm(z+L)μnm(x,y),g∈B¯ (80)

has a nontrivial (not necessarily unique) solution U˜(g,k¯j)={E˜(g,k¯j),H˜(g,k¯j)} be called the point spectrum of the grating. It is obvious that these solutions characterize the so-called free oscillations, whose field pattern, structure of their spatial harmonics and behavior of these harmonics for large |z| and t are determined by the value of k¯j=Rek¯j+iImk¯j and by a position of the point k¯j (the eigen frequency associated with a free oscillation U˜(g,k¯j) ) on the surface K [4,6,7]. By continuing analytically the problems (53), (54) and (60), (61) together with their solutions U˜(g,k)={E˜(g,k),H˜(g,k)} into the domain K of the complex-valued k , we detect poles of the function U˜(g,k) at the points k=k¯j . In the vicinity of these poles, the desired solutions can be represented by the Loran series in terms of the local on K variable κ [17]. The analytical findings of this kind may form the basis for detailed study of physical features of resonant wave scattering by one-dimensionally and two-dimensionally periodic structures [4,6,7,18,19].

Derive now the conditions that constrain existence of nontrivial solutions of the problem (79), (80). These conditions can be considered as uniqueness theorems for the problems (53), (54) and (60), (61) formulated for different domains of the surface K . Notice that the study of the uniqueness allows one to estimate roughly a domain where elements of the set Ωk are localized and simplify substantially the subsequent numerical solution of spectral problems owing to reduction of a search zone of the eigen frequencies. The uniqueness theorems serve also as a basis for application of the ‘meromorphic’ Fredholm theorem [20] when constructing well grounded algorithms for solving diffraction problems as well as when studying qualitative characteristics of gratings’ spectra [4,7].

Assume that grating scattering elements are nondispersive ( ε(g,k)=ε(g) , μ(g,k)=μ(g) , and σ(g,k)=σ(g) ). In this case, the analytical continuation of the spectral problem (79), (80) into the domain of complex-valued k are simplified considerably. From the complex power theorem in the integral form formulated for the nontrivial solutions U˜(g,k¯j) like

 ∮SL([E→˜×H→˜∗]⋅ds→)=∫QLdiv[E→˜×H→˜∗]dg=ikη0∫QLμ|H→˜|2dg−ik∗η0∫QLε|E→˜|2dg−∫QLσ|E→˜|2dg (81)

the following relations result:

 ∑n,m=−∞∞1λnm2{(ReΓnmRek+ImΓnmImk)(ImΓnmRek−ReΓnmImk)}[(|Anm(E)|2+|Bnm(E)|2)±η02(|Anm(H)|2+|Bnm(H)|2)]=1ε0{−Imk(V3+V2)−V1Rek(V3−V2)} (82)

Notation: k=k¯j , E˜=E˜(g,k¯j) , Γnm=Γnm(k¯j) , Anm(E)=Anm(E)(k¯j) , etc., and

V1=ε0η0QLσ|E˜|2dg , V2=QLε0ε|E˜|2dg , V3=QLμ0μ|H˜|2dg .

No free oscillations exist whose amplitudes do not satisfy equations (82). From this general statement, several important consequences follow. Below some of them are formulated for gratings with ε(g)>0 , μ(g)>0 , and σ(g)0 .

• There are no free oscillations whose eigen frequencies k¯j are located on the upper half-plane ( Imk0 ) of the first sheet of the surface K . This can be verified by taking into account the upper relation in (82), the function Γnm(k) on Ck , and the inequalities V10 , V2>0 , V3>0 .

• If σ(g)0 (the grating is non-absorptive), no free oscillations exist whose eigen frequencies k¯j are located on the bottom half-plane ( Imk<0 ) of the sheet Ck between the cuts corresponding to the least absolute values of knm± . In Figure 6, this region of the first sheet of K and the above-mentioned domain are shaded by horizontal lines.

• If σ(g)>0 on some set of zero-measure points gQL , then there are no elements k¯j of grating’s point spectrum Ωk that are located on the real axis of the plane Ck .

Investigation of the entire spectrum of a grating, i.e. a set of the points kK , for which the diffraction problems given by (53), (54) and (60), (61) are not uniquely solvable, is a complicated challenge. Therefore below we do no more than indicate basic stages for obtaining well grounded results. The first stage is associated with regularization of the boundary value problem describing excitation of a metal-dielectric grating by the currents J˜(g,k)J(g,t) located in the domain QL :

 {η0rotH→˜(g,k)=−ikε¯(g,k)E→˜(g,k)+J→˜(g,k),rotE→˜(g,k)=ikη0μ(g,k)H→˜(g,k),g∈QLD[E→˜(H→˜)](lx,y)=e2πiΦxD[E→˜(H→˜)](0,y),0≤y≤ly,|z|
 {E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Anm(E)(k)Anm(H)(k)}eiΓnm(z−L)μnm(x,y),g∈A¯,{E˜z(g,k)H˜z(g,k)}=∑n,m=−∞∞{Bnm(E)(k)Bnm(H)(k)}e−iΓnm(z+L)μnm(x,y),g∈B¯. (84)

By regularization is meant (see, for example, [7]) a reduction of the boundary value electrodynamic problem to the equivalent operator equation of the second kind

 [E+B(G˜0,S,ε¯,μ,k)]X=Y,EX=X (85)

with a compact (in some space W of vector fields) finite-meromorphic (in local on K variables κ ) operator-function B(G˜0,S,ε¯,μ,k) [20,21]. If the problem given by (83), (84) is considered separately for metal gratings ( intS¯ and S are sufficiently smooth surfaces; ε¯(g,k)=μ(g,k)1 ) and dielectric gratings ( intS¯= , ε¯(g,k)=ε(g) and μ(g,k)=μ(g) are sufficiently smooth functions), then its regularization can be performed by applying the potential theory methods [4,7,22].

In the second stage, the following statements should be proved: (i) the resolvent [E+B(k)]1 ( kK ) of the problem in (85) is a finite-meromorphic operator-function; (ii) its poles are located at the points k=k¯j ( j=1,2,3,... ); (iii) the entire spectrum coincides with its point spectrum Ωk ; (iv) Ωk is nothing more than a countable set without finite accumulation points. All these statements are corollaries of the ‘meromorphic’ Fredholm theorem [4,20,21] and the uniqueness theorem proved previously.

By inverting homogeneous operator equation (85), we can construct a numerical solution of the spectral problem given by (79), (80) [4,6], in other words, calculate the complex-valued eigen frequencies k¯j and associated eigen waves U˜(g,k¯j)={E˜(g,k¯j),H˜(g,k¯j)} or free oscillations of an open two-dimensionnaly periodic resonator. Commonly, this operation is reduced to an approximate solution of the characteristic equation like:

 det[C(k)]=0. (86)

Here C(k) is some infinite matrix-function; the compactness of the operator B(k) ensures (i) existence of the determinant det[C(k)] and (ii) the possibility to approximate the solutions k¯ of equation (86) by the solutions k¯(N) of the equation det[C(k,N)] with the matrix C(k,N) reduced to dimension N×N .

Let k¯ be a root of characteristic equation (86) that do not coincide with any pole of the operator-function B(k) . The multiplicity of this root determines the multiplicity of the eigen value k¯ of homogeneous operator equation (85), i.e. the value M=M(1)+M(2)+...+M(Q) [21]. Here, Q is the number of linearly-independent eigen functions U˜(q)(g,k¯) ; q=1,2,...,Q (the number of free oscillations) corresponding to the eigen value (eigenfrequency) k¯ , while M(q)1 is the number of the associated functions U˜(m)(q)(g,k¯) ; m=1,2,...,M(q)1 . The order of pole of the resolvent [E+B(k)]1 (and of the Green function G˜(g,p,k) of the problem in (83), (84)) for k=k¯ is determined by a maximal value of M(q) .

## 9. Conclusion

The analytical results presented in the chapter are of much interest in the development of rigorous theory of two-dimensionally periodic gratings as well as in numerical solution of the associated initial boundary value problems. We derived exact absorbing boundary conditions truncating the unbounded computational space of the initial boundary value problem for two-dimensionally periodic structures to a bounded part of the Floquet channel. Some important features of transient and steady-state fields in rectangular parts of the Floquet channel were discussed. The technique for calculating electrodynamic characteristics of multi-layered structure consisting of two-dimensionally periodic gratings was developed by introducing the transformation operators similar to generalized scattering matrices in the frequency domain. In the last section, the elements of spectral theory for two-dimensionally periodic gratings were discussed.

## References

1 - V. P. Shestopalov, L. N. Litvinenko, S. A. Masalov, V. G. Sologub, 1973 Wave Diffraction by Gratings. Kharkov: Kharkov State Univ. Press.
2 - R. Petit, editor 1980 Electromagnetic Theory of Gratings. New York: Springer.
3 - V. P. Shestopalov, A. A. Kirilenko, S. A. Masalov, Y. K. Sirenko, 1986 Resonance Wave Scattering. Vol. 1. Diffraction Gratings. Kiev: Naukova Dumka; (in Russian).
4 - V. P. Shestopalov, Y. K. Sirenko, 1989 Dynamic Theory of Gratings. Kiev: Naukova Dumka; (in Russian).
5 - M. Neviere, E. Popov, 2003 Light Propagation in Periodic Media: Differential Theory and Design New York: Marcel Dekker.
6 - Y. K. Sirenko, S. Strom, N. P. Yashina, 2007 Modeling and Analysis of Transient Processes in Open Resonant Structures. New Methods and Techniques New York: Springer.
7 - Y. K. Sirenko, S. Strom, editors 2010 Modern Theory of Gratings. Resonant Scattering: Analysis Techniques and Phenomena New York: Springer.
8 - O. A. Ladyzhenskaya, 1985 The Boundary Value Problems of Mathematical Physics. New York: Springer-Verlag.
9 - A. Taflove, S. C. Hagness, 2000 Computational Electrodynamics: the Finite-Difference Time-Domain Method. Boston: Artech House.
10 - K. Sirenko, V. Pazynin, Y. Sirenko, H. Bagci, 2011 An FFT-Accelerated FDTD Scheme with Exact Absorbing Conditions for Characterizing Axially Symmetric Resonant Structures Progress in Electromagnetics Research 111 331 64
11 - E. J. Rothwell, M. J. Cloud, 2001 Electromagnetics. New York: CRC Press.
12 - V. S. Vladimirov, 1971 Equations of Mathematical Physics. New York: Dekker.
13 - K. Y. Sirenko, Y. K. Sirenko, 2005 Exact’Absorbing’ Conditions in the Initial Boundary-Value Problems of the Theory of Open Waveguide Resonators Computational Mathematics and Mathematical Physics 45 490 506
14 - V. F. Kravchenko, Y. K. Sirenko, K. Y. Sirenko, 2011 Electromagnetic Wave Transformation and Radiation by the Open Resonant Structures. Modelling and Analysis of Transient and Steady-State Processes. Moscow: FizMathLit; (in Russian).
15 - E. Titchmarsh, 1948 Introduction to the Theory of Fourier Integrals Oxford: Clarendon Press.
16 - I. S. Gradshteyn, I. M. Ryzhik, 1994 Table of Integrals, Series, and Products. New York: Academic.
17 - A. von Hurwitz, R. von Courant, 1964 Allgemeine Funktionentheorie und Elliptische Funktionen: Geometrische Funktionentheorie. Berlin: Springer-Verlag; (in German).
18 - Y. K. Sirenko, L. G. Velychko, F. Erden, Time-Domain and Frequency-Domain Methods Combined in the Study of Open Resonance Structures of Complex Geometry Progress in Electromagnetics Research 2004 44 57 79
19 - L. G. Velychko, Y. K. Sirenko, O. S. Shafalyuk, 2006 Time-Domain Analysis of Open Resonators. Analytical Grounds Progress in Electromagnetics Research 61 1 26
20 - M. Reed, B. Simon, 1978 Methods of Modern Mathematical Physics IV: Analysis of Operators. New York: Academic Press.
21 - I. Z. Hokhberg, Y. I. Seagul, 1971 Operator Generalization of the Theorem about Logarithmic Residue and the Rouche Theorem Matematicheskii Sbornik 84 607 629 (in Russian)
22 - D. Colton, R. Kress, 1983 Integral Equation Methods in Scattering Theory New York: Wiley-Interscience.