Methods

Discretization of Maxwell’s equations

The structure under our consideration is a PC system consisting of periodic dielectric rods with finite thickness in a two-dimensional slab. Figure ¹ a shows a PC slab sandwiched between two high-dielectric-constant layers parallel to the x -axis direction. Also, due to total reflection, the electric field vanishes at the boundaries. Although the slab is limited in width, the x -axis direction is unlimited. A point defect can be easily exhibited by inserting an air rod instead of one of the dielectric lattice sites (Figure ¹ b) while the nearest neighboring dielectric rods are numbered. We begin our investigation by following Maxwell’s equations for electric field E and magnetic field H :

1ε(r)∇×{∇×E(r)}=ω2c2E(r),

**A PC slab.** A PC slab consisting of high-dielectric rods between low-dielectric background which is sandwiched in two high-dielectric-constant layers parallel to the x -axis direction, ( a ) the perfect system and ( b ) the imperfect system, while the central rod is removed. The neighboring sites are not removed but are indicated by numbers.

∇×1ε(r)∇×H(r)=ω2c2H(r),

where ω is the light frequency, c is the light velocity, and ε _r ( r ) is the dielectric constant. By assuming that the dielectric constant ε _r ( r ) is independent of the z coordinate, we rewrite the first equation for a two-dimensional slab. We just study the TM mode since the mentioned PC structure with dielectric rods is appropriate for TM modes [ ¹⁹ , ²⁰ ]. This formalism is general, and for TE modes the same procedures have to be reviewed. Extending the method implemented in [ ¹² ], we discretize Equation ¹ since for TM modes the main components of fields are ( E _z , H _x , H _y ); therefore, the electric field and the other two components can be solved by the following differential equations:

1ε(r)∂2Ez∂x2+∂2Ez∂y2+ω2c2Ez=0,Hx=−iωμ0∂Ez∂y,Hy=iωμ0∂Ez∂x.

The corresponding equation for E _z in the form of operator equation with a Hermitian differential operator could be defined as [ ¹² ]

LTMf=ωc2f,

where L and the function f for TM modes are defined by

LTM=−1√ε(r)∂2Ez∂x2+∂2Ez∂y21√ε(r),f=√εr(r)Ez.

This equation has to be discretized in the x direction, and the position of the center of the l th site in the x direction is indicated by x _l = l a ₀ and in the y direction by y ^α = α a ₀ , where a ₀ is the lattice constant and l and α are integers. While l is an infinite number, α is finite. By using the following discretization of Equation ⁵ ,

a02∂2∂x2ξ(x)f(x)→ξl+1αfl+1α−2ξlαflα+ξl−1αfl−1α,a02∂2∂y2ξ(y)f(y)→ξlα+1flα+1−2ξlαflα+ξlα−1flα−1,

where flα=√ξr,lαEz,lα and ξlα=1√ξr,lα , then a finite differential equation is obtained:

vlαflα−tl,l+1α,αfl+1α−tl,l−1α,αfl−1α−tl,lα,α+1flα+1−tl,lα,α−1flα−1=ωa0c2flα,

where the coefficients are defined as

vlα=4(ξlα)2,tl,l+1α,α=ξlαξl+1α,tl,l−1α,α=ξl−1αξlα,tl,lα,α+1=ξlα+1ξlα,tl,lα,α−1=ξlαξlα−1.

General form of tight binding model for a PC slab

For mapping of light propagation in a PC to a tight binding model, we introduce photon creation and annihilation operators al† and alα for the site ( l a ₀ , α a ₀ ) with the following bosonic particle commutation relations:

[alα,alα′†]=δl,l′α,α′†,[al†,alα′†]=0,[alα,al′α′]=0.

When the Hermitian operator L̂ acts on the photonic field stateket ∣ f >, we have the following eigenvalue-eigenfunction equation:

L̂∣f>=ωa0c2∣f>,

where L̂ and ∣ f > are defined by

L̂=∑l;α(vlαal†alα)−∑l,m;α,β(tl,mα,βal†alα),

∣f>=∑l;αflαal†∣0>,

where ∣0> is the photonic vacuum state. The result is equivalent to the finite differential equation ( ⁷ ). To distinguish defected and perfect systems, index zero is used for the perfect system. Equation ¹¹ , in terms of perfect system operators and onsite and hopping energies, is given by

L̂=L̂0+∑l;α(vl′α−vlα)al†alα−∑l,m;α,β(tl,m′α,β−tl,mα,β)al†amβ,

where prime vl′α=4(ξl′α)2 and tl,m′α,β=ξlαξm′β . For convenience, two parentheses on the right-hand side of Equation ¹³ are substituted with εlα and δtl,mα,β , respectively:

L̂=L̂0+∑l;αεlαal†alα−∑l,m;α,βδtl,mα,βal†amβ.

The matrix element of Equation ¹⁴ is Ll,mα,β=L0l,mα,β+εlαδl,mα,β−δtl,mα,β . This equation could be written in the following operator form:

L̂=L̂0+ν̂,

where ν̂ is defined by

ν̂=ε̂−δt̂.

Now an efficient propagator, the Green function [ ²¹ ], is defined:

ωa0c2−L̂G(ω)=1̂G0(ω)=ωa0c2−L̂0−1.

In the following part, we determine the Green functions for the described system to calculate the defected slab electric field.

PC slab Green function

According to Equation ¹⁰ , we act the operator L̂0 on stateket ∣ f k 0>, and since the eigenvalue equals ωka0c2 , we proceed a quantum mechanical treatment to obtain the PC slab Green function:

(E−L̂0)∣fk0>=ωa0c2−ωka0c2∣fk0>,

where E=(ωa0c)2

and k characterizes the field mode. Including all modes, we have

(E−L̂0)∑k1ωa0c2−ωka0c2∣fk0><fk0∣=∑k∣fk0><fk0∣,

then the Green function is defined:

G0(k)=∑k1(ωa0c)2−(ωka0c)2∣fk0><fk0∣.

The field stateket ∣ f > or the electric field ket can be determined according to [ ²⁰ ], but since the electric field vanishes at the limited boundaries, the following relation is adequate for the electric field stateket:

∣fkx,ν(xlα,ylα)>=∑l,νCkx,νexp(ikxxlα)×sinπνwylα∣ukx,ν(xlα)>,

where Ckx,ν is the normalization factor, ν is an integer indicating discrete modes along y axes, and w is the slab width. By defining r→ as a one-dimensional (1D) Bravias lattice site position and d→lα as a sub-site position with respect to the 1D Bravias lattice site in such a description as Figure ² which shows mapping of a PC slab to a 1D Bravis lattice site, the electric field ket Equation ¹⁹ is also written as

∣fkx,ν(r→)>=∑l,ν,αCkx,νexp(ikxxlα)×sinπνwd→lα∣ukx,ν(r→−r→l−d→lα)>;

**A 1D Bravias lattice site.** Mapping of a PC slab to a 1D Bravias lattice site is shown schematically. The horizontal axis is infinite, while the vertical axis is finite.

therefore, a general relation for the perfect slab Green function component is

Gl,m0α,β=∑kx,ν1ωa0c2−ωka0c2∣Ckx,ν∣2×<ukx,ν(r→−r→l−d→lα)∣fkx,ν(xlγ,ylγ)>×<fkx,ν(xmθ,ymθ)∣ukx,ν(r→−r→l−d→lβ)>;

however, for the case of the electric ket Equation ¹⁹ , we have a simpler form of this component:

Gl,m0α,β=∑kx,ν1ωa0c2−ωka0c2∣Ckx,ν∣2exp[ikx(xlmαβ)]×sinπνwylαsinπνwymβ.

This Green function, is used to relate that defected system field to the perfect system field by the Lippmann-Schwinger formalism.

Lippmann-Schwinger equation

The electric field component in the z -axis direction, ψ(r→lα) , is determined by multiplying <ukx,ν(r→−r→l−d→lα)∣ from the left to the stateket ∣fkx,ν(xlα,dlα)> :

ψ(r→lα)=<ukx,ν(r→−r→l−d→lα)∣fkx,ν(xlα,dlα)>;

inserting Equation ¹⁹ , we have

ψ(r→lα)=Ckx,νexp(ikxxlα)sinπνwdlα.

Perturbation modifies the electric field, but with first-order approximation of the Lippmann-Schwinger formalism, we relate the defected system electric field to perfect system’s. For this purpose, we write

(G0−1−ν)ψl=ωa0c2−ωka0c2ψl.

Using ψl≅ψl0 and ωk≅ωk0 , we have L0ψl≅(ωka0c)2ψl and

ωa0c2−ωka0c2=ωa0c2−L0=G0−1,

where the right-hand side of Equation ²⁵ is equal to G ⁰ ⁻¹ ψ _l ; then, all ends to the Lippmann-Schwinger equation:

ψ̂l=ψ̂l0+G0νψ̂l0.

For calculating G ⁰ , we should have the proper dispersion relation of the slab that is illustrated in Figure ¹ . Therefore, by Equations ⁷ and ¹⁹ and some mathematics, we obtain the desired relation:

ωkx,ν=ca0vlα−2tl,l±1α,αcoskxa0−2tl,lα,α±1cosπνwa012,

which is based on tight binding method and is applicable only for the described slab.

Application of the method to a PC slab

In order to relate the wave function of a defected system to a perfect system, Equation ²⁶ is rewritten while a perturbation term is included:

ψ̂lα=ψ0̂lα+∑l′,l′′;β,γG0̂l,l′α,β(εl′βδl′,l′′β,γ−δtl′,l′′β,γ)ψ0̂l′′γ.

We apply this model for a quasi-one-dimensional PC slab consisting of periodic dielectric rods in which the dielectric constant ε _r =8.9 and the central defect site is an air cylinder with a circle cross section. The slab is infinite along the x direction, but due to layering high-dielectric-constant material sheets on its borders parallel to the x direction, it has a finite y direction. Although we have considered the first neighboring sites, the difference of the dielectric constant term for neighboring sites, δtl′,l′′β,γ , can be generalized in such a way that other neighboring sites can be included. For numerical calculation of G0̂l,l′α,β , Equation ²⁰ or 22 is used in which due to the action of the photon creation and annihilation operators al† and alα on the center of sites, dlα is the distance between the center of the sites. The quantized parameter ν for ky=πνw is summed over all discrete modes of separate sites on the y direction. Figure ³ a,b shows plots of dispersion relation in terms of k _x presented for two different numerical integers ν =0,1,2,3,4 and ν =1,2,3,4,5, respectively. Mode numbers start upward to higher values. There is no degeneracy for modes in Figure ³ a, but since the third and fifth modes are degenerate, in Figure ³ b four modes are plotted. We calculated the electric field of the defected system at different sub-sites for a slab with five sub-sites. Figure ⁴ shows the ratio of electric field, and Figure ⁵ shows the ratio of electric field intensity with respect to its clean system at defect site l =0 in sub-site α =3 (central site). In Figure ⁴ at normalized frequency 0.819, the peak’s height is maximum 10.2, while it is 103 for the ratio of field intensity which is presented in Figure ⁵ , where these peaks and other peaks in the following figures correspond to the first Brillouin zone (BZ) edges. The ratio of field intensity for site l =±1 and sub-site α =3 (right side or left side of the defect site) is presented in Figure ⁶ . The peak’s height is maximum 10 at normalized frequency 0.909. The following figures are presented for the second case with ν =1,2,3,4,5. Figure ⁷ shows the ratio of field intensity for the defect site; the peak’s height is maximum 4,782 where the ratio of z-component of electric field is −69.2 at normalized frequency 0.909. These are a considerable amount for amplification of incoming pulses at a certain frequency and demonstrate bosonic localization because of the defect and the limited boundaries. The defect disturbs the electric field of other sites which depends on frequency, so a comparison between the ratio of electric field and field intensity of the defect site with the neighboring sites is shown in the following figures where k _x takes a constant value corresponding to the first BZ edge. This comparison in Figures ⁸ and ⁹ are for the neighboring sites that are indicated by numbers 1 and 3 in Figure ¹ b which are laid on the x direction having the same situation for the ratio of electric field and field intensity. According to Figure ⁸ , an important result is achieved, showing equality of electric field sign of the sites with defect sites for lower frequencies but difference with opposite sign for higher frequencies. Figures ¹⁰ and ¹¹ show this comparison for the neighboring sites that are indicated by numbers 2 and 4 in Figure ¹ b and are laid in the y direction. But different from the first neighboring site on the x direction, Figure ¹⁰ shows that the electrical field of these sites have opposite sign with the electric field of the defect site for both the lower and higher frequency range. In Figures ¹² and ¹³ , electrical field and field intensity of the second neighboring site on each corner of the defect site are compared with those of the two mentioned first neighboring sites. As it is evident in Figure ¹² , these sites have different electrical field with opposite sign with the first neighboring site in the x direction (no. 1 or 3 sites) but with equality sign for lower frequencies and opposite sign for higher frequencies with the first neighboring site in the y direction (no. 2 or 4 sites).