Some time ago some research led me to the following equality: \begin{equation} \frac{1}{\sin^2 \phi }\frac{\partial^2 }{\partial \theta^2} +\frac{\partial^2 }{\partial \phi^2} +\cot \phi \frac{\partial }{\partial \phi} = \frac{1}{\sin^2 \phi}(\frac{\partial}{\partial \theta} - i \sin \phi \frac{\partial}{\partial \phi})(\frac{\partial}{\partial \theta} + i \sin \phi\frac{\partial}{\partial \phi}) \end{equation} and it struck me as curious that it seemed to establish some specific connection between harmonic functions in the sphere and some sort of complex analytic functions also in the sphere. The factorization itself is elementary but I didn't remember seeing it before that. Is there any reference for this factorization?

Addendum: I will clarify what I mean by 'some sort of complex analytic functions'. If we focus on one of the factors in this factorization: \begin{equation} \frac{\partial}{\partial \theta} + i \sin \phi\frac{\partial}{\partial \phi} \end{equation} The first question to comes to mind is wether there is a non-trivial null-solution to this system. If we look for complex-valued functions that satisfy: \begin{equation} (\frac{\partial}{\partial \theta} + i \sin \phi\frac{\partial}{\partial \phi})(u + iv) = 0 \end{equation} Given the similarity of this system to Cauchy-Riemann (save for the sine factor) we can make the ansatz that \begin{equation} u = \theta \end{equation} And then comparing the real part of the above system to zero yields: \begin{equation} \frac{\partial v}{\partial \phi} = \frac{1}{\sin \phi} \end{equation} from which we can make \begin{equation} v = \ln \tan \frac{\phi}{2} \end{equation} So thus we can find a non-constant solution: \begin{equation} \theta + i \ln \tan \frac{\phi}{2} \end{equation} On the other side this solution suffers from discontinuities and singularities so using this as a starting point we can try to obtain a better kind of solution. For the first, one can check that solutions to this 'angular Cauchy-Riemann' are also closed under complex multiplication (so given two solutions their complex product is again a solution) and naturally addition. Moreover, since solutions are complex-valued then is possible to make complex functions act on these solutions by composition. One can straight-forwardly prove that these functions are closed under the composition by a complex analytical functions (in some suitable open set where of the composition makes sense, of course). In other words if $g$ is a solution to this 'angular Cauchy-Riemann' and $f$ is complex analytical then $f \circ g$ is again a solution to the angular Cauchy-Riemann. This is very convenient to improve the above non-trivial solution by first multiplying by $-i$ and then apply the complex exponential. We end up with the much better-looking: \begin{equation} \tan \frac{\phi}{2} e^{-i \theta} \end{equation} Which is (weakly) singular in the north pole but behaves fine everywhere else. I don't know if these are the 'best' solutions one can find but despite my efforts I could never find 'better' than this one. To conclude this addendum and in view of all the above: given any non-constant complex analytical function $f$ there is a 'angular analytic function': \begin{equation} f(\tan \frac{\phi}{2} e^{-i \theta}) \end{equation} whose real and imaginary parts are harmonic and that is defined in possibly all the sphere without one point. I couldn't make any substantial progress beyond this point.

Addendum 2:

This Cauchy-Riemann type factorization allowed us above to build some family of harmonic functions that depend only on parameters $\theta$ and $\phi$. I will try to draw a comparison with *irregular solid harmonics*. In order to consider them one has to allow singularities at the origin. An irregular solid harmonic is commonly written as:
\begin{equation}
I_{{\ell }}^{m}({\mathbf {r}})\equiv {\sqrt {{\frac {4\pi }{2\ell +1}}}}\;{\frac {Y_{{\ell }}^{m}(\theta ,\phi )}{r^{{\ell +1}}}}
\end{equation}
If we wanted to pick from these a solution whose radial part decreases like
\begin{equation}
\frac{1}{r}
\end{equation}
then $l=0$ leaves us with $Y_{{0 }}^{0}(\theta ,\phi )$, which is a constant. This is the best we get if we require our harmonic function to be defined in all the unit sphere. If we allow to remove one point of the sphere then we get as a counterpart harmonic functions whose radial part still decrease as $r^{-1}$ but are not constant as a function of the sphere minus one point removed. In fact the function:
\begin{equation}
\frac{\tan \frac{\phi}{2} e^{-i \theta}}{r}
\end{equation}
is harmonic in all its domain:
\begin{equation}
\Delta ( \frac{\tan \frac{\phi}{2} e^{-i \theta}}{r} ) = 0
\end{equation}