Is there a closed-form expression for a monomial $x^m$ in terms of a sum of Legendre polynomials $P_n(x)$?
$$ x^m = \sum_n a_n P_n(x) $$
How can I determine the coefficients $a_n$ in general?
According to this question (and here), the answer seems to be something along the lines of
$$ a_n = \frac{2n + 1}{2} \int_{-1}^1 x^m P_n(x) \,\mathrm{d}x $$
How can I use orthogonality to prove this? Can the integral be expressed in closed form?
Start from the definition of Legendre polynomials through its generating function:
$$\frac{1}{\sqrt{1-2tx+t^2}} = \sum_{n=0}^\infty P_n(x) t^n$$
Introduce variable $s$ such that
$$1 - ts = \sqrt{1-2tx + t^2} \iff x = \frac{1+t^2 - (1-ts)^2}{2t} = s + \frac{t}{2}(1-s^2) $$ We have $\displaystyle\;\frac{dx}{\sqrt{1-2tx+t^2}} = \frac{1+t(-s)}{1-ts} ds = ds\;$ and hence
$$\sum_{n=0}^\infty t^n \int_{-1}^1 x^m P_n(x) dx = \int_{-1}^1 \frac{x^m dx}{\sqrt{1-2tx+t^2}} = \int_{-1}^1 \left(s + \frac{t}{2}(1-s^2)\right)^m ds $$ By comparing the coefficients of $t^n$ on both sides, we find the coefficients defined by $$a_{m,n} \stackrel{def}{=} \int_{-1}^1 x^m P_n(x) dx$$ vanishes whenever $n > m$ and $m, n$ has different parity. When $n = m - 2k$ for some integer $k$, we have
$$\begin{align} a_{m,m-2k} &= 2^{2k-m}\binom{m}{2k}\int_{-1}^1 s^{2k} (1-s^2)^{m-2k} ds = 2^{2k-m}\binom{m}{2k}\int_0^1 y^{k-\frac12} (1-y)^{m-2k} dy\\ &= 2^{2k-m}\binom{m}{2k}\frac{\Gamma(k+\frac12)\Gamma(m-2k+1)}{\Gamma(m-k+\frac32)} = \frac{m!}{2^{k-1}k!(2m-2k+1)!!} \end{align} $$ Together with the orthogonality relation $$\frac{2\ell+1}{2}\int_{-1}^{1} P_\ell(x) P_\ell'(x) = \begin{cases}1,& \ell = \ell'\\0, & \text{ otherwise }\end{cases}$$
We find for integer $m$,
$$x^m = \sum_{k=0}^{\lfloor m/2\rfloor} \left(\frac{2m-4k+1}{2}\right) a_{m,m-2k} P_{m-2k}(x) = \sum_{k=0}^{\lfloor m/2\rfloor}\frac{m!(2m-4k+1)}{2^k k!(2m-2k+1)!!}P_{m-2k}(x) $$
