Definite integral of Gaussian divided by hyperbolic cosine

Question

Question: is there a nice formula for $\int_{-\infty}^{+\infty} \frac{e^{-rx^2}}{e^x + e^{-x}} dx$ ? For a real parameter $r > 0$.

Maybe useful: consider the bilateral Laplace transform $$J(z) = \int_{-\infty}^{+\infty} \frac{e^{-rx^2}}{e^x + e^{-x}} e^{zx} \, \, dx $$ $J$ is an entire even function. We want to compute $J(0)$.

$$ \partial^n J(z) = \int_{-\infty}^{+\infty} \frac{e^{-rx^2}}{e^x + e^{-x}} x^n e^{zx} \, dx $$

$$ (e^{\partial} + e^{-\partial})J(z) = \int_{-\infty}^{+\infty} e^{-rx^2} e^{zx} \, dx $$

$$J(z-1) + J(z+1) = \sqrt{\frac{\pi}{r}} e^{\frac{z^2}{4r}} $$

Since $J(1)=J(-1)$, we can compute $J(1)$. The functional equation now provides a formula for $J(n)$ where $n$ is an odd integer.

Answer 1

The integral $$I(r)=\int_0^{\infty} \frac{e^{-rx^2}}{\cosh x}\,dx$$ has no closed-form expression, however there is an exact identity relating small and large values of $r$: $$I(r)=\frac{1}{\sqrt{\pi r}}I\left(\frac{1}{\pi^2 r}\right),$$ with $I(0)=\pi/2$, $\lim_{r\rightarrow\infty}\sqrt{r}I(r)=\tfrac{1}{2}\sqrt{\pi}$ as special cases.

This follows from the fact that the cosine transforms ${F}(k)=\sqrt{2/\pi}\int_0^\infty f(x)\cos kx\,dx$ of $f(x)=e^{-x^2/2}$ and $g(x)=1/\cosh(\sqrt{\pi/2}x)$ are the same functions: ${F}(k)=f(k)$ and ${G}(k)=g(k)$. Substitution of the transform and exchange of the order of integration then gives $$\int_0^\infty f(x)g(\alpha x)\,dx=\int_0^\infty f(x){G}(\alpha x)\,dx=\int_0^\infty f(x)\left(\sqrt{2/\pi}\int_0^\infty g(y)\cos (\alpha xy)\,dy\right)\,dx$$ $$=\int_0^\infty F(\alpha y)g(y)\,dy=\int_0^\infty f(\alpha x)g(x)\,dx,$$ from which the identity follows.

Ramanujan gave this result in a letter to Hardy, in the symmetric form $$\sqrt{\alpha} \int_{0}^{\infty} \frac{e^{-x^2}}{\cosh\alpha x} dx=\sqrt{\beta} \int_{0}^{\infty} \frac{e^{-x^2}}{\cosh\beta x} dx,\;\;\text{if}\;\;\alpha\beta=\pi.$$ Hardy wrote back that he knew of the result and had published it.

Answer 2

This is a method for approximating the integral that I consider interesting. It makes use of the integrals $J(n)$ for $n$ odd defined at the question.

Consider the functions

\begin{equation} a_k(x) = e^{-rx^2 + (2k - 1/2)x} \in L^2(\mathbb R), \mbox{ for } k \in \mathbb Z \\ b_k(x) = \frac{a_k}{\|a_k\|_2} \\ h(x) = e^{-rx^2 + x/2} \\ f(x) = \frac{e^{-x/2}}{e^x + e^{-x}} \end{equation}

With these definitions, the integral we want to compute takes the form $\langle f, h \rangle$ in $L^2(\mathbb R)$.

Call $S$ the closure of the subspace generated by $\{b_k \}_{k \in \mathbb Z}$, and $P_S$ the orthogonal projection to $S$. After computing the norms $\|a_k\|_2$ we reach $$A_{ij} := \langle b_i, b_j \rangle = e^{-(i-j)^2 / (2r)}$$ We can compute $P_S(h)$ since the coordinates $\langle h, b_k \rangle$ are definite integrals of Gaussian functions. We can also compute $P_S(f)$ because each $\langle f, b_k \rangle$ is a constant times $J(2k-1)$.

This means that we can accurately and efficiently compute $\langle P_S(f), h \rangle$ by solving a linear system. I've done this, and found out that this is a very good approximation of $\langle f, h \rangle$. By comparing with an online integration tool, the error is in the order of $10^{-6}$ for $r=1$. It gets better for larger $r$ and worse for lower $r$. Therefore it seems that $f, h \notin S$. There are variants for the definitions of the auxiliary functions leading to similar results. Possibly there are choices that reduce the error. This method can also be applied to small values of $r$ thanks to the formula at Carlo Beenakker's answer.