Propagation of a flat-top laser beam in free space I have spent some time looking at how to design a field mapping system that turns a Gaussian beam profile into a top-hat beam. However, it would also be very useful to have a very sharp super-Gaussian profile for a different experiment. 
I have found a page (the RP Photonics article on flat-top beams) which argues that since the top-hat is not a free space mode of the Helmholtz equation, the profile changes as it propagates. From the simulated model (Figure 1) and the description given by the author it looks to me like the top-hat might initially collapse to something resembling a super-Gaussian but its a bit hand-wavy to deduce that from an intensity plot.

My question now is: Does anyone have a suggestion as to how to replicate this simulation or at least find a way to solve for the specific intensity profile I am looking for. 
Which equations do I use and how can I apply them for a computational solution (I'm guessing the solutions will have to be numerically evaluated)?
Or alternatively, is my reasoning just plain wrong and the profile does not contract into a super Gaussian as the beam propagates? 
 A: The initial flat-top function is often assumed to be a super-Gaussian function
$$ \exp\left(-\frac{r^n}{w^n}\right) . $$
Unfortunately, super-Gaussians are a bit unfriendly to use in analytic calculations. So it is probably better to due it numerically. The simplies why to do this is to use, what is called beam-propagation:


*

*Take the 2D Fourier transform of the function (numerically: 2D FFT)

*Multiply with the propagation kernel:
$ \exp\left(i\pi\lambda z r^2\right) $ where $r$ is the radial coordinate in the Fourier domain, $z$ is the propagation distance and $\lambda$ is the wavelength.

*Do the inverse Fourier transform to get the field at the new $z$-location.
There are a few subtleties with such numerical implementations, but I trust you'll figure them out.
A: Just like any beam profile can be expressed as a sum of plane-waves, it can also be decomposed in Hermite functions. 
To do so, consider a beam propagating in the z-direction. The electric field can be expressed as the product of an envelop $u$ and a carrier $e^{ikz}$ :
$$
E_{m,n}(x,y,z) = A\times u_{m,n}(x,y,z) \exp\left( ikz \right) 
$$
and the propagation of the envelop can be expressed as
$$
u_{m,n}(x,y,z)= \sqrt{\frac{2}{w(z)}} \phi_n(\frac{x\sqrt{2}}{w(z)}) \phi_m(\frac{y\sqrt{2}}{w(z)}) \exp\left(ik\frac{x^2 + y^2}{2R(z)} \right)\exp(-i \psi_{m,n}(z)) \quad (1)
$$
where
$$
\psi_{m,n}(z) = (n+m+1) \arctan\frac{z}{z_R} \\
w(z)=w_0 \sqrt{1+\frac{z^2}{z_R ^2}} \\
R(z) = z + \frac{z_R ^2}{z} \\
z_R = \frac{\pi w_0 ^2}{\lambda}
$$
and $\phi_n$ is the Hermite function defined as 
$$
\phi_n(x)=(-1)^n e^{x^2}\frac{\mathrm{d}^n}{\mathrm{d}x^n}e^{-x^2}
$$
Typically, for $n=m=0$, you recover the standard Gaussian propagation.
All this come from Maxwell equation. Now, there is a theorem that any function $v(x,y)$ can be expressed as a sum of functions $u_{n,m}$ :
$$
v(x,y)=\sum_{n,m} \, C_{n,m} u_{n,m}(x,y,0)
$$
with
$$
C_{n,m} = \int u^*_{n,m}(x,y,0)\times v(x,y) \,dxdy
$$
So you should take v(x,y) as a top hat and calculate $C_{n,m}$ - this will tell you how gaussian beams superpose to give a top hat at $z=0$. Then, the propagation equation (1) will tell you how this superposition looks like away from $z=0$.
