SphericalHarmonicY does not seem to be an eigenfunction of the spherical harmonic equation

The operator in the question can be simplified by first applying it to an arbitrary function f:

Simplify[(1/Sin[θ]) D[Sin[θ] D[f[θ,ϕ],θ],θ]+(1/Sin[θ]^2) D[f[θ,ϕ],{ϕ,2}]]

$$f^{(2,0)}(\theta ,\phi )+\cot (\theta ) f^{(1,0)}(\theta ,\phi )+\csc ^2(\theta ) f^{(0,2)}(\theta ,\phi )$$

This can be recognized as the negative of the squared angular momentum operator, $-L^2$, see below (and see also my answer here). It doesn't seem possible to do the simplification exactly as asked for in the question, but you can put together a Mathematica based proof of the desired eigenvalue property $$L^2 Y^{m}_{\ell} = \ell (\ell +1) Y^{m}_{\ell}$$ in several steps. In doing so, the advantage is that you can see more clearly how this equation is related to the properties of the angular-momentum operator (the brute-force approach in the question would be shorter, but less illuminating, even if it were possible).

The spherical harmonics are simultaneous eigenfunctions of the squared angular momentum and its z component. Here is how you can verify this by asking Mathematica to derive the two corresponding eigenvalues:

First I define the two operators to make the notation in the following steps clearer. Then I apply these operators in the eigenvalue equation with the spherical harmonic as the assumed solution. The proof consists in the fact that Mathematica does indeed find an eigenvalue that makes each relation true.

Clear[θ, ϕ, f, m, ℓ]
lzOp = Function[f, (-I)*D[f, {ϕ}]];
l2Op = Function[f, -D[f, {θ, 2}] - Csc[θ]^2*D[f, {ϕ, 2}] - Cot[θ]*D[f, {θ}]];

Simplify[
 lzOp[SphericalHarmonicY[ℓ, m, θ, ϕ]]]

m SphericalHarmonicY[ℓ, m, θ, ϕ]

Simplify[
 l2Op[SphericalHarmonicY[ℓ, ℓ, θ, ϕ]]]

ℓ (1 + ℓ) SphericalHarmonicY[ℓ, ℓ, θ, ϕ]

So both lzOp (the z-angular momentum) and l2Op (the squared angular momentum) reproduce the spherical harmonic with a prefactor that corresponds to the correct eigenvalue.

However, the last eigenvalue equation with l2Op (which is what the question focuses on) only works for m = ℓ. For the remaining values of m, I would proceed by using the ladder operators which turn $Y^m_{\ell}$ into $Y^{m\pm 1}_{\ell}$ and also commute with l2Op. These operations will then necessarily leave the eigenvalue ℓ (1 + ℓ) unchanged for l2Op.

To execute this idea, I have to define the lowering operator, lMinus:

lMinus = Function[f, Exp[-I ϕ] (Cot[θ] D[f, ϕ] + I D[f, θ])];

First I now show that lMinus indeed has the effect of lowering the z eigenvalue of any spherical harmonic:

Simplify[
 First[Solve[
   lzOp[lMinus@SphericalHarmonicY[ℓ, m, θ, ϕ]] ==
     eigenvalue lMinus@SphericalHarmonicY[ℓ, m, θ, ϕ], 
   eigenvalue]
 ]]

{eigenvalue -> -1 + m}

Here I had to use Solve to extract the eigenvalue property by defining the unknown eigenvalue as a prefactor in front of the expected functional form lMinus@SphericalHarmonicY. Without this trick, the functions are too complicated, and FullSimplify can't bring them into the same form on both sides. Note that I have only cheated a little here, because I didn't tell Mathematica what the new eigenvalue is expected to be. It found that out by itself.

The last eigenvalue equation means that if the original spherical harmonic had index m, then the result of lMinus has index m -1. We could now apply this operator to the "maximal state" $Y^{\ell}_{\ell}$ for which I had already shown that the eigenvalue equation $L^2 Y^{\ell}_{\ell} = \ell (\ell +1) Y^{\ell}_{\ell}$ is satisfied. Then it remains to be shown that the same eigenvalue of $L^2$ is obtained with the function created by repeated application of lMinus (one writes it $L_{-}$). The crucial fact that is needed for that proof is this commutation relation:

$$L^2 L_{-} = L_{-} L^2$$

which I show here by applying both orders to an arbitrary test function $f(\theta, \phi)$:

Simplify[
 l2Op[lMinus[f[θ, ϕ]]] == lMinus[l2Op[f[θ, ϕ]]]]

True

With this, I have all the building blocks to deduce that

$$L^2 Y^{m}_{\ell} = \ell (\ell +1) Y^{m}_{\ell}$$

holds. You start by setting $m=\ell$, for which the above is known to be true from earlier. Then apply lMinus on both sides to get the same equation for $m=\ell -1$, and repeat another $2\ell$ times to span the entire angular-momentum ladder of allowed m for the given ℓ.