[proofplan]
We first pass from invariance of the [homogeneous polynomial](/page/Homogeneous%20Polynomial) on the diagonal to invariance of its polarized symmetric multilinear form. Then we insert the one-parameter subgroup $\exp(tY)$ into this multilinear invariance identity and differentiate at $t=0$. The product rule contributes one derivative from each slot, and the differential of the adjoint action along $\exp(tY)$ is the Lie bracket $[Y,X_j]$.
[/proofplan]
[step:Polarize the invariant polynomial to obtain an invariant multilinear form]
Define \begin{align*} p:\mathfrak g\to \mathbb R \end{align*} to be the given homogeneous polynomial, and let
\begin{align*} P:\mathfrak g^k\to \mathbb R \end{align*} be its polarization. Thus $P$ is symmetric and $k$-linear, and $p(X)=P(X,\dots,X)$ for every $X\in\mathfrak g$.
We claim that $P$ is invariant under the diagonal adjoint action of $G$: \begin{align*} P(\operatorname{Ad}_g X_1,\dots,\operatorname{Ad}_g X_k)=P(X_1,\dots,X_k) \end{align*} for every $g\in G$ and every $X_1,\dots,X_k\in\mathfrak g$.
Indeed, the polarization identity expresses $P(X_1,\dots,X_k)$ as a fixed linear combination of values of $p$ at sums of the vectors $X_i$. Since $\operatorname{Ad}_g:\mathfrak g\to\mathfrak g$ is linear, each such sum is carried to the corresponding sum of the vectors $\operatorname{Ad}_g X_i$. Since $p$ is $\operatorname{Ad}$-invariant, every value of $p$ in the polarization formula is unchanged. Hence the polarized form $P$ is $\operatorname{Ad}$-invariant.
[guided]
We need invariance not only for $p(X)$ on the diagonal, but for the full multilinear expression $P(X_1,\dots,X_k)$. The bridge is the polarization identity: because $p$ is a homogeneous polynomial of degree $k$, its associated symmetric $k$-linear form $P$ is obtained as a fixed linear combination of values of $p$ evaluated at finite sums of the vectors $X_1,\dots,X_k$.
Fix $g\in G$ and $X_1,\dots,X_k\in\mathfrak g$. The adjoint map
\begin{align*}
\operatorname{Ad}_g:\mathfrak g\to\mathfrak g
\end{align*}
is linear, so every sum appearing in the polarization identity is carried to the corresponding sum of the adjoint-transformed vectors. For example, a sum $X_{i_1}+\cdots+X_{i_m}$ is sent to
\begin{align*}
\operatorname{Ad}_g(X_{i_1}+\cdots+X_{i_m})=\operatorname{Ad}_gX_{i_1}+\cdots+\operatorname{Ad}_gX_{i_m}.
\end{align*}
Since the theorem statement assumes $p(\operatorname{Ad}_g Z)=p(Z)$ for every $Z\in\mathfrak g$, every value of $p$ occurring in the polarization formula is unchanged after applying $\operatorname{Ad}_g$ to all slots. Therefore
\begin{align*}
P(\operatorname{Ad}_g X_1,\dots,\operatorname{Ad}_g X_k)=P(X_1,\dots,X_k).
\end{align*}
Thus the polarized form $P$ is invariant under the diagonal adjoint action of $G$.
[/guided]
[/step]
[step:Differentiate the multilinear invariance identity along $\exp(tY)$]
Let \begin{align*} \exp:\mathfrak g\to G \end{align*} denote the Lie group exponential map. Fix $Y,X_1,\dots,X_k\in\mathfrak g$. For each $j\in\{1,\dots,k\}$, define the smooth curve \begin{align*} \gamma_j:\mathbb R\to \mathfrak g \end{align*} by
\begin{align*}
\gamma_j(t)=\operatorname{Ad}_{\exp(tY)}X_j.
\end{align*}
By the definition of the infinitesimal adjoint action,
\begin{align*}
\gamma_j'(0)=[Y,X_j].
\end{align*}
Using the $\operatorname{Ad}$-invariance of $P$ proved above, define
\begin{align*}
F:\mathbb R\to \mathbb R
\end{align*}
by
\begin{align*}
F(t)=P(\gamma_1(t),\dots,\gamma_k(t)).
\end{align*}
For every $t\in\mathbb R$,
\begin{align*}
F(t)=P(X_1,\dots,X_k),
\end{align*}
so $F'(0)=0$. Since $P$ is multilinear, differentiating $F$ at $0$ gives
\begin{align*}
F'(0)=\sum_{j=1}^k P(\gamma_1(0),\dots,\gamma_{j-1}(0),\gamma_j'(0),\gamma_{j+1}(0),\dots,\gamma_k(0)).
\end{align*}
Substituting $\gamma_i(0)=X_i$ and $\gamma_j'(0)=[Y,X_j]$ yields
\begin{align*}
0=\sum_{j=1}^k P(X_1,\dots,X_{j-1},[Y,X_j],X_{j+1},\dots,X_k).
\end{align*}
[guided]
The goal is to convert global $\operatorname{Ad}$-invariance into its infinitesimal form. We do this by moving in $G$ along the one-parameter subgroup determined by $Y$. For each slot, define the smooth curve
\begin{align*}
\gamma_j:\mathbb R\to \mathfrak g
\end{align*}
by
\begin{align*}
\gamma_j(t)=\operatorname{Ad}_{\exp(tY)}X_j.
\end{align*}
The curve starts at $X_j$, because $\exp(0Y)=e$ and $\operatorname{Ad}_e$ is the identity map on $\mathfrak g$. Thus
\begin{align*}
\gamma_j(0)=X_j.
\end{align*}
Its derivative at $0$ is the infinitesimal adjoint action:
\begin{align*}
\gamma_j'(0)=[Y,X_j].
\end{align*}
Now package the invariant expression into a single real-valued function. Define
\begin{align*}
F:\mathbb R\to \mathbb R
\end{align*}
by
\begin{align*}
F(t)=P(\gamma_1(t),\dots,\gamma_k(t)).
\end{align*}
Because the previous step proved that the polarized form $P$ is invariant under the simultaneous adjoint action on all slots, we have
\begin{align*}
F(t)=P(\operatorname{Ad}_{\exp(tY)}X_1,\dots,\operatorname{Ad}_{\exp(tY)}X_k)=P(X_1,\dots,X_k)
\end{align*}
for every $t\in\mathbb R$. Hence $F$ is constant, and therefore
\begin{align*}
F'(0)=0.
\end{align*}
It remains to compute this same derivative using multilinearity. Since $P$ is $k$-linear, the derivative of $P(\gamma_1(t),\dots,\gamma_k(t))$ is obtained by differentiating one slot at a time while holding the other slots fixed. At $t=0$ this gives
\begin{align*}
F'(0)=\sum_{j=1}^k P(\gamma_1(0),\dots,\gamma_{j-1}(0),\gamma_j'(0),\gamma_{j+1}(0),\dots,\gamma_k(0)).
\end{align*}
Substituting the two identities $\gamma_i(0)=X_i$ and $\gamma_j'(0)=[Y,X_j]$ gives
\begin{align*}
F'(0)=\sum_{j=1}^k P(X_1,\dots,X_{j-1},[Y,X_j],X_{j+1},\dots,X_k).
\end{align*}
Since $F'(0)=0$, this is exactly
\begin{align*}
\sum_{j=1}^k P(X_1,\dots,X_{j-1},[Y,X_j],X_{j+1},\dots,X_k)=0.
\end{align*}
[/guided]
[/step]
[step:Conclude the infinitesimal invariance identity]
The elements $Y,X_1,\dots,X_k\in\mathfrak g$ were arbitrary. Therefore the displayed identity holds for every $Y,X_1,\dots,X_k\in\mathfrak g$, completing the proof.
[/step]