Optical beams diverge with propagation in free space due to their wave nature. With diffraction of the plane wavefront at an aperture, the linear divergence is determined by the corresponding angle in the far-field zone, proportional to the ratio of wavelength to aperture size. As a result, the product of aperture size by divergence angle divided by wavelength remains the same for a particular aperture shape, e.g., a cylindrical one, independent of aperture size. Similarly, one can characterize the diffraction properties of propagating localized aperture-free beams by the product of the minimal observed beam size times the divergence angle divided by the wavelength. This dimensionless parameter should depend only on the inner structure of the beam amplitude profile. One such propagation-invariant parameter introduced by A. Siegman describing the beam divergence in one transverse $x$-direction is defined as
$$M_x^2=2k\sqrt{\langle \Delta x^2_{\text{min}}\rangle} \sqrt{\langle\Delta \theta_x^2\rangle}.$$
Here $k=2\pi n_0/\lambda$ is the wave vector with vacuum wavelength $\lambda$ in linear lossless medium with refractive index (RI)
$n_0$; in the air, RI is taken as one. Second-order averages of the coordinate $x$ and angle $\theta_x$ are calculated from moments via
$\langle \Delta x^2\rangle = \langle x^2\rangle - \langle x\rangle ^2$.
The parameter $M_x^2$, also known as the
beam propagation ratio or
beam quality factor, is widely used for beam characterization in the laser science and industry.
The Helmholtz equation for electric field amplitude $E(\mathbf r =\{x,y,z\})$ in the case of a propagating coherent laser beam along $z$-direction is reduced to the scalar paraxial wave equation (PWE) for slow varying complex amplitude $U(\mathbf r)$:
$$E(\mathbf r)=U(\mathbf r)\cdot e^{ikz},\quad \frac{\partial U}{\partial z}-\frac{i}{2k} \left(\frac{\partial^2 U}{\partial x^2}+\frac{\partial^2 U}{\partial y^2}\right)=0.$$
If the transverse profile $U$ is known at some propagation position $z=z_1$, then this boundary condition determines the solution of PWE at
$z$ the Huygens integral in the Fresnel approximation. In the case of a normalized Gaussian amplitude profile
$U_\text{G}(x,y,0)=\sqrt{2/\pi}/w_0 \cdot \exp(-(x^2+y^2)/w_0^2)$, this well-known propagation solution of the PWE is:
$$ U_\text{G}(\mathbf r)=\tfrac{-ik}{2\pi z} \int{U_G(x',y',0) e^{ik\frac{(x-x')^2+(y-y')^2}{2z}}dx'dy'}=$$
$$=\tfrac{\sqrt{2/\pi}}{w_0(1+iz/z_\text{R})} e^{-\frac{x^2+y^2}{w_0^2(1+iz/z_\text{R})}},\quad \quad z_\text{R}=\tfrac{1}{2}k w_0^2, $$
$$ \langle x^2\rangle_\text{G}(z)=\frac{\int{x^2 |U_\text{G}|^2 dxdy}}{\int{|U_\text{G}|^2 dxdy=1}}=\tfrac{1}{4}w_0^2(1+z^2/z_\text{R}^2).$$
Here $z_\text{R}$ is the so-called Rayleigh length, where the square of the Gaussian beam size increases by two times.
The following propagation equations for the moments of coordinate and angle can be derived by applying the PWE and integrating by parts:
$$ \frac{d}{dz}\langle x\rangle=\langle \theta_x\rangle,\quad \frac{d}{dz}\langle \theta_x\rangle=0,$$
$$ \frac{d}{dz}\langle x^2\rangle=2\langle x\theta_x\rangle,\quad \frac{d}{dz}\langle x\theta_x\rangle=\langle \theta_x^2\rangle,\quad \frac{d}{dz}\langle \theta_x^2\rangle=0.$$
Explicit expressions for moments are presented in [1].
The propagation equations for the moments are leading to the propagation invariance of the beam quality parameter $M_x^2$ expressed as:
$$ M_x^2=2k\sqrt{\langle \Delta x^2\rangle \langle\Delta \theta_x^2\rangle - \langle\Delta x \Delta\theta_x\rangle^2 }.$$
This invariant can be calculated from the complex amplitude of the propagating beam $U(x,y,z)$ at any $z$;
here $\langle\Delta x \Delta\theta_x\rangle=\langle x \theta_x\rangle - \langle x \rangle \langle \theta_x\rangle$.
The propagation solution for $\langle\Delta x^2 \rangle (z)$ can be obtained in terms of the moments calculated at some particular $z=z_1$. Then the position $z_{\text{min}}$ of the beam waist with the minimal value of $\langle\Delta x^2 \rangle$ can be introduced correspondingly:
$$ \langle \Delta x^2\rangle(z)=\langle \Delta x^2\rangle_{\text{min}} + \langle \Delta \theta_x^2\rangle (z-z_{\text{min}})^2,$$
$$ z_{\text{min}}=z_1- \frac{\langle \Delta x \Delta\theta_x\rangle_1}{\langle \Delta \theta_x^2\rangle},\quad
\langle \Delta \theta_x^2\rangle = \frac{M_x^4}{4k^2\langle \Delta x^2\rangle_{\text{min}}}.$$
The beam has its waist at $z$-position where $\langle \Delta x \Delta\theta_x\rangle=0$. In particular, this occurs when the amplitude profile at this position is a real function.
Quadratic fitting of the experimental dependence $\langle \Delta x^2\rangle(z)$ gives $\langle \Delta x^2\rangle_{\text{min}}$ and $\langle \Delta \theta_x^2\rangle $, thus determining $M_x^2$.
In the case of a Gaussian beam $\langle \Delta x^2\rangle_{\text{min}}=w_0^2/4$, $\langle \Delta \theta_x^2\rangle =k^{-2}w_0^{-2}$ and $M_x^2=1$.
If, for an arbitrary beam, we consider a corresponding "underlying" Gaussian beam with the same size $\langle \Delta x^2\rangle_{\text{min}}$, then according to the equations above we obtain another equivalent definition of $M_x^2$ as the ratio of angular divergencies of the studied beam and a Gaussian beam of equal waist sizes:
$$ M_x^2=\sqrt{\frac{\langle \Delta \theta_x^2\rangle}{\langle \Delta \theta_x^2\rangle_\text{G}}}\quad\quad \text{at} \quad \langle \Delta x^2\rangle_{\text{min}}=\langle \Delta x^2\rangle_{\text{G,min}}.$$
Reference: [1] S. Mokhov, L. Glebov, and B. Zeldovich, "Quality deterioration of self-phase-modulated Gaussian beams," Laser Phys. Lett. 12, 015004−6 (2015). www.doi.org/10.1088/1612-2011/12/1/015004
Let us consider the beam propagating along the $z$-direction and diverging only at one transverse $x$-coordinate; and at some $z$-position, the beam amplitude $U(x)$ is represented as a sum of four lowest-order orthonormal Hermite-Gaussian functions $u(x)$ with complex coefficients $c_p=\alpha_p+i\beta_p$ (here $\alpha_p$ and $\beta_p$ are real):
$$ U(x)=\sum_{p=0}^3 c_p u_p(x),\quad u_p(x)=\tfrac{(2/\pi)^{1/4}}{\sqrt{2^p p!}}H_p(\sqrt{2}x)e^{-x^2},$$
$$ \int{u_p u_q dx}=\delta_{pq}:\quad u_0={\left(\tfrac{2}{\pi}\right)}^{\frac{1}{4}}e^{-x^2},\quad u_1={\left(\tfrac{2}{\pi}\right)}^{\frac{1}{4}}2x e^{-x^2},$$
$$ u_2={\left(\tfrac{2}{\pi}\right)}^{\frac{1}{4}}\tfrac{1}{\sqrt{2}}(4x^2-1)e^{-x^2},\quad u_3={\left(\tfrac{2}{\pi}\right)}^{\frac{1}{4}}\sqrt{\tfrac{2}{3}}(4x^3-3x)e^{-x^2}.$$
The computation script below calculates the beam propagation ratio $M^2$ based on inserted complex coefficients $c_p=\alpha_p+i\beta_p$. The plot shows the beam intensity profile $|U(x)|^2$. One can see that for one function $u_p$ alone, the $M^2=2p+1$.