Following Robert Israel's answer, we also scale everything to $[-1,1]$ (thus multiplying the result by 2). As he mentions, the optimal polynomial is always a square of some other polynomial, $p_{2n}=p_{2n+1}=q_n^2$, and $q_n$ is either even or odd (see Lemma below). So we are left to find the minimal $L_2[-1,1]$-norm of an odd/even polynomial $q_n$ such that $\deg q_n\leq n$ and $q_n(\pm1)=(\pm1)^n$. In other words (recall the division by 2 in the first line!), $a_{2n}=a_{2n+1}=d_n^2/2$, where $d_n$ is the distance from $0$ to the hyperplane defined by $q(1)=1$ in the space of all odd/even polynomials of degree at most $n$ with $L_2[-1,1]$-norm.

Now, this hyperplane is the affine hull of the Legendre polynomials $ P_i(x)=\frac1{2^ii!}\frac{d^i}{dx^i}(x^2-1)^i$, where $i$ has the same parity as $n$ (since we know that $ P_i(1)=1$ and $P_i(-1)=(-1)^i$). Next, by $\| P_i\|^2=\frac{2}{2i+1}$, our distance is
$$
\left(\sum_{j\leq n/2}\| P_{n-2j}\|^{-2}\right)^{-1/2}
=\left(\sum_{j\leq n/2}\frac{2(n-2j)+1}{2}\right)^{-1/2}=\sqrt{\frac4{(n+1)(n+2)}},
$$
attained at
$$
q_n(x)=\left(\sum_{i\leq n/2}\frac{2}{2(n-2i)+1}\right)^{-1}
\sum_{i\leq n/2}\frac{2P_{n-2i}(x)}{2(n-2i)+1}.
$$
Thus the answer for the initial question is $a_{2n}=a_{2n+1}=\frac2{(n+1)(n+2)}$ and $p_{2n}=p_{2n+1}=q_n^2$.

**Lemma.** For every $n$, one of optimal polynomials on $[-1,1]$ is a square of an odd or even polynomial.

*Proof.* Let $r(x)$ be any polynomial.which is nonnegative on $[-1,1]$ with $r(\pm 1)=1$. We will replace it by some other polynomial which has the required form, has the same (or less) degree and the same values at $\pm 1$, is also nonnegative on $[-1,1]$, and is not worse in the integral sense. Due to the compactness argument, this yields the required result.

Firstly, replacing $r(x)$ by $\frac12(r(x)+r(-x))$, we may assume that $r$ is even (and thus has an even degree $2n$). Let $\pm c_1,\dots,\pm c_{n}$ be all complex roots of $r$ (regarding multiplicities); then
$$
r(x)=\prod_{j=1}^n\frac{x^2-c_j^2}{1-c_j^2},
$$
due to $r(\pm 1)=1$.

Now, for all $c_j\notin[-1,1]$ we simultaneously perform the following procedure.

(a) If $c_j$ is real, then we replace $\pm c_j$ by $\pm x_j=0$. Notice that
$$
\frac{|x^2-c_j^2|}{|1-c_j^2|}\geq 1\geq \frac{|x^2-0^2|}{|1-0^2|}.
$$
for all $x\in[-1,1]$.

(b) If $c_j$ is non-real, then we choose $x_j\in[-1,1]$ such that $\frac{|c_j-1|}{|c_j+1|}=\frac{|x_j-1|}{|x_j+1|}$. Notice that all complex $z$ with $\frac{|c_j-z|}{|x_j-z|}=\frac{|c_j-1|}{|x_j-1|}$ form a circle passing through $-1$ and $1$, and the segment $[-1,1]$ is inside this circle. Therefore, for every $x\in[-1,1]$ we have
$$
\frac{|c_j-x|}{|x_j-x|}\geq\frac{|c_j-1|}{|x_j-1|},
$$
thus
$$
\frac{|x^2-c_j^2|}{|1-c_j^2|}
=\frac{|c_j-x|}{|c_j-1|}\cdot\frac{|c_j+x|}{|c_j+1|}
\geq \frac{|x_j-x|}{|x_j-1|}\cdot\frac{|x_j+x|}{|x_j+1|}
=\frac{|x^2-x_j^2|}{|1-x_j^2|}.
$$
So, we replace $\pm c_j$ and $\pm\bar c_j$ by $\pm x_j$ and $\pm x_j$ (or simply $\pm c_j$ by $\pm x_j$ if $c_j$ is purely imaginary).

After this procedure has been applied, we obtain a new polynomial whose roots are in $[-1,1]$ and have even multiplicities, and its values at $\pm1$ are equal to $1$. So it is a square of some polynomial which is even/odd (since the roots are still split into pairs of opposite numbers). On the other hand, its values at every $x\in[-1,1]$ do not exceed the values of $r$ at the same points, as was showed above. So the obtained polynomial is not worse in the integral sense, as required. The lemma is proved.