Is a Gaussian beam the outcome of the Fresnel diffraction of a Gaussian distribution?

Question

Rigorous diffraction theory states that if you know a scalar function that satisfies the helmholtz equation, $u(\vec{x})$, at a plane (for example the $xy$-Plane), given as $u(\vec{x}) = u_0(x,y)$, then you can calculate it at every other point by convoluting it with the function $$\tag{1}\partial_z \frac{e^{ik(|\vec{x}|)}}{|\vec{x}|}$$ (that would be Rayleigh-Sommerfeld-Diffraction Theory).

This would be equivalent to a Fourier-decomposition of $u(\vec{x})$ into plane waves with wave-number $k$, and propagating these through space. In this framework, you would get the solution by calculating the Fourier transform $\tilde{u}(\vec{k})$: $$\tilde{u}(\vec{k}) = e^{iz\sqrt{k^2-k_x^2 - k_y^2}} \tilde{u}_0(k_x, k_y) . $$

In this case, the function in (1) is exactly the Fourier transform of $$e^{iz\sqrt{k^2-k_x^2 - k_y^2}} . $$

The paraxial approximation assumes that there are only components of the solution with $k_x \ll k_z$ and $k_y \ll k_z$. If this is the case, we can calculate the Fourier transform as: $$\tilde{u}(\vec{k}) = e^{iz(k - \frac{k_x^2+k_y^2}{2k})} \tilde{u}_0(k_x, k_y)$$

Equivalently, obtain the solution $u(\vec{x})$ by convolving with the function $$h(x, y, z) = \frac{-i k}{z} e^{i(kz + \frac{k(x^2 + y^2)}{2 z})} \\ u(x,y,z) = \int dx' \int dy' u_0(x',y') h(x - x', y - y', z) $$ which is exactly the way to calculate diffraction patterns in Fresnel diffraction.

Question 1: Does that mean, Fresnel diffraction is equivalent to the assumption that only paraxial waves contribute to a function?

Question 2: Does that mean, the outcome of Fresnel diffraction $u(x,y,z)$ suffices the paraxial wave equation?

Question 3 (most important): Does that mean, Fresnel diffraction of the function $$u_0(x,y) = A_0 e^{-\frac{(x^2 + y^2)}{w_0^2}}$$ will yield the exact formula of a Gaussian beam? I tried to calculate it, but I failed: When I convolve $$u_0(x,y) = A_0 e^{-\frac{(x^2 + y^2)}{w_0^2}}$$ with $$h(x, y, z) = \frac{-i k}{z} e^{i(kz + \frac{k(x^2 + y^2)}{2 z})} , $$ I don't arrive at the formula for a Gaussian beam. I still think that this should be possible. Does somebody know where this is explicitly calculated?

Edit: Siegman states in his Book "Lasers" at the end of chapter 16: " This gaussian-spherical wave solution given by Equations 16.29 - 16.35. is still an exact mathematical solution to either the paraxial wave equation or the Huaygens-Fresnel- integral". That means Siegman states the answer to my question is yes. But he doesn't give a calculation!

Answer 1

I finally was able to calculate something satisfying: I assume a gaussian beam to be expressed as a complex parabolic wave: $$ u_\tilde{z}(\tilde{x}, \tilde{y})=\frac{u_0}{i q_0 + z} e^{-i(q_0 i + z)k - i k \frac{\tilde{x}^2 + \tilde{y}^2}{2(q_0 i + z) }} $$ Fresnel Diffraction (I will denote the outcome of the diffraction by $\bar{u}:$ $$ \bar{u}(x, y, z + \tilde{z})=\int d\tilde{x} \int d\tilde{y} u_\tilde{z}(\tilde{x}, \tilde{y}) \frac{i}{\lambda z} e^{-ikz-ik\frac{(x-\tilde{x})^2 + (y - \tilde{y})^2}{2 z} } $$ Plugging everything in will result in: $$ \frac{u_0}{i q_0 + \tilde{z}} e^{-i(q_0 i + \tilde{z})} \frac{i}{\lambda z} e^{-ikz} \int d \tilde{x} e^{-ik \frac{x^2}{2z}} e^{-(\frac{ik}{2(q_0 + \tilde{z})} + \frac{ik}{2z})\tilde{x}^2 - (-\frac{ikx}{z})\tilde{x} } \int d\tilde{y} \text{the same stuff} $$ Using $\int e^{-ax^2-bx} dx = \sqrt{\frac {\pi}{ a}} e^\frac{b^2}{a}$ you arrive at: $$ \frac{u_0}{i q_0 + \tilde{z}} e^{-i(q_0 i + \tilde{z})} \frac{i}{\lambda z} e^{-ikz} \frac{\pi}{ik(\frac{1}{2(\tilde{z} + q_0 i)}+ \frac{1}{2 z})} e^{\frac{-ikx^2}{2z} + ik \frac{x^2}{2z^2(\frac{1}{\tilde{z} + q_0 i}+ \frac{1}{ z})}} e^{ \text{the same stuff with y}} $$ Looking at the argument of the exponential function with the $x$-dependency: $$ -\frac{ikx^2}{2z} + \frac{ikx^2}{z^2}\frac{1}{\frac{1}{(\tilde{z} + q_0 i)}+ \frac{1}{ z}} = \frac{-ikx^2}{2}(\frac{1}{z} - \frac{1}{z^2(\frac{1}{\tilde{z} + q_0 i} + \frac{1}{z})}) \\ = \frac{-ikx^2}{2}\frac{z^2(\frac{1}{\tilde{z} + q_0 i} + \frac{1}{z})-z}{z^3(\frac{1}{\tilde{z} + q_0 i} + \frac{1}{z})} = \frac{-ikx^2}{2} \frac{\frac{z^2}{\tilde{z}+q_0 i}}{z^3 (\frac{1}{\tilde{z} + q_0 i}+ \frac{1}{z})} \\ = \frac{-i k x^2 }{2} \frac{1}{z(1 + \frac{\tilde{z} + q_0 i}{z} ) } = \frac{ik x^2}{2} \frac{1}{z + \tilde{z} + q_0 i} $$ That's the argument of the usual complex parabolic ray. In the same manner the factor infront becomes: $$ \frac{u_0}{z + \tilde{z} + q_0 i} $$

And all in all, we can conclude $\bar{u}(x, y, z+ \tilde{z}) = u(x, y, z + \tilde{z})$

I can conclude that fresnel diffraction won't change a gaussian ray.