I have to solve the following integral
$$ \int_{-\infty}^{\infty} dx \frac{e^x-1}{x} e^{-\frac{(x+y)^2}{2 \sigma^2}}$$
This integral converges, becaue $\frac{e^x-1}{x}$ goes to 1 for $x \to 0$. When I feed the previous integral to Mathematica, however, I get logarithms with negative arguments. For the case $y=0$, for example, the result of the integral is $-i \pi$. How can I have a complex argument from such a real integral? Is mathematica doing something wrong?
Mathematica is very wrong. If you use NIntegrate, it gets itself right. I am not sure why this is happening yet, but I can direct you toward the right answer for the simple case $y=0$ and $2 \sigma^2=1$.
Note that
$$\frac{e^x-1}{x} = \int_0^1 dt \, e^{x t} $$
Then the integral we want is
$$\int_{-\infty}^{\infty} dx \, e^{-x^2} \, \int_0^1 dt \, e^{x t} $$
We can reverse the order of integration and complete the square the exponential to get
$$\int_0^1 dt \, e^{t^2/4} \int_{-\infty}^{\infty} dx \, e^{-(x-t/2)^2} = \pi \operatorname{erfi}{\left (\frac12 \right )}$$
which agrees with the numerical results in Mathematica.