Thanks to @metamorphy and his insights about similiarity to Parseval's theorem, the answer is achieved, in essence, by expanding $\sigma$ as a Fourier-Legendre series, which I detail below. (Not in its full detail, though, mostly recollecting the high points.)
$\newcommand{\dd}{\mathrm{d}}$
$\newcommand{\para}[1]{\left( #1 \right)}$
$\newcommand{\encla}[1]{\langle #1 \rangle}$
We begin by first expanding $\sigma$ as a Fourier-Legendre series; recall this takes the form
$$
\sigma(x) = \sum_{n=0}^\infty c_n P_n(x) \text{ where } c_n = \frac{2n+1}{2} \int_{-1}^1 \sigma(x)P_n(x) \, \dd x
$$
We begin calculating the constants $c_n$. Note that, due to the piecewise nature of $\sigma$,
$$
c_n = \frac{2n+1}{2} \para{ \int_{-1}^0 -P_n(x) \, \dd x + \int_0^1 P_n(x) \, \dd x }
$$
We utilize that $P_n$ is even if $n$ is even here, and similar for $n$ odd. This will lead us to conclude $c_n = 0$ for $n$ even, and for $n$ odd,
$$
c_{\text{n, odd}} = (2n+1) \int_0^1 P_n(x) \, \dd x
$$
To calculate what remains, consider the recursion
$$
\frac{P_{n+1}' (x) - P_{n-1}'(x)}{2n+1} = P_n(x)
$$
Integrate both sides of this over $x \in (0,1)$. Then it readily follows that
$$
\int_0^1 P_n(x) \, \dd x =\left. \frac{1}{2n+1} \Big( P_{n+1}(x) - P_{n-1}(x) \Big) \right|_{x=0}^1
$$
Due to the normalization process, $P_n(1) = 1$ always. Meanwhile,
$$
P_n(0) = \begin{cases} (-1)^{n/2} \frac{(n-1)!!}{n!!} & n \equiv 0 \pmod 2 \\ 0 & n \equiv 1 \pmod 2 \end{cases}
$$
Using this, the assumption $n$ is odd, and a multitude of algebraic manipulations (mainly factoring and properties of double factorials) we conclude that, if $n = 2k+1$,
$$
\int_0^1 P_{2k+1} \, \dd x = (-1)^k \frac{(2k-1)!!}{(2k+2)!!}
$$
Thus, if $n$ is even, $c_n = 0$; if $n = 2k+1$ is odd, then
$$
c_{2k+1} = (4k+3) (-1)^k \frac{(2k-1)!!}{(2k+2)!!}
$$
This means that we can simply sum over the odd indices in our series for $\sigma$, switch our dummy variable from $k$ to $n$, and conclude
$$
\sigma(x) = \sum_{n=0}^\infty (4n+3) (-1)^n \frac{(2n-1)!!}{(2n+2)!!} P_{2n+1}(x)
$$
We then choose to square this representation of $\sigma$, using the Cauchy product:
$$
\sigma^2(x) = \left( \sum_{n=0}^\infty c_n P_n(x) \right)^2 = \sum_{n=0}^\infty \sum_{m=0}^n c_{2n+1} c_{2(n-m)+1} P_{2n+1}(x) P_{2(n-m)+1}(x)
$$
We now integrate throughout with respect to $x \in (-1,1)$, and assume there is no issue in interchanging summation and integration here. Then we see
$$
\int_{-1}^1 \sigma^2(x) \, \dd x = \sum_{n=0}^\infty \sum_{m=0}^n c_{2n+1} c_{2(n-m)+1} \int_{-1}^1 P_{2n+1}(x) P_{2(n-m)+1}(x) \, \dd x
$$
We recall the orthogonality relation, $\encla{P_p,P_q} = \frac{2}{2p+1} \delta_{p,q}$. This causes much simplification, yielding
$$
\int_{-1}^1 \sigma^2(x) \, \dd x = \sum_{n=0}^\infty c_{2n+1}^2 \frac{2}{4n+3}
$$
We bring back in our values for the $c_{2n+1}$, and square them. A factor of $4n+3$ cancels in this process. We can then bring the $2$ outside of the sum, to conclude with our desired result:
$$
\int_{-1}^1 \sigma^2(x) \; \dd x = 2 \sum_{n=0}^\infty (4n+3) \para{ \frac{(2n-1)!!}{(2n+2)!!} }^2
$$
Best Answer
Note that the spherical Bessel functions admit the integral representation $$j_n(z)=\frac{(-\mathrm i)^n}{2}\int_0^\pi \mathrm e^{\mathrm iz\cos\theta}P_n(\cos\theta)\sin\theta~\mathrm d\theta$$ Which, by the Leibniz rule, means $$j_n^{(n)}(z)=\frac{1}{2}\int_0^\pi \mathrm e^{\mathrm iz\cos\theta}\cos(\theta)^nP_n(\cos\theta)\sin\theta~\mathrm d\theta \\ j_n^{(n)}(0)=\frac{1}{2}\int_0^\pi \cos(\theta)^n P_n(\cos\theta)\sin\theta~\mathrm d\theta$$ Which, by a change of variable $\cos\theta=x$ yields $$j_n^{(n)}(0)=\frac{1}{2}\int_{-1}^1 x^nP_n(x)\mathrm dx$$
Which you already know!