[proofplan]
We rewrite the empirical-process vector as the normalized sum of independent identically distributed centred random vectors in $\mathbb R^k$. The finite second-moment assumptions imply that each coordinate is integrable and that the covariance matrix is finite. The [multivariate central limit theorem](/theorems/1854) then gives convergence to a centred Gaussian vector, and a direct computation identifies its covariance matrix with the displayed formula.
[/proofplan]
[step:Define the centred summands whose normalized sum is the empirical-process vector]
The theorem statement supplies the measurable space $(S,\mathcal A)$, the probability measure $P$ on $(S,\mathcal A)$, and the ambient probability space $(\Omega,\mathcal F,\mathbb P)$ on which the independent identically distributed random variables $X_m:(\Omega,\mathcal F)\to(S,\mathcal A)$ are defined. Let $\mathcal B(\mathbb R)$ denote the Borel $\sigma$-algebra on $\mathbb R$, and let $\mathcal B(\mathbb R^k)$ denote the Borel $\sigma$-algebra on $\mathbb R^k$. For any integrable real-valued [random variable](/page/Random%20Variable) $H:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$, write $\mathbb E[H]:=\int_\Omega H(\omega)\,d\mathbb P(\omega)$. For each $j\in\{1,\dots,k\}$, the assumption $P f_j^2<\infty$ implies $P|f_j|<\infty$ because $|a|\le 1+a^2$ for every $a\in\mathbb R$ and $P(S)=1$. Hence $P f_j$ is finite.
For each $m \in \mathbb N$, define the random vector $Y_m:(\Omega,\mathcal F)\to(\mathbb R^k,\mathcal B(\mathbb R^k))$ by $Y_m(\omega)=\bigl(f_1(X_m(\omega))-P f_1,\dots,f_k(X_m(\omega))-P f_k\bigr)$. Since each $f_j:(S,\mathcal A)\to(\mathbb R,\mathcal B(\mathbb R))$ is measurable and $X_m$ is measurable, each coordinate $f_j\circ X_m$ is real-valued measurable, so $Y_m$ is a Borel-measurable $\mathbb R^k$-valued random vector.
The random variables $X_1,X_2,\dots$ are independent and identically distributed with common law $P$. Therefore the random vectors $Y_1,Y_2,\dots$ are independent and identically distributed. If $Y_{m,j}:\Omega\to\mathbb R$ denotes the $j$-th coordinate of $Y_m$, then
\begin{align*}
\mathbb E[Y_{m,j}]=\int_\Omega \bigl(f_j(X_m(\omega))-P f_j\bigr)\,d\mathbb P(\omega)=\int_S f_j(x)\,dP(x)-P f_j=0.
\end{align*}
Finally, for every $n\in\mathbb N$, the identity of $\mathbb R^k$-valued random vectors is
\begin{align*}
\frac{1}{\sqrt n}\sum_{m=1}^n Y_m=\left(\sqrt n\left(\frac{1}{n}\sum_{m=1}^n f_1(X_m)-P f_1\right),\dots,\sqrt n\left(\frac{1}{n}\sum_{m=1}^n f_k(X_m)-P f_k\right)\right)=(\alpha_n(f_1),\dots,\alpha_n(f_k)).
\end{align*}
[guided]
The goal is to put the empirical-process vector into the exact form required by the multivariate [central limit theorem](/theorems/521): a normalized sum of independent identically distributed centred random vectors. The theorem statement supplies the measurable space $(S,\mathcal A)$, the probability measure $P$, and the ambient probability space $(\Omega,\mathcal F,\mathbb P)$. Let $\mathcal B(\mathbb R)$ denote the Borel $\sigma$-algebra on $\mathbb R$, and let $\mathcal B(\mathbb R^k)$ denote the Borel $\sigma$-algebra on $\mathbb R^k$. For any integrable real-valued random variable $H:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$, write $\mathbb E[H]:=\int_\Omega H(\omega)\,d\mathbb P(\omega)$. Thus every occurrence of $\mathbb E[H]$ below is integration with respect to $\mathbb P$, while every expression of the form $P f$ means integration over $S$ with respect to $P$.
Before subtracting $P f_j$, we verify that this number is finite. For each $j\in\{1,\dots,k\}$, the assumption $P f_j^2<\infty$ means
\begin{align*}
\int_S f_j(x)^2\,dP(x)<\infty.
\end{align*}
Since $|a|\le 1+a^2$ for every $a\in\mathbb R$ and $P(S)=1$, we have
\begin{align*}
\int_S |f_j(x)|\,dP(x)\le \int_S \bigl(1+f_j(x)^2\bigr)\,dP(x)<\infty.
\end{align*}
Therefore $P f_j=\int_S f_j(x)\,dP(x)$ is finite.
For each $m \in \mathbb N$, define $Y_m:(\Omega,\mathcal F)\to(\mathbb R^k,\mathcal B(\mathbb R^k))$ by $Y_m(\omega)=\bigl(f_1(X_m(\omega))-P f_1,\dots,f_k(X_m(\omega))-P f_k\bigr)$. This definition subtracts the population mean from each coordinate. Since $X_m:(\Omega,\mathcal F)\to(S,\mathcal A)$ is measurable and each $f_j:(S,\mathcal A)\to(\mathbb R,\mathcal B(\mathbb R))$ is measurable, the composition $f_j\circ X_m:\Omega\to\mathbb R$ is measurable. Hence every coordinate of $Y_m$ is measurable, and therefore $Y_m$ is a Borel-measurable $\mathbb R^k$-valued random vector.
Because $X_1,X_2,\dots$ are independent and identically distributed with law $P$, applying the same measurable map $x\mapsto\bigl(f_1(x)-P f_1,\dots,f_k(x)-P f_k\bigr)$ to each $X_m$ preserves independence and identical distribution. Thus $Y_1,Y_2,\dots$ are independent and identically distributed.
The centering is exact. If $Y_{m,j}:\Omega\to\mathbb R$ denotes the $j$-th coordinate of $Y_m$, then the law of $X_m$ is $P$, so
\begin{align*}
\mathbb E[Y_{m,j}]=\int_\Omega \bigl(f_j(X_m(\omega))-P f_j\bigr)\,d\mathbb P(\omega)=\int_S f_j(x)\,dP(x)-P f_j=0.
\end{align*}
Now compute the normalized sum coordinate by coordinate:
\begin{align*}
\frac{1}{\sqrt n}\sum_{m=1}^n Y_m=\left(\sqrt n\left(\frac{1}{n}\sum_{m=1}^n f_1(X_m)-P f_1\right),\dots,\sqrt n\left(\frac{1}{n}\sum_{m=1}^n f_k(X_m)-P f_k\right)\right)=(\alpha_n(f_1),\dots,\alpha_n(f_k)).
\end{align*}
This identity is the bridge from empirical-process notation to the finite-dimensional [central limit theorem](/theorems/1848).
[/guided]
[/step]
[step:Verify the finite covariance hypotheses]
For each $j \in \{1,\dots,k\}$, the first step has already shown that $P f_j$ is finite. For $i,j \in \{1,\dots,k\}$, the elementary inequality $2|ab|\le a^2+b^2$ gives
\begin{align*}
\int_S |f_i(x)f_j(x)|\,dP(x)\le\frac{1}{2}\int_S f_i(x)^2\,dP(x)+\frac{1}{2}\int_S f_j(x)^2\,dP(x)<\infty.
\end{align*}
Thus every coordinate product $Y_{1,i}Y_{1,j}$ is integrable, and the covariance matrix of $Y_1$ is finite.
[/step]
[step:Compute the covariance matrix of one summand]
Define the matrix $\Sigma \in \mathbb R^{k \times k}$ by
\begin{align*}
\Sigma_{ij}:=\mathbb E[Y_{1,i}Y_{1,j}], \qquad i,j \in \{1,\dots,k\}.
\end{align*}
Since $\mathbb E[Y_{1,i}]=\mathbb E[Y_{1,j}]=0$, this is the covariance matrix of $Y_1$. Using the law of $X_1$, for all $i,j \in \{1,\dots,k\}$ we obtain
\begin{align*}
\Sigma_{ij}=\int_\Omega \bigl(f_i(X_1(\omega))-P f_i\bigr)\bigl(f_j(X_1(\omega))-P f_j\bigr)\,d\mathbb P(\omega)=\int_S \bigl(f_i(x)-P f_i\bigr)\bigl(f_j(x)-P f_j\bigr)\,dP(x).
\end{align*}
Expanding the product and using $\int_S1\,dP(x)=1$ because $P$ is a probability measure gives
\begin{align*}
\Sigma_{ij}=\int_S f_i(x)f_j(x)\,dP(x)-(P f_i)\int_S f_j(x)\,dP(x)-(P f_j)\int_S f_i(x)\,dP(x)+(P f_i)(P f_j)=\int_S f_i(x)f_j(x)\,dP(x)-(P f_i)(P f_j).
\end{align*}
[/step]
[step:Apply the multivariate central limit theorem and identify the limit]
The sequence $Y_1,Y_2,\dots$ consists of independent identically distributed centred $\mathbb R^k$-valued random vectors with finite covariance matrix $\Sigma$. Therefore, by the multivariate central limit theorem applied to the independent identically distributed centred $\mathbb R^k$-valued random vectors $Y_m$ with finite covariance matrix $\Sigma$,
\begin{align*}
\frac{1}{\sqrt n}\sum_{m=1}^n Y_m \xrightarrow{d} Z,
\end{align*}
where $Z$ is a centred Gaussian random vector in $\mathbb R^k$ with covariance matrix $\Sigma$.
By the identity proved above,
\begin{align*}
\frac{1}{\sqrt n}\sum_{m=1}^n Y_m
=
(\alpha_n(f_1),\dots,\alpha_n(f_k)).
\end{align*}
Writing
\begin{align*}
Z=(G(f_1),\dots,G(f_k)),
\end{align*}
we obtain
\begin{align*}
(\alpha_n(f_1),\dots,\alpha_n(f_k)) \xrightarrow{d} (G(f_1),\dots,G(f_k)).
\end{align*}
The covariance computation gives, for every $i,j \in \{1,\dots,k\}$,
\begin{align*}
\operatorname{Cov}(G(f_i),G(f_j))
=
\Sigma_{ij}
=
\int_S f_i(x)f_j(x)\,dP(x)
-
\left(\int_S f_i(x)\,dP(x)\right)
\left(\int_S f_j(x)\,dP(x)\right).
\end{align*}
This is exactly the asserted finite-dimensional empirical-process limit.
[/step]