$\def\si{\sigma}$
$\def\cF{\mathcal{F}}$
$\def\ol{\overline}$
Here is an alternative derivation. Let $B_t=(\mu/\si)W_{\si^2t/\mu^2}$, $Y_t=t+B_t$, and $b=a\mu/\si^2$. Then $X_t=a$ if and only if $Y_{\mu^2t/\si^2}=b$. Hence, if $T$ is the last time $X$ hits $a$ and $\tau$ is the last time $Y$ hits $b$, then $T=\si^2\tau/\mu^2$.
By the Markov property,
\begin{align}
P(\tau < t) &= P(Y_t > b, Y_s > b \text{ for all $s\ge t$})\\
&= E[P(Y_t > b, Y_s > b \text{ for all $s\ge t$} \mid \cF_t)]\\
&= E[1_{\{Y_t > b\}}P^{Y_t}(\tau_b = \infty)],
\end{align}
where $\tau_b$ is the first hitting time of $b$. For $y>b$,
$$
P^y(\tau_b = \infty) = P(\tau_{b-y} = \infty) = 1 - e^{-2(y-b)}.
$$
Thus,
\begin{align}
P(\tau < t) &= P(Y_t > b)
- \frac1{\sqrt{2\pi t}}\int_b^\infty\exp\left({
-2(y - b) - \frac{(y - t)^2}{2t}
}\right)\,dy\\
&= P(Y_t > b)
- \frac1{\sqrt{2\pi t}}\int_b^\infty\exp\left({
2b - \frac{(y + t)^2}{2t}
}\right)\,dy\\
&= P(t + B_t > b) - e^{2b}P(-t + B_t > b)\\
&= \ol\Phi\left({\frac{b-t}{\sqrt t}}\right)
- e^{2b}\ol\Phi\left({\frac{b+t}{\sqrt t}}\right),
\end{align}
where $\ol\Phi=1-\Phi$, and $\Phi$ is the standard normal distribution function. Differentiating and doing a little algebra gives the density of $\tau$,
$$
f_\tau(t) = \frac1{\sqrt{2\pi t}}\exp\left({-\frac{(b-t)^2}{2t}}\right),
$$
and from here, we get the density of $T$,
$$
f_T(t) = \frac{\mu}{\si\sqrt{2\pi t}}\exp\left({
-\frac{(\si^2b-\mu^2t)^2}{2\mu^2\si^2t}
}\right).
$$
Not that you needed another derivation to tell you this, but there must be something wrong with your numerical simulations.