I think I found a way to apply the method of stationary phase (please give feedback if you think that the answer is not sound).
As Jon noted via partial integration the integral can be brought into the form
$$I(\alpha) = \alpha \int_0^\infty\!dx\,x^2 K_1(x) \sin[2 \alpha K_0(x)].$$
Next I perform the substitution $y=2K_0(x)$ (note that $K_0$ is a monotonously falling function such that the inverse $K_0^{-1}$ is uniquely defined). I get
$$ I(\alpha) = \frac{\alpha}2 \int_0^\infty\!dy\,\sin(\alpha y)[K_0^{-1}(y/2)]^2; $$
note that $K_1$ cancelled with the derivative of $K_0^{-1}$ via inverse function theorem and the fact that $K_0'(x)=-K_1(x)$.
This integral is almost ready for a stationary phase analysis, we replace $\sin(\alpha y) = \operatorname{Im} e^{i\alpha y}$ and note that the path of stationary phase starting from $y=0$ is along the imaginary axis. Thus, we substitute $y=i \zeta$ and get
$$I(\alpha)=\frac{\alpha}2 \int_0^\infty \!d\zeta\, e^{-\alpha \zeta} \operatorname{Re}[K_0^{-1}(i\zeta/2)]^2 ;$$
here, we have analytically continued $K_0^{-1}$ into the complex plane.
The integral is dominated at $\zeta\approx 0$. Thus we can Taylor-expand $K_0^{-1}(i\zeta)$. We have $K_0(x) \sim e^{-x}$ thus $K_0^{-1}(y) \sim -\log y$. Using this asymptotic relation, we obtain
$$I(\alpha)\sim \frac{\alpha}2 \int_0^\infty \!d\zeta\, e^{-\alpha \zeta} \overbrace{\operatorname{Re}[\log(i\zeta/2)]^2}^{\sim \log^2 (\zeta/2)}
\sim \frac{\alpha}2 \int_0^\infty \!d\zeta\, e^{-\alpha \zeta}\log^2 (\zeta/2)
\sim \tfrac12 \log^2\alpha. $$
I guess for the next term one has to work a bit harder and use $K_0^{-1}(y) \sim - \log y - \tfrac12 \log(-\log y)$. However, in principle it should be possible to extend this approach to next to leading terms.
One can try to do the following derivations.
$$\mathrm{Int}=\frac{1}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{(-\frac{x}{\alpha}-\frac{1}{2w^2}(x-\hat{x})^2)} \, I_0\left(\frac{x}{\beta}\right)\,dx$$
You can simplify the power of the exponent:
$$-\left(\frac{x}{\alpha}+\frac{1}{2w^2}(x-\hat{x})^2\right)=-\left(\frac{(x-\mu)^2}{2w^2}+\gamma \right),$$
where $\gamma=\left(\frac{\hat{x}}{\alpha}-\frac{w^2}{2\alpha^2}\right)$, $\mu=\hat{x}-\frac{w^2}{\alpha^2}$. So:
$$\mathrm{Int}=\frac{1}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{(-\frac{x}{\alpha}-\frac{1}{2w^2}(x-\hat{x})^2)} \, I_0\left(\frac{x}{\beta}\right)\,dx=\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{-\frac{(x-\mu)^2}{2w^2}} \, I_0\left(\frac{x}{\beta}\right)\,dx$$
Then you can change the variable of integration $y=x-\mu$:
$$\mathrm{Int}=\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, (y+\mu) \, e^{-\frac{y^2}{2w^2}} \, I_0\left(\frac{y+\mu}{\beta}\right)\,dx$$
and make use of the Neumann’s addition theorem:
$$\mathop{I_{{\nu}}}\nolimits\!\left(u\pm v\right)=\sum _{{k=-\infty}}^{\infty}(\pm 1)^{k}\mathop{I_{{\nu+k}}}\nolimits\!\left(u\right)\mathop{I_{{k}}}\nolimits\!\left(v\right)$$
setting $\nu=0$, assuming that $y,\mu,\beta\in\mathbb{R}$ and using connection formulas one will get:
$$\mathop{I_{0}}\nolimits\!\left(\frac{y+\mu}{\beta}\right)=\sum _{{k=-\infty}}^{\infty}\mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{k}}}\nolimits\!\left(\frac{\mu}{\beta}\right)=\mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{0}}}\nolimits\!\left(\frac{\mu}{\beta}\right)+2\sum _{{k=1}}^{\infty}\mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{k}}}\nolimits\!\left(\frac{\mu}{\beta}\right)$$
$$\mathrm{Int}=\!
\frac{\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)e^{-\gamma}}{\sqrt{2\pi w^2}}\!\int_{-\infty}^{+\infty} \, (y+\mu)\! e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx\!
+\!\frac{2 e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty}\! (y+\mu) \, e^{-\frac{y^2}{2w^2}} \! \sum _{k=1}^{\infty}\!\mathop{I_k}\!\!\left(\frac{y}{\beta}\right)\!\!\mathop{I_k}\nolimits\!\left(\frac{\mu}{\beta}\right) \,dx$$
Or reorganising terms:
$\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)$
$$\mathrm{Int}=\!
\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\!\left(\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)\int_{-\infty}^{+\infty} \, (y\!+\!\mu)\! e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx\!
+\!2\!\sum _{k=1}^{\infty}\!\mathop{I_k}\!\left(\frac{\mu}{\beta}\right)\!\int_{-\infty}^{+\infty}\! (y\!+\!\mu) \, e^{-\frac{y^2}{2w^2}} \! \!\mathop{I_k}\!\left(\frac{y}{\beta}\right)\!\!\,dx \right)$$
Let's set
$$\mathrm{Int}_k(\mu,\alpha,\beta)\!=\!\int_{-\infty}^{+\infty} \, (y\!+\!\mu) e^{-\frac{y^2}{2w^2}}\! \mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx$$
$$\mathrm{Int}_0(\mu,\alpha,\beta)\!=\!\int_{-\infty}^{+\infty} \, (y\!+\!\mu) e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx$$
Here one can play with the evenness/unevenness of the functions under integrals (as I did at first), or just use CAS and get:
$$\mathrm{Int}_k(\mu,\alpha,\beta)\!=\!\sqrt{\frac{\pi \alpha ^2}{8 \beta ^2}}\! e^{\frac{\alpha ^2}{4 \beta ^2}} \left(\!\!2 \beta \left(\!\!(-1)^k+1\!\right)\! \mu I_{\frac{k}{2}}\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)+\alpha ^2 \left(\!1-(-1)^k\!\right)\!\! \left(\!\!I_{\frac{k-1}{2}}\!\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)+I_{\frac{k+1}{2}}\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)\!\!\right)\!\!\right)$$
$$\mathrm{Int}_0(\mu,\alpha,\beta)\!=\! \sqrt{2 \pi } \alpha \mu e^{\frac{\alpha ^2}{4 \beta ^2}} I_0\left(\frac{\alpha ^2}{4 \beta ^2}\right)$$
So the original integral can be represented in the following way:
$$\mathrm{Int}(\mu,\alpha,\beta)\!=\! \frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\left(\mathrm{Int}_0(\mu,\alpha,\beta)+
2\!\sum _{k=1}^{\infty}\mathop{I_k}\left(\frac{\mu}{\beta}\right)
\mathrm{Int}_k(\mu,\alpha,\beta)
\right)$$
And I am not shure that there is any chance to simplify it futher (meaning finding the sum of the series).
Best Answer
Use integration by parts and the fact that $\int x K_0(x)dx = -x \frac{d}{dx}K_0(x)=-xK'_0(x)$
$$ \int x\ln(x)K_0(x)\,dx = -x\ln(x)K'_0(x) - \int (-x K'_0(x))(\frac{1}{x}) \, dx =\dots. $$