## 2 Solution of the Paraxial Equation in Elliptic Coordinates  
The paraxial Helmholtz equation is given by [3]:  
\[\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+ 2 ik\frac{\partial}{\partial z}\right)E(\mathbf{r},z)=0, \tag{1}\]  
where \((x,y,z)\) are the Cartesian coordinates; \(\mathbf{r}=(x,y)=(r\cos\varphi,r\sin\varphi)\) is a 2D vector; \(z\) is the coordinate along the optical axis; \(k\) is the wavenumber of light; and \(E(\mathbf{r},z)\) is the complex amplitude of a monochromatic light field. Below, we use the following dimensionless coordinates: \(x:=x/w_{0}\), \(y:=y/w_{0}\), \(z:=2z/(kw_{0}^{2})\), with \(w_{0}\) being the waist radius of the Gaussian beam. Then, in dimensionless variables, Equation (1) reads as follows:  
\[\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}} +4i\frac{\partial}{\partial z}\right)E(\mathbf{r},z)=0. \tag{2}\]  
Its simplest solution with a finite energy is a conventional Gaussian beam:  
\[G(\mathbf{r},z)=\frac{1}{1+iz}\exp\left(-\frac{r^{2}}{1+iz}\right). \tag{3}\]  
It is known that separation of variables in Equation (2) into Cartesian and polar coordinates makes it possible to obtain two solution families--HG and LG beams, respectively [17]:  
\[\begin{split}&\mathcal{H}_{n,m}(\mathbf{r},z)=G(\mathbf{r},z) \left(\frac{1-iz}{1+iz}\right)^{\frac{n+m}{2}}H_{n}\left(\frac{x\sqrt{2}}{ \sqrt{1+z^{2}}}\right)H_{m}\left(\frac{y\sqrt{2}}{\sqrt{1+z^{2}}}\right),\\ &\mathcal{L}_{n,\pm m}(\mathbf{r},z)=G(\mathbf{r},  
Both families are examples of structurally stable light fields, i.e., at any propagation distance \(z\), the transverse intensity pattern of such fields aligns--up to scale--with that in the initial plane \(z=0\):  
\[ \begin{array}{l}\mathcal{H}_{n,m}(\mathbf{r},z)=\frac{1}{\sqrt{1+z^{2}}} \mathcal{H}_{n,m}\Big{(}\frac{\mathbf{r}}{\sqrt{1+z^{2}}},0\Big{)}\exp\Big{(} \frac{izr^{2}}{1+z^{2}}-i(n+m+1)\text{arctan}z\Big{)},\\ \mathcal{L}_{n,\pm m}(\mathbf{r},z)=\frac{1}{\sqrt{1+z^{2}}}\mathcal{L}_{n,\pm m }\Big{(}\frac{\mathbf{r}}{\sqrt{1+z^{2}}},0\Big{)}\exp\Big{(}\frac{izr^{2}}{1+ z^{2}}-i(2n+m+1)\text{arctan}z\Big{)}.\end{array}  \tag{5}\]  
Thus, it will suffice to study both families at \(z=0\), omitting the last argument for brevity: \(\mathcal{H}_{n,m}(\mathbf{r},0)=\mathcal{H}_{n,m}(\mathbf{r})\) and \(\mathcal{L}_{n,\pm m}(\mathbf{r},0)=\mathcal{L}_{n,\pm m}(\mathbf{r})\). We also introduce normalized variants of each solutions family:  
\[ \begin{array}{l}\mathbf{HG}_{n,m}(\mathbf{r})=\frac{\mathcal{H}_{n,m}( \mathbf{r})}{\|\mathcal{H}_{n,m}(\mathbf{r})\|},\|\mathcal{H}_{n,m}(\mathbf{r })\|=\sqrt{\frac{\pi}{2}\cdot 2^{n+m}n!m!},\\ \mathbf{LG}_{n,\pm m}(\mathbf{r})=\frac{\mathcal{L}_{n,\pm m}(\mathbf{r})}{\| \mathcal{L}_{n,\pm m}(\mathbf{r})\|},\|\mathcal{L}_{n,\pm m}(\mathbf{r})\|= \sqrt{\frac{\pi}{2}\cdot\frac{(n+m)!}{2^{2}n!}}.\end{array}  \tag{6}\]  
When Equation (2) is solved, the variables are separated in the elliptic coordinates \((\xi,\eta)\)[3]; then, a family of IG beams, which are also structurally stable, is obtained:  
\[\mathcal{I}(\mathbf{r},z)=G(\mathbf{r},z)\bigg{(}\frac{1-iz}{1+iz}\bigg{)}^{p/ 2}P(\xi)Q(\eta). \tag{7}\]  
Here,  
\[\frac{x}{\sigma\sqrt{1+z^{2}}}=\cosh\xi\cos\eta,\frac{y}{\sigma\sqrt{1+z^{2}} }=\sinh\xi\sin\eta, \tag{8}\]  
and the parameter \(\sigma\) is supposed to be arbitrary. Substitution of Equation (7) into Equation (2) yields equations for the functions \(P(\xi)\) and \(Q(\eta)\):  
\[ \begin{array}{l}\Big{(}\frac{\partial^{2}}{\partial t^{2}}-2\sigma^{2}\sinh 2 \xi\frac{\partial}{\partial\xi}-(\lambda-2p\sigma^{2}\cosh 2\xi)\Big{)}P(\xi)=0,\\ \Big{(}\frac{\partial^{2}}{\partial t^{2}}+2\sigma^{2}\sin 2\eta\frac{ \partial}{\partial\eta}+(\lambda-2p\sigma^{2}\cos 2\eta)\Big{)}Q(\eta)=0.\end{array}