Is $-\det\big[\big(\frac{i^2+j^2}p\big)\big]_{1\le i,j\le (p-1)/2}$ always a square for each prime $p\equiv 3\pmod 4$? Let $p$ be an odd prime and let $S_p$ denote the determinant
$$\det\left[\left(\frac{i^2+j^2}p\right)\right]_{1\le i,j\le (p-1)/2}$$
with $(\frac{\cdot}p)$ the Legendre symbol. By Theorem 1.2 of my paper arXiv:1308.2900 available from http://arxiv.org/abs/1308.2900, $-S_p$ is a quadratic residue modulo $p$. Here I ask a further question.
QUESTION. Is it true that for each prime $p\equiv3\pmod4$ the number $-S_p$ is always a positive square divisible by $2^{(p-3)/2}$? 
Define $a_p=\sqrt{-S_p}/2^{(p-3)/4}$ for any prime $p\equiv3\pmod4$. Then
\begin{gather*}a_3=a_7=a_{11}=1,\ a_{19}=2,\ a_{23}=1,\ a_{31}=29,\ a_{43}=254,
\\a_{47}=367,\ a_{59}=9743,\ a_{67}=305092,\ a_{71}=29,\ a_{79}=1916927.
\end{gather*}
I have computed the values of $a_p$ for all primes $p\equiv3\pmod4$ with $p<2000$. Based on the numerical data, I conjecture that the above question has an affirmative answer but I'm unable to prove this. 
Any ideas towards the solution? Your comments are welcome!
 A: This is not an answer, but a reduction to supposedly simpler problem. (edited with more details)
Using quadratic Gauss sums, we can express Legendre symbol $\left(\frac{i^2+j^2}p\right)$
as
\begin{split}
\left(\frac{i^2+j^2}p\right) &= \frac{1}{I\sqrt{p}}\sum_{n=0}^{p-1} \zeta^{(i^2+j^2)n^2} \\
&= \frac{1}{I\sqrt{p}}\left(1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2}\right),\qquad(\star)
\end{split}
where $\zeta := \exp(\frac{2\pi I}p)$ and $I$ is the imaginary unit.
From $(\star)$ it follows that
$$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{1}{I\sqrt{p}}(J + 2B^2),$$
where $J$ is the matrix of all ones and $B:=[\zeta^{i^2j^2}]_{i,j=1}^{(p-1)/2}$. This enables application of pavl0's approach in evaluating the determinant in order to show that it is a square.

Alternatively, $(\star)$ allows us to extract the square root from the determinant without computing its value. To do so, we need to embed "$1+$" into the quadratic sum in $(\star)$ by introducing a parameter $\alpha$:
$$1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2} = 2\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha).$$
Since $2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2} = I\sqrt{p} - 1$ for any $i\not\equiv 0\pmod{p}$, we need $\alpha$ to satisfy the quadratic equation:
$$1 = 2(I\sqrt{p} - 1)\alpha + (p-1)\alpha^2.$$
So, we can set $\alpha := ((p^{1/4}-1)(1+Ip^{1/4}))^{-1}$, which is a root of this equation. Then $(\star)$ turns into
$$\left(\frac{i^2+j^2}p\right) = \frac{2}{I\sqrt{p}}\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha),$$
implying that
$$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{2}{I\sqrt{p}} A^2,$$
where matrix $A := \left[ \zeta^{i^2j^2}+\alpha \right]_{i,j=1}^{(p-1)/2}$. It follows that
$$-S_p = T_p^2,\quad \text{where}\quad T_p := \det\left(\frac{1+I}{p^{1/4}}A\right).$$
The original question reduces to showing that $T_p$ is an integer.
A: This is a partial answer, as to why the determinant is always divisible by $2^{(p-3)/2}.$ We can write matrix $M_p:=\left[\left(\dfrac{i^2+j^2}{p}\right)\right]_{1\le i,j\le (p-1)/2}$ as $J+2A$, where $J$ is a matrix of ones, and $A$ is some symmetric matrix with integer entries. If $A$ is invertible, then we can apply the matrix determinant lemma, in particular, we find that $\det(M_p)=\det(2A)+1^t\mathrm{adj}(2A)1=2^{(p-1)/2}\det(A)+2^{(p-3)/2}1^t\mathrm{adj}(A)1.$ 
A: It can be seen that $S_p$ is not divisible by $p$. The argument about the decomposition of the matrix as $\frac{2}{i\sqrt{p}}A^2$ suggested above implies that $-S_p$ is a square in $\mathbb{Q}[\zeta_p][\sqrt{\lambda_p}]$, where $\lambda_p=-2i\sqrt{p}$. Writing $-S_p=(a\sqrt{\lambda_p}+b)^2$ with $a,b\in\mathbb{Q}[\zeta_p]$, we have $ab=0$ (since $S_p^2\in\mathbb{Z}\subset\mathbb{Q}[\zeta_p]$ and $\sqrt{\lambda_p}\not\in\mathbb{Q}[\zeta_p]$). The case $b=0$ is impossible as then $S_p$ is divisible by $p$, which follows from the fact that the norm of $\lambda_p$ is divisible by an odd power of $p$. Thus $-S_p=b^2$ with $b\in\mathbb{Q}[\zeta_p]$, which implies $-S_p$ is a square since $S_p$ is not divisible by $p$.   
