Yet, another numerical variant of the Vandermonde matrix In my earlier (soft) MO post, an elementary response was given by Ofir Gorodetsky in regard to the determinant of the symbolic counterpart to the numerical matrix $\mathbf{M}_n=(i^j-j^i)_{i,j}^{1,n}$.
I'm now looking at yet another variation. Let $\mathbf{A}_n=(j^j-i^j)_{i,j}^{1,n}$ be an $n\times n$-matrix. Denote the (signed) Stirling numbers of the first kind by $s(n,k)$. For clarity, here are some contrasting examples:
$$\mathbf{M}_3=\begin{pmatrix} 0&-1&-2 \\ 1&0&-1 \\2&1&0 \end{pmatrix}
\qquad \text{and} \qquad \mathbf{A}_3=\begin{pmatrix}0&3&26 \\-1&0&19 \\ -2&-5&0 \end{pmatrix}.$$
I would like to ask:

QUESTION. Is this true? With the convention that $0^0=1$, we have
$$\det\mathbf{A}_n=\prod_{j=1}^{n-1}j!\cdot\sum_{k=0}^n\,s(n+1,k+1)\cdot k^k.$$

Remark. If it helps, since the factor $\prod_{j=1}^{n-1}j!$ is exactly the determinant of the Vandermonde matrix, $\mathbf{V}_n=(i^{j-1})_{i,j}^{1,n}$, we can say
$$\frac{\det\mathbf{A}_n}{\det \mathbf{V}_n}=\sum_{k=0}^n\,s(n+1,k+1)\cdot k^k.$$ This formalism goes in the "spirit" of the (specialization) $s_{\lambda}(1,\cdots,1)$ of the Schur polynomials.
 A: This is true and this is not about numbers $j^j$ at all. Consider a more genral matrix $(c_j-i^j)_{1\leqslant i,j\leqslant n}$. Denote $f_j(x)=c_j-x^j$, and find the numbers $\beta$, $\alpha_1,\ldots,\alpha_{n-1}$ such that $f_n(x)+\sum_{j=1}^{n-1}\alpha_j f_j(x)=\beta$ for all $x\in \{1,2,\ldots,n\}$. Then the matrix $(f_j(i))_{i,j}$ has the same determinant as if we replace the column with $f_n$ by a column of $\beta$'s, that is seen from operations with columns. Now by operations with columns we may replace each $f_j$ to $-x^j$, so we get almost a Vandermonde matrix (up to $n-1$ changes of signs in the columns and the cyclic shift, and $\beta$'s instead of 1's), so its determinant equals to the determinant of Vandermonde matrix for the numbers $1,\ldots,n$ times $(-1)^{n-1}(-1)^{n-1}\beta=\beta$. Thus, it appears to find $\beta$. For this, write
$$
f_n+\sum_{j=1}^{n-1}\alpha_jf_j(x)=\beta-(x-1)(x-2)\ldots(x-n)=
\beta-\sum_{k=0}^n s(n+1,k+1)x^k.
$$
Equalising the coefficients of $x^k$ we get $\alpha_k=s(n+1,k+1)$ for all $k=1,\ldots,n$ (where $\alpha_n=1$ by agreement), thus $\beta-s(n+1,1)=\sum_{k=1}^n c_k s(n+1,k+1)$, $\beta=\sum_{k=0}^n c_k s(n+1,k+1)$, where $c_0=1$ by agreement.
