Let $X_n$ be a sequence of Gaussian random vectors that converges in distribution to some random vector $X$. Is $X$ Gaussian?
Partial solution
If we can show that $\mu_n := E\left(X_n\right)$ and $\Sigma_n := V\left(X_n\right)$ converge (to $\mu$ and $\Sigma$, respectively), we're done, because then the characteric functions $\phi_n = \phi_{X_n}$ converge to
$$
\phi\left(x\right) = \phi_X\left(x\right) = \exp\left(-\frac{1}{2}x^T\Sigma x + i\ x^T\mu\right)
$$
which implies that $X$ is Gaussian with $E(X) = \mu$ and $V(X) = \Sigma$.
Now, since $\exp\left(-\frac{1}{2}x^T\Sigma_n x\right) = \left|\phi_n(x)\right| \rightarrow \left|\phi(x)\right|$, then whenever $\phi(x) \neq 0$ we have $\nu(x):=\lim_{n \rightarrow \infty}x^T\Sigma_n x \in \mathbb{R}$. But since $\phi(0) = 1$ and since $\phi$ is continuous (as are all characteristic functions), there is some neighborhood of $0$, where $x^T\Sigma_n X$ converges. If $x$ is outside this neighborhood, it can be represented as $x = \alpha y$ for some $y$ in this neighborhood, and then $x^T\Sigma_n x \rightarrow \alpha^2 \nu(y)$.
It is now possible by a judicious selection of $x$'s to show that $\Sigma_n$ converges. All that's left is to show that $\mu_n$ converges.
Best Answer
Firstly, we establish some notation. We will denote the $n$th vector of the sequence thus: $X^n$, rather than thus: $X_n$ as in the original post. We will denote by $N$ the common dimension of $X$ and the $X^n$'s, i.e. $X, X^1, X^2, \dots \in \mathbb{R}_{N \times 1}$. For every $k \in \left\{1, \dots, N\right\}$, $X_k$ will denote the $k$th component of $X$, i.e. $X = \left[X_1, \dots, X_N\right]^T$, and $X^n_k$ will denote the $k$th component of $X^n$, i.e. $X^n = \left[X^n_1, \dots, X^n_N\right]^T$.
With this notation, we would like to prove that the sequence $\left(\mu_n\right)_{n = 1}^\infty$ converges, i.e. that the sequence $\left(E\left(X^n\right)\right)_{n = 1}^\infty$ converges, i.e. that for every $k \in \left\{1, \dots, N\right\}$, the sequence $\left(E\left(X^n_k\right)\right)_{n = 1}^\infty$ converges.
Since by assumption $X^n \overset{D}{\longrightarrow} X$, according to the Cramér–Wold theorem for every $k \in \left\{1, \dots, N\right\}$ we have $X^n_k \overset{D}{\longrightarrow}X_k$. But for every $n \in \mathbb{N}_1$ and every $k \in \left\{1, \dots, N\right\}$, $X^n_k$ is (one-dimensional) Gaussian.
Hence, we have reduced the problem of proving that the sequence $\left(\mu_n\right)_{n = 1}^\infty$ converges to the following one-dimensional case.
Theorem
Let $\left(Y_n\right)_{n = 1}^\infty$ be a sequence of normally distributed random variables (not necessarily defined over a single probability space) that converges in distribution to the random variable $Z$. Then the sequence $\left(E\left(Y_n\right)\right)_{n = 1}^\infty$ converges.
Lemma 1
$(Y_n)_{n = 1}^\infty$ is tight, i.e. for each $\varepsilon \in \left(0, \infty\right)$, there is $T \in \mathbb{R}$, such that for each $n \in \mathbb{N}_1$, $P\left(\left|Y_n\right| > T\right) < \varepsilon$.
Proof (Davide Giraudo)
By the continuous mapping theorem, $\left|Y_n\right| \overset{D}{\longrightarrow} \left|Z\right|$. Denote $Z$'s cdf by $F$ and, for each $n \in \mathbb{N}_1$, denote $\left|Y_n\right|$'s cdf by $F_n$. By Froda's theorem, $F$ has a continuity point $t_1 \in \mathbb{R}$, such that $F(t_1)\gt 1-\varepsilon$. Then there is some $N \in \mathbb{N}_1$, such that for all $n \in \left\{N+1, N+2, \dots\right\}$, $F_n(t_1)\gt 1-\varepsilon$, i.e. $P\left(\left|Y_n\right| > t_1\right) < \varepsilon$. Since every finite collection of random variables is finite, there is some $t_2 \in \mathbb{R}$, such that for all $n \in \left\{1, \dots, N\right\}$, $P\left(\left|Y_n\right| > t_2\right) < \varepsilon$. Define $T := \max\left(t_1, t_2\right)$. Q.E.D.
Lemma 2
$\left\{E\left(Y_n\right)\ :\mid\ n \in \mathbb{N}_1\right\}$ is bounded.
Proof (D. Thomine)
For every $n \in \mathbb{N}_1$, define $\nu_n := \left|E\left(Y_n\right)\right|$. Suppose $\sup \left\{\nu_n\ :\mid\ n \in \mathbb{N}_1\right\} = \infty$ and for every $T \in \mathbb{R}$ let $n_T \in \mathbb{N}_1$ be such that $\nu_{n_T} \geq T$. Since a normal distribution is symmetric, $P\left(\left|Y_n\right| > \nu_n\right) \geq \frac{1}{2}$ for all $n \in \mathbb{N}_1$. Therefore, for any $T \in \mathbb{R}$, $P\left(\left|Y_{n_T}\right| > T\right) \geq P\left(\left|Y_{n_T}\right| > n_T\right) \geq \frac{1}{2}$, which contradicts tightness. Q.E.D.
Proof of the theorem
For each $n \in \mathbb{N}_1$, define $m_n := E\left(Y_n\right)$ and $s^2_n := V\left(Y_n\right)$ and denote $Y_n$'s characteristic function by $\varphi_n$. Denote $Z$'s characteristic function by $\varphi$. So $\lim_{n \rightarrow \infty} \varphi_n = \varphi$ and for each $n \in \mathbb{N}_1$, $\varphi_n = t \in \mathbb{R} \mapsto \exp\left(-\frac{1}{2}s^2_n t^2 + i m_n t\right)$. As OP (= me) showed in the original post, $\left(s_n\right)_{n = 1}^\infty$ converges. Define $s := \lim_{n \rightarrow \infty} s_n$. Then for every $t \in \mathbb{R}$, $\lim_{n \rightarrow \infty} e^{i m_n t} = \frac{\varphi\left(t\right)}{\exp\left(-\frac{1}{2}s^2t^2\right)} =: \rho\left(t\right) \in \mathbb{C}$.
We wish to show that $\left(m_n\right)_{n = 1}^\infty$ converges. Define $M := 1 + \sup \left\{\nu_n\ :\mid\ n \in \mathbb{N}_1\right\}$. Then with $t := \frac{\pi}{M}$ and $\alpha_n := tm_n$ we have $\lim_{n \rightarrow \infty} e^{i \alpha_n} = \rho\left(\frac{\pi}{M}\right)$. In particular, $\rho\left(\frac{\pi}{M}\right) \neq 0$. But since $\alpha_n \in \left(-\pi, \pi\right)$ for every $n \in \mathbb{N}_1$, and since the complex exponential function is injective over the domain $\left\{a + i b\ :\mid\ a \in \mathbb{R}, b \in \left(-\pi, \pi\right)\right\}$ with inverse $\textrm{Log}$ (the principal branch of the complex logarithm; cf. Definition 3.1.3, p. 40 in Ash, Robert B. and Novinger, W. P. "Complex Variables", 2nd edition, Dover Publications, 2004), we have $\lim_{n \rightarrow \infty} m_n = \frac{M}{i \pi}\textrm{Log}\left(\rho\left(\frac{\pi}{M}\right)\right)$.
Q.E.D.