I want to find an asymptotic approximation when $a\to 0^+$ for the integral $$I(a):=\int_0^\infty \int_0^{a/x}e^{-x-y}\ dy\ dx.$$ I found the following approximation: $$C_1\, a\, \mathrm{ln}(1/a) < I(a)< C_2 \sqrt{a},$$ where $C_1$ and $C_2$ are constant terms.
I would like to have a better approximation, something of the type $I(a)=\Theta(\sqrt{a})$ or $I(a)=\Theta(a\, \mathrm{ln}(1/a))$ when $a\to 0^+$.
In the folowing I explain how I found those bounds. First we split our integral: $$I(a):=\int_0^\sqrt{a} \int_0^\sqrt{a} e^{-x-y}\ dy\ dx + \int_0^\sqrt{a} \int_\sqrt{a}^{a/x}e^{-x-y}\ dy\ dx +\int_\sqrt{a}^\infty \int_0^{a/x}e^{-x-y}\ dy\ dx.$$ But the two last integral are equal (clear if you draw the picture of the domain of integration and observe the symmetry), so we have: $I(a)=I_1(a)+2\, I_2(a)$ with $$I_1(a):=\int_0^\sqrt{a} \int_0^\sqrt{a} e^{-x-y}\ dy\ dx\quad \text{ and }\quad I_2(a):=\int_\sqrt{a}^\infty \int_0^{a/x}e^{-x-y}\ dy\ dx.$$ We consider that $a$ is small, so we have $0<a<1$ and in the domain of integration of $I_1(a)$, we have $e^{-2}<e^{-x-y}<1$. Hence $e^{-2} a <I_1(a)< a$. That is $$I_1(a)=\Theta(a).$$ This first integral will not play any role for our bounds, because the second integral will be much bigger (for small values of $a$).
Now we observe that if $x>\sqrt(a)$, then $a/x<a/\sqrt{a}=\sqrt{a}$. Hence $$I_2(a)<\int_\sqrt{a}^\infty \int_0^\sqrt{a}e^{-x-y}\ dy\ dx <\int_\sqrt{a}^\infty e^{-x} \int_0^\sqrt{a}1\ dy\ dx =\sqrt{a}\,e^{-\sqrt{a}}<\sqrt{a}.$$ So we have proved the upper bound, and will show the above bound now. We have for $x>\sqrt{a}>a$ and $a<1$, $$\int_0^{a/x}e^{-y}\ dy> \frac{a}{x} e^{-a/x} >\frac{a}{x} e^{-1}.$$ So $$I_2(a)> a\,e^{-1} \int_\sqrt{a}^\infty \frac{e^{-x}}{x}\ dx >a\,e^{-1} \int_\sqrt{a}^1 \frac{e^{-x}}{x}\ dx >a\,e^{-2} \int_\sqrt{a}^1 \frac{1}{x}\ dx =\frac{e^{-2}}{2}\,a\,\textrm{ln}(1/a).$$ This gives us the lower bound.
Your method will definitely work, we just need to keep track of the errors along the way.
You showed that
$$ \begin{align} I(a) &= 2I_2(a) + \Theta(a) \\ &= 2\int_{\sqrt{a}}^\infty e^{-x} - e^{-a/x-x}\,dx + \Theta(a). \tag{1} \end{align} $$
As you observed, the tail of the integral is insignificant. In fact we have
$$ \begin{align} 0 &\leq \int_{1}^\infty e^{-x} - e^{-a/x-x}\,dx \\ &\leq \int_{1}^\infty e^{-x} - e^{-a/1-x}\,dx \\ &= (1-e^{-a}) \int_1^\infty e^{-x}\,dx, \end{align} $$
so that
$$ \int_{1}^\infty e^{-x} - e^{-a/x-x}\,dx = O(a). \tag{2} $$
In the remaining interval $x \in [\sqrt{a},1]$ we have $e^{-a/x} = 1 - a/x + O(a/x)^2$. Substituting this into the integral yields
$$ \begin{align} \int_\sqrt{a}^1 e^{-x} - e^{-a/x-x}\,dx &= \int_{\sqrt{a}}^1 e^{-x} - e^{-x}\Bigl[1 - a/x + O(a/x)^2\Bigr]\,dx \\ &= a \int_\sqrt{a}^1 \frac{e^{-x}}{x}\,dx + O\left(a^2\int_\sqrt{a}^1 \frac{e^{-x}}{x^2}\,dx\right). \tag{3} \end{align} $$
In the integrals we substitute $e^{-x} = 1 + O(x)$ to get
$$ \int_\sqrt{a}^1 \frac{e^{-x}}{x}\,dx = \int_\sqrt{a}^1 \frac{dx}{x} + O(1) = \frac{\log(1/a)}{2} + O(1) $$
and
$$ \int_\sqrt{a}^1 \frac{e^{-x}}{x^2}\,dx = \frac{1}{\sqrt{a}} + O(\log a). $$
Thus $(3)$ becomes
$$ \int_\sqrt{a}^1 e^{-x} - e^{-a/x-x}\,dx = \frac{a\log(1/a)}{2} + O(a). $$
Combining this with $(1)$ and $(2)$, we conclude that
$$ I(a) = a\log\left(\frac{1}{a}\right) + O(a) $$
as $a \to 0^+$.