Here's a sketch of a solution.
First, ignore the scale factor of $\sigma = 2^{-n+1}$ because it can be restored at the end. Because the characteristic function of a Gamma$(1/2)$ variable like $X$ or $Y$ is
$$\psi_X(t) = \psi_Y(t) = (1 - it)^{-1/2},$$
the characteristic function of $X-Y$ (a version of the Fourier Transform of its density) is
$$\psi_{X-Y}(t) = \psi_X(t)\psi_Y(-t) = (1+t^2)^{-1/2}.$$
The density function of $X - Y$ therefore is recovered by applying the inverse Fourier transform, giving
$$f_{X-Y}(x) = \widehat{\psi_{X-Y}}(t) = \frac{1}{\pi}K_0(|x|),$$
a multiple of the Bessel $K_0$ function.
Because all events of the form $|X-Y|\le t$ for $t\ge 0$ correspond to events $-t \le X-Y \le t$ and the symmetry (and continuity) of $f_{X-Y}$ show the events $-t \le X-Y\le 0$ and $0 \le X-Y \le t$ must have equal probabilities, the density of $|X-Y|$ is twice the density of $X-Y$ for all positive values (and otherwise is zero, of course). Thus
$$f_{|X-Y|}(x) = \frac{2}{\pi}K_0(|x|).$$
The mean therefore is
$$E\left[\,|X-Y|\,\right] = \int_0^\infty x f_{|X-Y|}(x)\,\mathrm{d}x = \int_0^\infty x \frac{2}{\pi}K_0(|x|)\,\mathrm{d}x = \frac{2}{\pi}.$$
The raw second moment similarly is
$$E\left[\,|X-Y|^2\,\right] = \int_0^\infty x^2 \frac{2}{\pi}K_0(|x|)\,\mathrm{d}x = 1.$$
(We could have obtained this earlier from $\psi_{X-Y}$ because the second moment of $X-Y$ is that of $|X-Y|$ and
$$\psi_{X-Y}(t) = (1+t^2)^{-1/2} = 1 + \binom{-1/2}{1}(t^2)^1 + \cdots = 1 - \frac{t^2}{2} + \cdots$$
and the second moment is, as always, $2!$ times the coefficient of $(it)^2.$)
Therefore
$$\operatorname{Var}(|X-Y|) = 1 - \left(\frac{2}{\pi}\right)^2 \approx 0.5947153.$$
Finally, when $X$ and $Y$ are multiplied by the same constant $\sigma,$ $|X-Y|$ is also multiplied by $\sigma,$ whence its mean is $(2/\pi)\sigma$ and its variance is $\left(1-(2/\pi)^2\right)\sigma^2.$
Results of a simulation of $|X-Y|$ compared to the theoretical density (shown in red):
This R code conducted the simulation and produced the graphic.
n <- 1e6
z <- abs(rgamma(n, 1/2) - rgamma(n, 1/2))
hist(z, xlim=c(0,5), breaks=200, freq=FALSE, col="Gray", main="Realizations of |X-Y|",
xlab=expression(abs(X-Y)))
curve(besselK(x, 0) * 2/pi, add=TRUE, col="Red", lwd=2, n=1501)
(cbind(Mean=c(Simulated=mean(z), Theoretical=2/pi),
Variance=c(var(z), 1 - (2/pi)^2)))
