I've recently come across the following integrals: \begin{equation} I_\pm(a)=\frac{1}{2}\int_0^1dx\int_0^1dy\frac{(y-2)x\pm2}{y}\exp\left[a\frac{((y-2)x+2)((y-2)x-2)}{xy}\right]. \end{equation} I'm really only interested in the asymptotic behavior as $a\rightarrow0^+$, but I'm not particularly familiar with asymptotic methods outside of a book by Bleistein and Handelsman I found at my library. Also, I've noticed that a simple series expansion doesn't work here, since we need the exponential to surpress the power divergence at 0.
Does anyone have any ideas for how to approach problems like these?
Out of curiosity, I numerically integrated the functions for various small positive values of $a$ and graphed the results:
There is a clear decay with respect to $a$ here, and this I what I expect physicsally ($a$ is related to the inverse energy scale in some scattering problem). I'm just hoping now to find some functional form for the behavior.
Thanks in advance for any suggestions!

That is a though one… I will concentrate on your interal $I_+(a):= I(a)$. We write:
$$ I(a)=I_1(a)+I_2(a)=\\ \frac{1}2\int_0^1dx\int_0^1dy\left(\frac{2-2x}{y}\right)e^{a g(x,y)}+\frac{1}2\int_0^1dx\int_0^1dyxe^{a g(x,y)}$$
where $g(x,y):=\frac{(y-2)^2x^2-4}{xy}$. It is now easy to show that $-\infty< g(x,y)< c_0$ inside the domain of integration (let us denote by $c_i, i\in \{0,1,2,..\}$ a set of constants) which implies that $I_2(a)<c_1$ independent of $a$. More informally speaking we can perform a simple Taylor expansion as $a\rightarrow0$. As OP already recognized, due to the $y^{-1}$-terms, things are not that easy for $I_1(a)$.
Let us further rewrite $g(x,y):=g_d(x,y)+g_f(x,y)$ where $g_d(x,y)=4 x/y-4/xy$ denotes the divergent part (in $y$) and $g_f(x,y)=yx-x^2y$ the corresponding finite part. It is clear, that to leading order we can put $e^{g_f(x,y)}=1$ since higher order Terms will always be multiplied by a small $a$. We therefore have (important: $1/x-x\geq 0$ for $x \in (0,1] $)
$$ 2 I_1(a)\sim\int_0^1(2-2x)dx\int_0^1dy\frac{1}{y}e^{a\tfrac4y(x-1/x)} $$
the inner integral might now be expressed as an incomplete Gammafunction
$$ 2 I_1(a)\sim \int_0^1 dx (2-2x)\Gamma\left(0,4a\left(\frac1x-x\right) \right) $$
Now we substitute $1/x-x=z$ which yields
$$ 2 I_1(a)\sim\int_0^{\infty}f(z)\Gamma(0,az) $$ with $f(z):=\tfrac{(\sqrt{z^2+4}-z)(\sqrt{z^2+4}-z-2)}{2 \sqrt{z^2+4}}$. This integral is now a perfect candidate for the method of asymptotic matching. Let us define $1<<\Delta<<1/a$. We have:
$$ 2 I_1(a)\sim \mathcal{I}_1(a)+\mathcal{I}_2(a)=\int_0^{\Delta}dzf(z)\Gamma(0,4az)+\int_{\Delta}^{\infty}dzf(z)\Gamma(0,4az) $$
For $\Delta\geq z\geq0$ we have $az\ll1$ so we can expand $\Gamma$ to leading order
$$ \mathcal{I}_1(a)\sim -\log(a)\int_0^{\Delta}f(z)dz=-\log(a)+O(\log(a)\Delta^{-1}) $$
For $\Delta\leq z\leq \infty$ we have $f(z)\sim -2/z^2$ so $$ \mathcal{I}_2(a)\sim-2\int_{\Delta}^{\infty}dz\frac{\Gamma(0,4az)}{z^2} $$
Using i.b.p with $\Gamma(0,q)'=e^{-q}/q,\quad \int dq q^{-2}=q^{-1}$ this integral can be done exactly: $$ \mathcal{I}_2(a)\sim -8a \Gamma(-1,4a \Delta)+2\frac{\Gamma(0,4a \Delta)}{\Delta} $$
or (recall $\Delta a\rightarrow 0$ as $a\rightarrow 0$) $$ \mathcal{I}_2(a)\sim \mathcal{O}(\log(a)\Delta^{-1}) $$
Choosing $\Delta = a^{-\epsilon}$ with $1 >\epsilon > 0$ we see that $\mathcal{I}_2(a)$ as well as the second term of $\mathcal{I}_1(a)$ are negligible compared with $\log(a)$ so we have $$ 2 I_1(a)\sim -\log(a) $$ and since furthermore $c_1$ is also small compared to $\log(a)$ we can also discard $I_2(a)$ which means:
PS: No Mathematica at Hand, so i can't say for sure i haven't made any mistakes PPS: The subleading term is most likely a constant