Solution (source code)

= Solution
{c}

Put $s=z-z_0$. In the source's <complex beam parameter> convention,
$$
\frac1q=\frac{s+iz_R}{s^2+z_R^2}
=\frac1{R(z)}+i\frac{2}{kw(z)^2}.
$$
Substitution separates the transverse <wave amplitude> and <phase>:
$$
U=\frac{U_0e^{i(kz-\omega t)}}{s-iz_R}
\exp\!\left[-\frac{\rho^2}{w(z)^2}\right]
\exp\!\left[i\frac{k\rho^2}{2R(z)}\right].
$$
Thus the <Gaussian beam> is circularly symmetric, has a <Gaussian function> transverse <light intensity>, and narrows to a <beam waist> at $z_0$ before spreading by <diffraction>. Its <wavefront> is planar at the <beam waist>, converging before it and diverging after it. The quadratic transverse <phase> identifies the signed <radius of curvature>:
$$
\boxed{R(z)=s+\frac{z_R^2}{s},\qquad R(z_0)=\infty.}
$$
Define spot size as the <wave amplitude> $1/e$ <radius>, equivalently the <light intensity> $1/e^2$ <radius>. Then
$$
\boxed{w(z)^2=\frac{2(s^2+z_R^2)}{kz_R}
=w_0^2\left(1+\frac{s^2}{z_R^2}\right),\qquad
w_0^2=\frac{2z_R}{k}=\frac{\lambda z_R}{\pi}.}
$$
The <Rayleigh range> is where $w=\sqrt2w_0$; the far-field <beam divergence> half-angle is $\lambda/(\pi w_0)$. The factor $1/q$ also produces a longitudinal <phase> shift relative to a <plane wave>, changing by $\pi$ across the <beam waist>; its sign here follows the given $e^{i(kz-\omega t)}$ convention. The on-axis <light intensity> decreases as $1/(s^2+z_R^2)$ while transverse integration leaves the <power> constant. These are scalar <paraxial approximation> results.