Properties of Vector Spherical Harmonics In section 5.3.2 of the book Advanced Classical Electromagnetism by Robert Wald, in deriving the multipole expansion for the retarded solution of electromagnetic field in presence of charge-current distribution, it was asserted that
$$\nabla \cdot\left(h_{\ell}^{(1)}\left(\frac{\omega r}{c}\right) \boldsymbol{r} \times \nabla Y_{\ell m}\right)=0$$
and
$$\left[\nabla^{2}+\frac{\omega^{2}}{c^{2}}\right]\left(h_{\ell}^{(1)}\left(\frac{\omega r}{c}\right) \boldsymbol{r} \times \nabla Y_{\ell m}\right)=0$$
where $h_{\ell}^{(1)}$ is the Hankel function of the first kind, and $\boldsymbol{r} \times \nabla Y_{\ell m}$ is vector spherical harmonics.
The first equation is a property of vector spherical harmonics, as I saw on Wikipedia, but I am not able to derive this property myself. I also tried to prove the second equation by brute force, using the identity $\nabla^2(\boldsymbol{r} \times \nabla Y_{\ell m})=-\frac{l(l+1)}{r^2}(\boldsymbol{r} \times \nabla Y_{\ell m})$, but the terms does not cancel at the end.
I wonder if someone can explain how can theses two identities be shown, or suggest some literature I can look into.
 A: A good overview over identities and definition can be found on en.wikipedia.org/wiki/Vector_spherical_harmonics. In the notation of wikipedia we are dealing with identities
$$
\begin{align}
\nabla\cdot(f(r)\, \mathbf{\Phi}_{lm}) &= 0,\tag{1}\\
\nabla^2(f(r)\, \mathbf{\Phi}_{lm}) &= \Delta(f(r)\, \mathbf{\Phi}_{lm})  =  \left(\frac{1}{r^2}\frac{\partial}{\partial r}r^2\frac{\partial f(r)}{\partial r}\right) \mathbf{\Phi}_{lm} +f(r)\Delta\mathbf{\Phi}_{lm}\\
&= \left(\frac{1}{r^2}\frac{\partial}{\partial r}r^2\frac{\partial f(r)}{\partial r}\right) \mathbf{\Phi}_{lm} - \frac{1}{r^2}l(l+1)f(r)\mathbf{\Phi}_{lm},\tag{2}
\end{align}
$$
for
$$
\mathbf{\Phi}_{lm}=\mathbf{r}\times\nabla Y_{lm}.
$$
The first identity is simple to prove using the definition
$$
Y_{lm}\equiv\sqrt{\frac{2 l+1}{4 \pi }} \sqrt{\frac{(l-m)!}{(l+m)!}} P_l^m(\cos \theta)\mathrm{e}^{\mathrm{i} m \phi }\equiv c_{lm} P_l^m(\cos \theta)\mathrm{e}^{\mkern1mu\mathrm{i}\mkern1mum\mkern1mu\phi }\tag{3}
$$
and the relations for the gradient and divergence in spherical coordinates:
$$
\begin{align}
\nabla\cdot(f(r)\, \mathbf{\Phi}_{lm})&=\nabla\cdot(f(r)\,\mathbf{r}\times\nabla Y_{lm}) \tag{4.1}\\[.75em]
&= \nabla\cdot(f(r)\,\mathbf{r}\times \frac{c_{lm}\mathrm{e}^{\mkern1mu\mathrm{i}\mkern1mum\mkern1mu\phi }}{r}(0,-\sin\theta (P_l^{m})'(\cos\theta),\mathrm{i}\,m\csc\theta P_l^m(\cos\theta))^T ) \tag{4.2}\\[.75em]
&=\nabla\cdot(c_{lm}f(r){\mathrm{e}^{\mkern1mu\mathrm{i}\mkern1mum\mkern1mu\phi }}(0,-\mathrm{i}\,m\csc\theta P_l^m(\cos\theta),-\sin\theta (P_l^{m})'(\cos\theta))^T ) \tag{4.3}\\[.75em]
&=0+  \mathrm{i}\,m \,c_{lm}f(r)(P_l^{m})'(\cos\theta)\mathrm{e}^{\mkern1mu\mathrm{i}\mkern1mum\mkern1mu\phi }(1-1)=0\qquad\square,\tag{4.4}
\end{align} 
$$
where we used the notation $(a,b,c)^T\equiv a \hat{r}+b \hat{\theta}+c \hat{\phi}$.
Proofing the second identity (2) is a bit more involved before we start with an auxillary computation of $\Delta\mathbf{\Phi}_{lm}$ we recall two properties of Associated Legendre polynomials:
$$
\begin{align}
(P_l^{m})'(\cos\theta)&=(l-m+1) \csc \theta  P_{l+1}^m(\cos \theta )-(l+1) \cot \theta P_l^m(\cos \theta),\tag{5.1}\\
(1 + l - m)P_{l+1}^m(\cos \theta)&=(2 l+1) \cos \theta  P_l^m(\cos \theta)-(l+m) P_{l-1}^m(\cos \theta). \tag{5.2}
\end{align} 
$$
We now first explicitly compute $\mathbf{\Phi}_{lm}$ using Eq. (5.1):
$$
\mathbf{\Phi}_{lm} = -c_{lm}\mathrm{e}^{\mkern1mu\mathrm{i}\mkern1mum\mkern1mu\phi }\csc\theta \,(0,\mathrm{i}\, m P_l^m(\cos \theta),(l+1) \cos \theta P_l^m(\cos \theta )-(l-m+1) P_{l+1}^m(\cos\theta ))^T. \tag{6}
$$
Using the Laplacian in spherical coordinates and Eq. (5.2) repeatedly for $l+3$ and $l+2$ we recover the identity
$$
\Delta \mathbf{\Phi}_{lm} = -\frac{1}{r^2}l(l+1) \mathbf{\Phi}_{lm},\tag{7}
$$
which proves Eq. (2) after using the chain rule $\square$.
To apply Eq. (2) to $f(r)=h_l^{(1)}(\frac{\omega}{c}r)$ we look up (see e.g. dlmf.nist.gov/10.6) an identity for the derivative of the Hankel function
$$
(h_l^{(1)})'(z)=\frac{l}{z}h_l^{(1)}(z)-h_{l+1}^{(1)}(z) \tag{8}
$$
and using the chain and product rule with Eq. (8) we compute
$$
\frac{1}{r^2}\frac{\partial}{\partial r}r^2\frac{\partial}{\partial r}h_l^{(1)}\left(\frac{\omega}{c}r\right)=\left(\frac{l (l+1)}{r^2}-\frac{\omega^2}{c^2}\right)h_l^{(1)}\left(\frac{\omega}{c}r\right), \tag{9}
$$
which together with Eq. (2) proves
$$
\left(\nabla^2+\frac{\omega^2}{c^2}\right)\left(h_l^{(1)}\left(\frac{\omega}{c}r\right)\, \mathbf{\Phi}_{lm}\right) = 0 \qquad \square.\tag{10}
$$
