Modeling the free space propagation of laser beams using Fourier transforms

I am trying to model the propagation of a laser beam in free space. I have an initial field $$E_{in}(x,z=0)$$ (a Gaussian beam) and need to find the fields at other points on the optical axis $$E(x,z=d)$$ for an arbitrary distance $$d$$.

By reading through a couple of texts, this is the approach that I have right now:

• Compute the Fourier transform of the initial field: $$\hat{E}(k_x) = \mathscr{F}[E_{in}(x,z=0)]$$
• Multiply $$\hat{E}$$ by the the free space transfer function $$e^{i k_z z0}$$ where $$k_z = \sqrt{k^2 - k_x^2}$$ to propagate it by a distance $$z0$$ along the optical axis.
• Inverse Fourier transform back to obtain $$E(x,z=z0)$$

This method makes sense to me. I think we are imagining the field as an infinite collection of plane waves and through the Fourier transform, we are essentially moving each of these plane waves by propagating through each of their respective wave numbers. I understand that the ABCD matrix method might be an easier technique, but I need a method that works for arbitrary beams and not just Gaussian beams.

I am implementing this on Mathematica at the moment and the resulting fields that I am getting do not match my expectations from Gaussian beam propagation (they do not follow the trends of spherical wave fronts). I would appreciate any help in figuring out if this is the right approach. I would also appreciate any help in finding other techniques that might be useful for this modeling.

• Have you tried a Fourier transform in time and a Laplace transform in z?
– webb
Aug 6, 2014 at 16:05

The problem with using the actual free space Fourier propagator is aliasing.

I learned this through trial and error as well, after a few wavelengths the numerical model really begins to behave poorly, probably due to aliasing and roundoff error.

The Fresnel approximation actually does a better job if you are anywhere further than the very near field region, it is numerically more stable...I'm sure you could write software that corrects the errors, but just use Fresnel, it is very accurate...

$$\newcommand{\Four}{\mbox{\mathcal{F}}} u_z(\mathbf{r}) \approx \frac{e^{ikz}}{i\lambda z}e^{i\pi\frac{\mathbf{r}^2}{\lambda z}}\Four\left\{ u_0(\mathbf{r}_0) e^{i\pi\mathbf{r}_0^2/\lambda z}\right\}_{\mathbf{\rho} = \frac{\mathbf{r}}{\lambda z}}$$

where the following constraint must be met : $$z^3 \gg \|\mathbf{r}-\mathbf{r}_0\|^4/\lambda$$ and $\mathbf{r}:=(x,y)$; the coordinates in the plane at your particular $z$, perpendicular to the $z$-axis, $\mathbf{r}_0:=(x_0,y_0)$; the coordinates in the plane at your initial field position $z_0 = 0$, perpendicular to the $z$-axis.

Fraunhofer (far-field) is valid when $$z \gg \|\mathbf{r}-\mathbf{r}_0\|^2/\lambda$$

• Thanks for the answer! I think I need to better understand my length scales here.I thought the Fresnel approximation was for the near-field. Why is it more stable at regions further from the near-field? The next steps of my problem involve diffraction by gratings. Would distances a few centimeters away from standard diffraction gratings be considered near field ? (i.e. Should I use Fresnel or Fraunhoffer diffraction?) Aug 11, 2014 at 16:19
• It is, but it looses accuracy in the very near field you could say...prior to 1 or 2 wavelengths from your initial field. Aug 11, 2014 at 18:46
• Also, if you write code for Fresnel, it will work in the far-field (Fraunhoffer) zone. I'll edit the above for the scales which are valid for each approximation. I believe that the Fresnel approximation is more stable numerically because some of the high frequency components of the actual free space transfer function are not well approximated when they are discretized. The Fresnel fixes this through a mathematical approximation...I didn't investigate fully however... Aug 11, 2014 at 18:52
• @daaxix Where did you get the constraints? Can you cite an article or book?
– DaP
May 24, 2017 at 15:41
• @DaP, Goodman and Gaskill both have it I believe, as does Barrett and Myers. Depending on your mood, and how you define the off-axis shape, the constants in these expressions can vary. I've also derived it myself directly from the Taylor series expansion used in the Fresnel approximation, my derivation gave me a multiple of $\pi$. May 25, 2017 at 0:01

I already wrote elsewhere that Gaussian beams are just some approximations to solutions of the Maxwell equations. For this reason, I derived some exact solutions of the Maxwell equations that are approximated by Gaussian beams extremely well ("asymptotically accurately") when the beam waist is much larger than the wavelength. Please see https://arxiv.org/abs/physics/0405091, Eq. 22. The solution describes a circularly polarized beam, but it is not difficult to derive a solution describing a linearly polarized beam.

As mentioned above, aliasing is an issue with Fresnel propagation using Fourier transforms, but it is very manageable. SPIE has published some good books on writing code for Fresnel propagation and how to minimize the effects of aliasing. These two cover the propagation math and have lots of Matlab code: https://doi.org/10.1117/3.858456

https://doi.org/10.1117/3.866274

This one is less about propagation and more about doing Fourier transforms in Mathematica: https://doi.org/10.1117/3.2574956

A lot of universities have institutional subscriptions to the SPIE Digital Library, so you might be able to download these books for free.