Plane wave expansion method: Difference between revisions

Revision as of 19:49, 4 March 2012

Latest revision as of 17:22, 9 October 2024

(33 intermediate revisions by 25 users not shown)

Line 1:

{{Short description|Technique in computational electromagnetism}}

{{Dubious|date=August 2009}}

{{Use American English|date=January 2019}}

'''Plane wave expansion method''' (PWE) refers to a computational technique in [[electromagnetics]] to solve the [[Maxwell's equations]] by formulating an [[eigenvalue]] problem out of the equation. This method is popular among the [[photonic crystal]] community as a method of solving for the [[band structure]] (dispersion relation) of specific photonic crystal geometries. PWE is traceable to the analytical formulations, and is useful in calculating modal solutions of Maxwell's equations over an inhomogeneous or periodic geometry.<ref>{{cite journal |last1=Andrianov |first1=Igor V. |last2=Danishevskyy |first2=Vladyslav V. |last3=Topol |first3=Heiko |last4=Rogerson |first4=Graham A. |title=Propagation of Floquet–Bloch shear waves in viscoelastic composites: analysis and comparison of interface/interphase models for imperfect bonding |journal=Acta Mechanica |date=25 November 2016 |volume=228 |page=1177-1196 |doi=10.1007/s00707-016-1765-4}}</ref> It is specifically tuned to solvingsolve problems in a time-harmonic forms, with [[Dispersion (optics)|non-dispersive]] media (a reformulation of the method named Inverse dispersion allows frequency-dependent refractive indices<ref >{{cite journal| title= Inverse dispersion method for calculation of complex photonic band diagram and PT symmetry |journal= Physical Review B|pages=165132 |year=2016|doi= 10.1103/PhysRevB.93.165132 | last1=Rybin |first1=Mikhail |last2=Limonov |first2=Mikhail |volume=93|issue= 16|arxiv=1707.02870 }}</ref>).

==Principles==

{{Dubious|date=August 2009}}

[[Plane wave]]s are solutions to the homogeneous [[Helmholtz equation]], and form a basis to represent fields in the periodic media. PWE as applied to photonic crystals as described is primarily sourced from Dr. Danner's tutorial.<ref>{{Cite web |last=Danner |first=Aaron J. |date=2011-01-31 |title=An introduction to the plane wave expansion method for calculating photonic crystal band diagrams |url=https://www.ece.nus.edu.sg/stfpage/eleadj/planewave.htm |url-status=live |archive-url=https://web.archive.org/web/20220615161702/https://www.ece.nus.edu.sg/stfpage/eleadj/planewave.htm |archive-date=2022-06-15 |access-date=2022-09-29 |website=Aaron Danner - NUS}}</ref>

from Dr. Danner's tutorial<ref>http://www.ece.nus.edu.sg/stfpage/eleadj/planewave.htm</ref>.

The electric or magnetic fields are expanded for each field component in terms of the [[fourier series]] components along the reciprocal lattice vector. Similarly, the dielectric permittivity (which is periodic along reciprocal lattice vector for photonic crystals) is also expanded through Fourier series components.

The electric or magnetic fields are expanded for each field component in terms of the [[fourierFourier series]] components along the reciprocal lattice vector. Similarly, the dielectric permittivity (which is periodic along reciprocal lattice vector for photonic crystals) is also expanded through Fourier series components.

:<math>

\frac{1}{\epsilon_r} = \sum_{m=-\infty}^{+\infty} K_m^{\epsilon_r} e^{-j \vec{G}.\vec{r}}

</math>

<math display="block">\frac{1}{\epsilon_r} = \sum_{m=-\infty}^{+\infty} K_m^{\epsilon_r} e^{-ji \vecmathbf{G}. \cdot \vecmathbf{r}}</math>

:<math>

<math display="block">E(\omega,\vecmathbf{r}) = \sum_{n=-\infty}^{+\infty} K_n^{E_y} e^{-ji\vecmathbf{G}. \cdot \vecmathbf{r}} e^{-ji\vecmathbf{k} \cdot \vecmathbf{r}}</math>

</math>

with the fourierFourier series coefficients being the K numbers subscripted by m, n respectively, and the reciprocal lattice [[Vector (mathematics and physics)|vector]] given by <math>\vecmathbf{G}</math>. In real modeling, the range of components considered will be reduced to just <math>\pm NmaxN_\max</math> instead of the ideal, infinite wave.

Using these expansions in any of the curl-curl relations like,

<math display="block">\frac{1}{\epsilon(\vecmathbf{r})} \nabla \times \nabla \times E(\vecmathbf{r},\omega) = \left( \frac{\omega}{c} \right)^2 E(\vecmathbf{r},\omega)</math>

:<math>

and simplifying under assumptions of a source free, linear, and non-dispersive region we obtain the [[eigen valueeigenvalue]] relations which can be solved.

\frac{1}{\epsilon(\vec{r})} \nabla \times \nabla \times E(\vec{r},\omega) = \left( \frac{\omega}{c} \right)^2 E(\vec{r},\omega)

</math>

and simplifying under assumptions of a source free, linear, and non-dispersive region we obtain the [[eigen value]] relations which can be solved.

==Example for 1D case==

For a y-polarized z-propagating electric wave, incident on a 1D-DBR periodic in only z-direction and homogeneous along x,y, with a lattice period of a. We then have the following simplified relations:

[[Image:Photonic Crystal 1D DBR aircore epsr12point25 DbyA0point8.png|thumb|right|Band structure of a 1D Photonic Crystal, DBR air-core calculated using plane wave expansion technique with 101 planewaves, for d/a=0.8, and dielectric contrast of 12.250.]]

For a y-polarized z-propagating electric wave, incident on a 1D-DBR periodic in only z-direction and homogeneous along x,y, with a lattice period of a. We then have the following simplified relations:

:<math>

<math display="block">\frac{1}{\epsilon_r} = \sum_{m=-\infty}^{+\infty} K_m^{\epsilon_r} e^{-ji \frac{2\pi m}{a}z}</math>

<math display="block">E(\omega,\vecmathbf{r}) = \sum_{n=-\infty}^{+\infty} K_n^{E_y} e^{-ji\frac{2\pi n}{a}z} e^{-ji \vecmathbf{k} \cdot \vecmathbf{r}}</math>

</math>

:<math>

E(\omega,\vec{r}) = \sum_{n=-\infty}^{+\infty} K_n^{E_y} e^{-j\frac{2\pi n}{a}z} e^{-j\vec{k}\vec{r}}

</math>

The constitutive eigeneigenvalue equation we finally have to solve becomes,

<math display="block">\sum_n{\left( \frac{2\pi n}{a} + k_z \right)\left( \frac{2\pi m}{a} + k_z \right) K_{m-n}^{\epsilon_r} K_{n}^{E_y}} = \frac{\omega^2}{c^2} K_{m}^{E_y}</math>

:<math>

\sum_n{\left( \frac{2\pi n}{a} + k_z \right)\left( \frac{2\pi m}{a} + k_z \right) K_{m-n}^{\epsilon_r} K_{n}^{E_y}} = \frac{\omega^2}{c^2}K_{m}^{E_y}

</math>

This can be solved by building a matrix for the terms in the left hand side, and finding its eigen valueeigenvalue and vectors. The eigen valueseigenvalues correspond to the modal solutions, while the corresponding magnetic or electric fields themselves can be plotted using the fourierFourier expansions. The [[coefficients]] of the field harmonics are obtained from the specific eigen vectorseigenvectors.

The resulting band-structure obtained through the eigen modeseigenmodes of this structure are shown to the right.

===Example code===

We can use the following code in [[MatlabMATLAB]] or [[GNU Octave]] to compute the same band structure,

<!--

% this code is released to Public Domain

% (c) 2007 Muthiah Annamalai

% Author: Muthiah Annamalai,

%

% This code was written based on Prof. Danner's resource on PWE:

Line 60 ⟶ 44:

-->

<syntaxhighlight lang="matlab">

<pre><nowiki>

%

% solve the DBR photonic band structure for a simple

Line 66 ⟶ 50:

% we assume an infinite stack of 1D alternating eps_r|air layers

% y-polarized, z-directed plane wave incident on the stack

% periodic in the z-direction;

%

% parameters

d = 8; % air gap

a = 10; % total periodicity

d_over_a = d / a;

eps_r = 12.2500; % dielectric constant, like GaAs,

% max F.S coefs for representing E field, and Eps(r), are

Mmax = 50;

% Q matrix is non-symmetric in this case, Qij != Qji

% Qmn = (2*pi*n + Kz)^2*Km-n

% Kn = delta_n / eps_r + (1 - 1/eps_r) (d/a) sinc(pi.n.d/a)

% here n runs from -Mmax to + Mmax,

freqs = [];

for Kz = - pi / a:pi / (10 * a): + pi / a

Q = zeros(2 * Mmax + 1);

for x = 1:2 * Mmax + 1

for y = 1:2 * Mmax + 1

X = x - Mmax;

Y = y - Mmax;

kn = (1 - 1 / eps_r) * d_over_a .* sinc((X - Y) .* d_over_a) + ((X - Y) == 0) * 1 / eps_r;

Q(x, y) = (2 * pi * (Y - 1) / a + Kz) .^ 2 * kn; % -Mmax<=(Y-1)<=Mmax

end

end

fprintf('Kz = %g\n', Kz)

omega_c = eig(Q);

omega_c = sort(sqrt(omega_c)); % important step.

freqs = [freqs; omega_c.'];

end

close()

fprintf('Kz = %g\n',Kz)

figure()

omega_c=eig(Q);

omega_c=sort(sqrt(omega_c));%important step.

freqs=[freqs; omega_c.'];

end

close()

figure()

hold on

idx = 1;

for idx = 1:length(- pi / a:pi / (10 * a): + pi / a)

plot(- pi / a:pi / (10 * a): + pi / a, freqs(:, idx), '.-')

end

hold off

xlabel('Kz')

ylabel('omega/c')

title(sprintf('PBG of 1D DBR with d/a=%g, Epsr=%g', d / a, eps_r))

</syntaxhighlight>

</nowiki></pre>

==Advantages==

PWE expansions are rigorous solutions. PWE is extremely well suited to the modal solution problem. Large size problems can be solved using iterative techniques like [[Conjugate gradient method]].

For both generalized and normal eigen valueeigenvalue problems, just a few band-index plots in the band-structure diagrams are required, usually lying on the [[brillouin zone edges]]. This corresponds to eigen modeseigenmodes solutions using iterative techniques, as opposed to diagonalization of the entire matrix.

The PWEM is highly efficient for calculating modes in periodic dielectric structures. Being a Fourier space method, it suffers from Gibb'sthe [[Gibbs phenomenon]] and slow convergence in some configuration when fast Fourier factorization is not used. It is the method of choice for calculating the band structure of photonic crystals. It is not easy to understand at first, but it is easy to implement.

==Disadvantages==

{{Dubious|date=August 2009}}

Sometimes spurious modes appear. Large problems scaled as ''O''(''n''³)'', with the number of the plane waves (''n'') used in the problem. This is both time consuming and complex in memory requirements.

Alternatives include Order-N spectral method, and methods using [[Finite-difference time-___domain]] (FDTD) which are simpler, and model transients.

If implemented correctly, spurious solutions are avoided. It is less efficient when index contrast is high or when metals are incorporated. It cannot be used for scattering analysis.

Being a Fourier-space method, Gibbs phenomenon affects the method's accuracy. This is particularly problematic for devices with high dielectric contrast.

==See also==

Line 143 ⟶ 129:

[[Category:Computational science]]

[[Category:Electrodynamics]]

[[Category:Computational electromagnetics]]

{{electromagnetism-stub}}

[[fr:Méthode de décomposition en onde plane]]