Say, we have a sequence of probability distributions $\mu_n$ on $\mathbb R$, that are uniformly subgaussian in the sense that $$\mu_n(\mathbb R\setminus[-R,R])\leq Ce^{-CR^2}$$ for some positive constant $C$. Assume that all moments of these distributions exist and that they converge to those of a standard normal distribution, i.e., $$\lim_{n\to\infty}\int x^k~d\mu_n(x)=\begin{cases}(k-1)!!&\text{if }k\text{ is even}\\ 0 &\text{else}\end{cases}.$$ Then $\mu_n$ converges weakly to a standard-normal distribution $\mu$, i.e., $$\lim_{n\to\infty} \int f~d\mu_n=\int f~d\mu$$ for all bounded continuous $f$.
The above claim is clearly true: The uniform subgaussian hypothesis implies that the $k$-th moments are uniformly bounded by $(Ck)^{k/2}$ and therefore the characteristic functions can be written as $$\int e^{itx}~d\mu_n(x)=\sum_{k=0}^\infty \frac{(it)^k}{k!}\int x^k~d\mu_n(x),$$ which converges point wise to the characteristic function of the normal distribution. The result now follows from Levy's Continuity Theorem.
Question: Is there a way to prove the above statement directly? (I am concerned with something related where the measures are random themselves and I can't argue with the characteristic function)
Thoughts: Fix $\epsilon>0$ and $f\in C_b(\mathbb R)$ and choose $\lambda<\infty$ such that $\mu_n(\mathbb R\setminus[-\lambda,\lambda])\leq \epsilon$. From the Weierstraß approximation theorem, we find a polynomial $p$ that approximates $f$ on the interval $[-\lambda,\lambda]$ up to an error of $\epsilon$. Thus we find $$\big\lvert\int f~d\mu_n-\int f~d\mu\big\lvert\leq \int_{[-\lambda,\lambda]}\lvert f-p\lvert~d\mu_n+\big\lvert\int_{[-\lambda,\lambda]^c} f-p~d\mu_n\big\lvert+\big\lvert\int p~d(\mu_n-\mu)\big\lvert+\int \lvert p-f\lvert~d\mu.$$ The first term is bounded by $\epsilon$ and the third term converges to zero by assumption. The fourth term can be bound by $\epsilon$ and the integral outside the finite interval. Since $f$ is bounded we can estimate $\lvert f(x)-p(x)\lvert\leq c x^{2k}$ outside the interval for some $c,k\in\mathbb N$. It would therefore remain to show that $$\int_{[-\lambda,\lambda]^c} x^{2k}~d(\mu_n+\mu)$$ is small. Does this somehow follow from the subgaussian assumption (which wasn't really used till now)?