Wave equation

We will need the Maxwell’s equations, which are shown below: \[\begin{eqnarray} \text{Gauss' Law for electric fields:}\qquad \mathbf{\nabla}\cdot\mathbf{E} &=& \frac{\rho}{\varepsilon_0}, \tag{1} \end{eqnarray}\] \[\begin{eqnarray} \text{Gauss' Law for magnetic fields:}\qquad \mathbf{\nabla}\cdot\mathbf{B} &=& 0 \tag{2}, \end{eqnarray}\] \[\begin{eqnarray} \text{Faraday's Law:}\qquad \mathbf{\nabla}\times\mathbf{E} &=& -\frac{\partial \mathbf{B}}{\partial t} \tag{3}, \end{eqnarray}\] \[\begin{eqnarray} \text{Ampere's Law:}\qquad \mathbf{\nabla}\times\mathbf{B} &=& \mu_0\mathbf{J}+\mu_0\varepsilon_0\frac{\partial\mathbf{E}}{\partial t} \tag{4}. \end{eqnarray}\]

Consider the Maxwell’s equations away from the sources ( \(\mathbf {J}=0\) and \(\rho=0\)). Lets take \(\mathbf{\nabla}\times\) of Eq. (3): \[\begin{eqnarray} \mathbf{\nabla}\times( \mathbf{\nabla}\times\mathbf{E}) &=& -\frac{\partial }{\partial t}\left( \mathbf{\nabla}\mathbf\times{B}\right). \tag{5} \end{eqnarray}\] We need the following vector identity: \[\begin{eqnarray} \left(\mathbf{\nabla}\times[ \mathbf{\nabla}\times\mathbf{E}] \right)^i=\epsilon^{ijk}\epsilon^{klm}\partial_j\partial_l E_m =\left(\delta^{il}\delta^{jm}-\delta^{im}\delta^{jl}\right)\partial_j\partial_l E_m=\mathbf{\nabla}^i\left(\mathbf{\nabla}\cdot \mathbf{E}\right)-\mathbf{\nabla}^2 E^i. \tag{6} \end{eqnarray}\] Using the Ampere’s law from Eq. (4) with \(\mathbf {J}=0\) , we get \[\begin{eqnarray} \left( \mu_0\varepsilon_0 \frac{\partial^2}{\partial t^2}-\mathbf{\nabla}^2 \right)\mathbf{E}=0. \tag{7} \end{eqnarray}\] This is a wave equation with speed \(c=\frac{1}{\sqrt{\mu_0\varepsilon_0 }}\), which is the speed of light. This shows that any change in the electric fields will travel at the speed of light. If we really wanted, we could do the identical derivation for magnetic fields.

We would like solve the wave equation: \[\begin{eqnarray} \left[\nabla ^2 -\partial^2_\tau\right]\mathbf{E}(\mathbf{r},t)=0 \tag{8}, \end{eqnarray}\] where \(\tau=c t\).

Figure 2: The spherical coordinates.

In the spherical coordinates, the Laplace operator reads: \[\begin{eqnarray} \nabla ^2 = \left[ \frac{1}{r^2}\frac{\partial }{\partial r}\left( {r^2\frac{\partial }{\partial r}} \right)+ \frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left( {\sin \theta \frac{\partial }{\partial \theta }} \right) +{ \frac{1}{r^2\sin ^2\theta }\frac{\partial ^2 }{\partial \phi ^2}} \right] \tag{9}. \end{eqnarray}\] However, the angle dependence drops out due to spherical symmetry. \[\begin{eqnarray} &&\left[\nabla ^2 -\partial^2_\tau\right]\mathbf{E}(\mathbf{r},t) =\left[ \frac{1}{r^2}\frac{\partial }{\partial r}\left( {r^2\frac{\partial }{\partial r}} \right) -\partial^2_\tau\right]\mathbf{E}(\mathbf{r},t)=\left[ \frac{2}{r}\frac{\partial }{\partial r}+\frac{\partial^2 }{\partial r^2} -\partial^2_\tau\right]\mathbf{E}(\mathbf{r},t) =\frac{\partial^2 }{\partial r^2}\left[r\mathbf{E}(\mathbf{r},t)\right] -\partial^2_\tau\left[r\mathbf{E}(\mathbf{r},t)\right]=0 \tag{10}, \end{eqnarray}\] where we used the Leibniz’s formula: \[\begin{eqnarray} \frac{d^m}{dx^m}\left[f(x)g(x)\right]=\sum_{k=0}^n {n \choose k} \frac{d^k f}{dx^k}\frac{d^{n-k}g}{dx^{n-k}} \tag{11}. \end{eqnarray}\] This is just like the one-dimensional wave equation, and we can find the solutions as incoming and outgoing waves: \[\begin{eqnarray} \left( \partial_r + \partial_\tau \right)\left( \partial_r - \partial_\tau \right) \left[r\mathbf{E}(\mathbf{r},t)\right] \tag{12}. \end{eqnarray}\]

The outgoing wave solution is \[\begin{eqnarray} \mathbf{E}(\mathbf{r},t)=\frac{\mathbf{\mathcal{E}}_0}{\lambda r}e^{-i(kr -\omega t+\delta)} \tag{13}, \end{eqnarray}\] where \(k\) is the wave-number, \(\lambda\) is the wavelength, \(\omega\) is the angular frequency, \(\mathbf{\mathcal{E}}_0\) is the polarization, and \(\delta\) is the initial phase. We defined the overall coefficient such that the expression actually represents the electric field surface density. The power radiated by the source is proportional to \(\mathbf{E}^2\) and it behaves as \(1/r^2\), which is the typical radius dependence.

Huygens-Fresnel principle

Figure 3: An illustration of the Huygens’ principle.

Let’s assume that there is an aperture located at \(z=0\) and a screen at \(z=z_0\). Per the Huygens–Fresnel principle, every differential area on the aperture can be considered as light source, and we just need to integrate over the aperture surface. We will assume that the distance to the screen is much larger than the aperture size. This enables us to be able to expand the radial distance from a point on the aperture \(\mathbf{r}'\) to a point on the screen \(\mathbf{r}\): \[\begin{eqnarray} |\mathbf{r}-\mathbf{r}'|&=&\sqrt{z_0^2+(x-x')^2+(y-y')^2}= z_0\sqrt{1+\frac{(x-x')^2}{z_0^2}+\frac{(y-y')^2}{z_0^2}}\simeq z+\frac{(x-x')^2}{2z_0}+\frac{(x-x')^2}{2z_0} \tag{14}. \end{eqnarray}\] Inserting this back in Eq. (13) yields: \[\begin{eqnarray} \mathbf{E}(\mathbf{r})\simeq \frac{ e^{-ikz_0}}{\lambda z_0} \iint dx'dy' \mathbf{\mathcal{E}}_0 (x',y') e^{-ik \frac{(x-x')^2}{2z_0}-ik\frac{(x-x')^2}{2z_0}} \tag{15}. \end{eqnarray}\] If the aperture is a rectangle, and the source is far away, i.e., \(\mathbf{\mathcal{E}}_0 (x',y')\simeq \text{constant}\) the integrals in Eq. (15) can be separated: \[\begin{eqnarray} \mathcal{I}&\equiv& \int_{a}^b dx' e^{-ik \frac{(x'-x)^2}{2z_0}} =\sqrt{\frac{z_0}{2 k}}\int_{\sqrt{\frac{k}{2 z_0}}(a-x)}^{\sqrt{\frac{k}{2 z_0}}(b-x)} du e^{-i u^2} =\sqrt{\frac{z_0}{2 k}}\left(\int_{0}^{\sqrt{\frac{k}{2 z_0}}(b-x)}e^{-i u^2} du- \int_0^{\sqrt{\frac{k}{2 z_0}}(a-x)} du e^{-i u^2}\right). \tag{16} \end{eqnarray}\] Now we need to study these interesting integrals.

Fresnel Integrals

The Fresnel integrals are defined as follows: \[\begin{eqnarray} S(t)&=& \int_0^t dr \sin r^2,\nonumber\\ C(t)&=& \int_0^t dr \cos r^2 \tag{17}. \end{eqnarray}\] For a general value of \(t\), the integrals need to be evaluated numerically. However, the asymptotic values \(C(t)\) and \(S(t)\) can be calculated via the closed contour integral below: \[\begin{eqnarray} I&=& \oint_C dz e^{- z^2} \tag{18}, \end{eqnarray}\] where the contour \(C\) is illustrated in Fig. 4.

Figure 4: The contour to evaluate the integral in Eq. (18). The return path, \(\gamma_1\), is chosen such that the integrand reduces to the regular Gaussian.

Let’s first evaluate the integral on \(\gamma_0\) in the limit \(R\to\infty\): \[\begin{eqnarray} I_{\gamma_0}&=& \lim_{R\to\infty} \int_0^R dr e^{- r^2}=\frac{\sqrt{\pi}}{2}, \tag{19} \end{eqnarray}\] where the details of the derivation can be found here. Now consider the (absolute value of the) integral on \(\gamma_R\) in the limit \(R\to\infty\): \[\begin{eqnarray} \left| I_{\gamma_R}\right|&=& \left| \lim_{R\to\infty} R \int_0^{\frac{\pi}{4}} d\theta e^{i\theta}e^{- R^2 (\cos^2\theta-\sin^2\theta+2i \cos\theta\sin\theta)}\right| = \left|\lim_{R\to\infty} R \int_0^{\frac{\pi}{4}} d\theta e^{i\theta+i\sin(2\theta)}e^{- R^2 \cos(2\theta)}\right|\nonumber\\ &\leq&\left|\lim_{R\to\infty} R \int_0^{\frac{\pi}{4}} d\theta e^{- R^2 \cos(2\theta)} \right|. \tag{20} \end{eqnarray}\] Let’s try to put a bound on \(\cos(2\theta)\) in the range \(0\leq\theta\leq\pi/4\). At \(\cos(2\theta)\vert_{\theta=0}=1\) and \(\cos(2\theta)\vert_{\theta=\pi/4}=0\). We can draw a line that connects these two points: \(1-\frac{4\theta}{\pi}\). Since \(\frac{d^2}{d\theta^2}\cos(2\theta)=-4 \cos(2\theta)<0\) for \(0<\theta<\pi/4\), we know that \(\cos(2\theta)<1-\frac{4\theta}{\pi}\) in this range. This observation is illustrated in Fig. 5.

Figure 5: \(\cos(2\theta)\) and a bound on it with the line \(1-\frac{4\theta}{\pi}\).

We can now go back to Eq. (22) make use of the bound:

\[\begin{eqnarray} \left| I_{\gamma_R}\right| &\leq&\left|\lim_{R\to\infty} R \int_0^{\frac{\pi}{4}} d\theta e^{- R^2 \cos(2\theta)} \right| \leq\left|\lim_{R\to\infty} R \int_0^{\frac{\pi}{4}} d\theta e^{- R^2 \left(1-\frac{4\theta}{\pi}\right)} \right| =\left|\lim_{R\to\infty} R e^{- R^2} \int_0^{\frac{\pi}{4}} d\theta e^{R^2\frac{4\theta}{\pi}} \right|\nonumber\\ &\leq& \left|\lim_{R\to\infty} R \frac{\pi}{4 R^2} \left(1- e^{- R^2}\right) \right| =\left|\lim_{R\to\infty} \frac{\pi}{4 R} \left(1- e^{- R^2}\right) \right|=0. \tag{21} \end{eqnarray}\]

Finally, let’s look at the integral on \(\gamma_1\) in the limit \(R\to\infty\): \[\begin{eqnarray} I_{\gamma_1}&=& \lim_{R\to\infty} R \int_R^{0} dr e^{\frac{i\pi}{4}}e^{- r^2 \frac{i\pi}{2}} =-\frac{1+i}{\sqrt{2}}\lim_{R\to\infty} \int_0^R dr \left( \cos r^2 -i\sin r^2\right) =-\frac{1+i}{\sqrt{2}} \int_0^\infty dr \left( \cos r^2 -i\sin r^2\right)\nonumber\\ &=&-\frac{1}{\sqrt{2}}\left( \int_0^\infty dr \cos r^2 +\int_0^\infty dr \sin r^2 +i\left[\int_0^\infty dr \cos r^2 -\int_0^\infty dr \sin r^2 \right] \right). \tag{22} \end{eqnarray}\] As we have computed individual pieces of the integral Eq.(18), we can assemble them and state that they need to add to \(0\) since \(e^{-z^2}\) is analytic everywhere. Therefore we have: \[\begin{eqnarray} \oint_C dz e^{- z^2}=0&=&I_{\gamma_0}+I_{\gamma_R}+I_{\gamma_1}\nonumber\\ &=&\frac{\sqrt{\pi}}{2}+0-\frac{1}{\sqrt{2}}\left( \int_0^\infty dr \cos r^2 +\int_0^\infty dr \sin r^2 +i\left[\int_0^\infty dr \cos r^2 -\int_0^\infty dr \sin r^2 \right] \right) \tag{23}. \end{eqnarray}\] Matching the real and imaginary parts, we get:

\[\begin{eqnarray} \frac{1}{\sqrt{2}}\left( \int_0^\infty dr \cos r^2 +\int_0^\infty dr \sin r^2 \right)&=&\frac{\sqrt{\pi}}{2},\nonumber\\ \frac{1}{\sqrt{2}}\left( \int_0^\infty dr \cos r^2 -\int_0^\infty dr \sin r^2 \right)&=&0, \tag{24} \end{eqnarray}\] from which we get \[\begin{eqnarray} \int_0^\infty dr \cos r^2 =\int_0^\infty dr \sin r^2 &=&\frac{1}{2}\sqrt{\frac{\pi}{2}}\simeq 0.626 \tag{25}. \end{eqnarray}\]

Now we can conclude with the plots of the Fresnel integrals in Fig. 6.

Figure 6: Left: Fresnel integrals as a function of their argument, Right: parametric plot of the integrals forming the Euler spiral.

Going back to Eq. (16), we see that we can express it in terms of the Fresnel integrals: \[\begin{eqnarray} \mathcal{I}&\equiv& \int_{a}^b dx' e^{-ik \frac{(x'-x)^2}{2z_0}} =\sqrt{\frac{z_0}{2 k}}\left(\int_{0}^{\sqrt{\frac{k}{2 z_0}}(b-x)}e^{-i u^2} du- \int_0^{\sqrt{\frac{k}{2 z_0}}(a-x)} du e^{-i u^2}\right)\nonumber\\ &=& \sqrt{\frac{z_0}{2 k}}\left\{ C(u) -C(v)+i[S(v)-S(u)] \right\}, \tag{26} \end{eqnarray}\] where \(u=\sqrt{\frac{k}{2 z_0}}(b-x)\) and \(v=\sqrt{\frac{k}{2 z_0}}(a-x)\). Remember that the aperture is between \(a\) and \(b\) with a span \(L=b-a\). We can shift the coordinate system so that the aperture is centered at \(0\) and it extends from \(a=-L/2\) to \(b=L/2\). Plugging this back in to Eq. (26) we get \[\begin{eqnarray} \mathcal{I}&=& \sqrt{\frac{z_0}{2 k}}\left\{ C(x^+) -C(x^-)+i [S(x^-)-S(x^+)] \right\} \equiv \sqrt{\frac{z_0}{2 k}} \mathcal{F}(x^-,x^+) \tag{27} \end{eqnarray}\] where \(x_\pm=x\pm L/2\) and \(\mathcal{F}(x^-,x^+)\equiv C(x^+) -C(x^-)+i [S(x^-)-S(x^+)]\). This is the \(x\) integral. We have an identical copy for \(y\) axis. Putting them together along with the overall coefficient in Eq. (15) we get: \[\begin{eqnarray} \mathbf{E}(\mathbf{r})&\simeq&\frac{\mathbf{\mathcal{E}}_0 e^{-ikz_0}}{2} \mathcal{F}(x^-,x^+)\mathcal{F}(y^-,y^+), \tag{28} \end{eqnarray}\] where \(y_\pm=y\pm L_y/2\).

Interactive Plots

We can now take E. (28) and plot the intensity in Fig. 7.

Figure 7: The plot of intensity for various input parameters.