Values of the determinants $\det[(j-k)^m+\delta_{jk}]_{1\le j,k\le n}\ (m=1,2,3,\ldots)$ For positive integers $m$ and $n$, let $D_m(n)$ denote the determinant $\det[(j-k)^m+\delta_{jk}]_{1\le j,k\le n}$, where the Kronecker delta $\delta_{jk}$ is $1$ or $0$ according as $j=k$ or not.
Recently, joint with my graduate student Han Wang, we proved in arXiv:2206.12317 the identity
$$D_1(n)=1+\frac{n^2(n^2-1)}{12}\tag{1}$$
via eigenvalues. (See also A079034 at OEIS.)
Based on my numerical computations, I conjectured that
$$\begin{align}D_2(n)&=1+\frac{n^2(n^2-1)}{1080}(n^5-5n^3-36n^2+4n+54)
\\&=\frac{(n^2-4)(n^2+2n+3)(n^5-2n^4-n^3-28n^2+60n-90)}{1080}\end{align}\tag{2}$$
and
$$\begin{align}D_3(n)=&1+\frac{n^2(n^2-1)}{672000}
\\&\times(n^{12}-19n^{10}+123n^8-337n^6+12376n^4-44144n^2+40000)\end{align}\tag{3}$$
(cf. A355175 and A355326 at OEIS). Note that the equation
$$x^5-2x^4-x^3-28x^2+60x-90=0$$
over $\mathbb Q$ is not solvable by radicals.
In view of $(1)-(3)$, I have the following general conjecture.
Conjecture. Let $m$ be any positive integer. Then $D_m(n)$ has the form
$1+n^2(n^2-1)P_m(n)$, where $P_m(n)$ is a polynomial in $n$ with rational number coefficients whose degree is $(m+1)^2-4$.
QUESTION. How to prove the formulas $(2)$ and $(3)$? Is the above general conjecture true?
Your comments are welcome!
 A: Edited to give a more complete (and cleaner) answer valid for all $m, n$. My original answer only applied for $n \geq m+1$, and only went through the proof for the case of $(x-y)^m$.
Say we have a polynomial $P(x, y) = \sum_{0 \leq a, b \leq m} c_{ab} x^a y^b$, and want to compute $\det(A + I)$ as a function of $n$, for $A = A(n) = (P(j, k))_{1 \leq j, k \leq n}$. Let $C = (c_{ab})_{0 \leq a, b \leq m}$, and let $V$ be the $n \times (m+1)$ matrix with columns $v_0, v_1, \dots, v_m$, where $v_d$ is the vector with $j$-th entry $j^d$ for $1 \leq j \leq n$. Then $v_a v_b^T$ is the $n \times n$ matrix with $(j, k)$-entry $j^a k^b$, so it follows that
$$A = \sum_{0 \leq a, b \leq m} c_{ab} v_a v_b^T = VCV^T$$
and thus by the Weinstein-Aronszajn identity
$$\det(A + I) = \det(VCV^T + I) = \det(CV^TV + I) = \det(CS + I),$$
where $S = V^T V$. Note that $C$, $S$, and hence also $CS+I$ are $(m+1) \times (m+1)$.
Now, for $0 \leq a, b \leq m$, the $(a, b)$-entry of $S$ is $v_a^T v_b = s_{a+b}(n)$, where $s_d(n) := \sum_{j=1}^n j^d$. Since each $s_d(n)$ can be expressed as a polynomial in $n$ of degree $d+1$ (with rational coefficients), all entries of $CS+I$ are fixed polynomials in $n$, and thus the equation above expresses $\det(A+I)$ as a polynomial in $n$ (with rational coefficients when all $c_{ab}$ are rational).
In the case of $P(x, y) = (x-y)^m$, the resulting formula is
$$D_m(n) = \det\left((-1)^{m-a} \tbinom{m}{a} s_{m-a+b}(n) + \delta_{ab}\right)_{0 \leq a, b \leq m}.$$
In this case, the degree of each term in the expansion of the determinant is $(m+1)^2$, so $D_m(n)$ has degree at most $(m+1)^2$. The calculated polynomial agrees with yours for $m=1$ and $m=2$, I haven't checked $m=3$ yet.
