How would you work out this integral as a series? The integral is:
$$f(a) = \int\limits_{-\infty}^\infty \frac{x e^{-a^2 x^2}}{\tanh(x)}dx$$
which seems to converge for all $a>0$. But I don't know how to get a sense of the function $f(a)$ such as writing it as a convergent series. The usual Taylor series has infinity for each term. Any ideas?
Edit:
I believe (using answers below) that when $a$ is close to zero we have (doing a substitution):
$$f(a) = \frac{1}{a^2}\int\limits_{-\infty}^\infty \frac{x e^{-x^2}}{\tanh (x/a)}dx$$
But using $\tanh(x/a)\rightarrow\operatorname{sign}(x)$ as $a\rightarrow 0^+$. So the above should become:
$$f(a) \approx \frac{2}{a^2}\int\limits_{0}^\infty x e^{-x^2}dx = \frac{1}{a^2}$$
So that gives the behaviour of $f(a)$ when $a$ is small. But I don't know how to give extra terms. (Also not sure if this is mathematically correct). From numerical calculations I find that near $0$,
$$
f(a)\approx \frac{1}{a^2} + \frac{\pi^6}{6} - \frac{\pi^4 }{60}a^2+\frac{\pi^6 }{252}a^4+... = \frac{1}{a^2}\sum\limits_{n=0}^\infty \frac{B_{2n} a^{2n}\pi^{2n}}{n!}
$$
although apparently this doesn't converge?
Comment
The answers below are two asymptotic series depending on whether $a$ is small or large. These give good approximations if we truncate the summation before begins to diverge. In the mid-range, when $a^2=1/\pi$, these two sums become term-by-term equal and the closest the truncated sum get to the true answer of $f(1/\sqrt{\pi})$ is to about 1% error.  Using both these sums, we can know any value to within about 1%-2% error, and if $a$ is small or large then much more accurately.
 A: A Taylor series exists in powers of $1/a$:
$$f(a) = \int\limits_{-\infty}^\infty \frac{x e^{-a^2 x^2}}{\tanh x}\,dx=a^{-2}\int\limits_{-\infty}^\infty x e^{-x^2}\,\text{cotanh}\, (x/a)\,dx$$
$$=\sum_{n=0}^{\infty}\frac{2^{2 n} B_{2 n}   \,\Gamma \left(n+\frac{1}{2}\right)}{(2 n)!a^{2 n+1}},$$
with $B_{2n}$ the Bernoulli number.

To assess whether this asymptotic series is useful for $a\gtrsim 1$, below I plot the sum $\sum_{n=0}^{10}$, so the first eleven terms, for $a=1,2,3\ldots 10$ (blue data points), and compare with a numerical evaluation of the integral (blue curve).

A: To obtain the series in $a$, separate off the leading term proportional to $1/a^2 $ and expand the Gaussian instead of the hyperbolic cotangent:
\begin{eqnarray}
f(a)&=& 2\int_{0}^{\infty } dx\, x \ e^{-a^2 x^2 } \left[ \coth x -1 +1 \right] \\
&=& 2\int_{0}^{\infty } dx\, x \ e^{-a^2 x^2 }+
2\sum_{n=0}^{\infty } \frac{(-a^2)^{n}}{n!} \int_{0}^{\infty} dx\, x^{2n+1} [\coth x -1] \\
&=& \frac{1}{a^2 } + 2\sum_{n=0}^{\infty } \frac{(-a^2)^{n}}{n!} \frac{ (-1)^n \pi^{2n+2} B_{2n+2} }{2n+2} \\
&=& \frac{1}{a^2} + \frac{1}{a^2} \sum_{n=1}^{\infty } \frac{(\pi a)^{2n} B_{2n} }{n!} \\
&=& \frac{1}{a^2} \sum_{n=0}^{\infty } \frac{(\pi a)^{2n} B_{2n} }{n!}
\end{eqnarray}
Also this is only an asymptotic expansion - the Bernoulli numbers diverge more rapidly than $n!$ at large $n$.
