Even if the formula I wrote has a typo
It has a typo, indeed. Bethe/Salpeter Q.M. of One- and Two-Electron Atoms, Plenum/Rosetta, 1977,p. 349, formula (A.37) reads
do you see the difference? And now the verification works in principle
In[13]:= FullSimplify[
Together[Sqrt[((l + m + 1) (l - m + 1))/((2 l + 3) (2 l +
1))] SphericalHarmonicY[l + 1,
m, \[Theta], \[CurlyPhi]] (f'[r] - l f[r]/r) +
Sqrt[((l + m) (l - m))/((2 l + 1) (2 l -
1))] SphericalHarmonicY[l - 1,
m, \[Theta], \[CurlyPhi]] (f'[
r] + (l + 1) f[r]/r) - (Cos[\[Theta]] D[#, r] -
Sin[\[Theta]]/r D[#, \[Theta]]) &[
f[r] SphericalHarmonicY[l, m, \[Theta], \[CurlyPhi]]]],
Element[l | m, Integers]]
Out[13]= $Aborted
In[15]:= FullSimplify[
Coefficient[
Sqrt[((l + m + 1) (l - m + 1))/((2 l + 3) (2 l +
1))] SphericalHarmonicY[l + 1,
m, \[Theta], \[CurlyPhi]] (f'[r] - l f[r]/r) +
Sqrt[((l + m) (l - m))/((2 l + 1) (2 l -
1))] SphericalHarmonicY[l - 1,
m, \[Theta], \[CurlyPhi]] (f'[
r] + (l + 1) f[r]/r) - (Cos[\[Theta]] D[#, r] -
Sin[\[Theta]]/r D[#, \[Theta]]) &[
f[r] SphericalHarmonicY[l, m, \[Theta], \[CurlyPhi]]], f[r]],
Element[l | m, Integers]]
Out[15]= (E^(-I \[CurlyPhi]) (Sqrt[Gamma[1 + l - m]] Sqrt[
Gamma[2 + l + m]]
Sin[\[Theta]] SphericalHarmonicY[l,
1 + m, \[Theta], \[CurlyPhi]] +
E^(I \[CurlyPhi]) Sqrt[Gamma[l - m]] Sqrt[
Gamma[1 + l +
m]] ((1 + l) Sqrt[((l - m) (l + m))/(-1 + 4 l^2)]
SphericalHarmonicY[-1 + l, m, \[Theta], \[CurlyPhi]] +
m Cos[\[Theta]] SphericalHarmonicY[l,
m, \[Theta], \[CurlyPhi]] -
l Sqrt[((1 + l - m) (1 + l + m))/(3 + 4 l (2 + l))]
SphericalHarmonicY[1 + l,
m, \[Theta], \[CurlyPhi]])))/(r Sqrt[Gamma[l - m]] Sqrt[
Gamma[1 + l + m]])
In[16]:= fsf = %15;
In[24]:= fsf /. {l -> RandomInteger[73],
m -> RandomInteger[99], \[Theta] ->
RandomReal[{0, \[Pi]}], \[CurlyPhi] -> RandomReal[{0, 2 \[Pi]}]}
Out[24]= (1.18791*10^-14 - 3.56613*10^-15 I)/r
In[17]:= FullSimplify[
Coefficient[
Sqrt[((l + m + 1) (l - m + 1))/((2 l + 3) (2 l +
1))] SphericalHarmonicY[l + 1,
m, \[Theta], \[CurlyPhi]] (f'[r] - l f[r]/r) +
Sqrt[((l + m) (l - m))/((2 l + 1) (2 l -
1))] SphericalHarmonicY[l - 1,
m, \[Theta], \[CurlyPhi]] (f'[
r] + (l + 1) f[r]/r) - (Cos[\[Theta]] D[#, r] -
Sin[\[Theta]]/r D[#, \[Theta]]) &[
f[r] SphericalHarmonicY[l, m, \[Theta], \[CurlyPhi]]], f'[r]],
Element[l | m, Integers]]
Out[17]= Sqrt[((l - m) (l + m))/(-1 + 4 l^2)]
SphericalHarmonicY[-1 + l, m, \[Theta], \[CurlyPhi]] -
Cos[\[Theta]] SphericalHarmonicY[l, m, \[Theta], \[CurlyPhi]] +
Sqrt[((1 + l - m) (1 + l + m))/(3 + 4 l (2 + l))]
SphericalHarmonicY[1 + l, m, \[Theta], \[CurlyPhi]]
In[18]:= fsfs = %17;
In[20]:= fsfs /. {l -> RandomInteger[73],
m -> RandomInteger[99], \[Theta] ->
RandomReal[{0, \[Pi]}], \[CurlyPhi] -> RandomReal[{0, 2 \[Pi]}]}
Out[20]= 0.
I say in principle because the fsfs term gave in all runs practically zero, but the fsf term might give also Indeterminate depending on the chosen quantum numbers and angles. That needs inspection. At least now there is a valuable candidate for the identity. Even more, there is also a possible candidate to solve p.o.'s original question because every summand on the r.h.s. now contains a SphericalHarmonicY guiding the multiple application of the rule. Go ahead!
