Assume we have a matrix $A$ of the following form
$$ 1\leq i,j \leq n \hspace{0.5cm} A = [a_{ij}] \in M_n(\{0,1\}) \hspace{0.4cm}a_{ij} = \begin{cases} 1 \hspace{0.3cm}i+j \in \mathbb{P} \\ 0 \hspace{0.3cm} i + j \notin \mathbb{P} \end{cases} $$
where $\mathbb{P}$ denotes the set of prime numbers. How can I prove the absolute value of its determinant is a perfect square?
So far, I've managed to figure out it's a symmetric matrix and has zeros on its diagonal except at $a_{11}$, which isn't much, and I haven't even made use of the prime property.
Here is a simplified approach to this problem inspired by the proof presented below.
First, we consider the case where $n$ is even. Rearrange both the rows and columns of $A$ via the permutation $$ [1,2,3,\dots,n] \mapsto [2,4,\dots,n,1,3,\dots,n-1]. $$ Rearranging the rows and columns simultaneously (i.e. applying a permutation similarity) does not affect the determinant. Because the $i,j$ row entry is zero whenever $i$ and $j$ are even, the resulting matrix has the form $$ B = \pmatrix{0 & P\\P^\top & Q} $$ for some matrices $P$ and $Q$ of size $n \times n$. We can rearrange the columns of this matrix (possibly changing the sign, but not the absolute value) to get the matrix $$ C = \pmatrix{P & 0\\Q& P^\top}. $$ Thus, we have $\det(A) = \det(B) = \pm \det(C) = \pm \det(P)\det(P^T) = \pm \det(P)^2$.
In the case where $n$ is odd, we can similarly apply the permutation $$ [1,2,3,\dots,n] \mapsto [1,3,\dots,n-2,n,2,4,\dots,n-1]. $$ Doing so gives us a matrix of the form $$ B = \pmatrix{ 1 & 0 & x^T\\ 0 & 0 & P\\ x & P^T & 0}. $$ We note that this has the same determinant as the matrix $$ C = \pmatrix{ 1 & 0 & 0 \\ 0&I&0\\ -x & 0 & I} B = \pmatrix{ 1 & 0 & x^T\\ 0 & 0 & P\\ 0 & P^T & -xx^T }. $$ Because $C$ is block upper-triangular, we have $$ \det(C) = 1 \cdot \det \pmatrix{0&P\\P^T & -xx^T}. $$ By our work in the previous case, we see that the second determinant must have an absolute value equal to a perfect square.
This is a slight expansion of the proof presented here, which is linked to in the comments to the OEIS sequence of these determinants here.
Proof: Consider the commutative ring $A=\Bbb Z[x_1,\dots,x_n]$ and let $M_n$ be the matrix of a quadratic form on the rational field of fractions of $A$. Let $e_0,\dots,e_n$ be the vectors for the canonical basis. Let $b$ be the associated bilinear form. By hypothesis, $b(e_i,e_j)=0$ if $i+j$ is even, except for $b(e_0,e_0)=1$.
Consider the odd case $n = 2k+1$. Let $$ f_0=e_n\\ f_1=e_n+e_{n-2}\\ f_2=e_n+e_{n-2}+e_{n-4}\\ \vdots\\ f_k=e_n+e_{n-2}+\cdots +e_1\\ f_{k+1}=e_{n-1}\\ f_{k+2}=e_{n-1}+e_{n-3}\\ \vdots \\ f_n=e_{n-1}+e_{n-3}+\cdots +e_0.\\ $$ Then the matrix of $b$ relative to the basis $\{f_0,f_1,\dots,f_n\}$ has the form $$ N = \pmatrix{0&P\\P^\top &Q}, $$ Where $Q$ is the matrix with a $1$ as its bottom-right entry and zeros elsewhere. Since $P$ is a square matrix, the determinant of $N$ is $\pm$ a square (sign depending on $n \bmod 4$). Indeed, we can circularly shift the columns of $N$ repeatedly to get the matrix $$ \pmatrix{P & 0\\Q & P^\top}, $$ which has determinant $\det(P)\det(P^\top) = \det(P)^2$. On the other hand, $N$ is obtained from $M$ via $N = X^\top MX$, where $X$ is the change of basis matrix from $\{e_0,\dots,e_n\}$ to $\{f_0,\dots,f_n\}$. Thus, we have $\det(M) = \det(N)/\det(X)^2$, which means that $\det(M)$ must also have an absolute value equal to a perfect square in $\Bbb Q(x_1,\dots,x_n)$, and this determinant is a polynomial in $\Bbb Z[x_1,\dots,x_n]$ which is integrally closed. Therefore, $\det(M)$ is $\pm$ a perfect square in $\Bbb Z[x_1,\dots,x_n]$.
For $n$ even , we have the same form for $N$. $P$ is not a square matrix, but the determinant of $N$ is (except sign) the same as for the matrix obtained after deleting last row and last column of $N$. From there, the steps from the case where $P$ is square apply.