Probablemente la mejor manera es regularizar la integral agregando un $x^s$, calcularlo en términos de funciones gamma incompletas, y luego dejar que $s \to 0$. Tenemos $$ \int_0^a x^{s-1}(e^{-x}-1)^2 \, dx = \int_0^a x^{s-1}(e^{-2x}-2e^{-x}+1) \, dx = \frac{a^s}{s} - 2\gamma(s,a) + 2^{-s}\gamma(s,2a), $$ en términos de la función gamma incompleta inferior. En términos de la función gamma incompleta superior, esto es $$ \frac{a^s}{s} + (2^{-s}-2)\Gamma(s) + 2\Gamma(s,a) - 2^{-s} \Gamma(s,2a). $$ Los dos últimos términos son analíticos en $s$, por lo que podemos tomar felizmente $s=0$ en ellos. Para el primero, tenemos $$ \frac{a^s}{s} = \frac{1}{s} + \log{a} + O(s), $$ usando la expansión en serie habitual para una potencia, y $$ (2^{-s}-2)\Gamma(s) = \left(-1-s\log{2}+O(s^2)\right) \left( \frac{1}{s} -\gamma + O(s) \right) = -\frac{1}{s} + \gamma - \log{2} + O(s), $$ usando la serie de Laurent de $\Gamma$ en $0$. Sumando estos, los divergentes $s^{-1}$ se cancelan como deberían, y obtenemos $$ \int_0^a x^{-1}(e^{-x}-1)^2 \, dx = \log{a}-\log{2}+\gamma + 2\Gamma(0,a) - \Gamma(0,2a). $$ (Los dos últimos términos también se pueden escribir usando integrales exponenciales si así se desea.)
Editado para agregar:
Recordemos una posible definición de $a^s$ es $a^s = e^{s\log{a}}$. Dado que sabemos cómo expandir la función exponencial como una serie de potencias, obtenemos la expansión $$ a^s = e^{s\log{a}} = \sum_{k=0}^{\infty} \frac{(\log{a})^k}{k!} s^k. $$ (alternativamente, considerar la definición de $e^x$ como el límite de $(1+x/n)^n$: reorganizar esto nos permite escribir $$ \log{y} = \lim_{n \to \infty} n(y^{1/n}-1) = \lim_{s \to 0} \frac{y^s-1}{s}, $$ que reconocemos como el cociente derivado de $y^s$ en $s=0$.)
$\gamma$ es la constante de Euler-Mascheroni, que básicamente se define como $-\Gamma'(1)$: para nuestros propósitos, proviene de la integral $$ \int_0^{\infty} e^{-t}\log{t} \, dt = -\gamma, $$ en la cual se reconoce la derivada de $x^{s-1} e^{-x} $ con respecto a $s$ en $s=1$.