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!

Thank you very much for the code. Indeed, the function shRubinD you have defined doesn't give the correct result when applied twice. It doesn't take into account that the indices 'l' and 'm' have changed.

My question isn't limited to SphericalHarmonics and can be refined as following. Assume we have the following recursive definition of a symbolic Operator which operates on a function H[y] Subscript[G, l][x]

Operator[H[y] Subscript[G, l][x]] = 
 l D[H[y], y] Subscript[G, l + 1][x] + 
(l + 1) D[H[y], y] Subscript[G, l - 1][x]

Here, 'x' and 'y' are two different arguments while 'l' is an index. How can one "teach" Mathematica to increase and decrease the index of the function G, as well as to take proper combination of the H function and yield a symbolic result due to applying the Operator arbitrary number of times? Particularly, this will enable to obtain the result of Operator^2.

Doesn't this tend to do the job?

In[47]:= Clear[shRubinD]
shRubinD[f_, x_, l_, m_, y_, z_] := 
  Sqrt[((l + m + 1) (l - m + 1))/((2 l + 3) (2 l + 
           1))] SphericalHarmonicY[l + 1, m, y, 
      z] (D[f, x] - l f/x) + 
    Sqrt[((l + m) (l - m))/((2 l + 1) (2 l - 1))] (D[f, 
        x] + (l + 1) f/x) /; 
   AtomQ[l] && AtomQ[m] && AtomQ[x] && AtomQ[y] && 
    AtomQ[z] && ! FreeQ[f, x];
shRubinD[f_, x_, n_Integer, l_, m_, y_, z_] := 
 Nest[shRubinD[#, x, l, m, y, z] &, f, n] /; n > 1

In[50]:= shRubinD[Sin[x], x, p, y, \[CurlyPhi], \[CurlyTheta]]

Out[50]= Sqrt[((p - y) (p + y))/((-1 + 2 p) (1 + 2 p))] (Cos[
     x] + ((1 + p) Sin[x])/x) + 
 Sqrt[((1 + p - y) (1 + p + y))/((1 + 2 p) (3 + 2 p))] (Cos[x] - (
    p Sin[x])/x) SphericalHarmonicY[1 + p, 
   y, \[CurlyPhi], \[CurlyTheta]]

In[52]:= shRubinD[Sin[x], x, 2, p, 
  y, \[CurlyPhi], \[CurlyTheta]] // FullSimplify

Out[52]= (1/(x^2))(((p - y) (p + 
     y) (2 (1 + p) x Cos[x] + (p + p^2 - x^2) Sin[x]))/(-1 + 4 p^2) - 
  Sqrt[((p - y) (p + y))/(-1 + 4 p^2)]
    Sqrt[((1 + p - y) (1 + p + y))/(
   3 + 4 p (2 + p))] (-2 x Cos[x] + (1 + 2 p (1 + p) + 2 x^2) Sin[
       x]) SphericalHarmonicY[1 + p, 
    y, \[CurlyPhi], \[CurlyTheta]] + ((1 + p - y) (1 + p + 
     y) (-2 p x Cos[x] + (p + p^2 - x^2) Sin[x]) SphericalHarmonicY[
    1 + p, y, \[CurlyPhi], \[CurlyTheta]]^2)/(3 + 4 p (2 + p)))

Well, it is extraordinarily ugly because it has SphericalHarmonicY as inmate. But you seem to intend that.

Thanks. Yet, the problem I like to solve seems to demand other solution.

As far as I see, the main difficulty is that the function needs to recognize the value of the first argument of the Spherical Harmonic function and then modify it accordingly. I don't know how to access the argument of the Spherical Harmonic function (or any other) in case it is multiplied by another function (which is relevant here according to the identity).

Thanks. Can you please share more details? It is hard for me to follow.

Which pattern matching function ti use? How do I access the argument of Spherical Harmonic function which is multiplied by another function?

Thanks for good intention, nevertheless I couldn't understand any of your answers.

Thanks. That was the first thing I tried. However, taking the second derivative of the Spherical Harmonic with respect to the variable 'z' directly (using the chain rule), yields a very long expression which isn't reduced upon using FullSimplify

Simplify[-Sin[\[Theta]]/
  r D[-Sin[\[Theta]]/
    r D[SphericalHarmonicY[l, m, \[Theta], phi], \[Theta]], \[Theta]]]

(E^(-2 I phi) (E^(I phi) Sqrt[
      Gamma[-1 + l - 
        m]] (E^(I phi) m (-1 + (1 + m) Cos[\[Theta]]^2) Sqrt[
         Gamma[l - m]] Sqrt[Gamma[1 + l + m]]
          SphericalHarmonicY[l, m, \[Theta], phi] + (1 + m) Sqrt[
         Gamma[1 + l - m]] Sqrt[Gamma[2 + l + m]]
          Sin[2 \[Theta]] SphericalHarmonicY[l, 1 + m, \[Theta], 
          phi]) + Sqrt[Gamma[l - m]] Sqrt[Gamma[1 + l - m]] Sqrt[
      Gamma[3 + l + m]]
       Sin[\[Theta]]^2 SphericalHarmonicY[l, 2 + m, \[Theta], 
       phi]))/(r^2 Sqrt[Gamma[-1 + l - m]] Sqrt[Gamma[l - m]] Sqrt[
   Gamma[1 + l + m]])

Unfortunately using FullSimplify doesn't use the internal relations between Spherical Harmonics and the compact relation for the derivative I wrote above isn't reproduced. Since I am interested in the simplest symbolic expression and not just on numeric value, it doesn't resolve the issue.

Just used the chain rule, yielding that the derivative with respect to the Cartesian 'z' coordinate in Spherical coordinates is given by

 Cos[\[Theta]] (\[PartialD]/\[PartialD]r)-(Sin[\[Theta]]/
   r) \[PartialD]/\[PartialD]\[Theta]

Next if interested in operating on 'r' independent function such as Spherical Harmonic just drop the partial derivative with respect to 'r'.

Probably better to define your function, rather than redefine the built-in derivative.