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.
As far as roots where the given function changes sign on both sides, they may be found using bracketing methods. In contrast, roots which do not have this property are local extrema.
For smooth functions, these cases fall into even and odd roots as you've mentioned, but not all functions are smooth (e.g. $\sqrt[3]x$ or $|x|$) and thus cannot be classified into either case.
After having read a fairly large volume of literature concerning the case in which roots may be bracketed, I have not found any consistent terminology for these two classes of roots.
If I had to give them names however, I would classify them as bracketable roots and local extrema roots due to the context in which I would differentiate the two, which would be in relation to bracketing methods (root-finding algorithms) versus optimization algorithms.
From a more numerical analysis standpoint, I've never heard the claim that roots which do not experience a sign change should not be classified as a root and certainly have seen discussions concerning the "local extrema roots".
Best Answer
To resolve the second question, note first that the Legendre polynomials are odd functions for odd order (0 then is one root of the polynomial), and even functions for even order. Thus, with regards to solubility in terms of radicals, you should be able to derive (possibly complicated!) radical expressions at least up until $P_9(x)$. To use that as an example, note that
$$\frac{P_9(\sqrt{x})}{\sqrt{x}}$$
is a quartic; thus, one can use the quartic formula to derive explicit expressions for its roots, and then you can easily derive the roots of $P_9(x)$ .
$P_{10}(x)$ is where your trouble starts. If we take a look at the polynomial
$$P_{10}(\sqrt{x})$$
we have a quintic to contend with. I'll skip the relatively tedious details, but you can verify that its Galois group is not a solvable group, and thus the solution cannot be expressed in terms of radicals (you can use theta or hypergeometric functions, though).
So, not much hope in the symbolic front. In the numeric front, things are much easier. The slickest way of getting accurate values of the roots of the Legendre polynomial is to use the Jacobi matrix in my previous answer. Since there exist stable and efficient algorithms (e.g. QR algorithm or divide-and-conquer) for the symmetric eigenproblem (in LAPACK, for instance), and things can be set such that only eigenvalues are returned, you have a good way of generating good approximate values of Legendre polynomial roots. (In the context of Gaussian quadrature, where the roots of orthogonal polynomials play a pivotal role, the scheme is referred to as the Golub-Welsch algorithm.)
Alternatively, as I mentioned in the comments, there exist asymptotic approximations for the roots, which can then be subsequently polished with a few applications of Newton-Raphson. One such asymptotic approximation is due to Francesco Tricomi. Letting $\xi_{n,k}$ be the $k$-th root of $P_n(x)$, ordered in decreasing order, we have
$$\xi_{n,k}\approx\left(1-\frac1{8n^2}+\frac1{8n^3}\right)\cos\left(\pi\frac{4k-1}{4n+2}\right)$$
and $O(n^{-4})$ and further terms are omitted. Other asymptotic approximations due to Luigi Gatteschi use roots of Bessel functions, but I won't say more about those.