Calculate the determinant of $n^\text{th}$ order

When after experimentation we find that the determinant is a product of some factors, it's a good idea to see if the matrix itself can be factored into a product of matrices. If the matrix $A = BC$, where the determinants of $B$ and $C$ are easily calculated, then we can recover $\det(A)$ using the identity $\det (BC) = \det (B) \det (C)$.

Here, as you noted, $\prod_{1 \leq i < j \leq n} (a_j - a_i)$ is a factor of the final determinant. It is well-known that this is the determinant of the Vandermonde matrix $$V = \begin{bmatrix} 1 & a_1 & a_1^2 & \cdots & a_1^{n-1} \\ 1 & a_2 & a_2^2 & \cdots & a_2^{n-1} \\ 1 & a_3 & a_3^2 & \cdots & a_3^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & a_n & a_n^2 & \cdots & a_n^{n-1} \end{bmatrix}$$ so this begs of whether we can factor the given matrix $A$ as $A = V B$ for some matrix $B$. And in fact, we can! In fact, we can show that with $$B = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1& (-1)^{n-1} e_n + 1 \\ 1 & 0 & 0 & \cdots & 0& (-1)^{n-2} e_{n-1} \\ 0 & 1 & 0 & \cdots & 0 & (-1)^{n-3} e_{n-2} \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots \\ \vdots & \vdots &\vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & (-1)^0 e_1\end{bmatrix}$$ the factorization $A = V B$ holds. Here the $e_i$ denotes the degree $i$ elementary symmetric polynomial in the variables $a_1, \cdots a_n$. The correctness of the first $n - 1$ columns of $A$ resulting from this product are easy to verify. What about the pesky last column? We basically want to verify that $$V \cdot \begin{bmatrix} (-1)^{n-1} e_n \\ (-1)^{n-2} e_{n-1} \\ \vdots \\ e_1 \end{bmatrix} = \begin{bmatrix} a_1^n \\ a_2^n \\ \vdots \\ a_n^n\end{bmatrix}$$ because if this is true, then the last column will also match, since by adding the extra $+1$ in the upper right entry of $B$ the last column will be $$V \cdot \left(~\begin{bmatrix} (-1)^{n-1} e_n \\ (-1)^{n-2} e_{n-1} \\ \vdots \\ e_1 \end{bmatrix} + \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}~\right) = V \cdot \begin{bmatrix} (-1)^{n-1} e_n \\ (-1)^{n-2} e_{n-1} \\ \vdots \\ e_1 \end{bmatrix} + \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} = \begin{bmatrix} 1 + a_1^n \\ 1+ a_2^n \\ \vdots \\ 1 + a_n^n \end{bmatrix}$$

Note that this is equivalent to asserting

Lemma: For all of the $a_i$, the following holds: $$a_i^n = \sum_{k = 0}^{n-1} (-1)^{n - k - 1} \cdot a_i^k \cdot e_{n - k}$$

Proof of Lemma: This fact follows because the elementary symmetric polynomials are the resulting coefficients of the monic polynomial $$p(\lambda) = \prod_{i = 1}^n (\lambda - a_i) = \lambda^n - e_1 \lambda^{n-1} + e_2 \lambda^{n-2} - \cdots + (-1)^n e_n$$ Plugging in $\lambda = a_i$, we have $$p(a_i) = \sum_{k = 0}^n (-1)^{n-k} \cdot a_i^k \cdot e_{n-k} = 0 \implies a_i^n = \sum_{k = 0}^{n-1} (-1)^{n-k - 1} \cdot a_i^k \cdot e_{n - k}$$

We've now verified our factorization is correct, so what remains is to calculate the determinant of $B$. Laplace expanding along the last column, we can see that \begin{align*} \det B = (-1)^{n+1} [(-1)^{n-1} \cdot e_n + e_0] &+ \sum_{k = 1}^{n-1} \left[(-1)^{n + k - 1} \cdot \det(B_{k, n}) \cdot (-1)^{n - k - 1} \cdot e_{n-k}\right] \\ \\ &= e_n + (-1)^{n-1} \cdot e_0 + \sum_{k = 1}^{n-1} \det(B_{k, n}) \cdot e_{n - k} \end{align*} where $B_{k, n}$ is the minor $B$ obtained by removing the $k$th row and $n$th column. It is not too difficult to verify that $\det(B_{k, n}) = (-1)^{k - 1}$ (although if it is not clear, I can post an addendum explaining why). Hence, the determinant of $B$ is $$e_n + (-1)^{n + 1} e_0 + \sum_{k = 1}^{n-1} \det(B_{k, n}) \cdot e_{n - k} = e_n + e_{n-1} - e_{n-2} + \cdots + (-1)^{n + 1} e_0$$ So it follows that $$\det A = \det(V) \det (B) = \left[\prod_{1 \leq i < j \leq n} (a_j - a_i)\right] \cdot (e_n + e_{n - 1} - e_{n-2} + \cdots + (-1)^{n+1} e_0)$$ verifying lhf's conjecture. (Note that the above expression also has the correct sign as well). $\square$