How to derive the equation relating the Laguerre Gaussian mode to the Hermite Gaussian mode of a laser beam 
Please help me derive the last equation, any tips are welcome. The $u^\text{LG}$ refers to Laguerre-Gaussian mode while $u^\text{HG}$ refers to the Hermite-Gaussian mode.
 A: Since the exponential factors dependent on $r$ and the $z$ dependent part of the profiles are easily matched up, it is enough to just check the part involving the Hermite and Laguerre polynomials.
For this keys step, Beijersbergen et. al. (1993) (which I assume is the paper where your expression comes from) cites a 1991 paper by Abramochkin and Volostnikov for a relation between the Hermite and Laguerre polynomials,
$$
\begin{multline}
\sum_{k=0}^{n+m}(2i)^k P_k^{(n-k,m-k)} (0) H_{n+m-k}(x)H_k(y) \\
= 2^{n+m}\times\begin{cases}
(-1)^m m! (x+iy)^{n-m} L_m^{n-m}(x^2+y^2) & \text{for } n\geq m\ ,\\
(-1)^n n! (x-iy)^{m-n} L_n^{m-n}(x^2+y^2) & \text{for } n \leq m\ ,
\end{cases}
\end{multline}
$$
where
$$
P_k^{(\mu,\nu)}(t) =\frac{(-1)^k}{2^k k!}(1-t)^{-\nu}(1-t)^{-\mu} \frac{d^k}{dt^k}[(1-t)^{k+\nu}(1+t)^{k+\mu}]\ .
$$
The main idea for proving of this formula is to cook up a generating function for both the Laguerre polynomials and the Hermite polynomials, then match up the two series to produce a relation between them.
A possibly cleaner presentation of this fact can be found in this paper. There the following generating functions are given for the two type of polynomials,
$$
\begin{gather}
\exp\left(2t_1 x + 2t_2 y - t_1^2 - t_2^2\right) = \sum_{n_1 = 0}^\infty\sum_{n_2 = 0}^\infty \frac{t_1^{n_1}t_2^{n_2}}{n_1!n_2!} H_{n_1}(x)H_{n_2}(y)\\
\begin{aligned}
\exp\left(s_1 z + s_2 \bar{z} - s_1 s_2\right) & = \sum_{m=0}^\infty\sum_{n=0}^\infty \frac{s_1^{m}s_2^{n}}{m!n!} (-1)^n n! z^{m-n} L_n^{m-n}(z \bar{z}) \\
 & = \sum_{m=0}^\infty\sum_{n=0}^\infty \frac{s_1^{m}s_2^{n}}{m!n!} (-1)^m m! \bar{z}^{n-m} L_m^{n-m}(z \bar{z})
\end{aligned}
\end{gather}
$$
If we set $t_1 = \frac{1}{2}(s_1 + s_2), t_2 = \frac{1}{2i}(s_2-s_1)$ and $z = x + i y, \bar{z} = x-iy$, then the two generating functions will become identical.
To match up the corresponding terms between the expansion in terms of $t_i$ and $s_i$, we can first write
$$
\begin{align}
t_1^{n_1} t_2^{n_2} & = \frac{1}{2^{n_1+n_2}i^{n_2}} (s_1 + s_2)^{n_1} (s_2 - s_1)^{n_2} \\
& = \frac{1}{2^{n_1+n_2}i^{n_2}} s_2^{n_1+n_2}\left(1 + \frac{s_1}{s_2}\right)^{n_1} \left(1 - \frac{s_1}{s_2}\right)^{n_2}.
\end{align}
$$
The last two factors dependent on $s_1/s_2\equiv t$ can be expanded as a Taylor series in $t$, so the $s_1^m s_2^n$ term in the generating function should have contributions from the $n$-th order term in the Taylor expansion for $t_1^{n_1} t_2^{n_2}$ where $n_1 + n_2 = m$.
You can now check that the $(2i)^{k}P_k^{(n-k,m-k)}(0)$ factor comes from the coefficient of the Taylor series and the factors of $2$ and $2i$ in the transformation between $s_i$ and $t_i$. Finally, the factorial prefactors in $b(n, m, k)$ are due to the factorials in the generating function as well as normalization factors for the polynomials. (Note that the Hermite and Laguerre polynomials defined above using generating functions have a different normalization from the conventional definition.)
The proof by Abramochkin and Volostnikov is essentially the same as the proof above, except that the generating function relations are rewritten term-wise using (essentially) Cauchy's integral formula.
