We wish to evaluate:
$$\int_0^{2\pi} P_l^m(\cos\theta)P_{l-1}^m(\cos\theta)\,\text{d}\theta $$
We use the transformation:
$$x = \cos\theta\quad \frac{dx}{d\theta} = -\sin\theta \quad d\theta = \frac{-1}{\sqrt{1-x^2}}dx$$ to obtain$$\int_1^1 \frac{-P_l^m(x)P_{l-1}^m(x)}{\sqrt{1-x^2}}dx = 0$$
Since you're integrating from 1 to 1.
Weird? Yes I thought so too. Seeing as the integral comes from spherical harmoics, I assume you are messing up your limits of integration (should be from $0$ to $\pi$, since that is how the zenith angle is normally defined).
Then, I think you should be working with:
$$\int_0^{\pi} P_l^m(\cos\theta)P_{l-1}^m(\cos\theta)\,\text{d}\theta $$
which after the transformation becomes:
$$I = \int_{-1}^1 \frac{P_l^m(x)P_{l-1}^m(x)}{\sqrt{1-x^2}}dx$$
The Legendre Polynomials satisfy the relationship:
$$P_l^m(-x) = (-1)^{l+m}P_l^m(x)$$
This means that if $l+m$ is even (odd), the Legendre Polynomial is also even (odd). Here we have the two numbers:
$l+m$ and $l+m-1$ One of these two has to be even and the other odd. So one of the polynomials is even and the other is odd. Since $\frac{1}{\sqrt{1-x^2}}$ is even, the integrand as a whole is odd. Thus the integral is also zero.
$$I = 0$$
Reference for parity relationship.
Best Answer
Let be $$\mathcal A_{k l}^m = \int_{-1}^1 P_k^m \left({x}\right) P_l^m \left({x}\right) \, \mathrm d x $$ where the associated Legendre functions are $$P^m_l \left({x}\right) = \left({1 - x^2}\right)^{m/2} \dfrac {\mathrm d^m P_l \left({x}\right)} {\mathrm d x^m}={\left({1 - x^2}\right)^m \frac {\mathrm d^{k + m} } {\mathrm d x^{k + m} } \left({x^2 - 1}\right)^k }\qquad 0 \le m \le l$$ Thus we have $$\mathcal A_{k l}^m = \frac 1 {2^{k + l} k! l!} \int_{-1}^1 \left({\left({1 - x^2}\right)^m \frac {\mathrm d^{k + m} } {\mathrm d x^{k + m} } \left({x^2 - 1}\right)^k}\right) \left({\frac{\mathrm d^{l + m} } {\mathrm d x^{l + m} } \left({x^2 - 1}\right)^l }\right) \, \mathrm d x$$ where $k$ and $l$ occur symmetrically. Let $l \ge k$.
We can integrate by parts $l + m$ times $$\int_{-1}^1 u v' \mathrm d x = \left.{u v}\right|_{-1}^1 - \int_{-1}^1 v u' \ \mathrm d x$$ where $u = \left({1 - x^2}\right)^m \frac {\mathrm d^{k + m} } {\mathrm d x^{k + m} } \left({x^2 - 1}\right)^k$ and $v' = \frac {\mathrm d^{l + m}} {\mathrm d x^{l + m}} \left({x^2 - 1}\right)^l$
For each of the first $m$ integrations by parts, $u$ in the $\left.{uv}\right|_{-1}^1$ term contains the factor $\left({1 - x^2}\right)$, so the term vanishes.
For each of the remaining $l$ integrations, $v$ in that term contains the factor $\left({x^2 - 1}\right)$ so the term also vanishes.
This means $$ \mathcal A_{k l}^m = \frac {\left({-1}\right)^{l + m} } {2^{k + l} k! l!} \int_{-1}^1 \left({x^2 - 1}\right)^l \frac {\mathrm d^{l + m} } {\mathrm d x^{l + m} } \left({\left({1 - x^2}\right)^m \frac {\mathrm d^{k + m} } {\mathrm d x^{k + m} } \left({x^2 - 1}\right)^k }\right) \, \mathrm d x$$
Expanding the second factor using Leibniz's Rule $$\frac {\mathrm d^{l + m} } {\mathrm d x^{l + m} } \left({1 - x^2}\right)^m \frac {\mathrm d^{k + m} } {\mathrm d x^{k + m} } \left({x^2 - 1}\right)^k= \sum_{r \mathop = 0}^{l + m} \binom {l + m} r \frac {\mathrm d^r} {\mathrm d x^r} \left({1 - x^2}\right)^m \frac {\mathrm d^{l + k + 2 m - r} } {\mathrm d x^{l + k + 2 m - r} } \left({x^2 - 1}\right)^k$$ the leftmost derivative in the sum is non-zero only when $r \le 2 m$ (remembering that $m \le l$) and the other derivative is non-zero only when $k + l + 2 m - r \le 2 k$, that is, when $r \ge 2 m + l - k$.
Because $l \ge k$, these two conditions imply that the only non-zero term in the sum occurs when $r = 2 m$ and $l = k$.
Thus $$\mathcal A_{k l}^m = \left({-1}\right)^l \delta_{k l} \frac {\left({-1}\right)^{l + m} } {2^{2 l} \left({l!}\right)^2} \binom {l + m} {2 m} \int_{-1}^1 \left({x^2 - 1}\right)^l \frac {\mathrm d^{2 m} } {\mathrm d x^{2 m} } \left({1 - x^2}\right)^m \frac {\mathrm d^{2 l} } {\mathrm d x^{2 l} } \left({1 - x^2}\right)^l \mathrm d x$$ where $\delta_{k l}$ is the Kronecker Delta.
The factor $(-1)^l$ at the front of $\mathcal A_{k l}^m$ comes from switching the sign of $x^2 - 1$ inside $\left({x^2 - 1}\right)^l$.
To evaluate the differentiated factors, expand $\left({1 - x^2}\right)^k$ using the Binomial Theorem $$\left({1 - x^2}\right)^k = \sum_{j \mathop = 0}^k \binom k j \left({-1}\right)^{k-j} x^{2 \left({k-j}\right)}$$
The only term that survives differentiation $2^k$ times is the $x^{2 k}$ term, which after differentiation gives $$\left({-1}\right)^k \binom k 0 2 k! = \left({-1}\right)^k \left({2k}\right)!$$
Therefore $$\mathcal A_{k l}^m = \left({-1}\right)^l \delta_{k l} \frac 1 {2^{2 l} \left({l!}\right)^2} \frac {\left({2 l}\right)! \left({l + m}\right)!} {\left({l - m}\right)!} \int_{-1}^1 \left({x^2 - 1}\right)^l \ \mathrm d x =\left({-1}\right)^l \delta_{k l} \frac 1 {2^{2 l} \left({l!}\right)^2} \frac {\left({2 l}\right)! \left({l + m}\right)!} {\left({l - m}\right)!} \mathcal B_l$$
The integral $$\mathcal B_l=\int_{-1}^1 \left({x^2 - 1}\right)^l \ \mathrm d x$$ can be evaluated by a change of variable $x = \cos \theta$
Thus $$\mathcal B_l= \left({-1}\right)^{l + 1} \int_\pi^0 \left({\sin \theta }\right)^{2 l + 1} \, \mathrm d \theta=\left({-1}\right)^{l} \int_0^\pi \left({\sin \theta }\right)^{2 l + 1} \, \mathrm d \theta$$
Integration of $$\frac {\mathrm d \left({\sin^{n - 1} \theta \cos \theta}\right)} {\mathrm d \theta} = \left({n-1}\right) \sin^{n-2} \theta - n \sin^n \theta$$ gives $$\int_0^\pi \sin^n \theta \, \mathrm d \theta = \frac {\left.{-\sin^{n - 1} \theta \cos \theta}\right|_0^\pi} n + \frac {\left({n - 1}\right)} n \int_0^\pi \sin^{n - 2} \theta \, \mathrm d \theta= \frac {\left({n - 1}\right)} n \int_0^\pi \sin ^{n - 2} \theta \, \mathrm d \theta$$
since $\displaystyle \left.{-\sin^{n-1} \theta \cos \theta}\right|_0^\pi = 0$ for $n > 1$. Applying this result to $\int_0^\pi \left({\sin \theta }\right)^{2 l + 1} \, \mathrm d \theta$ and changing the variable back to $x$ yields:
$$\int_{-1}^1 \left({x^2 - 1}\right)^l \, \mathrm d x = - \frac {2 l} {2 l + 1} \int_{-1}^1 \left({x^2 - 1}\right)^{l - 1} \, \mathrm d x\quad\text{for}\;l \ge 1$$
Using this recursively $$\displaystyle \int_{-1}^1 \left({x^2 - 1}\right)^l \, \mathrm d x = \left({-1}\right)^l \left({\frac {2 l} {2 l + 1} \frac {2 \left({l - 1}\right)} {2 l - 1} \frac {2 \left({l - 2}\right)} {2 l - 3} \cdots \frac 2 3}\right) \int_{-1}^1 \, \mathrm d x$$
and noting that
$$ \frac {2 l} {2 l + 1} \frac {2 \left({l-1}\right) } {2 l - 1} \frac {2 \left({l - 2}\right)} {2 l - 3} \cdots \frac 2 3= \frac {2^l l!} {\left({2 l + 1}\right) \left({2 l - 1}\right) \left({2 l - 3}\right) \cdots 3} = \frac{2^l l!} {\frac {\left({2 l + 1}\right)!} {2^l l!} } = \frac {2^{2 l} \left({l!}\right)^2} {\left({2 l + 1}\right)!} $$
it follows that $$\mathcal B_l=\int_{-1}^1 \left({x^2 - 1}\right)^l \mathrm d x = \ \left({-1}\right)^l \ \frac{2^{2l+1} \left({l!}\right)^2} {\left({2l+1}\right) !}$$
Therefore we have $$\mathcal A_{k l}^m =\left({-1}\right)^l \delta_{k l} \frac 1 {2^{2 l} \left({l!}\right)^2} \frac {\left({2 l}\right)! \left({l + m}\right)!} {\left({l - m}\right)!} \mathcal B_l= \delta _{k l} \frac 2 {2 l + 1} \frac {\left({l + m}\right) !} {\left({l - m}\right)!}$$ that is