It looks to me like $ \int_{-\infty}^\infty f_Z(x)^2\ dx = -{\frac {{\sigma}}{2{a}^{2}\sqrt {\pi }}}+\frac{1}{2a} {{\rm erf}\left({\frac {a}{{\sigma}}}\right)}+ \frac{\sigma}{2a^2 \sqrt{\pi}} {{\rm e}^{-a^2/\sigma^2}} $
EDIT: OK, here's the proof.
For convenience, scale distances so that $\sigma = 1$. We're looking at
$J =\frac{1}{8 \pi a^2} \int_{-\infty}^\infty dx \int_{x-a}^{x+a} ds \int_{x-a}^{x+a} dt\ e^{-(s^2+t^2)/2} $
Interchange the order of integration so this becomes
$\eqalign{ J &= \frac{1}{8 \pi a^2} \int_{-\infty}^\infty ds \int_{s-2a}^{s+2a} dt \int_{\max(s,t)-a}^{\min(s,t)+a} dx\ e^{-(s^2+t^2)/2}\cr &= \frac{1}{8 \pi a^2} \int_{-\infty}^\infty ds \int_{s-2a}^{s+2a} dt\ (2a + \min(s,t) -\max(s,t)) e^{-(s^2+t^2)/2} \cr}$
Break this up into two pieces, one where s and the other where $s>t$.
$ \eqalign{J_1 &= \frac{1}{8 \pi a^2} \int_{-\infty}^\infty ds \int_s^{s+2a} dt\ (2a + s-t) e^{-(s^2+t^2)/2}\cr J_2 &= \frac{1}{8 \pi a^2} \int_{-\infty}^\infty ds \int_{s-2a}^{s} dt\ (2a + t-s) e^{-(s^2+t^2)/2}\cr}$
In $J_1$, take $t = s+u$ (so that $-(s^2+t^2)/2 = -s^2 - us - t^2/2$); in $J_2$, take $t = s - u$ (so that $-(s^2+t^2)/2 = -s^2 +us - t^2/2$). With these changes of variables we recombine the integrals:
$ J = \frac{1}{8 \pi a^2} \int_{-\infty}^\infty ds \int_0^{2a} du\ (2a - u) e^{-s^2} (e^{-us} + e^{us}) e^{-u^2/2} $
Interchange the order of integration again
$ \eqalign{J &= \frac{1}{8 \pi a^2} \int_0^{2a} du \int_{-\infty}^\infty ds\ (2a-u) e^{-s^2} (e^{-us}+e^{us}) e^{-u^2/2} \cr &= \frac{1}{2\sqrt{\pi} a^2} \int_0^{2a} du\ (a-u/2) e^{-u^2/4}\cr &= \frac{1}{2a} \text{erf}(a) + \frac{1}{2 \sqrt{\pi} a^2} e^{-a^2} - \frac{1}{2 \sqrt{\pi} a^2}\cr}$