[proofplan]
We compute the characteristic function directly from the density. First translate $x$ by $\mu$, which separates the deterministic phase factor $e^{i t^\top \mu}$. Then diagonalise the positive definite covariance matrix and use the linear change of variables that reduces the quadratic form to $|z|^2$. The remaining integral factors into one-dimensional standard normal Fourier transforms, which are computed by a short integration-by-parts argument and then reassembled into $e^{-t^\top \Sigma t/2}$.
[/proofplan]
[step:Translate the density integral to center the Gaussian at the origin]
Fix $t \in \mathbb{R}^p$. By the density formula for $X$, the characteristic function satisfies
\begin{align*}
\phi_X(t)
&= \int_{\mathbb{R}^p} e^{i t^\top x} f_X(x)\, d\mathcal{L}^p(x) \\
&= \frac{1}{(2\pi)^{p/2}(\det \Sigma)^{1/2}}
\int_{\mathbb{R}^p}
e^{i t^\top x}
\exp\left(-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)\right)
\, d\mathcal{L}^p(x).
\end{align*}
Apply the translation change of variables $y=x-\mu$, so that $x=y+\mu$ and $d\mathcal{L}^p(x)=d\mathcal{L}^p(y)$. The domain $\mathbb{R}^p$ is invariant under translation. Hence
\begin{align*}
\phi_X(t)
&= e^{i t^\top \mu}
\frac{1}{(2\pi)^{p/2}(\det \Sigma)^{1/2}}
\int_{\mathbb{R}^p}
e^{i t^\top y}
\exp\left(-\frac{1}{2}y^\top \Sigma^{-1}y\right)
\, d\mathcal{L}^p(y).
\end{align*}
[guided]
Fix $t \in \mathbb{R}^p$. The characteristic function is an expectation of the complex-valued random variable $\omega \mapsto e^{i t^\top X(\omega)}$. Since $X$ has density $f_X$ with respect to $\mathcal{L}^p$, this expectation is the [Lebesgue integral](/page/Lebesgue%20Integral)
\begin{align*}
\phi_X(t)
&= \int_{\mathbb{R}^p} e^{i t^\top x} f_X(x)\, d\mathcal{L}^p(x) \\
&= \frac{1}{(2\pi)^{p/2}(\det \Sigma)^{1/2}}
\int_{\mathbb{R}^p}
e^{i t^\top x}
\exp\left(-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)\right)
\, d\mathcal{L}^p(x).
\end{align*}
The purpose of the first change of variables is to remove the mean from the quadratic term. Define the translated variable $y=x-\mu$, equivalently $x=y+\mu$. Translation preserves $p$-dimensional Lebesgue measure, so $d\mathcal{L}^p(x)=d\mathcal{L}^p(y)$, and it maps $\mathbb{R}^p$ onto $\mathbb{R}^p$. Under this substitution,
\begin{align*}
t^\top x = t^\top(y+\mu)=t^\top y+t^\top\mu.
\end{align*}
Therefore the factor depending only on $\mu$ separates from the integral:
\begin{align*}
\phi_X(t)
&= e^{i t^\top \mu}
\frac{1}{(2\pi)^{p/2}(\det \Sigma)^{1/2}}
\int_{\mathbb{R}^p}
e^{i t^\top y}
\exp\left(-\frac{1}{2}y^\top \Sigma^{-1}y\right)
\, d\mathcal{L}^p(y).
\end{align*}
This reduces the problem to the centered [Gaussian integral](/theorems/1140).
[/guided]
[/step]
[step:Diagonalize the covariance matrix and reduce the integral to standard normal form]
Since $\Sigma$ is real symmetric and positive definite, the spectral theorem for real symmetric matrices applies (citing a result not yet in the wiki: Spectral Theorem for Real Symmetric Matrices). Thus there exist an orthogonal matrix $Q \in \mathbb{R}^{p \times p}$ and positive numbers $\lambda_1,\dots,\lambda_p \in (0,\infty)$ such that
\begin{align*}
\Sigma = Q\Lambda Q^\top,
\qquad
\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_p).
\end{align*}
Define the positive diagonal square root
\begin{align*}
\Lambda^{1/2}=\operatorname{diag}(\lambda_1^{1/2},\dots,\lambda_p^{1/2})
\end{align*}
and the invertible [linear map](/page/Linear%20Map)
\begin{align*}
A: \mathbb{R}^p &\to \mathbb{R}^p \\
z &\mapsto Q\Lambda^{1/2}z.
\end{align*}
Apply the multivariable change-of-variables theorem (citing a result not yet in the wiki: Change of Variables Theorem) with $y=A z$. Since $|\det A|=(\det \Sigma)^{1/2}$, $d\mathcal{L}^p(y)=(\det \Sigma)^{1/2}\,d\mathcal{L}^p(z)$. Moreover,
\begin{align*}
y^\top \Sigma^{-1}y
&= z^\top (\Lambda^{1/2})^\top Q^\top Q\Lambda^{-1}Q^\top Q\Lambda^{1/2}z \\
&= z^\top z,
\end{align*}
and, defining
\begin{align*}
s := \Lambda^{1/2}Q^\top t \in \mathbb{R}^p,
\end{align*}
we have $t^\top y=s^\top z$. Hence
\begin{align*}
\phi_X(t)
&= e^{i t^\top \mu}
\frac{1}{(2\pi)^{p/2}}
\int_{\mathbb{R}^p} e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)\, d\mathcal{L}^p(z).
\end{align*}
[guided]
The centered integral still contains the covariance matrix $\Sigma$. The goal is to choose coordinates in which the quadratic form $y^\top\Sigma^{-1}y$ becomes the Euclidean square $|z|^2$.
Because $\Sigma$ is real, symmetric, and positive definite, the spectral theorem for real symmetric matrices applies (citing a result not yet in the wiki: Spectral Theorem for Real Symmetric Matrices). It gives an orthogonal matrix $Q \in \mathbb{R}^{p \times p}$ and positive eigenvalues $\lambda_1,\dots,\lambda_p \in (0,\infty)$ such that
\begin{align*}
\Sigma = Q\Lambda Q^\top,
\qquad
\Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_p).
\end{align*}
The positivity of the eigenvalues is exactly where positive definiteness is used. Define
\begin{align*}
\Lambda^{1/2}=\operatorname{diag}(\lambda_1^{1/2},\dots,\lambda_p^{1/2})
\end{align*}
and define the linear map
\begin{align*}
A: \mathbb{R}^p &\to \mathbb{R}^p \\
z &\mapsto Q\Lambda^{1/2}z.
\end{align*}
This map is invertible because each $\lambda_j^{1/2}$ is positive and $Q$ is invertible.
Now apply the multivariable change-of-variables theorem (citing a result not yet in the wiki: Change of Variables Theorem) using $y=A z=Q\Lambda^{1/2}z$. The domain $\mathbb{R}^p$ maps onto $\mathbb{R}^p$ because $A$ is invertible. The Jacobian factor is
\begin{align*}
|\det A|
&= |\det Q|\,\det(\Lambda^{1/2})
= \prod_{j=1}^p \lambda_j^{1/2}
= (\det \Sigma)^{1/2},
\end{align*}
since $|\det Q|=1$ for an orthogonal matrix and $\det \Sigma=\prod_{j=1}^p \lambda_j$. Therefore
\begin{align*}
d\mathcal{L}^p(y)=(\det \Sigma)^{1/2}\,d\mathcal{L}^p(z),
\end{align*}
so the Jacobian cancels the normalising factor $(\det\Sigma)^{-1/2}$ in the density.
The quadratic form becomes
\begin{align*}
y^\top \Sigma^{-1}y
&= (Q\Lambda^{1/2}z)^\top (Q\Lambda^{-1}Q^\top)(Q\Lambda^{1/2}z) \\
&= z^\top(\Lambda^{1/2})^\top Q^\top Q\Lambda^{-1}Q^\top Q\Lambda^{1/2}z \\
&= z^\top \Lambda^{1/2}\Lambda^{-1}\Lambda^{1/2}z \\
&= z^\top z.
\end{align*}
For the oscillatory factor, define
\begin{align*}
s := \Lambda^{1/2}Q^\top t \in \mathbb{R}^p.
\end{align*}
Then
\begin{align*}
t^\top y=t^\top Q\Lambda^{1/2}z=(\Lambda^{1/2}Q^\top t)^\top z=s^\top z.
\end{align*}
Substituting these identities into the centered integral gives
\begin{align*}
\phi_X(t)
&= e^{i t^\top \mu}
\frac{1}{(2\pi)^{p/2}}
\int_{\mathbb{R}^p} e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)\, d\mathcal{L}^p(z).
\end{align*}
[/guided]
[/step]
[step:Factor the standard normal integral into one-dimensional Fourier transforms]
Write $z=(z_1,\dots,z_p)$ and $s=(s_1,\dots,s_p)$. Since
\begin{align*}
s^\top z=\sum_{j=1}^p s_j z_j,
\qquad
z^\top z=\sum_{j=1}^p z_j^2,
\end{align*}
the integrand factors as
\begin{align*}
\frac{1}{(2\pi)^{p/2}}e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)
=
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}e^{i s_j z_j}\exp\left(-\frac{1}{2}z_j^2\right)
\right].
\end{align*}
The absolute value of this product is the product of the one-dimensional standard normal densities, whose integral over $\mathbb{R}^p$ is $1$. Therefore [Fubini's theorem](/theorems/2961) applies (citing a result not yet in the wiki: Fubini's Theorem), and
\begin{align*}
\frac{1}{(2\pi)^{p/2}}
\int_{\mathbb{R}^p} e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)\, d\mathcal{L}^p(z)
=
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}
\int_{\mathbb{R}} e^{i s_j u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
\right].
\end{align*}
[guided]
After the linear change of variables, the covariance matrix has disappeared. What remains is a product structure. Write
\begin{align*}
z=(z_1,\dots,z_p),
\qquad
s=(s_1,\dots,s_p).
\end{align*}
Then the dot product and Euclidean square split coordinatewise:
\begin{align*}
s^\top z=\sum_{j=1}^p s_j z_j,
\qquad
z^\top z=\sum_{j=1}^p z_j^2.
\end{align*}
Consequently,
\begin{align*}
e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)
&=
\prod_{j=1}^p e^{i s_j z_j}\exp\left(-\frac{1}{2}z_j^2\right),
\end{align*}
and including the normalising factor gives
\begin{align*}
\frac{1}{(2\pi)^{p/2}}e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)
=
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}e^{i s_j z_j}\exp\left(-\frac{1}{2}z_j^2\right)
\right].
\end{align*}
To justify turning the $p$-dimensional integral into a product of one-dimensional integrals, we check absolute integrability. The absolute value of the integrand is
\begin{align*}
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}\exp\left(-\frac{1}{2}z_j^2\right)
\right],
\end{align*}
which is the product of $p$ one-dimensional standard normal densities. Its integral over $\mathbb{R}^p$ is $1$. Thus Fubini's theorem applies (citing a result not yet in the wiki: Fubini's Theorem), and the integral factors:
\begin{align*}
\frac{1}{(2\pi)^{p/2}}
\int_{\mathbb{R}^p} e^{i s^\top z}\exp\left(-\frac{1}{2}z^\top z\right)\, d\mathcal{L}^p(z)
=
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}
\int_{\mathbb{R}} e^{i s_j u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
\right].
\end{align*}
[/guided]
[/step]
[step:Compute the one-dimensional standard normal Fourier transform]
For $a \in \mathbb{R}$, define
\begin{align*}
G: \mathbb{R} &\to \mathbb{C} \\
a &\mapsto \frac{1}{(2\pi)^{1/2}}\int_{\mathbb{R}} e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u).
\end{align*}
The function $G$ is differentiable, with differentiation under the integral justified by domination by the integrable function $u \mapsto (2\pi)^{-1/2}|u|\exp(-u^2/2)$. Thus
\begin{align*}
G'(a)
&= \frac{1}{(2\pi)^{1/2}}\int_{\mathbb{R}} i u e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u).
\end{align*}
Since
\begin{align*}
\frac{d}{du}\left[\exp\left(-\frac{1}{2}u^2\right)\right]
=
-u\exp\left(-\frac{1}{2}u^2\right),
\end{align*}
[integration by parts](/theorems/2098) on $\mathbb{R}$ gives
\begin{align*}
G'(a)
&= -\frac{i}{(2\pi)^{1/2}}
\int_{\mathbb{R}} e^{i a u}
\frac{d}{du}\left[\exp\left(-\frac{1}{2}u^2\right)\right]
\, d\mathcal{L}^1(u) \\
&= -\frac{i}{(2\pi)^{1/2}}
\left(
0-\int_{\mathbb{R}} i a e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
\right) \\
&= -aG(a).
\end{align*}
Also $G(0)=1$ by the normalisation of the one-dimensional standard normal density. Therefore $G$ solves the initial value problem $G'(a)=-aG(a)$ and $G(0)=1$, so
\begin{align*}
G(a)=\exp\left(-\frac{1}{2}a^2\right)
\end{align*}
for every $a \in \mathbb{R}$.
[guided]
We now compute the one-dimensional integral that appears in each coordinate. For $a \in \mathbb{R}$, define
\begin{align*}
G: \mathbb{R} &\to \mathbb{C} \\
a &\mapsto \frac{1}{(2\pi)^{1/2}}\int_{\mathbb{R}} e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u).
\end{align*}
The strategy is to show that $G$ satisfies a first-order differential equation. Differentiating the factor $e^{iau}$ with respect to $a$ introduces $iu$, and that factor can be converted into a derivative of the Gaussian density.
First justify differentiating under the integral. For each fixed compact interval $K \subset \mathbb{R}$ of $a$-values,
\begin{align*}
\left|\frac{\partial}{\partial a}\left(e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\right)\right|
=
|u|\exp\left(-\frac{1}{2}u^2\right),
\end{align*}
and the function $u \mapsto |u|\exp(-u^2/2)$ is integrable over $\mathbb{R}$ with respect to $\mathcal{L}^1$. Thus differentiation under the integral is valid, and
\begin{align*}
G'(a)
&= \frac{1}{(2\pi)^{1/2}}\int_{\mathbb{R}} i u e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u).
\end{align*}
The key identity is
\begin{align*}
\frac{d}{du}\left[\exp\left(-\frac{1}{2}u^2\right)\right]
=
-u\exp\left(-\frac{1}{2}u^2\right).
\end{align*}
Therefore
\begin{align*}
i u\exp\left(-\frac{1}{2}u^2\right)
=
-i\frac{d}{du}\left[\exp\left(-\frac{1}{2}u^2\right)\right].
\end{align*}
Substituting this into the formula for $G'(a)$ yields
\begin{align*}
G'(a)
&= -\frac{i}{(2\pi)^{1/2}}
\int_{\mathbb{R}} e^{i a u}
\frac{d}{du}\left[\exp\left(-\frac{1}{2}u^2\right)\right]
\, d\mathcal{L}^1(u).
\end{align*}
Integrating by parts over $\mathbb{R}$ is justified by first integrating over $[-R,R]$ and then letting $R \to \infty$; the boundary term $e^{iau}\exp(-u^2/2)$ tends to $0$ at both endpoints because the Gaussian factor decays to $0$. Hence
\begin{align*}
G'(a)
&= -\frac{i}{(2\pi)^{1/2}}
\left(
0-\int_{\mathbb{R}} i a e^{i a u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
\right) \\
&= -aG(a).
\end{align*}
At $a=0$,
\begin{align*}
G(0)
=
\frac{1}{(2\pi)^{1/2}}\int_{\mathbb{R}}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
=
1,
\end{align*}
by the normalisation of the one-dimensional standard normal density. Thus $G$ satisfies the initial value problem
\begin{align*}
G'(a)=-aG(a),
\qquad
G(0)=1.
\end{align*}
Solving this scalar differential equation gives
\begin{align*}
G(a)=\exp\left(-\frac{1}{2}a^2\right)
\end{align*}
for every $a \in \mathbb{R}$.
[/guided]
[/step]
[step:Reassemble the diagonalized variables to obtain the characteristic function]
Applying the previous step with $a=s_j$ for each $j \in \{1,\dots,p\}$ gives
\begin{align*}
\prod_{j=1}^p
\left[
\frac{1}{(2\pi)^{1/2}}
\int_{\mathbb{R}} e^{i s_j u}\exp\left(-\frac{1}{2}u^2\right)\, d\mathcal{L}^1(u)
\right]
&=
\prod_{j=1}^p \exp\left(-\frac{1}{2}s_j^2\right) \\
&= \exp\left(-\frac{1}{2}s^\top s\right).
\end{align*}
Because $s=\Lambda^{1/2}Q^\top t$,
\begin{align*}
s^\top s
&= t^\top Q\Lambda Q^\top t \\
&= t^\top \Sigma t.
\end{align*}
Substituting into the expression for $\phi_X(t)$ yields
\begin{align*}
\phi_X(t)
&= e^{i t^\top \mu}\exp\left(-\frac{1}{2}t^\top\Sigma t\right) \\
&= \exp\left(i 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]
