A binomial determinant formula: a new variant

Question

In a previous MO question, the OP asks a proof for $\det_{1\leq i,j\leq n}\left(\binom{i}{2j}+\binom{-i}{2j}\right)=1$. Subsequently, Gjergji Zaimi generalized the problem to $$\det_{1\le i,j\le n}\left( \binom{x_i}{2j}+ \binom{-x_i}{2j}\right)=\prod_{i=1}^n x_i^2 \prod_{i<j} (x_j^2-x_i^2) \prod_{j=1}^n \frac{2}{(2j)!}.$$ Let's define $$\Psi_n:=\prod_{i=1}^n x_i^2 \prod_{i<j} (x_j^2-x_i^2) \prod_{j=1}^{\mathbf{n+1}} \frac{2}{(2j)!}.$$ Recall the elementary symmetric functions $e_k=e_k(x_1^2,\dots,x_n^2)$. Now, I would like ask:

Question 1. Is it true that there exist positive integers $a_0, a_1, \dots, a_n$ such that $$\det_{1\le i,j\le n}\left( \binom{x_i}{2j+2}+ \binom{-x_i}{2j+2}\right)= (a_0+a_1e_1+\cdots+a_{n-1}e_{n-1}+a_ne_n)\cdot\Psi_n.$$ Question 2. What are these coefficients $a_0, a_1, \dots, a_n$?

NOTE. For any $n$, we observe that $a_n=1, a_{n-1}=11, a_{n-2}=661, a_{n-3}=151451$. These sequence does not appear on OEIS.

Answer 1

Let $M$ be the matrix in question. The entry $M_{ij}$ is of the form $\frac{2x_i^2}{(2j+2)!}p_j(x_i)$ for some even polynomial $p_j$ of degree $2j$. After factoring out the $2x_i^2$ terms from each row and the $\frac{1}{(2j+2)!}$ terms from each column, we are left with the matrix $(p_j(x_i))_{i,j=1}^n$. This can be column reduced to a matrix $P$ with rows $(x_i^2+a_1,x_i^4-a_2,x_i^6+a_3,\dots)$, where $a_1=11$, $a_2=661$, $a_3=151451$, etc. It remains to show that $\det P$ is $(a_ne_0+a_{n-1}e_1+\cdots+a_1e_{n-1}+e_n)\prod_{i<j}(x_j^2-x_i^2)$.

$P$ breaks into the sum of a Vandermonde matrix $V$ in the indeterminates $x_i^2$ and a rank-$1$ matrix $uv^\top$, where $u=(1,\dots,1)^\top$ and $v=(a_1,-a_2,a_3,\dots,(-1)^{n-1}a_n)^\top$. By the matrix determinant lemma, $\det P=\det V(1+v^\top V^{-1} u)$. Since $\det V=e_n\prod_{i<j}(x_j^2-x_i^2)$, it suffices to show that $V^{-1}u=(e_{n-1},-e_{n-2},\dots,(-1)^{n-1} e_0)^\top/e_n$. Indeed, the $i$th entry of $V(e_{n-1},-e_{n-2},\dots,(-1)^{n-1} e_0)^\top/e_n$ is $$(e_n-e_n+x_i^2e_{n-1}-\cdots\pm x_i^{2n}e_0)/e_n=(e_n\pm(x_i^2-x_1^2)\cdots(x_i^2-x_n^2))/e_n=1.$$