In this post, I seek to enlighten the reader on a particularly niche but useful semi-analytical technique used in Computational Electromagnetism(CEM). While this is my first post on CEM, I have decided to cover this method first as learning this was particularly memorable for me. Before continuing, I strongly encourage the reader to familiarise themselves with maxwell’s laws and fourier transforms(esp fourier series). If you want to go further and explore more stuff with this method, I highly suggest checking out the youtube channel EMpossible where I learnt most of my CEM techniques. In this particular post I do not just seek to regurgitate what I have been taught but to share my insights on an intuitive derivation, specifically for the 1D, and 2D PWEM methods. The motivation of the PWEM is to solve for the propagation of light or EM waves in a periodic medium. Usually such materials are referred to as photonic crystals for optical frequencies and if metallic structures are used for microwave to radiowave frequencies, metamaterial like behaviour is noticeable. Due to the periodic nature of the medium, this creates an analogous effect to periodic potentials in semiconductors – electron waves under the influence of a periodic potential split into electron bands and likewise, EM waves also split into optical bands when propagating through a periodic medium.

1. Bloch States

Bloch states are the solutions to the problem. Consider an ordinary scalar wave equation for a wave with an angular frequency \(\omega\):

\(\nabla^2 \psi = -(\frac{\omega}{c})^2\psi\)

And now if a medium is periodic, the phase speed of the wave, \(c\) varies periodically throughout the medium. The ansatz to this is:

\(\psi(\vec{r}) = u(\vec{r})e^{i\vec{k}\bullet\vec{r}}\)

In this equation, the function \(u(\vec{r})\) has the same periodicity as the medium and this entire wave propagates through the medium with the wavevector \(\vec{k}\). Using Fourier’s theorm, any periodic function can be represented as a fourier series:

\(u(\vec{r}) = \sum_{\vec{G}}a_{\vec{G}}e^{i\vec{G}\bullet \vec{r}}\)

\(\vec{G}\) refers to the valid reciprocal lattice vectors for the system. For example, in a simple 1D case where the system has a spatial period of a, the reciprocal lattice vectors are integer multiples of \(\frac{2\pi}{a}\). For a 2D system where the axes are orthogonal, the lattice vectors are also straightforward. If the periods along the x and y axes are a and b respectively, then the reciprocal lattice vectors are: \(\frac{2\pi}{a}\), \(\frac{2\pi}{b}\). For the more general case where the system is 3D or the lattice vectors aren’t orthogonal, then given any 3 lattice vectors (\(\vec{a_1},\vec{a_2},\vec{a_3}\)), the corresponding reciprocal lattice vectors are:

\(\vec{b_1} = \frac{2\pi (\vec{a_2}\times\vec{a_3})}{\vec{a_1}\bullet(\vec{a_2}\times\vec{a_3})}\)

\(\vec{b_1} = \frac{2\pi (\vec{a_3}\times\vec{a_1})}{\vec{a_2}\bullet(\vec{a_3}\times\vec{a_1})}\)

\(\vec{b_3} = \frac{2\pi (\vec{a_1}\times\vec{a_2})}{\vec{a_3}\bullet(\vec{a_1}\times\vec{a_2})}\)

This construction ensures that \(\vec{b_i}\bullet\vec{a_j} = 2\pi\delta_{ij}\) holds and this also applies to orthogonal coordinates. For 2D materials, this method is also valid and \(\vec{a_3}\) is set to be \(\hat{z}\). To find the fourier coefficients \(a_{\vec{G}}\) and inner product is done over one unit cell:

\(a_{\vec{G}} = \frac{1}{\Omega}\int_{unit cell}u(\vec{r})e^{-i\vec{G}\bullet \vec{r}}d\Omega\)

For a simple example, a 2D square lattice with sides a will have \(\Omega = a^2\) and \(d\Omega = dxdy\). For a 1D case, \(\Omega = a\) and \(d\Omega = dx\). Now that I have explained periodicity, it is time to apply the concepts to the EM bloch waves.

2. Maxwell’s equations and 1D case

I will be mainly using the last 2 maxwell equations to derive the necessary wave equations and I will be covering the derivation for the PWEM for 2D and 1D systems. The last 2 maxwell equations read:

\(\nabla\times\vec{E} = -\mu_0\partial_t\vec{H}\)

\(\nabla\times\vec{H} = \epsilon_0\epsilon_r\partial_t\vec{E}\)

I have assumed non-magnetic media and this is particularly applicable for most photonic applications and the materials are all linear. Consider the simple 1D case where there is only periodicity in 1 direction and the material is infinitely wide along the other 2 axes. The dielectric constant has its own fourier expansion given by:

\(\epsilon(z) = \sum_{m=-\infty}^{m=\infty}\xi_me^{i\frac{2\pi m}{a}z}\)

Applying a curl to the first equation and using the results of the second equation gives:

\(\nabla\times(\nabla\times\vec{E}) = (\frac{\omega}{c})^2\epsilon_r(\vec{r})\vec{E}\)

If the polarisation of the material is always perpendicular to the propagation axis, then the E field only has one component and it can be written as:

\(E(z) = e^{ikz}\sum_{m=-\infty}^{m=\infty}E_{m}e^{i\frac{2\pi m}{a}z}\)

Putting this in the equation gives:

\((k+2\pi m/a)^2E_{m} = (\frac{\omega}{c})^2\sum_{m’}\xi_{m-m’}E_{m’}\)

This is equivalent to solving the generalized eigenvalue problem:

\(D\vec{v}=(\frac{\omega}{c})^2M\vec{v}\)

Where \(D = diag(k + 2\pi m/a)^2\) for all values of m used in the truncation and \(\vec{v} = (E_{-M}…,E_{-1}, E_0,E_1,….,E_{M})^T\), the vector of all fourier coefficients of the(only) electric field component. The M matrix is a mixing matrix with toeplitz structure: \(M_{ij} = \xi_{i-j}\) and this holds if the fourier components are arranged from -M to M where M is the truncation limit of the expansion. Such a method yields 2M + 1 different bands and likewise different values of \(\frac{\omega}{c}\). Notice that the matrix D changes with the wave vector k. Hence, to find the full dispersion relation this had to be done for all values of \(k\in [\pi/a,\pi/a]\)(this is the first brillounin zone). Note that for every point in k space you choose to include, you are doing one additional eigenvalue solve of complexity ~ (2M+1)^3 which is pretty massive. Fortunately the matrices are sparse so please use sparse solvers whenever possible to speed things up. An example of a calculation is shown below:

Notice that for some ranges of \(\omega/c\) there seem to be no valid k values. This is commonly referred to as an optical bandgap. Just like how an electronic band gap indicates that there is a gap in the energy states and electrons cannot have those energy levels, the optical/photonic band gap indicates that EM waves with those range of frequencies cannot propagate inside the material. For this 1D case, a thick distributed bragg reflector at normal incidence behaves almost like a 1D photonic crystal and it relfects very strongly and uniformly for those frequencies within the bandgap.

3. The 2D PWEM

For this situation, I will be assuming a square or rectangular cell for convenience of calculation. Infact, even for a triangular/hexagonal lattice, the unit cell can be chosen to be a conventional unit cell. The result for this case is a bit more special and there are 2 possible solutions – one with the E field vector perpendicular to the plane(where the periodicities are) and the other being H. Thus, for this case, polarisations matter. Applying the curl to one maxwell equation and then substituiting in the othermaxwell equation gives 2 equations:

\(\nabla\times(\nabla\times\vec{E}) = (\frac{\omega}{c})^2\epsilon_r(\vec{r})\vec{E}\)

\(\nabla\times(\frac{1}{\epsilon_r(\vec{r})}\nabla\times\vec{E}) = (\frac{\omega}{c})^2\vec{E}\)

You can first consider the case that is the medium is periodic in the x-y plane and the EM waves only propagate in the plane. The consider the 2 cases H//z and E//z. Working through the equations will show that the first equation describes E//z and the second the case when H//z. The following expansions need to be calculated using the inner product earlier:

\(\epsilon_r(x,y) = \sum_{m,n=-\infty}^{m,n=\infty}\xi_{m,n}e^{i\frac{2\pi m}{a}x}e^{i\frac{2\pi n}{b}y}\)

\(\frac{1}{\epsilon_r}(x,y) = \sum_{m,n=-\infty}^{m,n=\infty}\eta_{m,n}e^{i\frac{2\pi m}{a}x}e^{i\frac{2\pi n}{b}y}\)

Now, only the relevant field components need to be thrown into the equation making it a ‘scalar’ wave equation which is much easier to solve. Following ansatz are used:

\(A(x,y) = \sum_{m,n=-\infty}^{m,n=\infty} a_{m,n}e^{i\frac{2\pi m}{a}x}e^{i\frac{2\pi n}{b}y}\)

where \(A, a_{m,n}\) are the spatial representations of the z component of E/H and the fourier representations respectively. Putting this ansatz into the relevant equation yields two different eigenvalue problems to be solved:

E-mode(generalized eigenvalue problem):

\([(k_x+2\pi m/a)^2 + (k_y+2\pi n/b)^2]E_{mn} = (\frac{\omega}{c})^2\sum_{m’,n’}\xi_{m-m’,n-n’}E_{m’,n’}\)

H-mode(eigenvalue problem):

\(\sum_{m’,n’}[(k_x+2\pi m/a)(k_x+2\pi m’/a) + (k_y+2\pi n/b)((k_y+2\pi n’/b))]\eta_{m-m’,n-n’}E_{m’n’}\)

\(=(\frac{\omega}{c})^2E_{m,n}\)

It is clear that the problems are not symmetric and hence, different polarisations will have dispersion diagrams. Typically, the dispersion diagrams for 2D and 3D systems will only go over a certain line in the brillounin zone giving diagrams such as:

I wont be going into the specifics of choosing the set of k values where the calculations are done and infact, for a 2D system, you could also choose to cover the entire 1st brillounin zone and your plot will just be sheets. But in the diagram above, I have chosen a set of points that form a line that cover the high symmetry points(\(\Gamma, X, K,M\)) These High symmstry points are often located at the edges or vertices of the Brillounin zone and thus shed some light on the properties of the material. For more information on the points, I suggest watching EMpossible or consulting a solid state physics textbook for the lattice points often used.

4. Ending

With that I conclude my post on the PWEM. I have managed to cover the intuitive quick derivations for the equations that can be used to give the necessary band diagrams for 1D and 2D cases. As for the additional stuff relating to the lattice structures, reciprocal lattices and high symmetry points in the brillounin zone, it is a pity that I didnt have time to cover those and they deserve an entire post dedicated to them.