Provide a proof to this interesting identity: $$\frac{\begin{vmatrix} 1^0 & 1^2 & 1^3 & \cdots & 1^n \\ 2^0 & 2^2 & 2^3 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^0 & n^2 & n^3 & \cdots & n^n \end{vmatrix}}{\begin{vmatrix} 1^1 & 1^2 & 1^3 & \cdots & 1^n \\ 2^1 & 2^2 & 2^3 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^1 & n^2 & n^3 & \cdots & n^n \end{vmatrix}} = H_n \tag{1}\label{harmonic}$$
where $H_n$ is the nth Harmonic Number.
I accidently stumbled upon this observation while Programming. The inverse of a $(n+1)\times (n+1)$ Vandermonde Matrix (say $V_n$) with $ [\alpha_1 , \alpha_2, \alpha_3, \cdots, \alpha_n , \alpha_{n+1} ] = [0,1,2,\cdots,n-1,n] $, where we treat $0^0$ as $1$, $$V_n = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 1^0 & 1^1 & 1^2 & \cdots & 1^n \\ 2^0 & 2^1 & 2^2 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^0 & n^1 & n^2 & \cdots & n^n \end{bmatrix} $$ has some interesting properties.
The first row of $V_n^{-1}$ is $ [ 1 \ 0 \ 0 \cdots \ 0]$ . We also have this recurrence relation:
$$V_n^{-1}[m][0] = - \sum_{k=1}^n \frac{V_n^{-1}[m-1][0]}{k} \tag{2}\label{recurrence}$$
Using \eqref{recurrence} we find: $V_n^{-1}[1][0] = -H_n$.
We can also find $V_n^{-1}[1][0]$ using $ \ adj(V_n)/det(V_n)$ which gives $V_n^{-1}[1][0] = C_{01}/det(V_n)$ where $C_{01}$ is the Cofactor of $V_n[0][1]$. But it is same as the -ve of LHS of \eqref{harmonic}. QCD
We can prove the same using Cramer's rule by Reframing the Arguments here a little. But None of these proofs are "nice" in a way of not using inverse of Vandermonde Matrix but only elementary arguments like Expansion of the determinants or Mathematical Inductions or something else. I've put exhaustive efforts but haven't succeded yet.
Beside an elementary proof, any other perspective on the Result would be helpful and appreciated.
As @achillehui suggested, we can consider the Vandermonde determinant $$ P(x) = \begin{vmatrix} x^0 & x^1 & x^2 & \cdots & x^n \\ 1^0 & 1^1 & 1^2 & \cdots & 1^n \\ 2^0 & 2^1 & 2^2 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^0 & n^1 & n^2 & \cdots & n^n \end{vmatrix} = \prod_{k=1}^n (x-k) \cdot \prod_{j<k} (j-k) \, . $$ $P$ is a polynomial, its constant term is exactly the denominator of your fraction, and its coefficient of $x^1$ is the negative of the numerator: $$ \frac{\begin{vmatrix} 1^0 & 1^2 & 1^3 & \cdots & 1^n \\ 2^0 & 2^2 & 2^3 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^0 & n^2 & n^3 & \cdots & n^n \end{vmatrix}}{\begin{vmatrix} 1^1 & 1^2 & 1^3 & \cdots & 1^n \\ 2^1 & 2^2 & 2^3 & \cdots & 2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^1 & n^2 & n^3 & \cdots & n^n \end{vmatrix}} = \frac{-P'(0)}{P(0)} \, . $$ The logarithmic derivative of a product is the sum of the logarithmic derivatives of each factor: $$ \frac{P'(x)}{P(x)} = \sum_{k=1}^n \frac{1}{x-k} $$ and this gives the desired result: $$ -\frac{P'(0)}{P(0)} = \sum_{k=1}^n \frac{1}{k} = H_n \, . $$
In the same way one can show that more generally $$ \frac{\begin{vmatrix} x_1^0 & x_1^2 & x_1^3 & \cdots & x_1^n \\ x_2^0 & x_2^2 & x_2^3 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_n^0 & x_n^2 & x_n^3 & \cdots & x_n^n \end{vmatrix}}{\begin{vmatrix} x_1^1 & x_1^2 & x_1^3 & \cdots & x_1^n \\ x_2^1 & x_2^2 & x_2^3 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_n^1 & x_n^2 & x_n^3 & \cdots & x_n^n \end{vmatrix}} = \sum_{k=1}^n \frac{1}{x_k} $$ for any non-zero distinct numbers $x_1, \ldots, x_n$.