As noted, you cannot use just any inner product.
Suppose $q(t)$ is a factor of $p_n$ such that $q(t)\ge0$ for $t$ in $[a,b]$. Let $r(t)=p(t)/q(t)$. Then the product $rp_n$ is non-negative on $[a,b]$ and is not identically zero, and therefore $\langle 1,rp_n\rangle>0$.
On the other hand, $p_n$ is orthogonal to all polynomials of degree less than $n$, and therefore if the degree of $q$ is greater than zero, $\langle r,p_n\rangle=0$. (Remark: this shows we cannot use just any inner product on polynomials, we need $\langle tp,q\rangle=\langle p,tq\rangle$
for all $p$ and $q$, that is, multiplication by $t$ must be self adjoint.)
Since $\langle1,rp_n\rangle = \langle r,p_n\rangle$, this is a contradiction. We conclude that $p_n$ has no non-constant non-negative factor. In particular
$(t-\lambda)^2$ cannot be a factor
and therefore $p_n$ has only simple zeros.
A minor variant of this argument also shows that the zeros of $p_n$ are all real.
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.
Best Answer
The following is an adapted translation of Wikipedia.fr, Suite_de_polynômes_orthogonaux#Relation_de_récurrence.
Let us prove that $$p_{n+1}=(x+b_n)p_n-c_np_{n-1}$$ where $$b_n=k_{n+1}- k_n,\qquad c_n=\frac{h_n}{h_{n-1}},$$ $k_j,h_j$ being given by $$p_j(x)=x^j+k_jx^{j-1}+\cdots,\quad h_j=\langle p_j,p_j \rangle.$$
(By convention, $c_0=p_{-1}=k_0=0.$)
The polynomial $(x+b_n)p_n–p_{n+1}$ has a degree $<n$ hence may be written $\sum_{j=0}^{n-1}\mu_{n,j}p_j,$ with $h_j\mu_{n,j}=\langle (x+b_n)p_n-p_{n+1},p_j\rangle=\langle xp_n,p_j\rangle$ (since for $j<n,$ $p_j$ is orthogonal to $p_n$ and to $p_{n+1}$).
Moreover, $\langle xp_n,p_j\rangle=\langle p_n,xp_j\rangle.$
We may conclude: $(x+b_n)p_n-p_{n+1}= c_np_{n-1},$ where $c_n=\mu_{n,n-1}=\frac{h_n}{h_{n-1}}.$