Here's the integral:
$$I(n) = \frac{1}{\alpha}\int_{-\infty}^{\infty}{\ \frac{dx}{\sqrt{1 + x^{4}}}\ \ x\ \ \sin \lbrack n\pi\phi(x) \rbrack \ }$$
where $n$ is an integer and
$$\phi(x) \equiv \frac{1}{\alpha} \int_{0}^{x}\frac{dy}{\sqrt{1 + y^{4}}}\ \ \ \ ;\ \ \ \ \alpha \equiv \int_{0}^{\infty}\frac{dy}{\sqrt{1 + y^{4}}} = \frac{\Gamma\left( 1\text{/}4 \right)^{2}}{4\sqrt{\pi}} \ \ . $$
The integral converges because $\ \phi(x\to\infty)\to 1-O(1/x)\ \Rightarrow\ \sin( n\pi \phi )\to O(1/x)$.
Update: As suggested below by Claude, one start is to use the hypergeometric series for $\phi(x)$, invert that series using e.g. mathematica to get $x(\phi)$, plug that into the more compact form below, then do the integral. The result is an infinite series over coefficients each of which is a product of a hypergeometric function and a series inverse expansion coefficient, all of which may or may not be tractable, and it's unclear (to me) if that series will converge for large $n$ (which is what is of interest), or how to extract the asymptotic behavior from that series.
Another is the elliptic function pole expansion suggested by reuns below.
So the question is which if any of these is useful for evaluating integrals of the form above, particularly in the large-$n$ limit (or would a steepest descent approach work?). Numerically, the matrix element below, which is the ultimate goal, shows interesting large-$n$ behavior: $B(n,n)\sim \ln{n}$, $\ B(m,n)\sim 1\text{/}|m-n|^{\beta}$, with $1\text{/}2 < \beta < 1 \ $.
$$B(m,n) = \frac{1}{\alpha}\int_{-\infty}^{\infty}{\ \frac{dx}{\sqrt{1 + x^{4}}}\ x\ \cos \lbrack (m+1\text{/}2)\pi\phi(x) \rbrack \ \sin \lbrack n\pi\phi(x) \rbrack \ }$$
The goal is to derive that behavior.
Note: The more compact form mentioned above is:
$$I(n) = \int_{- 1}^{1}{d\phi}\ \ x(\phi)\ \sin\lbrack n\pi \phi \rbrack$$
which requires the inverse $x(\phi)$.
Side question: Simplest derivation of the identity for $\alpha$? – see comment by J.G. below.
Best Answer
The asymptotics of $I(n)=\int_{-1}^1 x(\phi)\sin n\pi\phi\,d\phi$ depends on the behavior of $x(\phi)$ at $\phi\to\pm1$ (simple poles, considering $x(\phi)$ analytically continued). Observe that $\phi(x)+\phi(1/x)=1$, hence $x(\phi)x(1-\phi)=1$. Thus the residues of $x(\phi)$ at $\phi\in\{\pm1\}$ are equal to $-1/\alpha$, as well as of $$y(\phi)=\frac{\pi}{2\alpha}\tan\frac{\pi\phi}{2},$$ so that $x(\phi)-y(\phi)$ is regular at $\phi\in\{\pm1\}$. This gives, as $n\to\infty$, $$I(n)\asymp\int_{-1}^1 y(\phi)\sin n\pi\phi\,d\phi=\frac\pi\alpha(-1)^{n-1}.$$
The behavior of $B(m,n)$ may be analysed similarly, using $$\frac\pi2\int_{-1}^1\tan\frac{\pi\phi}{2}\,\cos\left(m+\frac12\right)\pi\phi\,\sin n\pi\phi\,d\phi=(-1)^{m+n-1}(S_{m,n}+S_{m,n-1}),\\S_{m,n}=\sum_{k=-n}^n\frac1{2m+2k+1}.$$
UPDATE. Diving into elliptic functions, one gets an exact formula $$x(\phi)=\frac{\pi}{2\alpha}\left(\tan\frac{\pi\phi}{2}+4\sum_{n=1}^\infty\frac{\sin n\pi\phi}{e^{n\pi}+(-1)^n}\right)$$ (should I elaborate this?) and thus an exact formula for $I(n)$.