Propagating higher-order Hermite-Gauss modes using the Complex Beam Parameter?

A Gaussian laser beam can be propagated through an optical system (consisting of free space, thin lenses, curved and flat interfaces, etc) by using the "ABCD" ray-transfer matrices, and the complex beam parameter $\tilde{q}$.

A higher-order Hermite-Gauss or Laguerre-Gauss laser beam will gain Gouy phase more quickly than the fundamental Gaussian mode. Is there a simple modification to the complex beam parameter propagation that will also work for these higher order modes?

The complex beam parameter $\tilde{q}$, otherwise known as the complex radius of curvature, describes the transformation of the fundamental Gaussian mode through an optical system. All of the parameters of the higher order modes can be related to this fundamental mode transformation.
In the case of Gouy phase; it can be calculated relative to the waist for the Hermite-Gauss modes by $$\eta(\tilde{q})=(m+n+1)\arctan\left(\frac{\Re(\tilde{q})}{\Im(\tilde{q})}\right),$$ where $m$ and $n$ describe the particular higher order mode of interest. A similar expression gives the Gouy phase of the Laguerre-Gauss modes. When working with Gouy phase in an optical system, the quantity of interest is the total, accumulated Gouy phase. To calculate this you need to calculate the Gouy phase picked up in each portion seperately and add them together.