This is a rather good challenge for numerical analysis. The good thing is that the parameter $a$ belongs in $\mathbb{R}^+$.
Well, a "not so basic" knowledge of Special Bessel Functions can help in getting the result. However we live in the world of programming hence why don't ask to technology for a help? As you already have gotten, W. Mathematical can compute that integral rather quickly, giving:
$$\boxed{\frac{2 (a-1) a K\left(\frac{a+1}{2 a+1}\right)+(8 a+4) E\left(\frac{a+1}{2 a+1}\right)}{3 \sqrt{a (2 a+1)}}}$$
Where $K$ and $E$ are the Elliptic integrals of, respectively, the first kind and second kind.
Knowing this result shall be "enough" for you, in order to understand the behaviour of this function. We can invoke series expansions and asymptotic analysis for a rather good study of the behaviour of this function, say we call it
$$A(a) = \frac{2 (a-1) a K\left(\frac{a+1}{2 a+1}\right)+(8 a+4) E\left(\frac{a+1}{2 a+1}\right)}{3 \sqrt{a (2 a+1)}}$$
First of all we observe that plotting the function already gives us a really good amount of information. I am going to show you two plots, one for $0 < a < 10$ and the other one for $0 < a < 2$:


Numerically, but also analytically, we deduce that
$$\lim_{a\to 0^+} A(a) = +\infty$$
$$\lim_{a\to +\infty}A(a) = +\infty$$
Yet we recognise the existence of a minimum for the function at a certain point $a^*$
$$0.3 < a^* < 0.7$$
The computation of that minimum is not really easy, but we can manage it via Mathematica: proceeding with the usual method for minimising a function, id est $A'(a) = 0$, this gives us
$$a^* \approx 0.671254$$
As we expected from the plot. (We did not expect THAT value, but it lies within the range we estimated).
In this way, the minimum reads
$$\text{min}A(a) \to A(a = a^*) \approx 2.82678$$
Which fits.
So I guess that the problem here is just to track down a behavioural analysis of the function, and for those situations plots help as well as a bit knowledge of programming tools (or W. Mathematica too).
For the Elliptic integrals are not easy to treat, it's not always a good idea to sink down in series expansions, or asymptotic expansion, although literature is "quite" full of Elliptic integral series expansions for small arguments, or asymptotic expansions for big arguments.
Here, yet, we have a constraint on $a$: $a > 0$. We clearly never deal with infinites, especially in physical applications (I wonder whether this integral came from a physical calculation, maybe in quantum cosmology?) but now you are aware of the fact that your integral exists as long as $a > 0$. Consider also that the function grows in a not really fast way. To say:
$$A(0.2) \to 3.5288$$
$$A(1) \to 2.9125$$
$$A(2) \to 3.5593$$
$$A(5) \to 6.0314$$
$$A(13) \to 12.9584$$
$$A(31) \to 28.6664$$
$$\ldots$$
$$A(101) \to 89.8355$$
$$A(1001) \to 876.448$$
$$A(10001) \to 8742. 62$$
$$\ldots$$
$$A(1000001) \to 874022$$
If we take a curious look at this we can observe a very interesting fact:
$$\text{as} ~~~ n ~~~ \text{goes to infinity},~~~ \dfrac{n}{A(n)} \approx 1.14$$
Pushed by the curiosity of this I computed the exact limit, finding:
$$\lim_{n\to +\infty} \dfrac{n}{A(n)} = \frac{6 \Gamma \left(\frac{3}{4}\right)}{\sqrt{\pi } \Gamma \left(\frac{1}{4}\right)} \approx 1.4414$$
Perhaps, due to the fact that it's more or less constant the ratio between $A(n)$ and $n$ for every real $n$ (or $a$, to turn back to our initial parameter), we are in good hands to find a sort of exponential asymptotic approximation for the function, like
$$A(a) \sim b^{g(a)} ~~~~~ a > 1$$
For some basis $1 < b < 2$.
Which works better and better the bigger is $a$.
I believe we could work a lot on this interesting integral. I do not know if I answered you in a proper way, maybe I just wrote down a lot of useless numerical information. I only hope that all of this, somehow, helps you!