Calculate the determinant of $n^\text{th}$ order: $$ \begin{vmatrix} 1 + a_1 & 1 + a_1^2 & \dots & 1 + a_1^n \\ 1 + a_2 & 1 + a_2^2 & \dots & 1 + a_2^n \\ \vdots & \vdots & \ddots & \vdots \\ 1 + a_n & 1 + a_n^2 & \dots & 1 + a_n^n \\ \end{vmatrix} $$
So whenever any two of the variables are equal the determinant becomes $0$. Therefore, it has $$\prod_{1 \le k < i \le n} (a_i - a_k)$$ as a factor. But I haven't been able to find the rest of the factors.
Any help is appreciated.
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
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$