Electromagnetic Waves
In the previous parts we have dealt with mechanics essentially. Even if we have described Brownian motion, this has been done by a particular type of Newton’s equation of motion, it is much like mechanics. Now we would like to have a look at some examples from electromagnetic waves. We will not solve the wave equation but look at some solution using the complex notion of the electric field. This shall train our use of complex numbers. The special solutions are the plane wave and the spherical wave and we will be able to simulate a number of things especially with the spherical waves as they are part of Huygens principle.
Plane waves
We will start with plane waves. Plane waves are solutions of the homogeneous wave equation and are the simplest solutions of the wave equation. They are also the basis for the description of more complicated waves. We will have a look at the electric field of a plane wave and its propagation in space and time.
Equations
A plane wave is a solution of the homogeneous wave equation and is given in its complex form by
\[\begin{equation} E=E_{0}e^{i\vec{k}\cdot \vec{r}}e^{-i\omega t} \end{equation}\]
where the two exponentials contain a spatial and a temporal phase. \(E_{0}\) denotes the amplitude of the plane wave. The plane is defined by the shape of the wavefront which is given by \(\vec{k}\cdot \vec{r}=const\), which is just the definition of a plane perpendicular to \(\vec{k}\).
A wave is a physical quantity which oscillates in space and time. Its energy current density is related to the square magnitude of the amplitude. We will include in the following the spatial and the temporal phase. For plotting just the spatial variation of the electric field, you may just use the spatial part of the equation
\[\begin{equation} E=E_{0}e^{i\vec{k}\cdot \vec{r}} \end{equation}\]
But since we also want to see the wave propagate, we will directly include also the temporal dependence in our function. In all of the examples below we set the amplitude of the wave \(E_{0}=1\).
The propagation of the wave is defined by wavevector \(\vec{k}\). In vacuum, the wavevector is just real valued
\[\begin{equation} \vec{k}_{0}= \begin{pmatrix} k_{0x} \\ k_{0y}\\ k_{0z}\\ \end{pmatrix} \end{equation}\]
The wavevector is providing the direction in which the wavefronts propagate. It is also proportional to the momentum of the wave, which will be important if we consider the refraction process a bit later. The magnitude of the wavevector is related to the wavelength \(\lambda\).
\[\begin{equation} k_{0}=\frac{2\pi}{\lambda_{0}}=\frac{\omega}{c_{0}} \end{equation}\]
At the same time, its magnitude is also given by the angular frequency divided by the speed of light. The latter is called a dispersion relation.
In a medium, the wavevector is by a factor of \(n\) longer, where n is the refractive index. Since the refractive index may be a complex number, e.g. \(n=\eta+i\kappa\), the wavevector can be complex as well. It is then given by
\[\begin{equation} \vec{k}=n\vec{k}_{0}= \begin{pmatrix} k_{x}^{\prime}+ik_{x}^{\prime\prime} \\ k_{y}^{\prime}+ik_{y}^{\prime\prime} \\ k_{z}^{\prime}+ik_{z}^{\prime\prime} \\ \end{pmatrix} \end{equation}\]
The complex refractive index means that there is some damping of the electromagnetic wave due to absorption, for example.
The wavelength is then related to
\[\begin{equation} \Re(k)=\eta \frac{2\pi}{\lambda_{0}} \end{equation}\]
and the imaginary part gives the damping
\[\begin{equation} \Im(k)=\kappa \frac{2\pi}{\lambda_{0}} \end{equation}\]
Electric field
Let’s have a look at waves and wave propagation. We want to create a wave, which has a wavelength of 532 nm in vacuum.
It shall propagate along the z-direction and we will have a look at the x-z plane.
We can plot the electric field in the x-z plane by defining a grid of points (x,z). This is done by the meshgrid function of numpy. The meshgrid returns a 2-dimensional array for each coordinate. Have a look at the values in the meshgrid.
In the last lines, we defined an array of X,0,Z, where X and Z are already 2-dimensional array. This finally gives an array 3D vectors, which we can use to calculate the electric field at any point in space. If we want to plot the electric field, we have to calculate the real part of the complex values, as the electric field is a physical quantity, which is always real. There is not much to see for a plane wave in the intensity plot, as the intensity of a plane wave is constant in space. Yet, if you want to plot it, you have to calculate the magnitude square of the electric field, e.g.
\[\begin{equation} I\propto |E|^{2} \end{equation}\]
Plane wave propagation
The above graph shows a static snapshot of the plane wave at a time \(t=0\). We know, however, that a plane wave is propagating in space and time. Since we know how to animate things, we may do that using the ipycanvas module.
x=np.linspace(-2.5e-6,2.5e-6,300)
z=np.linspace(0,5e-6,300)
X,Z=np.meshgrid(x,z)
r=np.array([X,0,Z],dtype=object)
canvas = Canvas(width=300, height=300,sync_image_data=False)
display(canvas)To do the animation I use a little trick to get the same color map as in the matplotlib plotting. The function below uses the matplotlib color map seismic and the corresponding mapping of values with a given minimum vmin and maximum vmax value. The mapping is done in the animation function with c=m.to_rgba(tmp).
#normalize the color map to a certain value range
norm = mpl.colors.Normalize(vmin=-1, vmax=1)
cmap = cm.seismic
# do the mapping of values to color values.
m = cm.ScalarMappable(norm=norm, cmap=cmap)This is our animation function, where I provide time and the wavevector as arguments, such that we may change both parameters easily.
def animate(k,time):
for t in time:
field=plane_wave(k,omega0,r,t)
tmp=np.real(field.transpose())
c=m.to_rgba(tmp)
with hold_canvas(canvas):
canvas.put_image_data(c[:,:,:3]*255,0,0)
#canvas.put_image_data(data*255,0,0)
sleep(0.02)With the call below, you may animate the wave now with different refractive indices.
eta=1.
kappa=0.0
n=eta+kappa*1j
k=n*k0*vec
time= np.linspace(0,5e-14,500)
animate(k,time)Imaginary wave vector
If we now create a material, which has an imaginary part of the refractive index, we see that the amplitude decays and the wave fades.
The above plots show the electric field amplitude in the x-z plane. We may also have a look at the field amplitude and intensity as a function of the z-position by choosing a single x-value. In the plot below, you may notice two things. The first is, that the wave decays exponentially with distance \(z\). Intensity and field decay with different decay lengths. The field decays with \(\exp(-\kappa k_{0}z)\) while the intensity of course decays twice as fast \(\exp(-2\kappa k_{0}z)\) due to the fact that the intensity is the square of the electric field.
Animation
Of course, we should not miss the animation.
display(canvas)k=n*k0*vec
time= np.linspace(0,5e-14,500)
animate(k,time)Interference of two plane waves
It is not very difficult to calculate from the definitions we did above now the interference of two plane waves, which have different directions of the wavevector. The total field in space is then just the sum of the two fields
\[\begin{equation} \vec{E}=\vec{E}_{1}+\vec{E}_{2} \end{equation}\]
The interesting thing is now to look at the intensity which
\[\begin{equation} I\propto |\vec{E}|^2=|\vec{E}_{1}|^2+|\vec{E}_{2}|^2 + \vec{E}_{1}^{*} \vec{E}_{2}+\vec{E}_{2}^{*}\vec{E}_{1} \end{equation}\]
While the field pattern still looks complicated, the intensity pattern is just a set of bright lines.
We want to go a bit further now and have a look at the wave at a boundary between vaccum and glass for example. At this boundary, the electromagnetic wave is reflected and refracted such that two new wavevectors arise. These are easily calculated by the law of reflection and the law of refraction. Besides that, also the amplitude of the waves change. To calculate the field we need the so-called Fresnel equations.
Fresnel equations
When electromagnetic waves hit a boundary, they will be reflected and refracted. The amplitude of the reflected and refracted wave is determined by the refractive index of the two materials, the angles and the polarizations. For the latter we differentiate between a polarization in the incident plane (the p-polarization) and perpendicular to the incident plane (s-polarization).
For each of the polarizations we in general obtain a coefficient for the reflection and one for the refraction. To make our calculation a bit simpler, we will assume only s-polarization. Then the two Fresnel coefficients are calculated as
\[\begin{equation} \left( \frac{E_{0t}}{E_{0e}} \right)_s = t_s =\frac{2n_1 \cos{\alpha}}{n_1\cos{\alpha}+n_2\cos{\beta}} \end{equation}\]
\[\begin{equation} \left( \frac{E_{0r}}{E_{0e}} \right)_s = r_s =\frac{n_1\cos{\alpha}-n_2\cos{\beta}}{n_1\cos{\alpha}+n_2\cos{\beta}} \end{equation}\]
where \(\alpha\) and \(\beta\) are the incident and refraction angles, respectively. Note that the Fresnel coefficients are for the amplitudes and can be negative to account for a phase jump by \(\pi\). To obtain the coefficients for the intensities, one has to square the Fresnel coefficients.
To bring everything correctly together, we therefore have to define a number of things. We will need a function calculating the outgoing angle from Snell’s law. And we need at least two functions calculating the reflection and transmission coefficient for one polarization. We use the s-polarization, where the electric field is always parallel to the interface.
With the definition of the Fresnel coefficients, we may now plot the reflection and the transmission coefficients. Note that the sum of reflection and transmission coefficients for the intensities have to add up to one if there is no absorption.
Incident wave
We want to study the electric fields and the intensities at various angles. The most interesting one is a case where we have total internal reflection. This happens if light is propagating from the higher refractive index to a lower refractive index. If we start in glass (\(n_1=1.5\)) and transmit to vacuum \(n_2=1\), then at all angles above \(\theta_{c}=\sin^{-1}(n_2/n_1)=41.810314895778596\) are total internally reflected.
We may now specify or calculate the corresponding wavevectors for an incident angle of \(45^{\circ}\). In general all waves (reflected, refracted) have to match with their phase at the boundary. If the boundary is along the x-direction, we therefore have
\[\begin{equation} k_{x,in}=k_{x,r}=k_{x,t} \end{equation}\]
This fixes one component of all wavevectors in the plane. What is then missing, is the z-component of the wavevectors. The incident wavevector is providing \(k_{z,in}\).
Reflected wave
For the reflected wave the z-component of the wavevector is just flipped in sign, e.g. \(k_{z,r}=-k_{z,in}\).
Refracted wave
The magnitude of the z-component of the transmitted wave can be obtained from the conservation of momentum. The momentum of the wave is proportional to the magnitude of the wavevector on both sides.
\[\begin{equation} k_{1}^2=k_{2}^{2} \end{equation}\]
which is, due to \(k=nk_{0}\) the same as
\[\begin{equation} n_{1}^2(k_{0x,in}^2+k_{0z,in}^{2})=n_2^2 (k_{0x,t}^{2}+k_{0z,t}^2) \end{equation}\]
from which we get
\[\begin{equation} k_{0z,t}=\pm \frac{1}{n_{1}}\sqrt{n_2^2 k_{0z,in}^2 -(n_{1}^2-n_{2}^2)k_{0x}^{2}} \end{equation}\]
If we go from a medium with high refrective index to a lower one, the second term in the root may surpass the first one and the whole solution will become imaginary. The wave in the lower refractive index medium \(n_{2}\) is then evanescent.
The total field thus contains three components. In medium 1, the field consists of the incident and the reflected wave. In medium 2, we just have the transmitted wave, with a possible evanescent solution.
The plots below show the electric field on the left side and the intensity on the right side. Interestingly, the intensity is that of a standing wave in medium 1, while it is just decaying in medium 2. Note that the electric field is oscillating along the interface in medium 2 but not at all in z-direction. This means that there is no energy transport along the z-direction anymore.
We will also have a look at the propagation of the wave by defining our animation.
canvas = Canvas(width=500, height=500,sync_image_data=True)
display(canvas)
def animate(k,time):
for t in time:
field=np.zeros([500,500],dtype=complex)
field1=plane_wave(k1,omega0,r1,t)
field2=plane_wave(k2,omega0,r1,t)
field3=plane_wave(k3,omega0,r2,t)
beta=snell(n1,n2,alpha)
r=rs(n1,n2,alpha,beta)
t=ts(n1,n2,alpha,beta)
field[0:250,:]=field1+r*field2
field[250:,:]=t*field3
tmp=np.real(field.transpose())
c=m.to_rgba(tmp)
with hold_canvas(canvas):
canvas.put_image_data(c[:,:,:3]*255,0,0)
sleep(0.02)
time= np.linspace(0,1e-14,100)
animate(k,time)As it is apparent from our simulation, the wave is longitudinal in medium 2 at this angle. Try to modify the incident angles yourself to see if the wave becomes propagating in medium 2.
In the last plot, we will have a look at the intensity in medium 1 and medium 2. What is nicely visible, is that the intensity decays in medium 2 with increasing distance. As compared to the absorbing case, there is not oscillation of the field in the z-direction, hence no energy transfer. Convince yourself that this is indeed an exponential decay by using the appropriate semilog plot.