Evaluating this integral using the Gamma function
$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\Li}[1]{\,{\rm Li}_{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ \begin{align}{\cal I}&=\ \overbrace{% \color{#66f}{\large\int_0^{\infty}t^{-1/2}\exp\pars{-a\bracks{t + t^{-1}}}\,\dd t}} ^{\ds{\color{#c00000}{t \equiv x^{2}}}}\ =\ 2\int_{0}^{\infty}\exp\pars{-a\bracks{x^{2} + {1 \over x^{2}}}}\,\dd x \end{align}
\begin{align} {\cal I}&=2\exp\pars{-2a}\int_{0}^{\infty}\exp\pars{-a\bracks{x - {1 \over x}}^{2}} \,\dd x\tag{1} \end{align}
With $\ds{x \equiv {1 \over u}}$ we'll get:
\begin{align} {\cal I}&=2\exp\pars{-2a}\int_{\infty}^{0}\exp\pars{-a\bracks{{1 \over u} - u}^{2}} \,\pars{-\,{\dd u \over u^{2}}} \end{align}
\begin{align} {\cal I}&=2\exp\pars{-2a}\int_{0}^{\infty}\exp\pars{-a\bracks{{1 \over u} - u}^{2}} {1 \over u^{2}}\,\dd u\tag{2} \end{align}
$\pars{1}$ and $\pars{2}$ lead to: \begin{align} 2{\cal I}&=2\exp\pars{-2a}\int_{0}^{\infty} \exp\pars{-a\bracks{x - {1 \over x}}^{2}}\pars{1 + {1 \over x^{2}}}\,\dd x \\[5mm]&=2\exp\pars{-2a}\int_{x\ =\ 0}^{x\ \to\ \infty} \exp\pars{-a\bracks{x - {1 \over x}}^{2}}\,\dd\bracks{x - {1 \over x}} \\[5mm]&=2\exp\pars{-2a}\,{1 \over \root{a}} \underbrace{\int_{-\infty}^{\infty}\exp\pars{-x^{2}}\,\dd x} _{\ds{\color{#c00000}{\root{\pi}}}}\ =\ {2\exp\pars{-2a} \over \root{a}}\,\root{\pi} \end{align}
$$ {\cal I}= \color{#66f}{\large\int_0^{\infty}t^{-1/2}\exp\pars{-a\bracks{t + t^{-1}}}\,\dd t ={\exp\pars{-2a} \over \root{a}}\,\root{\pi}} $$
$\ds{{\tt\mbox{The Gamma function}}\ \Gamma\pars{z}\ \mbox{appears in the evaluation of the Gaussian integral !!!}}$
Let $t=u^2$, and the integral becomes
$$2 \int_0^{\infty} du \, e^{-a \left (u^2 + \frac1{u^2} \right)} = 2 e^{2 a} \int_0^{\infty} du \, e^{-a \left ( u+\frac1{u} \right )^2}$$
Let $v=u+1/u$, then
$$u = \frac12 \left (v \pm \sqrt{v^2-4} \right ) $$
$$du = \frac12 \left (1 \pm \frac{v}{\sqrt{v^2-4}} \right ) dv $$
Now, it should be understood that as $u$ traverses from $0$ to $\infty$, $v$ traverses from $\infty$ down to a min of $2$ (corresponding to $u \in [0,1]$), then from $2$ back to $\infty$ (corresponding to $u \in [1,\infty)$). Therefore the integral is
$$e^{2 a} \int_{\infty}^{2} dv \left (1 - \frac{v}{\sqrt{v^2-4}}\right ) e^{-a v^2} + e^{2 a} \int_{2}^{\infty} du \left (1 + \frac{v}{\sqrt{v^2-4}}\right ) e^{-a v^2} $$
which is
$$\begin{align}2 e^{2 a}\int_2^{\infty} dv \frac{v}{\sqrt{v^2-4}} e^{-a v^2} &= e^{2 a} \int_4^{\infty} \frac{dy}{\sqrt{y-4}} e^{-a y}\\ &= e^{-2 a} \int_0^{\infty} dq \, q^{-1/2} \, e^{-a q} \end{align}$$
I guess the gamma function comes from this integral, but I find it easier to refer to gaussian integrals.
Alternatively, denote the integral as $I$ and use substitution $t=x^2$, the integral becomes $$I=2\int_0^\infty e^{-\Large a\left(x^2+\frac{1}{x^2}\right)}\,dx=2e^{-2\large a}\int_0^\infty e^{-\Large a\left(x-\frac{1}{x}\right)^2}\,dx$$ Since the integrand is an even function, then it can be rewritten as $$\begin{align}I&=e^{-2\large a}\int_{-\infty}^\infty e^{-\Large a\left(x-\frac{1}{x}\right)^2}\,dx\\&=e^{-2\large a}\int_{-\infty}^0 e^{-\Large a\left(x-\frac{1}{x}\right)^2}\,dx+e^{-2\large a}\int_{0}^\infty e^{-\Large a\left(x-\frac{1}{x}\right)^2}\,dx\\&=I_1+I_2\end{align}$$ Making substitution $x=-e^{-\theta}$ for $I_1$ and $x=e^{\theta}$ for $I_2$ then using $\sinh\theta=\dfrac{e^{\theta}-e^{-\theta}}{2}$ and $\cosh\theta=\dfrac{e^{\theta}+e^{-\theta}}{2}$, we have $$\begin{align}I&=e^{-2\large a}\int_{-\infty}^\infty e^{-2a\sinh^2\theta}e^{-\theta}\,d\theta+e^{-2\large a}\int_{-\infty}^\infty e^{-2a\sinh^2\theta}e^{\theta}\,d\theta\\&=2e^{-2a\large a}\int_{-\infty}^\infty e^{-2a\sinh^2\theta}\cosh\theta\,d\theta\\&=2e^{-2\large a}\int_{-\infty}^\infty e^{-2a\sinh^2\theta}\,d(\sinh\theta)\\&=e^{-2\large a}\int_{0}^\infty e^{-2az^2}\,dz\end{align}$$ We can get a gamma function by using substitution $y=2az^2$.