Let $X_{1},...,X_{N}$ be nondegenerate independent Gaussian random vectors from a probability space $(\Omega, \mathcal{F},P)$ to $\mathbb{R}^{n}$. Assume that each $X_{j}$ has its associate Gaussian distribution $\mu_{C_{j}}$, where: \begin{eqnarray} d\mu_{C_{j}}(x) = \frac{1}{\sqrt{(2\pi)^{n}\det C}}e^{-\frac{1}{2}\langle x, C_{j}^{-1}x\rangle}d\lambda^{n}(x) \tag{1}\label{1} \end{eqnarray} where $\lambda^{n}$ is the Lebesgue measure on $\mathbb{R}^{n}$ and each $C_{j}$ is a positive-definite matrix. I know for a fact that if we define a random vector $\tilde{X} := X_{1}+\cdots + X_{N}$, then $\tilde{X}$ will be Gaussian with distribution $\mu_{C}$, where $C = C_{1}+\cdots + C_{N}$, but I don't know how to prove it. I tried using the definition of distribution and using Characteristic functions but I gor nowhere. Can anyone help me?
EDIT: In the comments, there is some ideas on how to prove the result. Let me add one more ingredient. It is easily proved that, with the above notation, if $f: \Omega \to \mathbb{R}$ is measurable and $\mu_{C}$-integrable then: \begin{eqnarray} \int d\mu_{C}(x)f(x) = \int d\mu_{C_{1}}(x_{1})\int d\mu_{C_{2}}(x_{2}) \cdots \int d\mu_{C_{N}}(x_{N})f(x_{1}+\cdots+x_{N}) \tag{2}\label{2} \end{eqnarray} And (\ref{2}) seems to be an indicative of the result of the question. Any ideas on how to use (\ref{2}) to prove it?