Solve $\int _{x=0}^{\infty }\int _{t=-\infty }^{\infty }\exp \left(\frac{-a t^2+i b t}{3 t^2+1}+i t x\right)\frac{x}{3 t^2+1}\mathrm{d}t\mathrm{d}x$ How to solve this double integral?
$$f(a,b)=\int _{x=0}^{\infty }\int _{t=-\infty }^{\infty }\exp \left(\frac{-a t^2+i b t}{3 t^2+1}+i t x\right)\frac{x}{3 t^2+1}\mathrm{d}t\mathrm{d}x$$
$$\text{with }a>0,b\in \mathbb{R},i^2=-1$$
Known special solution for ${\bf b=0}$
$$f(a,0)=\frac{\pi}{\sqrt{3}\, {\rm exp}\left(\frac{a}{6}\right)}\left[(a+3) I_0\left(\frac{a}{6}\right)+a I_1\left(\frac{a}{6}\right)\right]$$
where $I_0,I_1$ are Bessel functions of order $0$ and $1$ (proof).
The difference of the special solution to the general solution $f(1,0)-f(1,b)$ is calculated numerically for $a=1$ and $-50<b<50$.

What I tried
I followed the first steps given here. Substitution of $t\rightarrow t \sqrt{3},x\rightarrow x/\sqrt{3}$ removes the factors in the denominator
$$f(a,b)=\sqrt{3}\int _{x=0}^{\infty }\int _{t=-\infty }^{\infty } \exp \left(\frac{1}{t^2+1}\left(\frac{i b t}{\sqrt{3}}-\frac{a t^2}{3}\right)+i t x\right)\frac{x}{t^2+1}\mathrm{d}t\mathrm{d}x$$
$$=\sqrt{3}\int _{x=0}^{\infty }{\rm d}x \frac{x}{\text{exp}(x)}\int _{t=-\infty }^{\infty} \exp \left(\frac{1}{t^2+1}\left(\frac{i b t}{\sqrt{3}}-\frac{a t^2}{3}\right)+i x(t-i) \right)\frac{1}{t^2+1}\mathrm{d}t$$
$$=\frac{\sqrt{3}}{{\rm exp}(a/3)}\int _{x=0}^{\infty }{\rm d}x \frac{x}{\text{exp}(x)}\underbrace{\int _{t=-\infty }^{\infty} \exp \left(\frac{1}{t^2+1}\left(\frac{i b t}{\sqrt{3}}+\frac{a}{3}\right)+i x(t-i) \right)\frac{1}{t^2+1}\mathrm{d}t}_{I(x)}.$$
Now $I(x)$ can be closed in the upper half plane since the contribution along the arc vanishes. Then this $t$-integral encloses the single essential singularity in the upper half plane at $t=i$. Hence we have
$$I(x)=2\pi i \, {\rm Res} \left(\exp \left(\frac{1}{t^2+1}\left(\frac{i b t}{\sqrt{3}}+\frac{a}{3}\right)+i x(t-i) \right)\frac{1}{t^2+1}\right)\Bigg|_{t=i} \, $$
where
$$\exp \left(\frac{1}{t^2+1}\left(\frac{i b t}{\sqrt{3}}+\frac{a}{3}\right)+i x(t-i) \right)\frac{1}{t^2+1}$$
can be written as the series
$$\sum_{n,m=0}^{\infty}\frac{(ix)^n}{n!}(t-i)^n\frac{1}{m!} \left(\frac{i b t}{\sqrt{3}}+\frac{a}{3}\right)^m\frac{1}{[(t+i)(t-i)]^{m+1}}$$
 A: I tried various approaches to get a result computationally - without success.
Then I decided to iterate (using two nested loops) both variables $a,b$, the first variable $a$ from $0.1$ to $10.1$ by steps of $0.5$ and the second variable $b$ from $-5$ to $5$ by steps of $0.5$. Whithin these loops, I am calculating the integral numerically and plot the surface of this function depending from both variables $a$ and $b$. To get a rough idea of that integral, here is my first picture:

Computation took hours and I need to optimize the code. After this I will be able to expand the iteration range. Hope this helps a bit.
Here is the code, if anyone would like to run it:
from mpmath import mp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

arr_x = np.arange(0.1,10.1,0.5)
arr_y = np.arange(-5,5,0.5)
X,Y = np.meshgrid(arr_x,arr_y)
Z = np.zeros([arr_x.size, arr_y.size])
I = mp.sqrt(-1)

for idx_a, a in enumerate(arr_x):
   for idx_b, b in enumerate(arr_y):
       f = lambda x, t: mp.exp((-a*t**2+I*b*t)/(3*t**2+1)+I*t*x)*(x/(3*t**2+1))
       m = mp.quad(f, [0, mp.inf], [-mp.inf, mp.inf])
       mreal = m.real
       mimag = m.imag
       Z[idx_a][idx_b] = mreal

fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(111, projection='3d')
mycmap = plt.get_cmap('gist_earth')
ax1.set_title('gist_earth color map')
surf1 = ax1.plot_surface(X, Y, Z, cmap=mycmap)
fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=5)
plt.show()

Thanks also to the Stack Overflow community who helped me.
Update (2021-12-31):
I ammended the code, performed a new run and updated the figure above that contains the plot.
A: Not a closed form, but maybe some progress. We write
$$f(3a,\sqrt{3}b)=\sqrt{3}\int _{x=0}^{\infty }\int _{t=-\infty }^{\infty } \exp \left(\frac{i b t-a t^2}{t^2+1}+i t x\right)\frac{x}{t^2+1} \, \mathrm{d}t \, \mathrm{d}x \equiv \sqrt{3} \, I(a)$$
and close the $t$-integral in the upper half-plane, since the integral along the semi-circle vanishes. Introducing a dampening factor $e^{-\epsilon x^2}$ for which we are ultimately interested in the case $\epsilon \rightarrow 0$, we can now interchange the integrals, since the $x$-integral converges properly. Using $$\int_0^\infty x\, e^{-\epsilon x^2 + itx} \, {\rm d}x = -\frac{1}{t^2} + O(\epsilon) \, ,$$
the remaining $t$-integral becomes
$$I(a)=-\oint _{-\infty }^{\infty } \frac{\exp \left(\frac{i b t-a t^2}{t^2+1}\right)}{t^2(t^2+1)} \, \mathrm{d}t \, .$$
The singularity at $t=0$ is not an issue, since the original integral (before interchanging the integrals) has the only singularity in the upper half plane at $t=i$ and the closed contour can be deformed at will without crossing $t=i$.
Therefore
$$I'(a)=\oint _{-\infty }^{\infty } \frac{\exp \left(\frac{i b t-a t^2}{t^2+1}\right)}{(t^2+1)^2} \, \mathrm{d}t = e^{-a} \frac{{\rm d}}{{\rm d}a} \oint _{-\infty }^{\infty } \frac{\exp \left(\frac{i b t+a}{t^2+1}\right)}{t^2+1} \, \mathrm{d}t = e^{-a} \frac{{\rm d}}{{\rm d}a} \int _{-\infty }^{\infty } \frac{\exp \left(\frac{i b t+a}{t^2+1}\right)}{t^2+1} \, \mathrm{d}t$$
where in the last step the vanishing integral over the semi-circle was removed again, giving an ordinary indefinite integral. Substituting $t=\tan(\phi/2)$ converts this into
$$I'(a) = \frac{e^{-a}}{2} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} \int _{-\pi }^{\pi} e^{i b/2 \, \sin(\phi) + a/2 \, \cos(\phi)} \, \mathrm{d}\phi \\
\stackrel{z=e^{i\phi}}{=} \frac{e^{-a}}{2} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} \oint_{|z|=1} \frac{{\rm d}z}{iz} \, e^{\frac{a+b}{4}z + \frac{a-b}{4z}} = \pi {e^{-a}} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} \, {\rm Res}\left( \frac{e^{\frac{a+b}{4} z}e^{\frac{a-b}{4z}}}{z} \right)\Bigg|_{z=0} \\
=\pi {e^{-a}} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} \, {\rm Res}\left( \sum_{n,m=0}^\infty \frac{\left(\frac{a+b}{4}\right)^n\left(\frac{a-b}{4}\right)^m}{n!m!} \, z^{n-m-1} \right)\Bigg|_{z=0} \\
\stackrel{n=m}{=} \pi {e^{-a}} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} \sum_{m=0}^\infty \frac{\left(\frac{\sqrt{a^2-b^2}}{4}\right)^{2m}}{m!^2} = \pi {e^{-a}} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} I_0\left( \frac{\sqrt{a^2-b^2}}{2} \right) \, .$$
Hence $$I(a)=\pi \int {\rm d}a \, {e^{-a}} \frac{{\rm d}}{{\rm d}a} \, e^{a/2} I_0\left( \frac{\sqrt{a^2-b^2}}{2} \right) \\
\stackrel{\text{P.I.}}{=} \pi {e^{-a/2}} I_0\left( \frac{\sqrt{a^2-b^2}}{2} \right) + \pi \int {\rm d}a \, {e^{-a/2}} I_0\left( \frac{\sqrt{a^2-b^2}}{2} \right) + {\rm const.} \, .$$
In the special case $b=0$, the integral has closed form giving the already known result. The constant is obtained from some limiting case e.g. $a=b$.

Notice that you can easily calculate $$\frac{f(3a,\sqrt{3}b)-f(3a,-\sqrt{3}b)}{\sqrt{3}} = -2\pi b$$
by exploiting the identity $$\int_{-\infty}^\infty e^{itx} \, {\rm d}x = 2\pi \delta(t) \, .$$
