Let $B_n(x)$ denote the Bernoulli polynomial of degree $n$. If we construct the symmetric matrix $$ \mathbf{B}_N(x) = \begin{pmatrix} B_0(x) & B_1(x) & \cdots & B_N(x) \\ B_1(x) & B_2(x) & \cdots & B_{N+1}(x) \\ \vdots & \vdots & \ddots & \vdots \\ B_N(x) & B_{N+1}(x) & \cdots & B_{2N}(x) \end{pmatrix} $$ it would seem reasonable to guess that the determinant $|\mathbf{B}_N(x)|$ is a polynomial in $x$ of degree $(2N)!!$ However, at least for small $N$, the determinant turns out to be independent of $x$. For example, \begin{align} |\mathbf{B}_0(x)| &= 1 \\ |\mathbf{B}_1(x)| &= -\frac{1}{12} \\ |\mathbf{B}_2(x)| &= -\frac{1}{540} \\ |\mathbf{B}_3(x)| &= \frac{1}{42000} \end{align} and so on (I've tried it up to $N=20$ in Mathematica). My questions are naturally how we can explain this independence. Does it hold for every $N$? Can we find a simple expression for the value of the determinant?
The independence of $x$ doesn't hold for other matrix properties e.g. the trace, the eigenvalues and the first minors (except of course the leading principal ones).
The matrix $\mathbf{B}_N(x)$ is a Hankel matrix and its determinant is a Hankel determinant. It is a known result that the Binomial transform of a sequence preserves the Hankel determinant of the sequence. The Bernoulli polynomials are explicitly the binomial transform of the Bernoulli numbers and the result follows.
In other words, suppose we have any sequence $\;b_0,b_1,b_2,\dots\;$ of numbers and define the sequence of polynomials $\;b_n(x) := \sum_{k=0}^n {n\choose k} b_k x^k,\;$ then $\;b_n(0) = b_n\;$ and the Hankel determinant of the sequence of polynomials is independent of $x$.
In this case, let $\;a_n := |\mathbf{B}_N(x)|.\;$ Then it safisfies $a_0 = 1,\; a_n = (-1)^n a_{n-1}n!^6/((2n)!(2n+1)!).$