Heuristics
Since $e^{-|x|^a}$ behaves like $1-|x|^a$ plus better behaved terms close to zero, and since formally the Fourier trasform of $|x|^a$ is proportional to $|\xi|^{-a-n}$, I would expect that the decay at infinity of your object is something like $|\xi|^{-a-n}$, which by the way gives the correct prediction for $a=1$ (in that case, the Fourier transform is explicit and should be $(1+|x|^2)^{-(n+1)/2}$, up to rescaling and multiplicative constants).
So, the answer should be that it lies in $L^1(\mathbb R^n)$ for all $a>0$ (and actually in any $L^p$ with $n/(n+a)<p\leq\infty$). You can most likely turn it into a proof, but I did not work out all the details.
An idea for a different proof of $L^p$ estimates.
In a paper by Weinstein et al. it is proved the estimate for the fractional chain rule:
$$ \|D^s F\circ u\|_{L^r}\lesssim \|G\circ u\|_{L^p}\|D^s u\|_{L^q} $$
for $0<s<1$, $1<p,q,r<\infty$, $1/p+1/q=1/r$, $F$ a continuous (maybe locally Lipschitz) function s.t. $F(0)=0$, and $G$ is a non-negative function such that
$$ |F(x)-F(y)|\lesssim (G(x)+G(y))|x-y| $$
(often, you can take $G=|F’|$).
You can also find the statement on some lecture notes of Harmonic Analysis by Monica Visan that are available online (chapter 22).
Then, you can write $\phi (x)e^{-|x|^a}$ as a composite function $F\circ u$, where
$$ F(t)=\exp(-1/t)\eta(t), $$
$$ u(x)=|x|^{-a}\chi(x), $$
where $\eta$ is a smooth cutoff function that equals $0$ for $t<2\epsilon$ and equals $1$ for $t>4\epsilon$, and $\chi$ is a cutoff function on $\mathbb R^n$ like $\phi$, that equals $1$ for $|x|<(\epsilon)^{-1/a}$.
Now you can try to use the above estimate with some $r\leq 2$ (the larger the better) and prove that $\|D^s\phi e^{-|x|^a}\|_{L^r}$ is finite. By Young inequality, this is implies a weighted estimate for the Fourier transform of your function:
$$ \||\xi|^{s}\widehat{\phi e^{-|x|^a}}\|_{L^{r’}}<\infty. $$
Then you can turn that weighted estimate into an estimate in $L^q$ with $q<r’$ by Hölder’s inequality :)
It could work, since computing fractional derivatives of $|x|^a$ is probably not so hard. I did not try to apply this argument because the calculations are a little messy, and I don’t know if you can go all the way to $L^1$, but if you have enough patience you could try to apply this method to a given example with fixed $n$ and $a$. There are surely better proofs for this simple case, but this method works for much more general compositions of functions, without symmetry or smoothness.