Proper and rigourous derivation of Gaussian beam? Gaussian beams are known solutions to the Paraxial Wave Equation:
$$ \frac{\partial^2 \Psi(x,y,t)}{\partial^2 x} + \frac{\partial^2 \Psi(x,y,t)}{\partial^2 y} = 2ik\frac{\partial \Psi(x,y,t)}{\partial t} $$
I am looking for a proof (or a reference to it) that the general solution to this equation is a Gaussian beam?
Since they usually just assume that a Gaussian beam is OK, they plug it in and show that it works.
 A: I don't think "the general solution to this equation is a Gaussian beam" - this is a linear equation, so a superposition of solutions is also a solution, and an arbitrary superposition of Gaussian beams is not always a Gaussian beam.
A: The general solution to the paraxial wave equation,
$$\left[ \frac{\partial^2 }{\partial^2 x} + \frac{\partial^2 }{\partial^2 y} - 2ik\frac{\partial}{\partial z} \right]\Psi(x,y,z)=0,$$
is not a gaussian beam. To see this, you should note that this equation is now of parabolic type, which essentially means that you can take any plane at $z=z_0=\text{const}$ and set the values of your scalar function there,
$$
\Psi(x,y,z_0)
$$
as your boundary conditions. This means that there is a one-to-one correspondence between the solutions of the paraxial wave equation and the (nice enough) scalar functions of two variables. The latter are much more varied than just gaussians, so that the solutions need to be more general as well.
Physically, this means that you can set essentially any optical mask at any point in your beam and it will still propagate. (Of course, you need to ensure that the spatial pattern of your mask does not include spatial frequencies higher than (or comparable to) $k$, or you will step outside of the paraxial regime.)
Another way to see this is to notice that the paraxial wave equation is exactly the same as the time-dependent Schrödinger equation for a free particle in two dimensions, which you've already noticed. The general solution to this is, of course, much more general than a gaussian.

So, why do people jump so quickly on the gaussian solution? The best answer to this is that you're not usually out to describe the paraxial wave equation; you're out to describe laser beams and that means that you only care about some specific subset of solutions, so it's sort of OK to not go for the general solution.
Nevertheless, the gaussian solution is indeed 'special'. In particular, it will minimize the divergence of the beam in the far field for a given fixed width at the beam waist. I'm not completely sure of how to prove this (or even the best rigorous formulation) but it is at least morally true.
If you want to go beyond the simple gaussian solution, there are several important families of solutions which are used to describe almost-gaussian beams, as well as some rather different solutions. The main two families are:


*

*Hermite-Gaussian modes, which are of the form
$$
\begin{align*}
\Psi(x,y,z)&=\sqrt{2\over\pi}\frac{1/w_0}{\sqrt{2^{n_x+n_y}n_x!n_y!}}\frac{q_0}{q(z)}
\left[\frac{q_0}{q_0^*}\frac{q^*(z)}{q(z)}\right]^\frac{n_x+n_y}{2}
\\ & \qquad\qquad \times
H_{n_x}\left(\frac{\sqrt{2}x}{w(z)}\right)H_{n_y}\left(\frac{\sqrt{2}y}{w(z)}\right)
\exp\left[-i\frac{kx^2}{2q(z)}\right],
\end{align*}
$$
and

*Laguerre-Gaussian modes, which go like
$$
\begin{align*}
\Psi(r,\phi,z)&=
\frac{C_{\ell p}}{w(z)}\left(\frac{r\sqrt 2}{w(z)}\right)^{|\ell |}e^{i\ell \phi}
\exp\left[-\frac{r^2}{w^2(z)}\right]
L_{p}^{|\ell |}\left(\frac{2r^2}{w^2(z)}\right)
\\ & \qquad\qquad \times
\exp\left[ik\frac{r^2}{2R(z)}\right]
e^{i(2p+|\ell |+1)\zeta(z)},
\end{align*}
$$
though other families exist for more exotic geometries.
In essence, these families provide convenient bases in which to treat a beam which is almost, but not quite, gaussian. Which one you use depends on what geometry is most relevant to your problem, but in fact they are very directly equivalent.
Alternatively, you can use these solutions on their own, and then they describe very special beams. For example, inserting a thin wire in the middle of your beam inside the laser cavity will make it lase into a first-order Hermite-Gaussian mode, with a node down the middle. To generate a Laguerre-Gaussian beam you use a spatial light modulator with a phase singularity, and the resulting beam will have orbital angular momentum which you can use to spin particles around an optical-tweezer trap. A large field of research is dedicated to making and using light beams in these and more exotic geometries.
One thing to note, though, is that in all such families, the (non-gaussian) beams of higher order always have some nontrivial polynomial of the coordinates over $w(z)$ multiplying the core gaussian dependence. This polynomial inherently makes them wider at their beam waist, but their angular divergence is the same as the gaussian. This is, partly, the basis of the claim I made above.
A: The other answer is correct that the gaussian beam isn't a general solution. What I'm showing is the motivation of the gaussian beam as a particular solution. The derivation can be made by starting with the spherical wave solution to the wave equation
$$U(\textbf{r})=\frac{A_0}{r}e^{-ikr}$$
and assuming that far from the origin, a spherical surface is approximately paraboloidal. We also assume that we're looking at points that are very close to the z axis, so that z is large, but $x^2+y^2$ is small.
$$r = \sqrt{x^2+y^2+z^2}=z\sqrt{\frac{x^2+y^2}{z^2}+1}$$
$$\theta ^2 :=\frac{x^2+y^2}{z^2}$$
By taylor expansion
$$r=z\sqrt{\theta^2+1}\approx z(1+\theta^2/2)=z+\frac{x^2+y^2}{2z}$$
We substitute this expression into the phase, and use the less accurate approximation $z\approx r$ for the amplitude. This gives us the expression
$$U(\textbf{r})=\frac{A_0}{z}\exp[-ik\frac{x^2+y^2}{2z}] e^{-ikz}$$
We notice that this takes the form of a slowly varying envelope $U(\textbf{r})=A(\textbf{r})e^{-ikz}$, and the paraxial wave equation is usually derived assuming a solution of this form. So taking this particular solution as a given, we notice the substitution $z \rightarrow z+iz_0$ also produces a solution. The final expression for the gaussian beam is then found by taking this equation
$$U(\textbf{r})=\frac{A_0}{z+iz_0}\exp[-ik\frac{x^2+y^2}{2(z+iz_0)}] e^{-ikz}$$
rearranging it to separate the complex phase from the amplitude, and choosing variables to make the expression look nicer.
