[proofplan] We represent the multivariate normal vector as an affine image of a standard normal vector. The main calculation is the moment generating function of a standard Gaussian vector, obtained by completing the square in the Gaussian density and using translation invariance of Lebesgue measure. Substituting the affine representation then converts the Euclidean norm $|A^\top t|^2$ into the quadratic form $t^\top \Sigma t$. [/proofplan] [step:Represent the Gaussian vector as an affine image of a standard normal vector] Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $\Sigma$ is symmetric and positive semidefinite, the [positive semidefinite square root construction](/page/Positive%20Semidefinite%20Square%20Root) gives a symmetric positive semidefinite matrix $A \in \mathbb{R}^{p \times p}$ such that $A A^\top=\Sigma$. Let \begin{align*} Z: (\Omega_0,\mathcal F_0,\mathbb P_0) &\to (\mathbb{R}^p,\mathcal B(\mathbb{R}^p)) \end{align*} be a standard normal random vector with distribution $\mathcal N_p(0,I_p)$. Let $\mathbb E_0$ denote expectation with respect to the probability measure $\mathbb P_0$. By the defining affine representation in the [multivariate normal distribution](/page/Multivariate%20Normal%20Distribution), the random vectors $X$ and $\mu+AZ$ have the same distribution. Therefore, for every $t \in \mathbb{R}^p$, \begin{align*} \mathbb E[\exp(t^\top X)] = \mathbb E_0[\exp(t^\top(\mu+AZ))], \end{align*} provided either side is finite. [guided] The covariance matrix $\Sigma$ may be singular, so we avoid an argument that depends on a non-degenerate density for $X$. Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $\Sigma$ is symmetric and positive semidefinite, the [positive semidefinite square root construction](/page/Positive%20Semidefinite%20Square%20Root) gives a symmetric positive semidefinite matrix $A \in \mathbb{R}^{p \times p}$ satisfying \begin{align*} A A^\top=\Sigma. \end{align*} This matrix is a square root of $\Sigma$. Let \begin{align*} Z: (\Omega_0,\mathcal F_0,\mathbb P_0) &\to (\mathbb{R}^p,\mathcal B(\mathbb{R}^p)) \end{align*} be a standard normal random vector, meaning $Z \sim \mathcal N_p(0,I_p)$. Let $\mathbb E_0$ denote expectation with respect to the probability measure $\mathbb P_0$. The defining affine construction in the [multivariate normal distribution](/page/Multivariate%20Normal%20Distribution) says that $\mu+AZ$ has distribution $\mathcal N_p(\mu,AA^\top)=\mathcal N_p(\mu,\Sigma)$. Since $X$ also has distribution $\mathcal N_p(\mu,\Sigma)$, the random vectors $X$ and $\mu+AZ$ have the same law. Hence every non-negative Borel function has the same expectation under the two laws. Applying this to the non-negative Borel function \begin{align*} g_t: \mathbb{R}^p &\to (0,\infty) \\ x &\mapsto \exp(t^\top x), \end{align*} we obtain \begin{align*} \mathbb E[\exp(t^\top X)] = \mathbb E_0[\exp(t^\top(\mu+AZ))], \end{align*} with finiteness to be established in the next step. [/guided] [/step] [step:Compute the moment generating function of a standard normal vector] Fix $a \in \mathbb{R}^p$. Let $\mathcal L^p$ denote $p$-dimensional Lebesgue measure on $\mathbb{R}^p$. The density of $Z$ with respect to $\mathcal L^p$ is \begin{align*} \phi_p: \mathbb{R}^p &\to (0,\infty) \\ z &\mapsto (2\pi)^{-p/2}\exp\left(-\frac{1}{2}|z|^2\right). \end{align*} Thus \begin{align*} \mathbb E_0[\exp(a^\top Z)] &= \int_{\mathbb{R}^p}\exp(a^\top z)\phi_p(z)\,d\mathcal L^p(z) \\ &= (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(a^\top z-\frac{1}{2}|z|^2\right) \,d\mathcal L^p(z). \end{align*} Completing the square gives \begin{align*} a^\top z-\frac{1}{2}|z|^2 = -\frac{1}{2}|z-a|^2+\frac{1}{2}|a|^2. \end{align*} Therefore \begin{align*} \mathbb E_0[\exp(a^\top Z)] &= \exp\left(\frac{1}{2}|a|^2\right) (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|z-a|^2\right) \,d\mathcal L^p(z). \end{align*} Using the translation $y=z-a$, for which $d\mathcal L^p(y)=d\mathcal L^p(z)$ and $\mathbb{R}^p-a=\mathbb{R}^p$, the remaining integral is the total mass of the standard normal density: \begin{align*} (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|z-a|^2\right) \,d\mathcal L^p(z) = (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|y|^2\right) \,d\mathcal L^p(y) = 1. \end{align*} Hence \begin{align*} \mathbb E_0[\exp(a^\top Z)] = \exp\left(\frac{1}{2}|a|^2\right)<\infty. \end{align*} [guided] We now compute the exponential moment of the standard Gaussian vector directly from its density. Fix a vector $a \in \mathbb{R}^p$. Since $Z \sim \mathcal N_p(0,I_p)$, its density with respect to $p$-dimensional Lebesgue measure $\mathcal L^p$ is the function \begin{align*} \phi_p: \mathbb{R}^p &\to (0,\infty) \\ z &\mapsto (2\pi)^{-p/2}\exp\left(-\frac{1}{2}|z|^2\right). \end{align*} Therefore the expectation is the [Lebesgue integral](/page/Lebesgue%20Integral) \begin{align*} \mathbb E_0[\exp(a^\top Z)] &= \int_{\mathbb{R}^p}\exp(a^\top z)\phi_p(z)\,d\mathcal L^p(z) \\ &= (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(a^\top z-\frac{1}{2}|z|^2\right) \,d\mathcal L^p(z). \end{align*} The reason this integral is finite for every real vector $a$ is that the linear term $a^\top z$ can be absorbed into the Gaussian quadratic term. Completing the square gives, for every $z \in \mathbb{R}^p$, \begin{align*} a^\top z-\frac{1}{2}|z|^2 = -\frac{1}{2}|z-a|^2+\frac{1}{2}|a|^2. \end{align*} Substituting this identity into the integral yields \begin{align*} \mathbb E_0[\exp(a^\top Z)] &= \exp\left(\frac{1}{2}|a|^2\right) (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|z-a|^2\right) \,d\mathcal L^p(z). \end{align*} Now apply the translation change of variables \begin{align*} y=z-a. \end{align*} Lebesgue measure is translation invariant, so $d\mathcal L^p(y)=d\mathcal L^p(z)$, and the domain $\mathbb{R}^p$ is mapped onto $\mathbb{R}^p$. Hence \begin{align*} (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|z-a|^2\right) \,d\mathcal L^p(z) &= (2\pi)^{-p/2} \int_{\mathbb{R}^p} \exp\left(-\frac{1}{2}|y|^2\right) \,d\mathcal L^p(y) \\ &=1, \end{align*} because $\phi_p$ is a probability density. Therefore \begin{align*} \mathbb E_0[\exp(a^\top Z)] = \exp\left(\frac{1}{2}|a|^2\right), \end{align*} and this number is finite for every $a \in \mathbb{R}^p$. [/guided] [/step] [step:Substitute the affine representation and simplify the quadratic form] Fix $t \in \mathbb{R}^p$ and define \begin{align*} a:=A^\top t \in \mathbb{R}^p. \end{align*} Using $t^\top(\mu+AZ)=t^\top\mu+(A^\top t)^\top Z$, the previous step gives \begin{align*} M_X(t) &= \mathbb E_0[\exp(t^\top(\mu+AZ))] \\ &= \exp(t^\top\mu)\mathbb E_0[\exp(a^\top Z)] \\ &= \exp(t^\top\mu)\exp\left(\frac{1}{2}|a|^2\right). \end{align*} Finally, \begin{align*} |a|^2 = |A^\top t|^2 = (A^\top t)^\top(A^\top t) = t^\top A A^\top t = t^\top\Sigma t. \end{align*} Thus \begin{align*} M_X(t) = \exp\left(t^\top\mu+\frac{1}{2}t^\top\Sigma t\right). \end{align*} Since $t \in \mathbb{R}^p$ was arbitrary, the formula holds for every $t \in \mathbb{R}^p$. [/step]