The Plane-Wave Expansion Method

After a decade of research in photonic-band-gap materials, several techniques have been pro-posed for solving Maxwell’s equations [John, S., et al. (1988); Pendry, J. B., et al. (1992); Qiu, M., et al. A (2000)]. Nevertheless, the plane-wave expansion method has become the standard de facto for computing the band structure of semiconductor-based photonic crystals. The method is based on the truncation of the sum in Eq. (1.11), reducing the master equation to a matrix eigenvalue problem. The coefficients c_n(k+G) and the eigenfrequencies ω_n(k) are obtained by standard numerical diagonalization of the resulting “hamiltonian”.

Since the basis functions f_k(r) are plane waves, the Eq. (1.11) is explicitly written as is the approximation imposed by the numerical method; in fact, it would be impossible to store infinite arrays in the computer memory. Eq. (1.20) is nothing but the Fourier expansion of the magnetic field truncated to a cut-off and c_σ(k+G) are its Fourier coefficients. The next step is to rewrite the master equation (1.5a) in the Fourier space by calculating the matrix elements of the

“hamiltonian” operator Oˆ on the plane-wave basis. The result is X

and the matrix [[η]] = [[²]]⁻¹ is the inverse of the dielectric function Fourier transform

²_G,G⁰ = 1 A_c

Z

Ac

²(r)e^{ı(G−G0)·r}dr, (1.23)

whereA_cis the space occupied by the unit cell. [[H]] is a square matrix with dimensions 2N×2N, [[η]] and [[²]] have dimensions N ×N instead; N is the number of G vectors below the cut-off.

Likewise the “hamiltonian” operator Oˆ associated to the master equation, H is hermitian with non-negative eigenvalues.

Eq. (1.22) is the eigenvalue problem for [[H]]. Standard numerical diagonalization of [[H]] yields the eigenfrequencies ω_n(k) and, optionally, the coefficients c_n(k+G). Notice that for each diagonal-ization, the routine outputs a set ω_n(k), withn= 1, . . . ,2N, corresponding to the energies of 2N bands for a fixed Bloch vectork. In order to calculate the whole band structure, the operation has to be repeated for a certain ensemble of kvectors in the B.Z., usually the edges of the irreducible B.Z..

The truncation of the sum in Eq. (1.20) is the trick that allowed numerical solution of the mas-ter equation. In fact, the exact Fourier transform of the operator Oˆ would be a matrix of infinite dimensions, whereas [[H]] is limited to 2N ×2N. For this reason, the matrix eigenvalue problem is an approximation of Maxwell’s equations and an error will occur in the calculated eigenfrequen-cies. Calling λ_n(k) the true eigenvalue and λ^(N)_n (k) the eigenvalue calculated with the numerical routine, the error committed in evaluating the band energies will be ∆^(N)_n (k) =|λ_n(k)−λ^(N_n ⁾(k)|.

As N → ∞, the “hamiltonian” [[H]]^(2N) should approach the true Fourier transform of Oˆ and, consequently, the band energies should converge to the true values, i.e. ∀k,∀n; ∆^(N)_n (k) → 0 as N → ∞. Now the issue is: is this statement true or not? Also, if true, how fast is the convergence?

The core of the problem consists in the Fourier transform of the dielectric function ²(r), which enters the operatorOˆ as 1/²(r). When the master equation is transformed into the Fourier space, there are two choices for the dielectric function: take the Fourier transform of 1/²(r)⇒[[η]] = [[1/²]]

or take the Fourier transform of ²(r) and invert the matrix ⇒[[η]] = [[²]]⁻¹. Contrary to intuition, as N → ∞, [[1/²]]^(N) 6= [[²]]^(N)−1. This result stems from the jump discontinuities of the dielec-tric function at the interfaces among the media of the photonic crystal. The inequality signifies that the convergence of the plane-wave expansion method is not a trivial problem at all. In a paper on the factorization of Fourier series of discontinuous periodic functions, Li, L. (1996) gave mathematical foundations that the correct choice for uniform convergence is [[η]] = [[²]]⁻¹, which is called the inverse rule. The fact is that [[²]]⁻¹ uniformly satisfies the boundary conditions for the electromagnetic field, while [[1/²]] does it only non uniformly. The idea of taking the inverse of the Fourier transform has been used since Ho, K. M.,et al. (1990), but without the rigorous motivation given by Li, L. (1996). That is why it is also known as Ho’s method. For instance, considering two-dimensional photonic crystals, the ∆^(N)n (k) caused by [[η]] = [[²]]⁻¹ is below 1% already withN of the order of 100. With [[η]] = [[1/²]], ∆^(N)_n (k) is still above 1% for more than 300 plane waves.

The plane-wave expansion method is thus able to output accurate eigenfrequencies with a moderate CPU time, provided the inverse rule is used. Notice that standard diagonalization and inversion routines areO(N³) operations; i.e. doubling the number of plane waves increases the CPU time by eight times. That is why the method is very efficiently applied to two-dimensional photonic crys-tals. For three-dimensional photonic crystals, where the number ofGvectors is of the order of one

thousand or more, the method start to approach its numerical limit as regards the CPU time and the accuracy of the diagonalization process. Nevertheless, improved algorithms of the plane-wave expansion method allow to use thousands of Gvectors for an extremely accurate determination of the eigenfrequencies [Meade, R. D., et al. (1993)].

It has been shown that the electromagnetic problem for a photonic crystal can be treated within an operator formalism by recasting Maxwell’s equations in a closed form for either the electric either the magnetic fields. The translational symmetry of the dielectric function implies that the solutions have to be Bloch waves. The frequency spectrum is organized in the band structure picture, with the classification of the energy levels in terms of Bloch vector k and band index n. The band structure is obtained by numerical solution of the master equation by means of the plane-wave expansion method. These concepts and tools can now be applied to the study of semiconductor-based two-dimensional photonic crystals.