[proofplan]
We separate the degenerate case $\sigma = 0$ from the non-degenerate case $\sigma > 0$. In the degenerate case, the normal distribution is the point mass at $\mu$, so the moment is immediate. In the non-degenerate case, we write the expectation using the Gaussian density, make the affine change of variables $x=\mu+\sigma z$, and reduce the problem to the finiteness of the standard Gaussian absolute moment. The final step proves this finiteness directly by showing that polynomial growth is dominated by the Gaussian exponential tail.
[/proofplan]
[step:Handle the degenerate Gaussian as a point mass]
Assume first that $\sigma = 0$. By the standard convention for the degenerate normal distribution, $X \sim \mathcal N(\mu,0)$ means
\begin{align*}
\mathbb P(X=\mu)=1.
\end{align*}
Let $k \ge 1$ be an integer. Since $|X|^k=|\mu|^k$ almost surely with respect to $\mathbb P$, the definition of expectation gives
\begin{align*}
\mathbb E[|X|^k]=|\mu|^k<\infty.
\end{align*}
Thus the desired conclusion holds when $\sigma=0$.
[/step]
[step:Reduce the moment to the standard Gaussian absolute moment]
Assume now that $\sigma>0$, and let $k \ge 1$ be an integer. Let $\mathcal L^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb R$. Since $X\sim \mathcal N(\mu,\sigma^2)$, define its density $f_X:\mathbb R\to [0,\infty)$ by
\begin{align*}
f_X(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
\end{align*}
This density satisfies, by the density formula for expectations of non-negative Borel [measurable functions](/page/Measurable%20Functions), for every non-negative Borel [measurable function](/page/Measurable%20Function) $g:\mathbb R\to[0,\infty]$,
\begin{align*}
\mathbb E[g(X)] = \int_{\mathbb R} g(x) f_X(x)\, d\mathcal L^1(x).
\end{align*}
Applying this with the Borel measurable function $g:\mathbb R\to [0,\infty)$ defined by $g(x)=|x|^k$ gives
\begin{align*}
\mathbb E[|X|^k]
=
\frac{1}{\sigma\sqrt{2\pi}}
\int_{\mathbb R}|x|^k
\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\, d\mathcal L^1(x).
\end{align*}
Use the affine substitution $x=\mu+\sigma z$. By the one-dimensional change-of-variables theorem, since $\sigma>0$, this substitution is a $C^1$ bijection from $\mathbb R$ to $\mathbb R$, and the one-dimensional Lebesgue measure transforms as $d\mathcal L^1(x)=\sigma\,d\mathcal L^1(z)$. Therefore
\begin{align*}
\mathbb E[|X|^k]
=
\frac{1}{\sqrt{2\pi}}
\int_{\mathbb R}|\mu+\sigma z|^k e^{-z^2/2}\, d\mathcal L^1(z).
\end{align*}
For all $a,b\in\mathbb R$ and all integers $k\ge1$, convexity of $t\mapsto t^k$ on $[0,\infty)$ gives
\begin{align*}
|a+b|^k \le 2^{k-1}(|a|^k+|b|^k).
\end{align*}
Applying this with $a=\mu$ and $b=\sigma z$ gives
\begin{align*}
|\mu+\sigma z|^k
\le
2^{k-1}\left(|\mu|^k+\sigma^k |z|^k\right).
\end{align*}
Hence
\begin{align*}
\mathbb E[|X|^k]
\le
\frac{2^{k-1}}{\sqrt{2\pi}}
\left(
|\mu|^k\int_{\mathbb R}e^{-z^2/2}\,d\mathcal L^1(z)
+
\sigma^k\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)
\right).
\end{align*}
Thus it remains only to prove that
\begin{align*}
\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)<\infty.
\end{align*}
[guided]
The goal is to isolate the only possible source of divergence: the polynomial factor in the tails. Let $\mathcal L^1$ denote one-dimensional Lebesgue measure on $\mathbb R$. Because $\sigma>0$, the normal law has the density $f_X:\mathbb R\to [0,\infty)$ defined by
\begin{align*}
f_X(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
\end{align*}
The function whose expectation we want is the Borel measurable function $g:\mathbb R\to [0,\infty)$ defined by $g(x)=|x|^k$.
It is Borel measurable and non-negative, so the density formula for expectations of non-negative Borel measurable functions applies to $g$ and gives
\begin{align*}
\mathbb E[|X|^k]
=
\frac{1}{\sigma\sqrt{2\pi}}
\int_{\mathbb R}|x|^k
\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\, d\mathcal L^1(x).
\end{align*}
Now we remove the location and scale parameters. Define the affine change of variables by $x=\mu+\sigma z$. The one-dimensional change-of-variables theorem applies because $\sigma>0$, so this map is a bijective $C^1$ change of variables from $\mathbb R$ onto $\mathbb R$, and its derivative has absolute value $\sigma$. Therefore one-dimensional Lebesgue measure transforms as $d\mathcal L^1(x)=\sigma\,d\mathcal L^1(z)$. The exponent becomes
\begin{align*}
\frac{(x-\mu)^2}{2\sigma^2}
=
\frac{z^2}{2},
\end{align*}
and so
\begin{align*}
\mathbb E[|X|^k]
=
\frac{1}{\sqrt{2\pi}}
\int_{\mathbb R}|\mu+\sigma z|^k e^{-z^2/2}\, d\mathcal L^1(z).
\end{align*}
The remaining issue is that $|\mu+\sigma z|^k$ is not exactly $|z|^k$. We bound it by separating the fixed shift $\mu$ from the scaled variable $\sigma z$. For all $a,b\in\mathbb R$, the triangle inequality gives $|a+b|\le |a|+|b|$, and convexity of $t\mapsto t^k$ on $[0,\infty)$ gives
\begin{align*}
(|a|+|b|)^k
\le
2^{k-1}(|a|^k+|b|^k).
\end{align*}
Using $a=\mu$ and $b=\sigma z$, we obtain
\begin{align*}
|\mu+\sigma z|^k
\le
2^{k-1}\left(|\mu|^k+\sigma^k |z|^k\right).
\end{align*}
Substituting this estimate into the moment integral yields
\begin{align*}
\mathbb E[|X|^k]
\le
\frac{2^{k-1}}{\sqrt{2\pi}}
\left(
|\mu|^k\int_{\mathbb R}e^{-z^2/2}\,d\mathcal L^1(z)
+
\sigma^k\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)
\right).
\end{align*}
The first integral is the total mass factor of the standard Gaussian density and is finite. Therefore the only remaining point is to prove
\begin{align*}
\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)<\infty.
\end{align*}
[/guided]
[/step]
[step:Dominate the polynomial factor by the Gaussian tail]
Let $m\in\mathbb N$ be an integer satisfying $2m>k$. For every $z\in\mathbb R$, the exponential series with non-negative terms gives
\begin{align*}
e^{z^2/4}\ge \frac{1}{m!}\left(\frac{z^2}{4}\right)^m.
\end{align*}
If $|z|\ge 1$, then
\begin{align*}
|z|^k e^{-z^2/2}
\le
4^m m! |z|^{k-2m} e^{-z^2/4}
\le
4^m m! e^{-z^2/4},
\end{align*}
because $k-2m<0$ and $|z|\ge1$. On the bounded interval $[-1,1]$,
\begin{align*}
|z|^k e^{-z^2/2}\le 1.
\end{align*}
Therefore
\begin{align*}
\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)
\le
2+
4^m m!
\int_{\{z\in\mathbb R: |z|\ge1\}} e^{-z^2/4}\,d\mathcal L^1(z).
\end{align*}
It remains to note that the last integral is finite. The interval $[1,\infty)$ is covered by the intervals $[n,n+1]$ with integers $n\ge1$. For $z\in[n,n+1]$ with an integer $n\ge1$, we have $z^2\ge n^2$, so
\begin{align*}
\int_n^{n+1} e^{-z^2/4}\,d\mathcal L^1(z)
\le
e^{-n^2/4}.
\end{align*}
Since $n^2\ge n$ for $n\in\mathbb N$, the geometric-series comparison gives
\begin{align*}
\sum_{n=1}^{\infty}e^{-n^2/4}
\le
\sum_{n=1}^{\infty}e^{-n/4}
<\infty.
\end{align*}
The same estimate applies on the negative half-line by the evenness of $z\mapsto e^{-z^2/4}$. Hence
\begin{align*}
\int_{\{z\in\mathbb R: |z|\ge1\}} e^{-z^2/4}\,d\mathcal L^1(z)<\infty,
\end{align*}
and consequently
\begin{align*}
\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)<\infty.
\end{align*}
[/step]
[step:Conclude the finiteness of the Gaussian polynomial moment]
From the preceding step,
\begin{align*}
\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)<\infty.
\end{align*}
Also,
\begin{align*}
\int_{\mathbb R}e^{-z^2/2}\,d\mathcal L^1(z)<\infty,
\end{align*}
which follows from the same tail argument with $k=0$. Substituting these two finite quantities into the estimate
\begin{align*}
\mathbb E[|X|^k]
\le
\frac{2^{k-1}}{\sqrt{2\pi}}
\left(
|\mu|^k\int_{\mathbb R}e^{-z^2/2}\,d\mathcal L^1(z)
+
\sigma^k\int_{\mathbb R}|z|^k e^{-z^2/2}\,d\mathcal L^1(z)
\right)
\end{align*}
shows that $\mathbb E[|X|^k]<\infty$. Since $k\ge1$ was arbitrary, every polynomial absolute moment of $X$ is finite.
[/step]