As same as all studies of the electromagnetism, analyses to the propagation of light in a photonic crystal also start with four macroscopic Maxwell’s equations. In cgs units, they are
πρ
where E and H are the macroscopic electric and magnetic fields, D and B are the electric displacement and magnetic induction fields, and ρ and J are free charge and current densities, respectively. Here we are concerned with the behavior of an electromagnetic wave in a source-free region where free charge ρ and free current J in Eq. (2.1) are both zero.
2-1 Introduction
In order to solve the wave equations derived from Maxwell’s equations, we need the constitution equations relating D to E and B to H. Since we do not deal with magnetic material, we assume the magnetic permeability µ is very close to unity and we may set
.
As for D and E, quite generally the components Di of the displacement field are related to the electric field components Ei by the following power series [1]:
.)
To simplify the question, we make four assumptions. First we usually assume the field strengths are small enough so that we are in the linear regime. It means χ and all higher order terms can be ignored. Second, we assume the material is macroscopic and isotropic, so that E(r,ω) and D(r,ω) are related by a scalar dielectric constant ε(r,ω). Third, any explicit frequency dependence of the dielectric constant are also been ignored. The last
assumption is that we focus only on low-loss dielectrics, which means ε(r) is treated as pure real. Hence, we have a brief expression relating D and E fields as
).
With four assumptions above, the Maxwell’s equations [Eq. (2.1)] become
0
The field functions E and H generally are both complicated functions of time and space, but thanks to the linearity of Maxwell's equation, it is convenient to look for solutions in form of harmonic fields:
Because there is no free charge and current, the electromagnetic waves considered to be transverse. By substituting Eq. (2.5) into Eq. (2.4) we can obtain the following equations:
)
Solving Eqs. (2.6) and (2.7) is to solve the eigen-value problems, and is a Hermitian operator. The eigenvectors H(r) and
ΘH of the harmonic modes, and the eigenvalues ( )2
c
ω are proportional to the square frequencies
of those modes.
The Maxwell’s equations are the most important kernel of following calculations (both PWE and FDTD) and analyses in the next chapter except only the tight-binding approximation by solid-state physics that we’ll discuss later.
2-2 Plane-wave expansion method
Photonic crystals is a periodically arranged structure (i.e., its dielectric constant is periodic distributed), so we assume that the dielectric constant is real, isotropic, perfectly periodic with the spatial coordinate rv , and does not depend on frequency. Hence we can write its dielectric function as
( )r (r i)
ε v =ε v+ av , i=1,2,3, (2.8) where are the primitive lattice vectors of the photonic crystal. Because of the spatial periodicity, we introduce the primitive reciprocal lattice vectors {b
{ }avi
i ; i=1,2,3} and the reciprocal lattice vector can be defined as {G}:
i⋅ j =2πδij
Because ε is a periodic function of the spatial coordinate r, we can apply Bloch’s theorem to Eqs. (2.6) and (2.7). and are thus characterized by a wave vector k in the first Brillouin zone and a band index n and expressed as
)
two Fourier expansions of the fields can be derived as the following form of the The expansion coefficients in reciprocal lattice space, i.e., and are denoted by the same symbols as the original ones in real space. Substituting Eq. (2.15), (2.16), (2.11) and (2.12) into (2.6) and (2.7), we obtain the following eigenvalue equations for the expansion coefficients { } and { }:
kn( ) electromagnetic field in the 2D photonic lattice can be decomposed into two independent polarization components, i.e., an E polarization (TM mode) for which the electric field is parallel to the rod axis (E
) (r
Ekn Hkn(r)
z only), and an H polarization (TE mode) for which the magnetic field is parallel to the rod axis (Hz only). In two-dimensional photonic crystals, Eq. (2.19) and Eq. (2.20) reduce to
where Eq. (2.21) is the master equation of TM mode. Similarly, the master equation of TE mode can be written as
1 2
For the photonic band calculating, the expansion coefficients { } in Eq. (2.10) is necessary to be calculated by the plane-wave expansion method. The inverse Fourier transform gives
where V is the volume of the unit cell of the photonic crystal. In general, this integral should be numerically evaluated by FFT method. However, if the shapes of the dielectric components in the unit cell are simple enough, we can calculate it analytically.
2-3 Finite-difference time domain method (FDTD) [29]
The finite-difference time-domain method is introduced by Yee in 1966 [30]. During the 1970s and 1980s, several defense agencies working in the areas motivated large-scale solutions of Maxwell’s equations. The entire field of computation electrodynamics is shifting rapidly in high-speed communications and computing. In 1990, engineers in the general electromagnetic community became aware of the modeling capabilities afforded by FDTD and related techniques, and the interest in this area has expanded well beyond defense technology. The main reason to introduce FDTD method to solve photonic crystal is that when the structure is too complex, it is hard to solve Maxwell’s equation in frequency domain.
FDTD provide a straight forward way to solve it in time domain. With this method, we can see the field distribution in photonic crystals. In addition, there are several advantages in FDTD method. First, FDTD is accurate and robust. The sources of error are well known.
Second, being a time domain technology, FDTD treats impulsive behavior and nonlinear behavior naturally. Third, FDTD uses no linear algebra. Being a fully explicit computation, FDTD avoids the difficulties with linear algebra that limit the size of frequency-domain integral-equation.
When the differential forms of Maxwell's equations are examined, it can be seen that the time derivative of the E field is related to the curl of the H field ( ). This can be simplified to state that the rate of the change in the E field (the time derivative) depends on the change in the H field across space (the curl). The results in the basic FDTD equations are that the new value of the E field is related to the old value of the E field (hence the
×H
∇
difference in time) and the difference of old values of the H fields on either side of the E field point in space. Naturally, this is a simplified description as illustrated in Fig. 2-1.
Z
Fig. 2-1 Interleaving of the E and H fields in space and time in the FDTD formulation.
2-3.1 FDTD method in One-dimensional case
Now we will start with simple one-dimensional differential equations. The time-dependent Maxwell’s curl equations in free space are
t H Here E and H are vectors in three dimensions. When we consider only in one dimension case, E and H simply have Ex and Hy components, so Eq. (2.24) and (2.25) become Above equations mean the electric field oriented in the x direction and the magnetic field
oriented in the y direction both traveling in the z direction. Taking the central difference approximation for both the temporal and spatial derivatives gives
z (2.29) assume that E and H fields are interleaved in both space and time. H uses the arguments and to indicate that the H field values are assumed to be located between the E field values. Similarly, superscript
2
that it occurs slightly after or before n, respectively. Eq. (2.28) and (2.29) can be rearranged as
The calculations are interleaved in both space and time. This is the fundamental paradigm of the finite-difference time-domain (FDTD) method. Eqs. (2.30) and (2.31) are very similar, but because ε0 and µ0 differ by several orders of magnitude. This is circumvented by making the following change of variables:
~ .
0 0 E
E µ
= ε (2.32)
Substituting (2.32) into (2.30) and (2.31) gives
)]
If the cell size ∆z is chosen, the time step ∆ can be determined by t to Eq. (2.30) related to the stability of the FDTD method. An electromagnetic wave propagating in free space cannot go faster than the speed of light. To propagate a distance of one cell ∆z needs a minimum time of
c0
t = ∆z
∆ . When we get to two-dimensional
simulation, we have to allow for the propagation in the diagonal direction, which brings the time requirement to
2c0
t= ∆z
∆ . Obviously, three-dimensional simulation requires
3c0
t ∆z
=
∆ . This is summarized by the well-known “Courant Condition” [31, 32]:
, where d is the dimension of the simulation. Hence we will determine in Eq. (2.37).
This is not necessarily the best formula! Therefore,
∆t
Substituting (2.35) into (2.33) and (2.34), those equations become
1/ 2 1/ 2 0.5 Besides the last two iterative equations, we still need to add incident wave source condition and absorbing boundary condition. It is a great subject in dealing with the wave source condition. For simplicity, we divide it into two categories in 1-D case: hard source and soft source. In a hard source, a propagation wave will see that value and be reflected, because the hard value of Ex looks like a metal wall to FDTD. However a soft source is added to Ex at a certain point and a propagating pulse will just pass through. In calculating
photonic crystal, we must consider the field scattering from the material. Therefore we use a In order to keep outgoing E and H fields from being reflected by the calculation boundary and back into the problem space, so the absorbing boundary conditions (ABC) are necessary to consider. The fields at the edge must be propagating outward. In one time step of the FDTD algorithm it travels
distance = This equation basically explains that it takes two time steps for a wave front to cross one cell.
So a common sense approach tells us that an ABC might be
2 steps before in Ex(0). Boundary conditions such as these have been implemented at both ends of the Ex array. Below are the examples of C computer code in one-dimensional absorbing boundary conditions. Additional parameters are used to store the boundary value for two time steps during the calculation loop.
ex[0] = ex_low_m2;
2-3.2 Two-dimensional photonic crystal formulation and perfectly matched layer (PML) boundary condition
We start again with the normalized Maxwell’s equations:
t H
six different fields. One is the transverse magnetic (TM) mode, which is composed of E~z, , and . Another is the transverse electric (TE) mode, which is composed of , The two-dimensional systemic interleaving of the calculated fields is more complex than one dimension. That is illustrated in Fig. 2-2 below.
Y
Hy
Hy
Put Eqs. (2.48), (2.50) and (2.51) into the finite difference scheme, and take equivalent incremental step in x and y direction, these equations become
))
We have briefly mentioned the issue of absorbing boundary conditions (ABCs) in discussion of one dimension. In the two-dimensional simulations, the program contains two-dimensional matrices for the values of all the fields (i.e. dz, ez, hx and hy). Assume we
X
Fig. 2-2 Interleaving of the E and H fields for the two-dimensional TM formulation.
are simulating a wave generated from a point source propagating in the free space. As the wave propagates outward, it will eventually come to the edge of the allowable space, which is dictated by how the matrices have been dimensioned in the program. If we had done nothing about this, reflections would be produced and would go back the problem space. Then we will have no way to differentiate between the real wave and the reflected wave. This is why the ABCs must exist. The most effective ABCs is the perfectly matched layer (PML) developed by Berenger [32]. How PML works can be easily understood by the following description. If a wave propagating in medium A and it impinges upon medium B, the amount of reflection can be determined by the intrinsic impedances of two media
B A
B
A η
η η η
+
= −
Γ , (2.55)
where the impedance is . ε
η = µ If µ changed with ε so η still remained a constant, Γ would
be zero and no reflection will occur. But this is still helpless to our problem, because waves will continue propagating in the new medium. We really want is a medium that is also lossy so the wave will decrease before it hits the boundary. Hence we mark both ε and µ of complex due to their imaginary parts causing decay.
To simulate a plane wave propagating in a 2D FDTD program, the space of problem will be divided into two regions, the total field and the scattered field (Fig. 2-3). There are two reasons for doing this: (1) The propagating plane wave should not interact with the absorbing boundary conditions; (2) the load on the absorbing boundary conditions should be minimized.
These boundary conditions are not perfect. By subtracting the incident field, the amount of the radiating field hitting the boundary is minimized, thereby reducing the calculation error.
PML
Incident plane wave is subtracted out here.
Total field
Calculation samples
Incident plane wave is subtracted out here.
Scattered field
Fig. 2-3 Total field/scattered field of the two-dimensional problem space.
2-4 Tight binding method in solid state physics
In the later discussion of the coupling between the photonic crystal waveguides we will apply the tight binding approximation to support our argument. So we here do a simple introduction of what is the tight binding method and its meaning in solid state physics [33].
In atoms the electrons are tightly bound to their nuclei. If the atoms are so close that their separations become comparable to the lattice constant in solids, their wave functions will overlap each other. If we consider only two atoms, their combined wave functions are
B
A ψ
ψ ± . The electron energy of state ψA+ψB is lower than one of state ψA−ψB. After they approach to each other, the Coulomb force between nucleuses and electrons can cause the energy level splitting and becomes energy band. The approximation method to obtain the energy band structure by calculating the free atomic wave functions is called tight binding approximation (TB) or linear combination of atomic orbitals (LCAO). In covalently bonded semiconductors the valence electrons are concentrated mainly in the bonds.
Therefore the wave functions of valence electrons should be very similar to bonding orbitals
found in molecules. In addition to being a good approximation for calculating the valence bond structure, the TB method has the advantage that the band structure can be defined in terms of a small number of overlap parameters. The overlap parameters have a simple physical interpretation as representing interactions between electrons of adjacent atoms.
Assume that an electron with ground state ϕ(r) exercises within a single atom’s potential U(r), where ϕ(r) denotes as the s state. It is too complex if we solve the energy band problem by using degenerated atomic energy levels. Therefore, we assume that the influence between two atoms is quite small, and then the wave function can be expanded as following： equation can be expressed as
rj
where T is the primitive vector connecting two lattice points. To calculate the 1st level energy by doing the Hamiltonian matrix diagonalization as follow:
, In Eq. (2.59), we do the integration to only an atom and other atoms nearby which are tied up by ρ . We can rewrite it as:
∫
dVϕ*(r)Hϕ(r)= a− ;∫
dVϕ*(r−ρ)Hϕ(r)=−γ (2.60) To set k k =1, the 1st level energy is∑
⋅ = The relation between overlapping energy γ and atomic spacing ρ in two hydrogen atoms which are both in the 1s state can be clearly calculated. Using the Rydberg-energy unit,, we have Considering to a simple cubic structure, the positions of the closest atoms are
); Other example likes the fcc structure which has twelve closest atoms and its band structure can be described as
2 )
Hence the tight-binding approximation method provides a very simple way to do the atomic energy band structure analysis. This way can also be applied to discussion of the small coupling effect inside a photonic crystal coupled-cavities waveguides (CCWs).