For $a,b>0$, Mathematica says that $$P \int_{-\infty}^\infty \frac{e^{-ax^2-bx}}{x} dx = -i\pi,$$ where $P$ denotes the principal-value integral. How can I derive the above?
Due to the factor $e^{-ax^2}$, the integrand grows exponentialy as $x \to \pm i \infty$. Hence, contour method with respect to the semicircle contour is not applicable.
My Mathematica 13 claims the answer is $i\pi$:
This is clearly wrong, though, since the principal value of $\int_{-\infty}^{\infty} f(x) \, \mathrm{d}x$ of a real-valued function $f$ on $\mathbb{R}$ cannot be imaginary.
1. What is the correct answer? We have
\begin{align*} \mathrm{PV}\!\int_{-\infty}^{\infty} \frac{e^{-ax^2-bx}}{x} \, \mathrm{d}x &= \frac{1}{2} \int_{-\infty}^{\infty} \frac{e^{-ax^2-bx} - e^{-ax^2+bx}}{x} \, \mathrm{d}x \\ &= -\frac{1}{2} \int_{-\infty}^{\infty} e^{-ax^2} \left( \int_{-b}^{b} e^{xs} \, \mathrm{d}s \right) \, \mathrm{d}x \\ &= -\frac{1}{2} \int_{-b}^{b} \left( \int_{-\infty}^{\infty} e^{-ax^2+sx} \, \mathrm{d}x \right) \, \mathrm{d}s \\ &= -\frac{1}{2} \int_{-b}^{b} \left( \sqrt{\frac{\pi}{a}} e^{s^2/4a} \right) \, \mathrm{d}s\\ &= -\pi \operatorname{erfi}\left(\frac{b}{2 \sqrt{a}}\right). \end{align*}
Here, $\operatorname{erfi}(\cdot)$ stands for the imaginary error function.
Addendum. The key idea in this solution is in fact an integral version of Feynman's trick. This means that we can adapt this logic to give a solution using Feynman's trick. Indeed, set
$$ f(b) = \mathrm{PV}\!\int_{-\infty}^{\infty} \frac{e^{-ax^2-bx}}{x} \, \mathrm{d}x. $$
Then $f(0) = 0$ since the integrand is an odd function. Moreover,
$$ f'(b) = - \mathrm{PV}\!\int_{-\infty}^{\infty} e^{-ax^2-bx} \, \mathrm{d}x = - \sqrt{\frac{\pi}{a}} e^{b^2/4a}. $$
Therefore,
$$ f(b) = f(0) + \int_{0}^{b} f'(s) \, \mathrm{d}s = - \int_{0}^{b} \left( \sqrt{\frac{\pi}{a}} e^{s^2/4a} \right) \, \mathrm{d}s = -\pi \operatorname{erfi}\left(\frac{b}{2 \sqrt{a}}\right). $$
2. What might have gone wrong? I suspect that what Mathematica computed is
\begin{align*} \int_{-\infty}^{\infty} \frac{e^{-ax^2-bx}}{x} \, \mathrm{d}x &= \lim_{\substack{\varepsilon & \to 0^+ \\ R & \to \infty}} \left( \int_{-R}^{-\varepsilon} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z + \int_{\varepsilon}^{R} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z \right) \\ &= \lim_{\substack{\varepsilon & \to 0^+ \\ R & \to \infty}} \left( \int_{\gamma_{\varepsilon}} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z - \int_{\gamma_{R}} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z \right), \end{align*}
where $\gamma_r$ denotes the counter-clockwise oriented, upper semicircular contour centered at $0$. Then
$$ \lim_{\varepsilon \to 0^+} \int_{\gamma_{\varepsilon}} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z = i\pi \, \underset{z=0}{\mathrm{Res}} \, \frac{e^{-az^2-bz}}{z} = i\pi, $$
but I suspect that Mathematica incorrectly assumed
$$ \color{red}{ \int_{\gamma_{R}} \frac{e^{-az^2-bz}}{z} \, \mathrm{d}z \quad \xrightarrow[?]{R\to\infty} \quad 0 }, $$
which is a fairly bold claim (and is indeed false in light of the above computation.)