[proofplan]
We cut out the image by the quadratic binomial equations expressing equality of products of degree-$d$ monomials with the same total exponent. The main point is local: on the chart where the pure coordinate $Y_{d e_i}$ is nonzero, the ratios $Y_{(d-1)e_i+e_j}/Y_{d e_i}$ recover the affine coordinates $x_j/x_i$. These recovered coordinates determine every Veronese coordinate by the binomial relations, so the closed binomial locus is locally exactly the Veronese image. The standard affine charts then give injectivity, closedness of the image, and an isomorphism from $\mathbb P_k^n$ onto the image.
[/proofplan]
[step:Handle the zero-dimensional projective space case]
If $n=0$, then
\begin{align*}
N=\binom{d}{d}-1=0.
\end{align*}
There is one degree-$d$ multi-index, namely $(d)$, and the map
\begin{align*}
\nu_d:\mathbb P_k^0\longrightarrow \mathbb P_k^0
\end{align*}
sends the unique point $[x_0]$ to $[x_0^d]$. Since $x_0\neq 0$ for a point of $\mathbb P_k^0$, this is the identity map on a one-point projective variety. Thus it is injective, its image is closed, and it is an isomorphism onto its image. Henceforth assume $n\geq 1$.
[/step]
[step:Verify that the projective formula is well-defined]
Let $(x_0,\dots,x_n)\in k^{n+1}\setminus\{0\}$ be an affine representative of a point of $\mathbb P_k^n$. Choose $i\in\{0,\dots,n\}$ with $x_i\neq 0$. Then the degree-$d$ coordinate indexed by $d e_i\in M_d$ is
\begin{align*}
x^{d e_i}=x_i^d\neq 0,
\end{align*}
so the target coordinate vector $(x^\alpha)_{\alpha\in M_d}$ is nonzero. If $\lambda\in k^\times$, then for every $\alpha\in M_d$,
\begin{align*}
(\lambda x)^\alpha=\lambda^{|\alpha|}x^\alpha=\lambda^d x^\alpha.
\end{align*}
Thus replacing the affine representative by $\lambda(x_0,\dots,x_n)$ rescales all target homogeneous coordinates by the same nonzero scalar $\lambda^d$. Therefore the formula defines a well-defined map $\nu_d:\mathbb P_k^n\to\mathbb P_k^N$, and since the coordinate functions $x^\alpha$ are homogeneous degree-$d$ polynomials with no common zero on $k^{n+1}\setminus\{0\}$, it is a morphism of projective varieties.
[/step]
[step:Define the closed binomial locus containing the Veronese image]
Let
\begin{align*}
M_d:=\{\alpha\in\mathbb N_0^{n+1}:|\alpha|=d\}
\end{align*}
be the set indexing the homogeneous coordinates of $\mathbb P_k^N$. For $\alpha,\beta,\gamma,\delta\in M_d$ satisfying
\begin{align*}
\alpha+\beta=\gamma+\delta,
\end{align*}
define the homogeneous quadratic polynomial
\begin{align*}
Q_{\alpha,\beta,\gamma,\delta}:=Y_\alpha Y_\beta-Y_\gamma Y_\delta\in k[Y_\eta:\eta\in M_d].
\end{align*}
For $\alpha\in M_d$ and $i\in\{0,\dots,n\}$, define the homogeneous degree-$d$ polynomial
\begin{align*}
R_{\alpha,i}:=Y_\alpha Y_{d e_i}^{d-1}-\prod_{j=0}^n Y_{(d-1)e_i+e_j}^{\alpha_j}\in k[Y_\eta:\eta\in M_d].
\end{align*}
For $\alpha\in M_d$, define the homogeneous degree-$d$ polynomial
\begin{align*}
S_\alpha:=Y_\alpha^d-\prod_{i=0}^n Y_{d e_i}^{\alpha_i}\in k[Y_\eta:\eta\in M_d].
\end{align*}
For a family $E$ of homogeneous polynomials in $k[Y_\eta:\eta\in M_d]$, let $V_+(E)\subseteq\mathbb P_k^N$ denote the set of projective points at which every polynomial in $E$ vanishes. Let $Z\subseteq\mathbb P_k^N$ be the common zero set of all these homogeneous polynomials:
\begin{align*}
Z:=V_+\bigl(Q_{\alpha,\beta,\gamma,\delta},R_{\alpha,i},S_\alpha:\alpha,\beta,\gamma,\delta\in M_d,\ \alpha+\beta=\gamma+\delta,\ i\in\{0,\dots,n\}\bigr).
\end{align*}
Since $Z$ is defined by homogeneous equations, it is a closed projective algebraic subset of $\mathbb P_k^N$.
For any point $x=[x_0:\cdots:x_n]\in\mathbb P_k^n$ and any $\alpha,\beta,\gamma,\delta\in M_d$ with $\alpha+\beta=\gamma+\delta$, the corresponding Veronese coordinates satisfy
\begin{align*}
x^\alpha x^\beta=x^{\alpha+\beta}=x^{\gamma+\delta}=x^\gamma x^\delta.
\end{align*}
Thus every polynomial $Q_{\alpha,\beta,\gamma,\delta}$ vanishes at $\nu_d(x)$. Also, for each $\alpha\in M_d$ and $i\in\{0,\dots,n\}$,
\begin{align*}
x^\alpha (x_i^d)^{d-1}=\prod_{j=0}^n (x_i^{d-1}x_j)^{\alpha_j},
\end{align*}
because $|\alpha|=d$. Hence every $R_{\alpha,i}$ vanishes at $\nu_d(x)$. Finally,
\begin{align*}
(x^\alpha)^d=\prod_{i=0}^n (x_i^d)^{\alpha_i},
\end{align*}
so every $S_\alpha$ vanishes at $\nu_d(x)$. Therefore
\begin{align*}
\nu_d(\mathbb P_k^n)\subseteq Z.
\end{align*}
[guided]
The target coordinates $Y_\alpha$ are meant to behave like the degree-$d$ monomials $x^\alpha$. If two products of such monomials have the same total exponent, then they are the same monomial in the $x_i$. This motivates imposing exactly those identities as equations on the target projective space.
Formally, define
\begin{align*}
M_d:=\{\alpha\in\mathbb N_0^{n+1}:|\alpha|=d\}.
\end{align*}
For $\alpha,\beta,\gamma,\delta\in M_d$ with $\alpha+\beta=\gamma+\delta$, define
\begin{align*}
Q_{\alpha,\beta,\gamma,\delta}:=Y_\alpha Y_\beta-Y_\gamma Y_\delta.
\end{align*}
For $\alpha\in M_d$ and $i\in\{0,\dots,n\}$, define
\begin{align*}
R_{\alpha,i}:=Y_\alpha Y_{d e_i}^{d-1}-\prod_{j=0}^n Y_{(d-1)e_i+e_j}^{\alpha_j}.
\end{align*}
For $\alpha\in M_d$, define
\begin{align*}
S_\alpha:=Y_\alpha^d-\prod_{i=0}^n Y_{d e_i}^{\alpha_i}.
\end{align*}
These are homogeneous polynomials, of degrees $2$, $d$, and $d$ respectively, so their common zero set
\begin{align*}
Z:=V_+\bigl(Q_{\alpha,\beta,\gamma,\delta},R_{\alpha,i},S_\alpha:\alpha,\beta,\gamma,\delta\in M_d,\ \alpha+\beta=\gamma+\delta,\ i\in\{0,\dots,n\}\bigr)
\end{align*}
is a closed projective algebraic subset of $\mathbb P_k^N$.
Now take a point
\begin{align*}
x=[x_0:\cdots:x_n]\in\mathbb P_k^n.
\end{align*}
Its Veronese image has $\alpha$-coordinate $x^\alpha$. If $\alpha+\beta=\gamma+\delta$, then multiplication of monomials gives
\begin{align*}
x^\alpha x^\beta=x^{\alpha+\beta}=x^{\gamma+\delta}=x^\gamma x^\delta.
\end{align*}
Therefore
\begin{align*}
Q_{\alpha,\beta,\gamma,\delta}(\nu_d(x))=0.
\end{align*}
The additional equations also vanish on the Veronese image. Indeed, for $\alpha\in M_d$ and $i\in\{0,\dots,n\}$,
\begin{align*}
x^\alpha (x_i^d)^{d-1}=\prod_{j=0}^n (x_i^{d-1}x_j)^{\alpha_j},
\end{align*}
since $|\alpha|=d$, and
\begin{align*}
(x^\alpha)^d=\prod_{i=0}^n (x_i^d)^{\alpha_i}.
\end{align*}
Thus all polynomials defining $Z$ vanish on $\nu_d(x)$, and we obtain
\begin{align*}
\nu_d(\mathbb P_k^n)\subseteq Z.
\end{align*}
[/guided]
[/step]
[step:Show the binomial locus is covered by the pure-power affine charts]
We prove that
\begin{align*}
Z\subseteq \bigcup_{i=0}^n W_i.
\end{align*}
Let $y=[y_\alpha]_{\alpha\in M_d}\in Z$. Since $y$ is a point of projective space, there exists $\alpha\in M_d$ such that $y_\alpha\neq 0$.
Since $y\in Z$, it satisfies the defining equation $S_\alpha=0$. Hence
\begin{align*}
y_\alpha^d=\prod_{i=0}^n y_{d e_i}^{\alpha_i}.
\end{align*}
Since $y_\alpha\neq 0$, the right-hand side is nonzero. Hence at least one factor $y_{d e_i}$ is nonzero, and therefore $y\in W_i$ for some $i$. Thus $Z$ is covered by the charts $Z\cap W_i$.
[/step]
[step:Recover affine coordinates on each pure-power chart]
Fix $i\in\{0,\dots,n\}$. Define the affine chart
\begin{align*}
U_i:=\{[x_0:\cdots:x_n]\in\mathbb P_k^n:x_i\neq 0\}
\end{align*}
and the affine chart
\begin{align*}
W_i:=\{[Y_\alpha]_{\alpha\in M_d}\in\mathbb P_k^N:Y_{d e_i}\neq 0\}.
\end{align*}
On $U_i$, define affine coordinates
\begin{align*}
t_j:U_i\longrightarrow k
\end{align*}
for $j\neq i$ by
\begin{align*}
[x_0:\cdots:x_n]\longmapsto \frac{x_j}{x_i}.
\end{align*}
On $Z\cap W_i$, define regular functions
\begin{align*}
s_j:Z\cap W_i\longrightarrow k
\end{align*}
for $j\neq i$ by
\begin{align*}
[y_\alpha]_{\alpha\in M_d}\longmapsto \frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}.
\end{align*}
For $x\in U_i$, the Veronese coordinates give
\begin{align*}
s_j(\nu_d(x))=\frac{x_i^{d-1}x_j}{x_i^d}=\frac{x_j}{x_i}=t_j(x).
\end{align*}
Thus the functions $s_j$ recover the affine coordinates of $x$ from $\nu_d(x)$.
Now let $y=[y_\alpha]_{\alpha\in M_d}\in Z\cap W_i$. Define a point $x(y)\in U_i$ by setting the homogeneous coordinates
\begin{align*}
x_i(y):=1
\end{align*}
and, for $j\neq i$,
\begin{align*}
x_j(y):=s_j(y)=\frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}.
\end{align*}
We claim that $\nu_d(x(y))=y$.
Let $\alpha\in M_d$. Since $y\in Z$, it satisfies the defining equation $R_{\alpha,i}=0$. Thus
\begin{align*}
y_\alpha y_{d e_i}^{d-1}=\prod_{j=0}^n y_{(d-1)e_i+e_j}^{\alpha_j}.
\end{align*}
Here the factor with $j=i$ is $y_{d e_i}^{\alpha_i}$. Dividing by $y_{d e_i}^d$, which is allowed because $y\in W_i$, gives
\begin{align*}
\frac{y_\alpha}{y_{d e_i}}=\prod_{j\neq i}\left(\frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}\right)^{\alpha_j}.
\end{align*}
The right-hand side is precisely $x(y)^\alpha$, because $x_i(y)=1$ and $x_j(y)=s_j(y)$ for $j\neq i$. Hence
\begin{align*}
[y_\alpha]_{\alpha\in M_d}=[x(y)^\alpha]_{\alpha\in M_d}.
\end{align*}
Therefore $y=\nu_d(x(y))$.
[guided]
The purpose of this step is to prove that the equations defining $Z$ do not introduce extra points. We do this on an affine chart, where we can divide by one nonzero coordinate and write explicit formulas.
Fix $i\in\{0,\dots,n\}$. On the source chart
\begin{align*}
U_i=\{[x_0:\cdots:x_n]\in\mathbb P_k^n:x_i\neq 0\},
\end{align*}
the ordinary affine coordinates are the ratios
\begin{align*}
t_j([x_0:\cdots:x_n])=\frac{x_j}{x_i}
\end{align*}
for $j\neq i$. On the target chart
\begin{align*}
W_i=\{[Y_\alpha]_{\alpha\in M_d}:Y_{d e_i}\neq 0\},
\end{align*}
the Veronese coordinates contain the monomial $x_i^{d-1}x_j$. Therefore the ratio
\begin{align*}
\frac{Y_{(d-1)e_i+e_j}}{Y_{d e_i}}
\end{align*}
should recover $x_j/x_i$. We make this precise by defining
\begin{align*}
s_j:Z\cap W_i\longrightarrow k
\end{align*}
for $j\neq i$ by
\begin{align*}
[y_\alpha]_{\alpha\in M_d}\longmapsto \frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}.
\end{align*}
This is a regular function on the affine chart because the denominator $y_{d e_i}$ is nonzero there. If $y=\nu_d(x)$ for some $x\in U_i$, then
\begin{align*}
s_j(\nu_d(x))=\frac{x_i^{d-1}x_j}{x_i^d}=\frac{x_j}{x_i}.
\end{align*}
So the functions $s_j$ recover the affine coordinates of a Veronese point.
Now take an arbitrary point
\begin{align*}
y=[y_\alpha]_{\alpha\in M_d}\in Z\cap W_i.
\end{align*}
We use the recovered coordinates to define a candidate preimage:
\begin{align*}
x_i(y):=1
\end{align*}
and, for $j\neq i$,
\begin{align*}
x_j(y):=\frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}.
\end{align*}
This gives a point $x(y)\in U_i$.
It remains to prove that every coordinate of $y$ agrees projectively with the corresponding Veronese monomial of $x(y)$. Let $\alpha\in M_d$. Since $y\in Z$, the defining equation $R_{\alpha,i}=0$ gives
\begin{align*}
y_\alpha y_{d e_i}^{d-1}=\prod_{j=0}^n y_{(d-1)e_i+e_j}^{\alpha_j}.
\end{align*}
This equation is precisely the monomial identity needed on the chart $W_i$. Since $y\in W_i$, we have $y_{d e_i}\neq 0$, so we may divide by $y_{d e_i}^d$:
\begin{align*}
\frac{y_\alpha}{y_{d e_i}}=\prod_{j\neq i}\left(\frac{y_{(d-1)e_i+e_j}}{y_{d e_i}}\right)^{\alpha_j}.
\end{align*}
The factor with $j=i$ contributes $1$, because $(d-1)e_i+e_i=d e_i$. By the definition of $x(y)$, the right-hand side is $x(y)^\alpha$. Thus
\begin{align*}
\frac{y_\alpha}{y_{d e_i}}=x(y)^\alpha.
\end{align*}
Equivalently,
\begin{align*}
[y_\alpha]_{\alpha\in M_d}=[x(y)^\alpha]_{\alpha\in M_d}.
\end{align*}
Therefore $y=\nu_d(x(y))$. This proves that every point of $Z\cap W_i$ is actually a Veronese point, with preimage given by the explicit coordinate recovery formulas.
[/guided]
[/step]
[step:Identify each affine chart with its image by explicit inverse maps]
The preceding step proves
\begin{align*}
Z\cap W_i=\nu_d(U_i)
\end{align*}
for every $i\in\{0,\dots,n\}$. The map
\begin{align*}
\nu_d|_{U_i}:U_i\longrightarrow Z\cap W_i
\end{align*}
is regular because, in the affine coordinates $t_j=x_j/x_i$ on $U_i$, every target affine coordinate
\begin{align*}
\frac{Y_\alpha}{Y_{d e_i}}
\end{align*}
is the polynomial
\begin{align*}
\prod_{j\neq i} t_j^{\alpha_j}.
\end{align*}
Its inverse
\begin{align*}
\rho_i:Z\cap W_i\longrightarrow U_i
\end{align*}
is the regular map whose affine coordinates are the functions $s_j$ defined by
\begin{align*}
s_j=\frac{Y_{(d-1)e_i+e_j}}{Y_{d e_i}}
\end{align*}
for $j\neq i$. The computations in the preceding step show
\begin{align*}
\rho_i\circ \nu_d|_{U_i}=\operatorname{id}_{U_i}
\end{align*}
and
\begin{align*}
\nu_d|_{U_i}\circ \rho_i=\operatorname{id}_{Z\cap W_i}.
\end{align*}
Thus $\nu_d|_{U_i}:U_i\to Z\cap W_i$ is an isomorphism of affine varieties. Consequently, it induces an isomorphism on affine coordinate rings.
[/step]
[step:Conclude closedness, injectivity, and the global isomorphism onto the image]
From the chart-covering step,
\begin{align*}
Z\subseteq \bigcup_{i=0}^n W_i.
\end{align*}
From the chart identification step,
\begin{align*}
Z\cap W_i=\nu_d(U_i)
\end{align*}
for every $i$. The standard affine charts $U_i$ cover $\mathbb P_k^n$: every projective point has at least one nonzero homogeneous coordinate. Therefore we obtain
\begin{align*}
Z=\bigcup_{i=0}^n (Z\cap W_i)=\bigcup_{i=0}^n \nu_d(U_i)=\nu_d(\mathbb P_k^n).
\end{align*}
Hence the image of $\nu_d$ is the closed projective algebraic subset $Z$.
The same chartwise inverse formulas also prove injectivity. If $\nu_d(x)=\nu_d(x')$, choose $i$ with $x_i\neq 0$. Then $\nu_d(x)\in W_i$, so $\nu_d(x')\in W_i$, and the inverse $\rho_i$ gives
\begin{align*}
x=\rho_i(\nu_d(x))=\rho_i(\nu_d(x'))=x'.
\end{align*}
Finally, the affine isomorphisms
\begin{align*}
\nu_d|_{U_i}:U_i\longrightarrow \nu_d(\mathbb P_k^n)\cap W_i
\end{align*}
agree on overlaps because they are restrictions of the same global map $\nu_d$, and their inverses agree on overlaps because they recover the same homogeneous point in $\mathbb P_k^n$. Therefore these affine chart isomorphisms glue to an isomorphism of projective varieties
\begin{align*}
\mathbb P_k^n\longrightarrow \nu_d(\mathbb P_k^n).
\end{align*}
This proves that the degree-$d$ Veronese map is injective, has closed image, and is an isomorphism onto its image.
[/step]