Proof of Mean-Square Limits of Gaussian Random Variables are Gaussian

Each $X_k$ is square-integrable by Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector.

Step 1 (Convergence of moments). Write $\mu_k=\mathbb{E}[X_k]$ and $\mu=\mathbb{E}[X]$. By the Cauchy-Schwarz inequality applied to $X_k-X$ and the constant random variable $1$,

By the triangle inequality of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm, $\bigl|\lVert X_k\rVert_{2}-\lVert X\rVert_{2}\bigr|\le\lVert X_k-X\rVert_{2}\to0$, so $\mathbb{E}[X_k^{2}]\to\mathbb{E}[X^{2}]$. Hence, with the variance identity $\operatorname{Var}(Y)=\mathbb{E}[Y^{2}]-\mathbb{E}[Y]^{2}$ recorded there,

This proves the two limit formulas; it remains to show $X$ is Gaussian.

Step 2 (Degenerate case). If $\sigma^{2}=0$, then $\mathbb{E}[(X-\mu)^{2}]=0$, so $P(X=\mu)=1$ by the null-equivalence clause of Square-Integrable Random Variables and the Mean-Square Inner Product. By Gaussian Random Vectors and Jointly Gaussian Random Variables with $m=0$ (empty sum), $X$ is Gaussian.

Step 3 (CDF identification in the nondegenerate case). Assume $\sigma^{2}>0$ and let $\sigma$ be its positive square root. Since $\operatorname{Var}(X_k)\to\sigma^{2}>0$, there is $K$ with $\sigma_k^{2}:=\operatorname{Var}(X_k)>0$ for $k\ge K$. By Markov's inequality applied to the nonnegative random variable $(X_k-X)^{2}$, for every $h>0$,

For $k\ge K$, Standardization and Cumulative Distribution Function of a Gaussian Random Variable gives $F_{X_k}(u)=\Phi\bigl((u-\mu_k)/\sigma_k\bigr)$ for the cumulative distribution function $F_{X_k}$, with $\Phi$ the standard normal cumulative distribution function, which is continuous at every point. Define $G(u)=\Phi\bigl((u-\mu)/\sigma\bigr)$; since $\mu_k\to\mu$, $\sigma_k\to\sigma>0$ (continuity of the square root at positive arguments, e.g. $|\sigma_k-\sigma|=|\sigma_k^{2}-\sigma^{2}|/(\sigma_k+\sigma)\le|\sigma_k^{2}-\sigma^{2}|/\sigma$), and $\Phi$ is continuous, we get $F_{X_k}(u)\to G(u)$ for every real $u$.

Now fix $u\in\mathbb{R}$ and $h>0$. From the inclusion $\{X\le u\}\subseteq\{X_k\le u+h\}\cup\{|X_k-X|>h\}$ and monotonicity and finite subadditivity of the measure $P$,

Similarly, from $\{X_k\le u-h\}\subseteq\{X\le u\}\cup\{|X_k-X|>h\}$ we get $F_{X_k}(u-h)-P(|X_k-X|>h)\le F_X(u)$, and letting $k\to\infty$ gives $G(u-h)\le F_X(u)$. Thus $G(u-h)\le F_X(u)\le G(u+h)$ for every $h>0$, and letting $h\downarrow0$ using the continuity of $G$ gives $F_X(u)=G(u)$ for every real $u$.

Step 4 (The standardized variable is standard normal). Set $Z=(X-\mu)/\sigma$, a random variable. For every real $u$,

Let $P_Z$ denote the distribution of $Z$ and $N$ the standard normal distribution; both are probability measures on the Borel $\sigma$-algebra. They agree on the family $\mathcal{P}=\{(-\infty,u]:u\in\mathbb{R}\}$, which is a $\pi$-system ($(-\infty,u]\cap(-\infty,v]=(-\infty,\min(u,v)]$). The collection $\Lambda=\{B:P_Z(B)=N(B)\}$ contains $\mathbb{R}$ (both measures have total mass $1$), is closed under proper differences (finite additivity) and under increasing countable unions (continuity from below of measures, a consequence of countable additivity), so it is a $\lambda$-system in the sense of Dynkin's Pi-Lambda Theorem. Moreover $\sigma(\mathcal{P})=\mathcal{B}(\mathbb{R})$: each $(-\infty,u]$ is Borel, and conversely every open interval satisfies $(a,b)=\bigcup_{n\ge1}\bigl((-\infty,b-(b-a)/(2n)]\setminus(-\infty,a]\bigr)$, every open set is a countable union of open intervals with rational endpoints (using the density of the rationals), and the open sets generate the Borel $\sigma$-algebra by Borel Sigma-Algebra on the Real Line. By Dynkin's theorem, $\Lambda\supseteq\mathcal{B}(\mathbb{R})$, so $P_Z=N$ and $Z$ is standard normal.

Step 5 (Gaussian representation). Pointwise $X=\mu+\sigma Z$, so $P(X=\mu+\sigma Z)=1$. This is a Gaussian representation of $(X)$ in the sense of Gaussian Random Vectors and Jointly Gaussian Random Variables with $d=1$, $m=1$, $\mu_1=\mu$, $a_{11}=\sigma$ (the independence requirement on a single standard normal variable is vacuous). Hence $X$ is Gaussian. $\square$