Assume all $X$ variables are standardized by the correlation transformation, like you mentioned, unit length scaled version of $\mathbf{X}$. The standardized model does not change the correlation between $X$ variables. $VIF$ can be calculated when standardized transformation of the original linear model is made. Let's denote the design matrix after standardized transformation as
\begin{align*}
\mathbf{X^*} = \begin{bmatrix} 1& X_{11}& \ldots &X_{1,p-1} \\ 1& X_{21}& \ldots &X_{2,p-1} \\ \vdots & \vdots & \vdots & \vdots \\ 1& X_{n1}& \ldots &X_{n,p-1} \\ \end{bmatrix}.
\end{align*}
Then
\begin{align*}
\mathbf{X^{*'}X^*} = \begin{bmatrix} n & \mathbf{0}' \\ \mathbf{0} & \mathbf{r}_{XX} \end{bmatrix},
\end{align*}
where $\mathbf{r}_{XX}$ is the correlation matrix of $X$ variables. We also know that
\begin{align*}
\sigma^2\{\hat{\beta}\}
& = \sigma^2 (\mathbf{X^{*'}X^*})^{-1}\\
& = \sigma^2 \begin{bmatrix} \frac{1}{n} & \mathbf{0}' \\ \mathbf{0} & \mathbf{r}^{-1}_{XX}. \end{bmatrix}\\
\end{align*}
$VIF_k$ for $k=1,2,\ldots,p-1$ is the $k$-th diagonal term of $\mathbf{r}^{-1}_{XX}$. We only need to prove this for $k = 1$ because you can permute the rows and columns of $r_{XX}$ to get the result for other $k$.
Let's define:
\begin{align*}
\mathbf{X}_{(-1)} = \begin{bmatrix} X_{12}&\ldots&X_{1,p-1} \\X_{22}&\ldots&X_{2,p-1}\\ \vdots & \vdots & \vdots \\ X_{n2}&\ldots&X_{n,p-1} \\ \end{bmatrix}, \mathbf{X}_1 = \begin{bmatrix} X_{11} \\ X_{21} \\ \vdots \\ X_{n1} \\ \end{bmatrix}.
\end{align*}
Note that both matrices are different from design matrices. Since we only care about the coefficients of $X$ variables, the $1$-vector of a design matrix can be ignored in our calculation. Hence, by using Schur's complement,
\begin{align*}
r^{-1}_{XX} (1,1)
& = (r_{11} - r_{1\mathbf{X}_{(-1)}} r^{-1}_{\mathbf{X}_{(-1)}\mathbf{X}_{(-1)}} r_{\mathbf{X}_{(-1)}1})^{-1} \\
& = (r_{11} - [r_{1\mathbf{X}_{(-1)}} r^{-1}_{\mathbf{X}_{(-1)}\mathbf{X}_{(-1)}}] r_{\mathbf{X}_{(-1)}\mathbf{X}_{(-1)}} [r^{-1}_{\mathbf{X}_{(-1)}\mathbf{X}_{(-1)}} r_{\mathbf{X}_{(-1)}1}])^{-1} \\
& = (1-\beta_{1\mathbf{X}_{(-1)}}' \mathbf{X}_{(-1)}' \mathbf{X}_{(-1)} \beta_{1\mathbf{X}_{(-1)}} )^{-1},
\end{align*}
where $\beta_{1\mathbf{X}_{(-1)}}$ is the regression coefficients of $X_1$ on $X_2, \ldots, X_{p-1}$ except the intercept. In fact, the intercept should be the origin, since all $X$ variables are standardized with mean zero.
On the other hand, (it would be more straightforward if we can write everything in explicit matrix form)
\begin{align*}
R_1^2
& = \frac{SSR}{SSTO} = \frac{\beta_{1\mathbf{X}_{(-1)}}' \mathbf{X}_{(-1)}' \mathbf{X}_{(-1)} \beta_{1\mathbf{X}_{(-1)}}}{1} \\
& = \beta_{1\mathbf{X}_{(-1)}}' \mathbf{X}_{(-1)}' \mathbf{X}_{(-1)} \beta_{1\mathbf{X}_{(-1)}}.
\end{align*}
Therefore
\begin{align*}
VIF_1 = r^{-1}_{XX} (1,1) = \frac{1}{1-R_1^2}.
\end{align*}
