[proofplan]
We center the [Gaussian random variable](/page/Gaussian%20Random%20Variable) by setting $Y=X-\mu$ and compare it with its reflection $Z=\mu-X=-Y$. The density of $Y$ is the centered Gaussian density with scale $\sigma$, and that density is an even function. Therefore integrating over an arbitrary Borel set gives the same probability for $Y$ and $Z$, which is exactly equality in distribution.
[/proofplan]
[step:Center the Gaussian random variable and record its density]
Let $Y:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$ and $Z:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$ be the real-valued random variables defined by
\begin{align*}
Y=X-\mu
\end{align*}
and
\begin{align*}
Z=\mu-X.
\end{align*}
Since translations and reflections of $\mathbb R$ are Borel measurable, both $Y$ and $Z$ are Borel-measurable random variables.
By the definition of $X\sim\mathcal N(\mu,\sigma^2)$ with $\sigma>0$, the law of $X$ has density $f_X:\mathbb R\to[0,\infty)$ with respect to one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) $\mathcal L^1$, where
\begin{align*}
f_X(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
\end{align*}
Thus for every Borel set $B\in\mathcal B(\mathbb R)$,
\begin{align*}
\mathbb P(X\in B)=\int_B f_X(x)\,d\mathcal L^1(x).
\end{align*}
Define the centered Gaussian density $\varphi_\sigma:\mathbb R\to[0,\infty)$ by
\begin{align*}
\varphi_\sigma(y)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{y^2}{2\sigma^2}\right).
\end{align*}
[/step]
[step:Compute the probability of the translated event]
Fix a Borel set $A\in\mathcal B(\mathbb R)$. Define the translation map $\tau:\mathbb R\to\mathbb R$ by $\tau(y)=\mu+y$, and write $\mu+A:=\tau(A)$. Since $\tau$ is a homeomorphism, $\mu+A$ is Borel. The event $\{Y\in A\}$ equals $\{X\in\mu+A\}$, so
\begin{align*}
\mathbb P(Y\in A)=\int_{\mu+A} f_X(x)\,d\mathcal L^1(x).
\end{align*}
Using the one-dimensional Lebesgue change-of-variables formula for the translation $x=\mu+y$ (citing a result not yet in the wiki: one-dimensional Lebesgue change-of-variables formula), whose absolute Jacobian is $1$, we obtain
\begin{align*}
\mathbb P(Y\in A)=\int_A f_X(\mu+y)\,d\mathcal L^1(y).
\end{align*}
For every $y\in\mathbb R$,
\begin{align*}
f_X(\mu+y)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{y^2}{2\sigma^2}\right)=\varphi_\sigma(y).
\end{align*}
Therefore
\begin{align*}
\mathbb P(Y\in A)=\int_A \varphi_\sigma(y)\,d\mathcal L^1(y).
\end{align*}
[guided]
We want to express the probability of the event $\{X-\mu\in A\}$ as an integral over the set $A$ itself. The reason for doing this is that the symmetry will become visible only after the Gaussian is centered at $0$.
Define the translation map $\tau:\mathbb R\to\mathbb R$ by $\tau(y)=\mu+y$, and define $\mu+A:=\tau(A)$. The map $\tau$ is a homeomorphism of $\mathbb R$, so it maps Borel sets to Borel sets. Hence $\mu+A\in\mathcal B(\mathbb R)$. For $\omega\in\Omega$, the condition $Y(\omega)\in A$ means $X(\omega)-\mu\in A$, which is equivalent to $X(\omega)\in\mu+A$. Thus the two events are equal:
\begin{align*}
\{Y\in A\}=\{X\in\mu+A\}.
\end{align*}
Since the law of $X$ has density $f_X$ with respect to $\mathcal L^1$, this gives
\begin{align*}
\mathbb P(Y\in A)=\int_{\mu+A} f_X(x)\,d\mathcal L^1(x).
\end{align*}
Now apply the one-dimensional Lebesgue change-of-variables formula to the translation $x=\mu+y$ (citing a result not yet in the wiki: one-dimensional Lebesgue change-of-variables formula). The map $y\mapsto\mu+y$ carries $A$ onto $\mu+A$, and its absolute Jacobian is $1$. Therefore
\begin{align*}
\mathbb P(Y\in A)=\int_A f_X(\mu+y)\,d\mathcal L^1(y).
\end{align*}
Substituting $x=\mu+y$ into the Gaussian density gives
\begin{align*}
f_X(\mu+y)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(\mu+y-\mu)^2}{2\sigma^2}\right).
\end{align*}
Since $\mu+y-\mu=y$, this becomes
\begin{align*}
f_X(\mu+y)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{y^2}{2\sigma^2}\right)=\varphi_\sigma(y).
\end{align*}
Hence
\begin{align*}
\mathbb P(Y\in A)=\int_A \varphi_\sigma(y)\,d\mathcal L^1(y).
\end{align*}
[/guided]
[/step]
[step:Compute the probability of the reflected event using evenness]
Define the reflection-translation map $\rho:\mathbb R\to\mathbb R$ by $\rho(y)=\mu-y$, and write $\mu-A:=\rho(A)$. Since $\rho$ is a homeomorphism, $\mu-A$ is Borel. The event $\{Z\in A\}$ equals $\{X\in\mu-A\}$, so
\begin{align*}
\mathbb P(Z\in A)=\int_{\mu-A} f_X(x)\,d\mathcal L^1(x).
\end{align*}
Using the one-dimensional Lebesgue change-of-variables formula for the reflection $x=\mu-y$ (citing a result not yet in the wiki: one-dimensional Lebesgue change-of-variables formula), whose absolute Jacobian is $1$, we obtain
\begin{align*}
\mathbb P(Z\in A)=\int_A f_X(\mu-y)\,d\mathcal L^1(y).
\end{align*}
For every $y\in\mathbb R$,
\begin{align*}
f_X(\mu-y)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(-y)^2}{2\sigma^2}\right)=\varphi_\sigma(y).
\end{align*}
Therefore
\begin{align*}
\mathbb P(Z\in A)=\int_A \varphi_\sigma(y)\,d\mathcal L^1(y).
\end{align*}
[/step]
[step:Conclude equality of probabilities and equality in distribution]
The preceding two steps show that for the fixed Borel set $A\in\mathcal B(\mathbb R)$,
\begin{align*}
\mathbb P(X-\mu\in A)=\mathbb P(Y\in A)=\int_A \varphi_\sigma(y)\,d\mathcal L^1(y).
\end{align*}
They also show
\begin{align*}
\mathbb P(\mu-X\in A)=\mathbb P(Z\in A)=\int_A \varphi_\sigma(y)\,d\mathcal L^1(y).
\end{align*}
Thus
\begin{align*}
\mathbb P(X-\mu\in A)=\mathbb P(\mu-X\in A).
\end{align*}
Since $A\in\mathcal B(\mathbb R)$ was arbitrary, the laws of $X-\mu$ and $\mu-X$ agree on every Borel set. By the definition of equality in distribution, $X-\mu$ and $\mu-X$ have the same distribution.
[/step]
