Bounds on integral $x^{-a} \int_{1}^x y^{a-1} \exp(-y a) dy$ Consider the function
$$
I(a,x) = x^{-a} \int_{1}^x y^{a-1} \exp(-y a) dy
$$
where $x \geq 1$, and $a \geq 0$.
I am not really interested in the parameter $x$, so define
$$
I(a) = \sup_{x \geq 1} I(a,x)
$$
Now what is the asymptotic behavior of $I(a)$ as $a \rightarrow 0$? 
I can show the bound $I(a) = \log(1/a) + O(1)$. This is by noting that $I(a,x) \leq \int_{1}^x \exp(-y a)/y \ dy \leq \int_{1}^{\infty} \exp(-y a)/y \ dy = \Gamma(0, a) = \log(1/a) + O(1)$.
The asymptotic behavior of $I(a)$ appears to be smaller than this, though, something like $I(a) \sim 0.9 \log(1/a)$
 A: In a first step, we use the variable transform $y=(y-1)/(x-1)$ and obtain
$$I(a,x)= \int_0^1\!dt\, x^{-a} (x-1) [1+ t(x-1)]^{a-1} e^{-a [1+ t (x-1)]}.$$
As we are interested in $a\to0$, we expand the integrand to first order in $a$ and obtain
$$\begin{align}I(a,x) &= \int_0^1\!dt\, \left[ \frac{x-1}{1+t(x-1)} + a \frac{(x-1)(\log[1+t (x-1)] - \log x -1 - t(x-1)}{1+ t (x-1)} + O(a^2) \right] \\
  &=\log x +a \left(1-x - \tfrac12 \log^2 x\right) + O(a^2)
\end{align}.$$
The function $I(a,x)$ assumes its maximum at $x^* = a^{-1} + O(\log a) $. Thus, we have $$I(a) = I(a,x^*) = \log a^{-1} -1 + o(a). $$
Edit:
As  oenamen pointed out the answer is not self consistent. In expanding the exponent to get from the first expression to the second, I assumed that $ax\ll1$. Then I found that $x \simeq a^{-1}$ which of course outside the scope of the first assumption. However, it is not difficult to convince oneself that the $\log a^{-1}$ scaling (with unit prefactor and not like the OP assumed with a different pre factor) is indeed the correct asymptotic expression. Thus $$I(a) = \log a^{-1} + O(1).$$
A: We must be careful when expanding in small $a$, we do not know a priori that $a x$ is small.
It is in fact necessary to treat $a$ as small, and $a x$ to be at least of order one.
While the maximum value of the integral goes like $\ln a^{-1}$,
there is an error of order one in the accepted answer.
Change variables.
Let $z = a y$.
The integral can then be written as
$$I(a,x) = (a x)^{-a}\left[\Gamma(a,a) - \Gamma(a,a x)\right],$$
where
$\Gamma(s,x) = \int_x^\infty dt \ t^{s-1} e^{-t}$
is the incomplete gamma function.
Look for an extremum, $X$.
(The reader can verify the extremum found below is a maximum.)
We find ${\partial I}/{\partial x}|_{x=X} = 0$ implies
\begin{equation*}
a\Gamma(a,a X) + (a X)^a e^{-a X} = a \Gamma(a,a) \tag{1}
\end{equation*}
so that
\begin{equation*}
I(a,X) = \frac{1}{a} e^{-a X}.
\end{equation*}
Integrating $a\Gamma(a,a X)$ in (1) by parts we find the condition on $X$ is
$$\frac{1}{a}\Gamma(a+1,a X) = \Gamma(a,a).$$
So far no expansion has been made.
Expanding the left hand side in small $a$ (but not small $a x$), we find
$$\begin{eqnarray}
\frac{1}{a}\Gamma(a+1,a X) &=& \frac{1}{a}e^{-a X} + \mathrm{h.o.} \\
&=& I(a,X) + \mathrm{h.o.}.
\end{eqnarray}$$
We find the higher order terms go like $e^{-a X} \ln a X$. 
We assume $a X$ is large enough so these terms are suppressed. 
We will find this assumption to be self consistent. 
Therefore, $I(a,X) = \Gamma(a,a) +  \mathrm{h.o.}$.
Expanding $\Gamma(a,a)$ in small $a$ we find
$$\begin{eqnarray}
I(a,X) &=& -\mathrm{Ei}(-a) + \mathrm{h.o.}  \\
&=& -\gamma + \ln a^{-1} + \mathrm{h.o.},
\end{eqnarray}$$
where $\mathrm{Ei}(x) = \int_{-\infty}^x d t \ t^{-1} e^t$ is the exponential integral and $\gamma$ is the Euler-Mascheroni constant.
The location of the maximum is
$$X = -\frac{1}{a} \ln\left[a(\ln a^{-1}-\gamma)\right] + \mathrm{h.o.}$$
For $a = e^{-k}$ we find $a X \approx k$ so the higher order terms discussed above go like $e^{-k}\ln k$.
Thus, since $k$ is large the expansion is valid. 
Below is a plot of $I(a,x)$ for $a=10^{-4}$.
The predicted maximum is $I(a,X) \approx 8.63$ at $X \approx 7.05\times 10^4$.
Clearly $\ln a^{-1} -1 \approx 8.21$ underestimates the maximum value of the integral.

