Because $\mathbb{Z}_p^*$ is cyclic, the subgroup $K$ of index $2$ is unique and consists of non-zero squares. We get the non-zero squares twice ($x^2=(-x)^2$), so
$$
2a+1=\sum_{k=0}^{p-1}\zeta^{k^2}.
$$
Therefore
$$
(2a+1)^2=\sum_{k=0}^{p-1}\sum_{\ell=0}^{p-1}\zeta^{k^2+\ell^2}.
$$
Let us try and determine how many times a given exponent $x\in\mathbb{F}_p$ appears here.
Denote that number by $N(x)$. It is the number of solutions of the equation
$k^2+\ell^2=x$ with $k,\ell$ ranging over all of $\mathbb{F}_p$.
It will turn out that the value of $N(x)$ depends on whether $p\equiv 1\pmod4$ or
$p\equiv 3\pmod4$. In the former case we have $N(0)=2p-1$, and $N(x)=p-1$ when $x\neq0$.
In the latter case we have $N(0)=1$ and $N(x)=p+1$. Furthermore, we know that
$$
S=\sum_{n=0}^{p-1}\zeta^n=0,
$$
as this is a full period of a geometric sum repeating after $p$ terms.
Therefore we get
$$
(2a+1)^2=\begin{cases}p\zeta^0+(p-1)S=p,&\text{if }\ p\equiv1\pmod4\cr
-p+(p+1)S=-p,&\text{if }\ p\equiv3\pmod4.\cr\end{cases}
$$
As $(2a+1)^2=4a^2+4a+1=4(a^2+a)+1$ your claim follows from this.
I'm fairly sure that the proof for the stated result about the numbers $N(x)$ is in e.g. Ireland & Rosen. Trying to recall how it goes. Anyway, this is a basic result about quadratic Gauss sums.
If $a\in\mathbb{F}_p^*$ we clearly have $N(x)=N(a^2x)$, because $(k,\ell)\mapsto (ak,a\ell)$ is a bijection from one set of solutions to the other. Therefore $N(x)=N(kx)$
for all $k\in K$, and we are left with three unknowns $N(0)$, $N(k), k\in K$ and $N(k'),k'\notin K$. We select $1$ to represent elements of $K$ and let $q\notin K$ be a fixed quadratic non-residue.
Assume first that $p\equiv1\pmod4$. Then we know that $-1\in K$, and $-1=i^2$ for some $i\in\mathbb{F}_p$. We have $k^2+\ell^2=0$ if and only if $k=\pm i\ell$. Therefore in this case we get $N(0)=2p-1$ - one solution is $(k,\ell)=(0,0)$ and the other $2(p-1)$ are $(\pm i\ell, \ell)$ with $\ell\neq0$. Let us then determine $N(1)$. If we replace $\ell$ with $i\ell$ we might as well look at the number of solutions of $k^2-\ell^2=1$ or $(k-\ell)(k+\ell)=1$. Here we can let $k-\ell=a$ be an arbitrary non-zero element $a\in\mathbb{F}_p^*$ and we can then solve $k+\ell=1/a$. This linear system determines the pair $(k,\ell)$ uniquely, so there are $p-1$ such pairs. Therefore $N(1)=p-1$, and $N(k)=p-1$ for all $k\in K$. As we have $p^2$ pairs $(k,\ell)$ altogether, we have
$$
p^2=N(0)+|K| N(1)+|K| N(q),
$$
we can solve for $N(q)$ and find that $N(q)=p-1$ as well.
In the case $p\equiv 3\pmod4$ we know that $-1\notin K$. Therefore $k^2+\ell^2=0$ only, if
$k=\ell=0$. So in this case $N(0)=1$. One of seeing that $N(x)=p+1$ for all $x\neq 0$ is to observe that in this case $k^2+\ell^2$ is the norm of the element $k+i\ell$ from $\mathbb{F}_{p^2}^*$. Here $i^2=-1$ and $\mathbb{F}_{p^2}=\mathbb{F}_p[i]$. The cyclicity of the group $\mathbb{F}_{p^2}^*$ together with the fact that the norm is just $(p+1)^{th}$ power then gives the claim. Another way would be to use the quadratic character. I'm fairly sure that there is a similar more elementary argument like the one in $p\equiv1\pmod4$ case, but I can't reproduce it now.