The solution of the Helmholtz equation

As we have seen in the previous lecture, the Helmholtz equation describes the propagation of the electric field in a anisotropic medium. Partially this equation can describe the propagation of the field in free space between two planes of the optical system. The solution of the Helmholtz equation in free space can be found using the plane wave decomposition, which leads to the Rayleigh-Sommerfeld diffraction integral for scalar fields. This integral allows us to calculate the field at any point in space given the field distribution on a plane, which is essential for understanding how light propagates through optical systems and how it diffracts around obstacles.

Green’s function

The Green’s function allows us to solve the differential equation with a source term:

where $L$ is a linear differential operator, $u(\boldsymbol{r})$ is the unknown function we want to find, and $f(\boldsymbol{r})$ is a known source term. The Green’s function, denoted as $G(\boldsymbol{r}, \boldsymbol{r}')$, is defined as the solution to the equation

where $\delta(\boldsymbol{r} - \boldsymbol{r}')$ is the Dirac delta function, which represents a point source located at $\boldsymbol{r}'$. The Green’s function can be used to express the solution to the original equation as an integral:

Physically, the Green’s function represents the response of the system to a point source, and the integral sums up the contributions from all sources in the volume $V$ to find the total solution at point $\boldsymbol{r}$.

Weyl representation of spherical waves

The Weyl representation of spherical waves is a mathematical technique used to express spherical wave solutions of the Helmholtz equation in terms of plane waves. This representation is particularly useful for explaining the propagation of light in free space using the Raleygh-Sommerfeld diffraction integral, which will be discussed further.

The spherical wave is given by the expression:

where $k$ is the wavenumber and $r=\sqrt{x^2 + y^2 + z^2}$ is the distance from the source. The spherical wave is the solution of the Helmholtz equation in free space with a point source, described by the Dirac delta function. It’s possible to decompose the spherical wave into the superposition of plane waves:

where

is the Fourier transform of the spherical wave.

The wave equation for the point source can be written as

Substituting the plane wave decomposition into this equation, we can derive the expression for the 3D-Fourier transform of the spherical wave. Implementing the Laplace operator in the Fourier domain, we get

where $|\boldsymbol{k}|^2 = k_x^2 + k_y^2 + k_z^2$ is the square of the magnitude of the wavevector.

Substituting this expression into the wave equation, we obtain

Further, using the representation of the Dirac delta function as an integral of plane waves, we can express the right-hand side as

and substitute it into the equation to get

Equating the integrands on both sides of the equation, we find the expression for the Fourier transform of the spherical wave:

In order to obtain the 2D Fourier transform of the spherical wave, we can integrate the 3D Fourier transform over the $k_z$ component:

Assuming that $z>0$ for waves propagating in the positive $z$ direction and denoting $\alpha = \sqrt{k^2 - k_x^2 - k_y^2}$, we can rewrite the integral as:

We can observe that the integrand has two simple poles at $k_z = \pm \alpha$. To evaluate the integral, we use the regularization technique by introducing a small imaginary part to the poles. The sign of the imaginary part is determined by the condition that the wave should decay as $z \to \infty$(Rayleigh-Sommerfeld condition). Therefore, we can write $\alpha$ as $\alpha + i\varepsilon$ with $\varepsilon > 0$. Substituting this into the integral, we get the poles

As $z>0$ and $\varepsilon>0$, the pole at $k_z = \alpha + i\varepsilon$ is located in the upper half of the complex plane, while the pole at $k_z = -\alpha - i\varepsilon$ is located in the lower half of the complex plane. To evaluate the integral, we can close the contour in the upper half-plane and apply the Cauchy residue theorem. Thus, the integral can be evaluated as

Substituting back $\alpha = \sqrt{k^2 - k_x^2 - k_y^2}$, we obtain the Weyl representation of the spherical wave:

or

where $\boldsymbol{k_{\perp}} = (k_x, k_y)$ and $\boldsymbol{r_{\perp}} = (x, y)$ are the transverse components of the wavevector and position vector, respectively.

Raleygh-Sommerfeld diffraction integral

Raleygh-Sommerfeld diffraction integral is a mathematical expression that describes the propagation of light in free space, particularly how it diffracts around obstacles. It is the strict solution of the Helmholtz equation in free space and can be derived using the plane wave decomposition. Thus, the general scalar equation describing the propagation of the field in free space can be written as

where $\boldsymbol{r} = (x, y, z>0)$. The 2D-plane wave decomposition of the field is given by

where $\boldsymbol{k_{\perp}} = (k_x, k_y)$ and $\boldsymbol{r_{\perp}} = (x, y)$ are the transverse components of the wavevector and position vector, respectively.

Substituting this equation into the Helmholtz equation, we get

Denoting $\alpha = \sqrt{k^2 - |k_\perp|^2}$, we can write the solution of this equation in spatial-frequency domain as

The choice of functions $A(\boldsymbol{k_{\perp}})$ and $B(\boldsymbol{k_{\perp}})$ is determined by the boundary conditions of the problem.

Firstly, the Fourier-transform of the field at $z=0$ gives us the relation

Secondly, the Rayleigh-Sommerfeld condition states that the wave should decay as $z \to \infty$, which implies that $B(\boldsymbol{k_{\perp}}) = 0$. Therefore, we can express the solution of the Helmholtz equation in spatial-frequency domain as

The important note is that the equation under the square root in the expression for $\alpha$ can be negative, which corresponds to evanescent waves that decay exponentially with distance. It is preferable to include these evanescent waves in the solution to ensure the completeness of the solution, especially when dealing with near-field effects or when the field distribution at $z=0$ contains high spatial frequencies.

Further, it’s necessary to apply the inverse Fourier transform to obtain the solution in spatial domain:

Notice that

Thus,

Finally, the solution of the Helmholtz equation in the spatial-frequency domain can be expressed as the Raleygh-Sommerfeld diffraction integral by using the convolution theorem:

where $u(\boldsymbol{r_0})$ is the field distribution at $z=0$ and the convolution is performed in the transverse coordinates, $R=\sqrt{(x-x')^2 + (y-y')^2 + z^2}$. This integral describes how the field propagates from the plane at $z=0$ to any point in space $z>0$, taking into account both the amplitude and phase changes due to diffraction. The physical interpretation of this integral is that the field at any point in space can be considered as a superposition of contributions from all points on the plane at $z=0$, where each contribution is weighted by the modified spherical wave

Longitudinal component of the field

In the previous sections, we have focused on the transverse components of the field, which are typically dominant in many optical systems. However, the longitudinal component of the field can also play a significant role, especially in tightly focused beams or near-field optics. The longitudinal component arises along the optical axis of the system(usually denotes $z$) due to the vector nature of the electromagnetic field and can be derived from the transverse components using divergence-free condition of the electric field in free space:

The Raleygh-Sommerfeld diffraction integral can be expressed in terms of a partial derivative with respect to $z$ as follows:

where $E_{\alpha}(\boldsymbol{r'})=E_{\alpha}(x',y',z=0)$ is the field distribution at $z=0$ and $R=\sqrt{(x-x')^2 + (y-y')^2 + z^2}$ is the distance from the source point $\boldsymbol{r'}$ to the observation point $\boldsymbol{r}$, $\alpha=x,y$.

Thus,

Taking into account that

we can express the longitudinal component of the field as

$E_z$-component of the field usually ignored in many problems due to its small magnitude compared to the transverse components, i.e $\frac{x-x'}{R}$ and $\frac{y-y'}{R}$ are much smaller than $\frac{z}{R}$. This regime is called the par-axial approximation.