The problem is to determine the expansion of the unit step function in terms of Legendre polynomials on the interval $[-1,1]$.
Here the Legendre polynomials are the family of orthogonal polynomials on the interval $[-1,1]$ with the orthogonality relation,
$$ \int_{-1}^1 P_l(x) P_{l'}(x) dx = \frac{2}{2l+1} \delta_{l,l'}$$
Where,
$$ P_0(x)=1, \quad P_1(x) = x, \quad P_2(x) = (3x^2-1)/2, \quad \dots$$
With the recurrence relation,
$$\qquad P_{l+1} = \frac{2l+1}{l+1}xP_l(x)-\frac{l}{l+1}P_{l-1}(x)$$
The unit step function is, $$\Theta(x) = \begin{cases} 1 \qquad x \geq 0 \\ 0 \qquad x <0 \end{cases}.$$
We wish to find coefficients $a_l$ such that,
$$\boxed{ \Theta(x) = \sum_{l=0}^\infty a_l P_l(x) }$$
Best Answer
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.