También es posible evaluar un determinante más general, donde la inteligente observación de @user1551 no está disponible. Sea $c_0,\dots, c_1$ sea arbitraria y que $X$ sea una variable polinómica. Entonces evaluaremos $$ \begin{vmatrix} c_0&1&\cdots & 1 \\ c_1 X&a_2&\cdots &a_n\\ c_2 X^2&a_2^{2}&\cdots&a_n^{2}\\ \vdots&\vdots&\ddots& \vdots\\ c_{n-1}X^{n-1}&a_2^{n-1}&\cdots&a_n^{n-1}\\ \end{vmatrix}. $$
Mediante una expansión de Lagrange en la primera columna obtenemos que es $$ \sum_{0}^{n-1} (-1)^{k} c_k X^k V_k \hspace{10em}\text{(*)} $$ donde $V_k$ es el determinante de la correspondiente $(n-1)\times(n-1)$ menor en $$ \begin{vmatrix} 1&\cdots & 1 \\ a_2&\cdots &a_n\\ a_2^{2}&\cdots&a_n^{2}\\ \vdots&\ddots& \vdots\\ a_2^{n-1}&\cdots&a_n^{n-1}\\ \end{vmatrix}. $$ Pero ahora consideremos el caso especial en el que todos $c_k=1$ es el determinante estándar de Vandermonde, por lo que tenemos $$ \sum_{0}^{n-1} (-1)^{k} X^k V_k = \prod_{r=2}^{n} (X-a_r)\prod_{2\leqslant s < t \leqslant n} (a_t - a_s). $$
Escribimos, de la forma habitual, $$ \prod_{r=2}^{n} (X-a_r)=\sum_{k=0}^{n-1}(-1)^{k} \sigma_{n-1-k} X^k; $$ no son más que las funciones simétricas habituales del $a_2,\dots, a_{n}$ tomando su suma, su suma de dos en dos, y así sucesivamente.
Los coeficientes $V_k$ entonces todos tienen un factor común $\prod_{2\leqslant s < t \leqslant n} (a_t - a_s)$ el Vandermonde de $a_2, \dots, a_n$ (que se ocupa de la simetría oblicua en estas variables); el otro factor en cada caso es simplemente la función simétrica $\sigma_{n-k}$ : $$ V_k= \sigma_{n-k}\prod_{2\leqslant s < t \leqslant n} (a_t - a_s).$$
En $c_0=1, c_1=\dots=c_n=1$ y sustituyendo $X=a_1$ en (*) da el valor del determinante de la pregunta original.