Bloch wave – MoM method

Bloch wave – MoM is a first principles technique for determining the photonic band structure of triply-periodic electromagnetic media such as photonic crystals. It is based on the 3-dimensional spectral domain method,[1] specialized to triply-periodic media. This technique uses the method of moments (MoM) in combination with a Bloch wave expansion of the electromagnetic field to yield a matrix eigenvalue equation for the propagation bands. The eigenvalue is the frequency (for a given propagation constant) and the eigenvector is the set of current amplitudes on the surface of the scatterers. Bloch wave - MoM is similar in principle to the plane wave expansion method, but since it additionally employs the method of moments to produce a surface integral equation, it is significantly more efficient both in terms of the number of unknowns and the number of plane waves needed for good convergence.

Bloch wave - MoM is the extension to 3 dimensions of the spectral domain MoM method commonly used for analyzing 2D periodic structures such as frequency selective surfaces (FSS). In both cases, the field is expanded as a set of eigenfunction modes (either a Bloch wave in 3D or a discrete plane wave - aka Floquet mode - spectrum in 2D), and an integral equation is enforced on the surface of the scatterers in each unit cell. In the FSS case, the unit cell is 2-dimensional and in the photonic crystal case, the unit cell is 3-dimensional.

Field equations for 3D PEC photonic crystal structures

The Bloch wave - MoM approach will be illustrated here for the case of perfectly electrically conducting (PEC) structures admitting only electric current sources, J. However, it can also be readily expanded to dielectric structures, using the well-known interior and exterior equivalent problems commonly employed in ordinary spatial domain method of moments formulations.[2] In dielectric problems, there are twice as many unknowns - J & M - and also twice as many equations to enforce - continuity of tangential E & H - at the dielectric interfaces.[3]

For PEC structures, the electric field E is related to the vector magnetic potential A via the well-known relation:

${\displaystyle \mathbf {E} ~=~-jk\eta \left[\mathbf {A} +{\frac {1}{k^{2}}}\nabla (\nabla \bullet \mathbf {A} )\right]~~~~~~~~~~~~~~(1.1)}$

and the vector magnetic potential is in turn related to the source currents via:

${\displaystyle \nabla ^{2}\mathbf {A} +k^{2}\mathbf {A} ~=~-\mathbf {J} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(1.2)}$

where

${\displaystyle k^{2}~=~\omega ^{2}(\mu \epsilon )~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(1.3)}$

Bloch wave expansion of the fields

To solve equations (1.1) and (1.2) within the infinite periodic volume, we may assume a Bloch wave expansion for all currents, fields and potentials:

${\displaystyle \mathbf {J} (x,y,z)~=~\sum _{mnp}~\mathbf {J} (\alpha _{m},\beta _{n},\gamma _{p})~e^{j(\alpha _{m}x+\beta _{n}y+\gamma _{p}z)}~~~~~(2.1a)}$

${\displaystyle \mathbf {E} (x,y,z)~=~\sum _{mnp}~\mathbf {E} (\alpha _{m},\beta _{n},\gamma _{p})~e^{j(\alpha _{m}x+\beta _{n}y+\gamma _{p}z)}~~~~(2.1b)}$

${\displaystyle \mathbf {A} (x,y,z)~=~\sum _{mnp}~\mathbf {A} (\alpha _{m},\beta _{n},\gamma _{p})~e^{j(\alpha _{m}x+\beta _{n}y+\gamma _{p}z)}~~~(2.1c)}$

where for simplicity, we assume an orthogonal lattice in which α only depends on m, β only depends on n and γ only depends on p. With this assumption,

${\displaystyle \alpha _{m}~=~k_{0}~\sin \theta _{0}~\cos \phi _{0}~+~{\frac {2m\pi }{l_{x}}}~~~~~~~~~~~(2.2a)}$

${\displaystyle \beta _{n}~=~k_{0}~\sin \theta _{0}~\sin \phi _{0}~+~{\frac {2n\pi }{l_{y}}}~~~~~~~~~~~~~(2.2b)}$

${\displaystyle \gamma _{p}~=~k_{0}~\cos \theta _{0}~+~{\frac {2p\pi }{l_{z}}}~~~~~~~~~~~~~~~~~~~~~~(2.2c)}$

and,

${\displaystyle k_{0}~=~2\pi /\lambda ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(2.3)}$

where l_x, l_y, l_z are the unit cell dimensions in the x,y,z directions respectively, λ is the effective wavelength in the crystal and θ₀, φ₀ are the directions of propagation in spherical coordinates.

The quantity k in equations (1.1) and (1.2) comes originally from the time derivative in Maxwell's equations and is the free space propagation constant (actually, the propagation constant of whatever dielectric medium the metallic scatterers are embedded in), proportional to frequency as in equation (1.3). On the other hand, k₀ in the equations above comes from the assumed Bloch wave solution given by equations (2.1) & (2.2). As a result, it represents the propagation constant inside the periodic medium, inversely proportional to wavelength. These two k's, i.e. the free space propagation constant (proportional to frequency) and the propagation constant of the Bloch wave (inversely proportional to wavelength), are in general different thereby allowing for dispersion in the solution. The band diagram is essentially a plot of k as a function of k₀.

The Bloch wave expansions in equations (2.1) are nothing more than exponential Fourier series multiplied by the cell-to-cell propagation factor: ${\displaystyle e^{j(\alpha _{0}x+\beta _{0}y+\gamma _{0}z)}.}$ The Bloch wave expansions are chosen since any field solution within an infinite periodic volume must have the same periodicity as the medium itself, or stated another way, the fields in neighboring cells must be identical up to a (real or complex) propagation factor. In passbands the propagation factor is an exponential function with purely imaginary argument and in the stop bands (or band gaps) it is a decaying exponential function whose argument has a real component.

The wave numbers α₀, β₀ and γ₀ satisfy the relations: ${\displaystyle ~~|\alpha _{0}|~\leq ~{\frac {\pi }{l_{x}}}~,~~~~|\beta _{0}|~\leq ~{\frac {\pi }{l_{y}}}~,~~~~|\gamma _{0}|~\leq ~{\frac {\pi }{l_{z}}}~~~~}$ and outside of these ranges, the bands are periodic.

The Bloch waves are periodic functions of space, with periods l_x, l_y, l_z and the bands are periodic functions of wavenumber, with periods: ${\displaystyle {\frac {2\pi }{l_{x}}}}$, ${\displaystyle {\frac {2\pi }{l_{y}}}}$ and ${\displaystyle {\frac {2\pi }{l_{z}}}}$

Integral equation for PEC media

Substituting equations (2.1) into (1.1) and (1.2) yields the spectral domain Greens function relating the radiated electric field to its source currents:

${\displaystyle \mathbf {E} (\alpha _{m},\beta _{n},\gamma _{p})~=~{\frac {jk\eta }{k^{2}-\alpha _{m}^{2}-\beta _{n}^{2}-\gamma _{p}^{2}}}~\mathbf {G} _{mnp}~\mathbf {J} (\alpha _{m},\beta _{n},\gamma _{p})~~~~~~~~~~~~~~~~~~~~~~~~~(3.1)}$

where,

${\displaystyle \mathbf {G} _{mnp}~=~\left({\begin{matrix}1-{\frac {\alpha _{m}^{2}}{k^{2}}}&-{\frac {\alpha _{m}\beta _{n}}{k^{2}}}&-{\frac {\alpha _{m}\gamma _{p}}{k^{2}}}\\-{\frac {\alpha _{m}\beta _{n}}{k^{2}}}&1-{\frac {\beta _{n}^{2}}{k^{2}}}&-{\frac {\beta _{n}\gamma _{p}}{k^{2}}}\\-{\frac {\alpha _{m}\gamma _{p}}{k^{2}}}&-{\frac {\beta _{n}\gamma _{p}}{k^{2}}}&1-{\frac {\gamma _{p}^{2}}{k^{2}}}\end{matrix}}\right)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(3.2)}$

is the tensor Green's function in the spectral domain. Note that the spatial domain convolution has been transformed into a simple multiplication in the spectral domain, consistent with the convolution theorem for Fourier transforms.

With this equation for the electric field, the electric field boundary condition (requiring that the total tangential electric field be zero on the surface of PEC scatterer) becomes:

${\displaystyle ~\sum _{mnp}~{\frac {1}{k^{2}-\alpha _{m}^{2}-\beta _{n}^{2}-\gamma _{p}^{2}}}~\mathbf {G} _{mnp}~\mathbf {J} (\alpha _{m},\beta _{n},\gamma _{p})~e^{j(\alpha _{m}x+\beta _{n}y+\gamma _{p}z)}~=~\mathbf {0} ~~~~~~~~~~~~~~~~(3.3)}$

Since we are seeking characteristic modes (eigenmodes) of the structure, there is no impressed E-field on the RHS of this electric field integral equation (EFIE). Equation (3.3) is not strictly correct however, since it is only the tangential components of electric field that are actually zero on the surface of the PEC scatterer. This inexactness will be resolved presently, when we test this equation with the electric current basis functions - defined as residing on the surface of the scatterer.

Method of moments (MoM) solution

As is usual in the method of moments, the source currents are now expanded in terms of a sum over some known set of basis functions with unknown weighting coefficients J_j :

${\displaystyle ~\mathbf {J} (x,y,z)~=~\sum _{j}~J_{j}~\mathbf {J} _{j}(x,y,z)~~~~~~~~~~~~~~~~~(4.1)}$

Different structures will have different sets of basis functions for representing the currents on the elements and as in the ordinary spatial domain method of moments, the solution (in this case, the band diagram) is a function of the set of basis functions used.

Substituting (4.1) into (3.3) and then testing the resulting equation with the i-th current basis function (i.e., dotting from the left and integrating over the domain of the i-th current basis function, thereby completing the quadratic form) produces the i-th row of the matrix eigenvalue equation for a 3-dimensional array of PEC scatterers as:

${\displaystyle ~\sum _{j}~J_{j}~\left[~\sum _{mnp}~{\frac {\mathbf {J} _{i}(-\alpha _{m},-\beta _{n},-\gamma _{p})~\mathbf {G} _{mnp}~\mathbf {J} _{j}(\alpha _{m},\beta _{n},\gamma _{p})}{k^{2}-\alpha _{m}^{2}-\beta _{n}^{2}-\gamma _{p}^{2}}}\right]~=~\mathbf {0} ~~~~~~~~~~~~~~~(4.2)}$

As in all MoM formulations, the reaction concept in electromagnetics[2][4] was used in obtaining this equation. The electric field boundary/continuity conditions are "tested" (or enforced) by being integrated against electric current basis functions (for dielectric structures, the magnetic field continuity conditions are additionally tested by being integrated against magnetic current basis functions), and this is how the electric (and magnetic) field boundary conditions are converted into a matrix equation via the method of moments. This process is wholly analogous to that used to decompose a periodic function into its Fourier sine & cosine components, the only difference being that in this case the basis functions are not necessarily orthogonal, merely linearly independent.

This matrix equation is easy to implement and requires only that the 3D Fourier transform (FT) of the basis functions be computed, preferably in closed form.[3] In fact, computing bands of a 3D photonic crystal with this method is no more difficult than computing reflection and transmission from a 2D periodic surface using the spectral domain method . This is because equation (4.2) is identical to the basic EFIE for a freestanding PEC FSS (see Frequency selective surface eq. (4.2)),[5] the only difference being the stronger singularity in 3D which significantly accelerates convergence of the triple sums, and of course the fact that the vectors are now 3-dimensional. As a result, an ordinary PC is sufficient to compute bands of many types of photonic crystals.

It's evident from (4.2) that the EFIE could become singular whenever the free space wavenumber is exactly equal to one of the wave numbers in any of the 3 periodic coordinate directions. This can happen for example when the free space wavelength exactly equals the lattice spacing. This is a statistically rare occurrence in computational practice and corresponds to a propagation anomaly similar to a Wood's reflection anomaly for gratings.

Computing bands

To compute bands of the crystal (i.e. k-k₀ diagrams), successive values of frequency (k) are tried - in conjunction with pre-selected values of propagation constant (k₀) and propagation direction (θ₀ & φ₀) - until a combination is found which drives the determinant of the matrix to zero. Equation (4.2) has been used to compute bands in various types of doped and undoped photonic crystals.[3][6] Not surprisingly, doping photonic crystals with defects provides a means for designing photonic passbands, in just the same way that doping semiconductors with chemical impurities provides a means for designing electronic passbands.

For many subsectional basis functions, such as those having a half-sine or triangular shape along a round wire, the FT of the basis function for negative wave numbers -α, -β, -γ is the complex conjugate of the basis function FT for positive wave numbers. As a result, the matrix in eqn. (4.2) is Hermitian. And as a result of that, only half the matrix needs to be computed. And a second result is that the determinant is a purely real function of the real-valued wavenumber k. Zeroes generally occur at zero-crossings (inflection points, where curvature is zero), so a simple root-finding algorithm such as Newton's method is usually sufficient to find the roots to a very high degree of accuracy. If may still be useful however to plot the determinant as a function of k, to observe its behavior near the zeros.

As a matter of computational convenience, whenever the matrix is larger than 2x2 it's vastly more efficient to compute the determinant either by reducing the matrix to upper triangular form using QR decomposition or to compute the determinant by reducing to echelon form using Gaussian elimination, rather than trying to actually compute the determinant of the matrix directly.

See also

• Bloch wave
• Photonic crystal
• Metamaterial
• Frequency selective surface
• Fourier optics
• Bragg diffraction

Notes

1. Kastner 1987.
2. Harrington 1961.
3. Scott 1998.
4. Rumsey 1954.
5. Scott 1989.
6. Scott 2002.

References

• Harrington, Roger F. (1961), Time-Harmonic Electromagnetic Fields, McGraw-Hill, pp. 106–118
• Kastner, Raphael (1987), On the Singularity of the Full Spectral Green's Dyad, IEEE Trans. on Antennas and Propagation, vol. AP-35, No. 11, pp. 1303–1305
• Rumsey, V. H. (1954), The Reaction Concept in Electromagnetic Theory, Physical Review 94(6):1483-1491, available on researchgate.net
• Scott, Craig (1989), The Spectral Domain Method in Electromagnetics, Artech House, ISBN 0-89006-349-4
• Scott, Craig (1998), Analysis, Design and Testing of Integrated Structural Radomes Built Using Photonic Bandgap Structures, 1998 IEEE Aerospace Conf. Aspen CO, pp. 463 - 479, available on researchgate.net
• Scott, Craig (2002), Spectral Domain Analysis of Doped Electromagnetic Crystal Radomes Using the Method of Moments, 2002 IEEE Aerospace Conf. Big Sky MT, paper #504, pp. 2-957 - 2-963, available on researchgate.net