From Gradshteyn & Ryzhik $3.692.6$ we know that $$\int_{-\infty }^{\infty } \left(\cos \left(\sqrt{x^2-1}\right)-\cos \left(\sqrt{x^2+1}\right)\right) \ dx=\pi (J_1(1)+I_1(1))$$ How can we establish it? Thanks for helping.
Update: I found a proof. One may start from the well-known formula $$\int_{-\infty }^{\infty } \frac{\sin \left(p \sqrt{a^2+x^2}\right)}{\sqrt{a^2+x^2}} \, dx=\pi J_0(a p)$$ Let $a=i a$, subtract it from the original result, differntiate with respect to $p$ yields: $$\int_{-\infty }^{\infty } \left(\cos \left(p \sqrt{x^2-a^2}\right)-\cos \left(p \sqrt{a^2+x^2}\right)\right) \, dx=\pi a (J_1(a p)+I_1(a p))$$ Now letting $p=a=1$ completes the proof. This also verified @skbmoore's generating function identity.
Furthermore, using the same technique in my answer of this post, a formula revealing beautiful symmetry is found: $$\int_{-\infty }^{\infty } \left(\cos \left(p \sqrt{x^2-a^2}\right)-\cos \left(p \sqrt{a^2+x^2}\right)\right) \, dx=\sum _{n=-\infty }^{\infty } \left(\cos \left(p \sqrt{n^2-a^2}\right)-\cos \left(p \sqrt{a^2+n^2}\right)\right)=\pi a (J_1(a p)+I_1(a p))$$
Here's a sketch of a proof. I intend to use Ramanujan's Master Theorem (RMT), which states that for a function $F(x)=\sum_{k=0}^\infty \phi(k) (-x)^k/k!,$ with $\phi(0) \neq 0$, then
$$ \int_0^\infty x^{n-1} F(x) dx = \Gamma(n) \phi(-n). $$ This is true for where the integral converges. A proof for non-integer n appears in L. Bougoffa, ArXiv 1902.01539v1, 5 Feb 2019. By splitting the integral at $x=0$ and scaling, we'll show the equivalent
$$ (1)\quad \int_0^\infty \Big( \cos{(\sqrt{2u-a^2}\ )} - \cos{(\sqrt{2u+a^2}\ )} \Big) \frac{du}{\sqrt{u}} = \frac{\pi \ a}{\sqrt{2} }\Big( J_1(a) + I_1(a) \Big) .$$ The formula of the OP is the special case of $a \to 1.$ Naturally we will let $n=1/2$ in RMT so we need the Taylor expansion of the function in the big parentheses. I worked out the first 40 terms (with a symbolic computer program) and discovered a pattern:
$$ (2) \quad \cos{(\sqrt{2u-a^2}\ )} - \cos{(\sqrt{2u+a^2}\ )} =\sqrt{\pi a/2} \sum_{k=0}^\infty \Big(-J_{k-1/2}(a)+I_{k-1/2}(a)\Big) (-u/a)^k/k! $$ which means $$ \phi(k) = \sqrt{\frac{\pi a}{2}} a^{-k} \Big(-J_{k-1/2}(a)+I_{k-1/2}(a)\Big). $$ Putting in $k=-1/2, \ \Gamma(1/2)=\sqrt{\pi}, $ and $ -J_{-1}(a) = J_1(a)$ gives the answer (1).
Of course I haven't proved (2), but it would surprise me if this expansion is not known. I have checked (1) numerically for many $0<a<\sqrt{2}$, using the PrincipalValue-> True argument in the numerical integration, with the point about which the principal value is taken at $u=a.$ Thus I think the suitable generalization is
$$ (3)\quad \int_{-\infty}^\infty \Big( \cos{(\sqrt{x^2-a^2}\ )} - \cos{(\sqrt{x^2+a^2}\ )} \Big) dx = \pi \ a\Big( J_1(a) + I_1(a) \Big), \quad 0<a\le 1 $$ as long as the integral is interpreted as a principal value.