Is a Gaussian beam the outcome of the Fresnel diffraction of a Gaussian distribution? 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!
 A: 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. 
