Paraxial Wave Equation

Table of Contents:

  1. Derivation
  2. General Solution
  3. Gaussian Solution
  4. Standard form of Gaussian Beam


In the last section, we started with a general solution (angular spectrum) to the Helmholtz equation:

\begin{equation} (\nabla^2+k^2)E(x,y,z) = 0\end{equation}

which we found specific solutions to by considering the propagation of a beam at small angles to the x-axis in the spatial frequency domain (Fresnel approximation). Another approach is to find a differential equation that approximates paraxial field propagation. In other words, as opposed to the last section when we found exact solutions to the Helmholtz equation using the angular spectrum that we then propagates through space using linear response theory, we are now making approximations to the Helmholtz equation itself by assuming paraxial propagation from the start such that we can rewrite the differential equation. We will then find solutions for this equation (in next part of this page, in fact!).

The field for a paraxial equation propagating primarily along the z-axis can be written as:

\begin{equation} E(x,y,z) = \varepsilon(x,y,z)e^{-ikz} \end{equation}

where we can note that what makes this “paraxial” is the fact that the phase terms only include propagation in the z-direction. $\varepsilon$ will be a slowly varying envelope function that modulates the carrier wave during propagation along the z-direction. Now, if we were to substitute this equation into the Helmholtz equation, we would first have via Chain’s Rule (hehe):

\begin{equation} \nabla^2E(x,y,z) = (\nabla^2\varepsilon – 2ik\hat{z}\cdot\nabla{\varepsilon}-k^2\varepsilon)e^{-ikz} \end{equation}

such that Helmholtz’s equation reads:

\begin{equation} \nabla^2\epsilon – 2ik\frac{\partial\varepsilon}{\partial{z}} = 0 \end{equation}

If we then consider that the envelope varies slowly, such that

\begin{equation} |\frac{\partial^2\varepsilon}{\partial{z}^2}| << 2k|\frac{\partial\epsilon}{\partial{z}}| \end{equation}

Then we have the paraxial wave equation:

\begin{equation} \frac{\partial^2\varepsilon}{\partial{x^2}}+\frac{\partial^2\varepsilon}{\partial{y^2}} – 2ik\frac{\partial\epsilon}{\partial{z}} = 0 \end{equation}

Here, we can see how the Fresnel and paraxial approximations are equivalent. Consider that the spherical wave $e^{-ikr}/r$ is an exact solution to the scalar Helmholtz equation. We can then write the radius “r” as:

\begin{equation} r = \sqrt{x^2+y^2+z^2} = z\sqrt{1+\frac{x^2+y^2}{z^2}} \end{equation}

so the field is then given in terms of “r” as:

\begin{equation} E(r) = E_0\frac{e^{-ikr}}{r} = \frac{E_0}{z\sqrt{1+\frac{x^2+y^2}{z^2}}}e^{-ikz\sqrt{1+\frac{x^2+y^2}{z^2}}} \end{equation}

Now here is the important part! If we make the Fresnel approximation such that $x^2+y^2 << z^2$, then:

\begin{equation} E(r) = \frac{E_0}{z}e^{-ikz(1+\frac{x^2+y^2}{2z^2})} \end{equation}

The result is a paraboloidal wave, given as:

\begin{equation} \varepsilon(x,y,z) = \frac{E_0}{z}e^{-ik(x^2+y^2)/2z} \end{equation}


Since the paraxial equation is just the Helmholtz equation with simplifying assumptions, we can use our basic solution to the Helmholtz, the angular spectrum, multiplied by a transfer function to find the field at an arbitrary distance z:

\begin{equation} E(x,y,z) = \iint^{\infty}_{-\infty}A(k_x,k_y;0)H(k_x,k_y;z)e^{-i(k_xx+k_yy)}dk_xdk_y \end{equation}

Here is where we make the situation specific for paraxial/Fresnel approximations. In these approximations, the transfer function of free space is:

\begin{equation} H(k_x,k_y;z) = e^{-ikz}e^{i(k_x^2+k_y^2)z/2k} \end{equation}

where $k = 2\pi/\lambda$ is the magnitude of the wavevector. Then, combining the last two equations:

\begin{equation} E(x,y,z) = e^{-ikz}\iint^{\infty}_{-\infty}A(k_x,k_y;0)e^{i(k_x^2+k_y^2)z/2k}e^{-i(k_xx+k_yy)}dk_xdk_y \end{equation}

so that obviously the field envelope is just the integral part:

\begin{equation} \varepsilon(x,y,z) = \iint^{\infty}_{-\infty}A(k_x,k_y;0)e^{i(k_x^2+k_y^2)z/2k}e^{-i(k_xx+k_yy)}dk_xdk_y \end{equation}

This is then the most general solution to the paraxial wave equation.


The parabolodal wave is an exact solution to the paraxial wave equation as we found above, so we can thus use this to find the Gaussian solution to the paraxial wave equation. The key mathematical insight is that the solution of a differential equation must be independent of origin. Thus, we can shift to the position:

\begin{equation} \varepsilon = \frac{E_0}{z-\zeta}e^{-ik(x^2+y^2)/2(z-\zeta)} \end{equation}

Since $\zeta$ is just a number, it can also be imaginary, so we can just try substituting $\zeta = -iz_R$ and find the envelope as:

\begin{equation} \varepsilon(x,y,z) = \frac{E_0}{q(z)}e^{-ik(x^2+y^2)/2q(z)} \end{equation}

Since we have $q(z)$ in the denominator of the exponent, we can break it into its real and imaginary parts as:

\begin{equation} \frac{1}{q} = \frac{1}{q_r} – i\frac{1}{q_i} \end{equation}

so that our envelope is:

\begin{equation} \varepsilon = \frac{E_0}{q(z)}e^{-k(x^2+y^2)/2q_i}e^{-ik(x^2+y^2)/2q_r} \end{equation}

We can choose to shift the origin of our paraboloidal equation, since the solution should be invariant upon translation. If we consider the shift to be $\eta = -iz_R$, then the envelope becomes:

\begin{equation} \varepsilon(x,y,z) = \frac{E_0}{q(z)}e^{-ik(x^2+y^2)/2q(z)} \end{equation}


We can define the origin as the position where the beam has its minimum beam radius, i.e.

\begin{equation} \omega(z=0) = \omega_0 \end{equation}

such that at this position, the wavefront radius R(z=0) = $\infty$. Then the q-parameter becomes:

\begin{equation} \frac{1}{q(z)} = \frac{1}{R(z)}-i\frac{\lambda}{\pi\omega^2(z)} \rightarrow 0 -i\frac{\lambda}{\pi\omega_0^2} = \frac{1}{q_0} \end{equation}

or rather

\begin{equation} q_0 = i\frac{\pi\omega_0^2}{\lambda} = iz_R \end{equation}

where as before we had the Rayleigh range defined as:

\begin{equation} z_R = \frac{\pi\omega_0^2}{\lambda} \end{equation}

Now we can write out our main three relations for a Gaussian beam:

\begin{equation} \omega(z) = \omega_0\sqrt{1+(\frac{z}{z_R})^2} \end{equation}

\begin{equation} R(z) = z[1+(\frac{z_R}{z})^2] \end{equation}

\begin{equation} \phi(z) = tan^{-1}(\frac{z}{z_R}) \end{equation}

Note that when $z=z_R$ then $\omega(z)=\sqrt{2}\omega_0$, so the beam radius has expanded by a factor of $\sqrt{2}$. This then gives us the physical intuition into what the Rayleigh range means: its a measure of how far the beam is approximately collimated (i.e. has flat wavefronts) before the divergence due to diffraction becomes significant.

A common nomenclature defines the confocal parameter of the beam as:

\begin{equation} b = 2z_R \end{equation}

Since $z_R = \pi\omega_0^2/\lambda$, a small beam waist corresponds to a short confocal parameter. As the beam propagates further into the far field, ie many Rayleigh ranges away, the beam expands essentially linearly with distance such that the shape approaches a cone with a divergence angle $\theta$. We can calculate the divergence angle by first finding the beam radius:

\begin{equation} \omega(z) = \omega_0\sqrt{1+(\frac{z}{z_R})^2} \approx \omega_0(\frac{z}{z_R}) \end{equation}

when z >> $z_R$. Then,

\begin{equation} \theta \approx tan\theta = \frac{\omega(z)}{z_R} = \frac{\omega_0}{\pi\omega_0^2/\lambda} = \frac{\lambda}{\pi\omega_0} \end{equation}

The angular spread of the Gaussian beam is then defined as:

\begin{equation} \theta = \frac{\lambda}{\pi\omega_0} \end{equation}

Thus, a small beam waist corresponds to a large diffraction angle, and vice versa.

Some other notes:

  1. For the wavefront radius, at z = $z_R$, the wavefront radius is at its minimum value. We can see this as the wavefront radius begins at infinity at the beam waist, then acquires curvature on diffraction, and then again looks like an infinite plane wave at infinity.
  2. The phase will asymptotically approach $\pi/2$ as z $\rightarrow \infty$.
  3. At z = $z_R$, the beam waist is $\sqrt{2}\omega_0$ and the beam diameter is $2\sqrt{2}\omega_0$.
%d bloggers like this: