Disclaimer since I am interested in calculating the answer efficiently there are many steps which I will not justify nor even worry about. For instance I will exchange the order of integration and summation without proving uniform convergence of the infinite sums. Later I will provide some numerical evidence to support the final answer but I will not provide proof.
If someone is interested in proving all the assertions which must be made to make the argument fully rigorous I would be very interested in seeing the result.
To find the unit step functions expansion we start with the fact that its derivative is the dirac delta function. Finding the Legendre polynomial expansion of $\delta(x)$ is quite simple,
$$ \delta(x) = \sum_{l'=0}^\infty d_{l'} P_{l'}(x) $$
$$ \int_{-1}^1 \delta(x) P_l(x) dx = \sum_{l'=0}^\infty \int_{-1}^1 d_{l'} P_{l'}(x) P_l(x) dx$$
$$ P_l(0) = \frac{2}{2l+1} d_l $$
$$ \Rightarrow \underline{\delta(x) = \sum_{l=0}^\infty \dfrac{2l+1}{2} P_l(0) P_l(x)} \qquad \textbf{(1)}$$
We will use the identity,
$$ \int_x^1 P_l(x') dx' = \frac{1-x^2}{l(l+1)} \frac{dP_l(x)}{dx} \qquad \textbf{(2)},$$
to help us integrate $\textbf{(1)}$ in order to get $\Theta(x)$.
$$ \int_x^1 \delta(x') dx' = \sum_{l=0}^\infty \dfrac{2l+1}{2} P_l(0) \int_x^1 P_l(x') dx'$$
$$ \int_x^1 \delta(x') dx' = \dfrac{1}{2} P_0(0) \int_x^1 P_0(x') dx' + \sum_{l=1}^\infty \dfrac{2l+1}{2} P_l(0) \int_x^1 P_l(x') dx'$$
$$ \Theta(1)-\Theta(x) = \left(\dfrac{1}{2} - \dfrac{x}{2}\right) +\sum_{l=1}^\infty \dfrac{2l+1}{2} P_l(0) \frac{1-x^2}{l(l+1)} \frac{dP_l(x)}{dx}$$
$$ \Theta(x) = \underline{ \dfrac{1}{2}+\dfrac{x}{2}-\sum_{l=0}^\infty \dfrac{2l+1}{2} P_l(0) \frac{1-x^2}{l(l+1)} \frac{dP_l(x)}{dx}} \qquad \textbf{(3)}$$
In order to proceed with $\textbf{(3)}$ we need the identity,
$$ (1-x^2)\frac{dP_l(x)}{dx} = -lxP_l(x)+lP_{l-1}(x) \qquad \textbf{(4)}, $$
but first combining $\textbf{(4)}$ with the recurrence relation gives us,
$$ (1-x^2)\frac{dP_l(x)}{dx} = -l \left( \frac{l+1}{2l+1}P_{l+1}(x)+\frac{l}{2l+1}P_{l-1}(x)\right)+lP_{l-1}(x) $$
$$ = \frac{-l(l+1)}{2l+1}P_{l+1}(x) +\frac{l(l+1)}{2l+1}P_{l-1}(x) \qquad \textbf{(5)}. $$
Substituting $\textbf{(5)}$ into $\textbf{(3)}$ we get,
$$ \Theta(x) = \dfrac{1}{2}+\dfrac{x}{2}-\sum_{l=1}^\infty \dfrac{2l+1}{2} P_l(0) \frac{1}{l(l+1)} \left( \frac{-l(l+1)}{2l+1}P_{l+1}(x) +\frac{l(l+1)}{2l+1}P_{l-1}(x)\right) $$
$$ = \dfrac{1}{2}+\dfrac{x}{2}-\sum_{l=1}^\infty \dfrac{1}{2} P_l(0) \left( -P_{l+1}(x) + P_{l-1}(x)\right)$$
$$ = \dfrac{1}{2}+\dfrac{x}{2} + \sum_{l=1}^\infty \dfrac{1}{2} P_l(0) P_{l+1}(x) - \sum_{l=1}^\infty \dfrac{1}{2} P_l(0) P_{l-1}(x) $$
$$ = \dfrac{1}{2}+\dfrac{x}{2} + \sum_{l=2}^\infty \dfrac{1}{2} P_{l-1}(0) P_{l}(x) - \sum_{l=1}^\infty \dfrac{1}{2} P_{l+1}(0) P_{l}(x) \qquad \textbf{(6)}$$
$$ = \dfrac{1}{2}+ \dfrac{1}{4} x+\sum_{l=2}^\infty \dfrac{1}{2}\left(-P_{l+1}(0) + P_{l-1}(0) \right) P_{l}(x) $$
$$ = \dfrac{1}{2}+ \sum_{l=1}^\infty \dfrac{1}{2}\left(-P_{l+1}(0) + P_{l-1}(0) \right) P_{l}(x) $$
So now we have evidence that, $$ \underline{\Theta(x) = 1/2 + \sum_{l=1}^\infty \dfrac{ P_{l-1}(0) - P_{l+1}(0)}{2} P_{l}(x)} \qquad \textbf{(7)}. $$
Because of the lack of rigor I feel that I need to present some numerical evidence that these are the correct expansion coefficients. Below are three plots, two are of the delta function using the $d_l$ coefficients. The last one is of the step funciton using our final result. In each case the curves are color coded with titles corresponding to the number of nonzero terms from the expansion that were used.
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
Basically you're asking about the Fourier expansion of $$P_n(\cos\theta)=\frac{a_{n,0}}{2}+\sum_{k=1}^n a_{n,k} T_k(\cos\theta)=\frac{a_{n,0}}{2}+\sum_{k=1}^n a_{n,k}\cos k\theta.$$
A relatively easy way is to consider the generating function $$\sum_{n=0}^\infty P_n(\cos\theta)t^n=(1-2t\cos\theta+t^2)^{-1/2}=(1-te^{i\theta})^{-1/2}(1-te^{-i\theta})^{-1/2}$$ and use the binomial series $$(1-z)^{-1/2}=\sum_{n=0}^\infty b_n z^n,\qquad b_n=\frac{(2n-1)!!}{(2n)!!}$$ (here $0!!=(-1)!!=1$). Computing the product, we find $$\sum_{n=0}^\infty P_n(\cos\theta)t^n=\sum_{j,k=0}^\infty b_j b_k t^{j+k}e^{i(j-k)\theta}=\sum_{n=0}^\infty t^n\sum_{k=0}^n b_k b_{n-k}e^{i(n-2k)\theta}.$$
Hence, "renaming" the indices, $$a_{n,k}=\begin{cases}\hfill 0,\hfill&n-k\text{ is odd}\\2b_{(n-k)/2}b_{(n+k)/2},&n-k\text{ is even}\end{cases}.$$