Intro

Since I graduated in Electrical Engineering, I always wanted to simulate what I think is THE most fundamental system of equations for us (Electrical Engineers): The Maxwell’s equations. A couple of years back I finally had time to do so, and I’ll share with you the work I did. Everything is on GitHub, so if you want to try yourself, please feel free to fork the repository and play around!

In this post I’ll assume that you have at least a vague idea of what Maxwells equations are and that they describe how electromagnetic fields behave. Also, I assume that you know that light is a “moving electromagnetic field” (called electromagnetic wave). So, if we are able to describe and simulate electromagnetic waves, we can actually describe and simulate light!

Before show the simulations, let me start with a bit of explanation about the electromagnetic theory. The most obvious place to start is with Maxwell’s equations themselves: $$ \begin{equation} \begin{gathered} \nabla \cdot \vec E = \frac{\rho }{{{\varepsilon _0}}} \hfill \\ \nabla \cdot \vec B = 0 \hfill \\ \nabla \times \vec E = – \frac{{\partial \vec B}}{{\partial t}} \hfill \\ {c^2}\nabla \times \vec B = \frac{{\vec J}}{{{\varepsilon _0}}} + \frac{{\partial \vec E}}{{\partial t}} \hfill \\ \end{gathered} \end{equation} $$ These equations together with a couple of technical details like continuity constraints and boundary conditions, simply describes ALL physics related to electromagnetic fields. That includes radio waves, light, electricity, magnetism and etc. It is not the goal of this post to talk about Maxwell’s equation themselves, but rather about its solutions. In particular, about numeric solutions. Also, I do not intend to go deep into the theory of the simulation method, but rather to show its potential and encourage you to try! This post is not a tutorial, but I’ll try to point you in the right direction if you want to do it yourself. Also, I have to assume that you know some Calculus (at some not so basic level, I’m afraid).

First of all, lets understand what a solution to Maxwell’s equation means. Any equation that describes a physical phenomenon is a kind of mathematical description of that phenomenon. For instance, the famous Newton’s equation \(F = m\,a\) tells us that if you have a mass of \(10Kg\) and it is accelerating with an acceleration of \(2m/s^2\), there should be a force of \(20N\) being applied to the mass. In the same way, Maxwells equation tells us that if a medium has electrical permittivity \(\varepsilon_0\), the charge density is \(\rho\), etc, etc, then the electric and magnetic fields MUST be such that the divergent of the electric field is \(\rho / \varepsilon_0 \) and etc. For instance, Maxwell’s equation directly states that, whatever the situation is, the divergent of the magnetic field is zero. And so on…

So, that means that if you know the environment, in another words, you know how its permittivity, permissively, charge density, etc are distributed, you can predict how the electromagnetic fields will behave on this environment. In order to do that, you solve Maxwell’s equations to find \({\vec E}\) and \({\vec B}\).

It turns out that Maxwell’s equations, although they have a “simple” look, they are kind of “extensive” in the sense that the electric and magnetic fields have each one 3 components x, y and z and each component varies in space and time. Also they are partial differential equations, which involve all the derivatives of the solution either in x, y, z and t! So, to solve the Maxwell’s equation mean to find 6 functions of 4 variables. In general that is not possible to do analytically except in special cases. Hence, we have to appeal to numeric solutions. There are tons of methods for numerically solve a partial differential equation. In this post we will talk about the FD-TD method (due to Kane S. Yee) or method of Finite Diferences in the Time Domain (sometimes called Yee method). The method itself is pretty straightforward and, as all numerical method, it begins by exchanging the functions themselves by a grid of discrete values. As I mention, I won’t go deep into the method itself, but I’ll try to give an idea of how it works. Another important thing: I implemented the 2D version of the FD-TD method. So, the simulation will consist in a fields that is the same in the z direction form minus infinity to infinity. Unfortunately, cool things like polarization of light won’t be possible here (by definition, a 2D wave is aways polarized). But don’t worry, I promise that you can still have A LOT of fun with this apparently frustrating limitation.

In summary up to now, we will solve Maxwell’s equation numerically using the FD-TD method. So, what does this mean in practice? That means that we are going to divide our environment in a grid (2D for our case) and set each point on the grid with objects like conductors, sources, walls and etc. After that, we compute the electric and magnetic field on each point on the grid at each instant of time. The result is that, hopefully, we are going to see an “animation” of the fields interacting with our small “world”. Figure 1 tries to show a possible setup for a simple simulation.

Figure 1: 2D grid with a source of electric field and some objects with different permittivities conductivities.

Finite-Difference Time-Domain method

There is still one “trick” that must be done in the numerical setup for the FD-TD method. Since se have relationships with the curl of the fields, we have to consider the resulting vector to the in the “center” of the diferences. So, in the TM approach, the electric field have to be in the center of the magnetic fields (and no magnetic field must be present). So, for the method to work, we have to “stagger” the grid, moving it “half step” in x and y for one of the fields. So, this produces a grid for the electric field that is halfway displayed from the magnetic field. This is important for the assignment of the spacial indices for the electric and magnetic fields. Figure 2 shows the staggered grid that represents the electric and magnetic fields in the TM configuration.

Figure 2: Staggered grid for the electric and magnetic fields in the TM scenario

Since we dealing with the 2D simulation, this simplifies the difference equations for the whole set of Maxwell’s equations. In fact we end up with this set of difference equations[3] bellow. Observe that the order of computation here matters… This is a technical detail that has to do with the stability of the method. Since we are integrating, if we do not compute things right, we might generate a systematic error in one of the fields, causing it to integrate more then necessary and blow up the values. $$ \begin{array}{l} {H_x}\left[ {m,n + \frac{1}{2},q + \frac{1}{2}} \right] = {H_x}\left[ {m,n + \frac{1}{2},q – \frac{1}{2}} \right] – \\ \frac{{\Delta t}}{{\mu \Delta y}}\left( {{E_z}\left[ {m,n + 1,q} \right] – {E_z}\left[ {m,n,q} \right]} \right) \end{array} $$ $$ \begin{array}{l} {H_y}\left[ {m + \frac{1}{2},n,q + \frac{1}{2}} \right] = {H_y}\left[ {m + \frac{1}{2},n,q – \frac{1}{2}} \right] + \\ \frac{{\Delta t}}{{\mu \Delta x}}\left( {{E_z}\left[ {m + 1,n,q} \right] – {E_z}\left[ {m,n,q} \right]} \right) \end{array} $$ $$ \begin{array}{l} {E_z}\left[ {m,n,q + 1} \right] = \frac{{2\varepsilon – \rho \Delta t}}{{2\varepsilon + \rho \Delta t}}{E_z}\left[ {m,n,q} \right] + \\ \frac{{2\varepsilon }}{{2\varepsilon + \rho \Delta t}}\left( \begin{array}{l} \frac{{\Delta t}}{{\varepsilon \Delta x}}\left( {{H_y}\left[ {m + \frac{1}{2},n,q + \frac{1}{2}} \right] – {H_y}\left[ {m – \frac{1}{2},n,q + \frac{1}{2}} \right]} \right)\\ – \frac{{\Delta t}}{{\varepsilon \Delta y}}\left( {{H_x}\left[ {m,n + \frac{1}{2},q + \frac{1}{2}} \right] – {H_x}\left[ {m,n – \frac{1}{2},q + \frac{1}{2}} \right]} \right) \end{array} \right) \end{array} $$ Now let me just comment a bit about the notation. \(E_z[m,n,q]\) is the electric field in z at index m and n in the grid at instant q in time (same for the other fields). \(\rho\), \(\epsilon\) and \(\mu\) are the conductivity, permittivity and permeability of the medium. Although not explicitly in the equations, each one also depend on the position on the grid m and n (\(\rho[m,n]\), \(\epsilon[m,n]\) and \(\mu[m,n]\)). The deltas are the spacing for the discretization in space and time.

Actually, the “algorithm” of FD-TD is pretty much only that! Its a loop that keeps updating the fields according to those equations. The heart of the implementation is the step that do that update. The code is as simple as this one

function [Ezo, Hxo, Hyo] = leapFrog2D(Ezi, Hxi, Hyi, parameters)

    Ezo = Ezi;
    Hxo = Hxi;
    Hyo = Hyi;


    c1 = parameters.c1;
    c2 = parameters.c2;
    c3 = parameters.c3;
    c4 = parameters.c4;
    dx = parameters.dx;
    dy = parameters.dy;
    boundary = parameters.boundary; 
    wire = parameters.wire;
    
   
    Nx = size(Ezi,1);
    Ny = size(Ezi,2);
    
    j_p_0_E = 1:Nx-1;
    j_p_0_E = 1:Ny-1;

    i_p_1_E = 2:Nx;
    j_p_1_E = 2:Ny;

    i_p_0_H = 1:Nx-1;
    j_p_0_H = 1:Ny-1;

    i_p_1_H = 2:Nx;
    j_p_1_H = 2:Ny;

   
    % transverse electric
  
    Ezo(i_p_1_E, j_p_1_E) = c1(i_p_1_E, j_p_1_E).*Ezi(i_p_1_E, j_p_1_E) + c2(i_p_1_E, j_p_1_E).*(  ( wire(i_p_1_E, j_p_1_E).*Hyi(i_p_1_E, j_p_1_E) - wire(i_p_0_E, j_p_1_E).*Hyi(i_p_0_E, j_p_1_E) )/dx - ( wire(i_p_1_E, j_p_1_E).*Hxi(i_p_1_E, j_p_1_E) - wire(i_p_1_E, j_p_0_E).*Hxi(i_p_1_E, j_p_0_E) )/dy   );

    Hxo(i_p_0_H, j_p_0_H) = Hxi(i_p_0_H, j_p_0_H) + c3(i_p_0_H, j_p_0_H).*(   ( wire(i_p_0_H, j_p_0_H).*Ezo(i_p_0_H, j_p_0_H) - wire(i_p_0_H, j_p_1_H).*Ezo(i_p_0_H, j_p_1_H) )/dy   );
    Hyo(i_p_0_H, j_p_0_H) = Hyi(i_p_0_H, j_p_0_H) + c4(i_p_0_H, j_p_0_H).*(   ( wire(i_p_1_H, j_p_0_H).*Ezo(i_p_1_H, j_p_0_H) - wire(i_p_0_H, j_p_0_H).*Ezo(i_p_0_H, j_p_0_H) )/dx   );
end

Having fun with it

A simple pulse of light

Single lens simulation

Another interesting thing that you can learn with this kind of simple example is that there is no such thing as “collimated” light! In another words, lasers are not perfect cylinders of light. They diverge! The difference is that the laser diverges as little as possible keeping the solutions for the wave propagation valid.

Actually, the subject of wave optics is a discipline in on itself. One cool thing that you can learn is that, at focal distance of the lens, instead of a “point” of light, we actually see the Fourier transform of whatever passes through the aperture as plane a wave!!!😱 For this example, since we have more or less a “plain” wave and the aperture is like a “square” (at least in 2D), what we see in the focal plane is the Fourier transform of the square pulse (a sinc function). If you look carefully at the focal distance of the lens, you can notice the first “side lobes” of the sinc forming along the y direction! This might become a new post entitles some thing like: “How to compute Fourier transforms with a piece of glass” %-).

Ha, I can’t forget to mention the “power density loss” that happens with electromagnetic waves. You can clearly see that, as the wave propagates from the source, the amplitude decays (the waves become fainter in the vídeo). Thats because the energy of the wave have to be distributed in a larger area. So, as the wave becomes “bigger” when it leaves the source, its amplitude have to decrease. Its like each “ring” of electric field has to have the same “volume”. So, small rings can be taller than the larger ones.  

FSS

So, in a sense, this vídeo is nothing less than the simulation of a 2D microwave oven 😇! In the simulation you first see a pulse of high frequency content. It passes through the FSS (dashed line in the middle), reflecting some part. Then a second pulse of much higher wavelength is shot. This time, you barely see any wave transmitted, as this configuration is a high-pass FSS filter.  

Antireflex coating and refraction

The second concept present in this simulation is the concept of anti-reflex coating. This is a process that optical engineers have to do when fabricating lens, or even sheets of glass, to avoid reflection on the surface when light hits it. In the simulation, the upper “glass” has no treatment and the lower does. Because of that, you can clearly see that, when the wave hits the object whiteout the coating, there is a reflected wave going up, while in the surface with the coating, no wave is reflected.  

Optical fiber

The video is practically self explanatory %-). We have two fibers. The cool thing about this is that, due to the 2D nature of the simulation, this fiber is not an actual “tube”. Rather, it is an infinite retangular region. But, so is our wave! Hence, the net effect is that we have the same behavior of internal reflection as an real optical fiber! Moreover, depending on the wavelength of our source, the propagation won’t be like a “reflecting” ray. Instead, if we have a larger wavelength source compared to the fiber diameter, the light will propagate like blobs of light. You can see that in the first fiber. There, the light seems to propagate in a straight path. That is the famous “mono-mode fiber”! The multimode fiber is the second one, where you have the reflection very well identifiable.  

Crystal photonics wave guide

Wavefront sensor

Optical circuit

Going a bit further with the fun

Even further…

Conclusions

As I said in the beginning of this post, this is one of my favorite “weekend project”. Not only because I always wanted to simulate Maxwell’s equations themselves. Actually, after I started to play with it, I learned a lot of what is considered very complicated concepts without even notice that I was learning it! One situation that I will aways remember is that, after I arrive in my second post-doc in the Genevra Observatory, I was very scary that I would not understand the optics part of the project. But for my surprise, as I started hearing terms like “diffraction limit”, “point spread function”, “anti-reflex coating”, “sub-aperture” and etc, I began to associate those terms to stuff that I actually saw in my situations, hence knowing exactly what they were!!! That is a feeling that I, as a professor, can’t have enough of.

Thank you for reading this post. I hope you liked reading it as much as I liked writing it! As aways, feel free to comment or share this post as much as you want!

References

[2] – CEM Lectures

[3] – John B. Schneider, Understanding the Finite-Difference Time-Domain Method

[4] – FD-TF GitHub

[5] – Electro-optical modulation in two-dimensional photonic crystal linear waveguides

[6] – Waveguide Bends in Photonic Crystals

Intro

The fun…

In the case of the meme in this post, the 4th equation is nothing like a linear combination of the values. In fact, it is some weird looking complex mathematical aberration! This is where the meme ends for most people. The funny in that meme is that nobody can solve it. So, you suppose to send that to your annoying friend that keeps sending regular linear combination meme to you, solving all of them to show that he is a math genius.

Well… It turns out that the 4th equation is solvable! Moreover, it involves the integral of a very famous function (\(sinc(x)\)) which has a very elegant solution for when the limits involves infinity and zero. And yet more, that is despite the fact that it have no closed form solution for arbitrary limits.

Translating the emojis to variables, we arrive at the following integral $$\int\limits_{2b – a}^\infty {\frac{{a\sin (x)}}{{c\,x}}{\text{d}}x}$$ which, with numbers (solving the trivial linear system of emojis), becomes $$ 5\int\limits_0^\infty {\frac{{\sin (x)}}{x}{\text{d}}x} $$ The goal of this post them, is to solve this integral using almost all steps as purely algebraic manipulations. We are also going to use rules of calculus too, obviously. Moreover, we will have to go though some complex numbers too (yes, we go that deep in a innocent internet meme). Lets now solve the sinc function integral.

We start by acknowledging that \(sin(x)/x\) is a function with an anti-derivative that is very hard to find. If only we had something to relate easier functions to integrate (like exponentials)… But wait, we have! The first exponential that we can exchange, is the sine function (as you probably already guessed). So lets use Euler’s identity: $$ \sin (s) = \frac{{{e^{sj}} – {e^{ – sj}}}}{{2j}} $$ Although this is not a pure algebraic step, I will save the Taylor expansion prove of this one, for it is a very, very known relation.

The second one is the “jump of the cat” as we say in my country. I found that in an awesome StackExchange answer[1]. It has to do with the Laplace transform of the step function $u(t)$. In short, it states that the Laplace of the step function is the inverse function. The relation has a simple one line of proof as $$ L\left\{ {u(t)} \right\} = \int\limits_0^\infty {{e^{ – st}}{\text{d}}t} = \left. {\frac{{{e^{ – st}} – {e^{ – st}}}}{{ – s}}} \right|_0^\infty = \frac{1}{s} $$ Hence $$ \int\limits_0^\infty {{e^{ – st}}{\text{d}}t} = \frac{1}{s} $$ Holly butter ball! Now we can relate the inverse function with the exponential 😱. I wish I could use emojis on technical papers. They convey so much of the authors wonder when some beautiful math trick is used… %-)

Now lets change the names of the variables to $$ \int\limits_0^\infty {\frac{{\sin (s)}}{s}{\text{d}}s} $$ which makes our original sinc integral to become $$ \int\limits_0^\infty {\int\limits_0^\infty {{e^{ – st}}\frac{{{e^{sj}} – {e^{ – sj}}}}{{2j}}{\text{d}}t} {\text{d}}s} $$ Ok, the \(s\) integral is easy and here are the steps to solve it (it’s an exponential integral) $$ \begin{gathered} \frac{1}{{2j}}\int\limits_0^\infty {\int\limits_0^\infty {{e^{ – s(t – j)}} – {e^{ – s(t + j)}}{\text{d}}s} {\text{d}}t} \\ \frac{1}{{2j}}\int\limits_0^\infty {\left. { – \frac{{{e^{ – s(t – j)}}}}{{t – j}} + \frac{{{e^{ – s(t + j)}}}}{{t + j}}} \right|_0^\infty {\text{d}}t} \\ \frac{1}{{2j}}\int\limits_0^\infty { – \frac{{{e^{ – \infty (t – j)}}}}{{t – j}} – \frac{{{e^{ – \infty (t + j)}}}}{{t + j}} + \frac{{{e^{ – 0(t – j)}}}}{{t – j}} – \frac{{{e^{ – 0(t + j)}}}}{{t + j}}{\text{d}}t} \\ \frac{1}{{2j}}\int\limits_0^\infty {\frac{1}{{t – j}} – \frac{1}{{t + j}}{\text{d}}t} \end{gathered} $$ Now we have a integral in \(t\). This one is a lot of fun. Lets start by separating the two integrals as $$ \begin{gathered} \frac{1}{{2j}}\int\limits_0^\infty {\frac{1}{{t – j}}{\text{d}}t} – \frac{1}{{2j}}\int\limits_0^\infty {\frac{1}{{t + j}}{\text{d}}t} \end{gathered} $$ and do the following substitution. $$ \begin{gathered} {u_1} = t – j \\ d{u_1} = dt \\ {u_2} = t + j \\ d{u_2} = dt \\ \end{gathered} $$ Now we foresee that we are going to have some trouble with this infinity and the logs that will appear… Hum… Lets be careful (to not say formal) and wrap that around a limit. We deal with the limit later. Ok, so we have the following inverse integral (with is also simple to solve inside the limit) $$ \begin{gathered} \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \int\limits_{ – j}^{a – j} {\frac{1}{{{u_1}}}{\text{d}}{u_1}} – \int\limits_j^{a + j} {\frac{1}{{{u_2}}}{\text{d}}{u_2}} \\ \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \left. {\ln ({u_1})} \right|_{ – j}^{a – j} – \left. {\ln ({u_1})} \right|_j^{a + j} \\ \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \ln (a – j) – \ln ( – j) – \ln (a + j) + \ln (j) \\ \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \ln (a – j) – \ln (a + j) + \ln (j) – \ln ( – j) \\ \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \ln (a – j) – \ln (a + j) + \ln \left( {\frac{j}{{ – j}}} \right) \\ \frac{{\ln \left( { – 1} \right)}}{{2j}} + \frac{1}{{2j}}\mathop {\lim }\limits_{a \to \infty } \ln (a – j) – \ln (a + j) \end{gathered} $$ The lets now investigate this peculiar limit (and begin the fun part). We start by acknowledging that it is a complex-number limit, hence not straight forward to evaluate. Then, we do the second “jump of the cat” (small when compared with the Laplace one before) and change the limit of infinity with a limit to zero as $$ \begin{gathered} \mathop {\lim }\limits_{a \to \infty } \ln (a – j) – \ln (a + j) \hfill \\ \mathop {\lim }\limits_{a \to 0} \ln \left( {\frac{1}{a} – j} \right) – \ln \left( {\frac{1}{a} + j} \right) \end{gathered} $$ Now we just do some algebraic steps as follows $$ \begin{gathered} \mathop {\lim }\limits_{a \to 0} \ln \left( {\frac{{1 – aj}}{a}} \right) – \ln \left( {\frac{{1 + aj}}{a}} \right) \hfill \\ \mathop {\lim }\limits_{a \to 0} \ln \left( {1 – aj} \right) – \ln (a) – \ln \left( {1 + aj} \right) + \ln (a) \hfill \\ \mathop {\lim }\limits_{a \to 0} \ln \left( {1 – aj} \right) – \ln \left( {1 + aj} \right) \end{gathered} $$ And voilá! Put \(a=0\) and we have out first part of the integral $$ \ln \left( 1 \right) – \ln \left( 1 \right) = 0 $$ Now we have to figure what the heck is this negative logarithm thing… $$ \begin{equation} \frac{{\ln \left( { – 1} \right)}}{{2j}} \end{equation} $$ To do that, lets make the cat jump again and write -1 as \(e^{j\,\pi}\)

now we have $$ \begin{equation} \frac{{\ln \left( { – 1} \right)}}{{2j}} = \frac{{\ln \left( {{e^{j\pi }}} \right)}}{{2j}} \end{equation} $$ which is just \(\frac{\pi }{2} \). %-)

Finally coming back to the original integral we have $$ 5\int\limits_0^\infty {\frac{{\sin (x)}}{x}{\text{d}}x} = \frac{5\,\pi}{2}$$ Which is the actual answer for the meme %-).

Conclusion

I hope you like this post! Feel free to share it or use any of this derivations on your own proofs! See you in the next post!

References: