You have to multiply $f_n$ by $(1-x^2)^2$ to get an orthogonal family of polynomials. In general, if $\{P_n\}_{n\ge 1}$ is the system of orthogonal polynomials with respect to the weight $W(x)$ then $k-$ the derivatives form an orthogonal system with respect to the weight $W^k(x)a(x)$ where $b(x)$ comes from the differential equation for orthogonal polynomials:
$$a(x)y''(x)+b(x)y'+\lambda_n y=0.$$
In this particular case, Legandre polynomials arise as solutions of the differential equation
$$(1-x^2)y''(x)-2xy'+n(n+1)y=0.$$
and $W=1$ which makes $4- $th derivatives to be orthogonal with respect to the weight $(1-x^2)^4.$ In other words $(1-x^2)^2f_n$ are going to be orthogonal with respect to the weight $1.$
One way to prove the statement above is to use the given differential equation together with integration by parts.
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 Legendre Polynomials are orthogonal with respect to the inner product $\langle f,g\rangle=\int_\mathbb R w(x)f(x) g(x)\, dx$ for a certain weight function $w(x).$