Evaluate $\sum _{n=1}^{\infty } \sin \left(\pi \sqrt{n^2+1}\right)$

Let's see: $$ \sin\left(\pi\sqrt{n^2+1}\right)=\sin\left(\pi n+\frac{\pi}{n+\sqrt{n^2+1}}\right)=(-1)^n\sin\left(\frac{\pi}{n+\sqrt{n^2+1}}\right) $$ hence the series is convergent by Leibniz' test. The collateral series $$ S(m)=\sum_{n\geq 1}(-1)^n\left(\frac{1}{n+\sqrt{n^2+1}}\right)^{2m+1} $$ are interesting objects, related to some series due to Ramanujan. By the (inverse) Laplace transform $$ S(0)=\sum_{n\geq 1}\frac{(-1)^n}{n+\sqrt{n^2+1}}=\int_{0}^{+\infty}\frac{J_1(s)}{s}\sum_{n\geq 1}(-1)^n e^{-ns}\,ds=-\int_{0}^{+\infty}\frac{J_1(s)}{s(e^s+1)}\,ds $$ and by the integral representation for the Bessel function $J_1$ $$ \frac{J_1(s)}{s}=\frac{1}{\pi}\int_{0}^{\pi}\cos(s\cos\theta)\sin^2(\theta)\,d\theta $$ such that $$ S(0) = -\frac{1}{2\pi}\cdot\Re\int_{0}^{\pi}\psi\left(\frac{2+i\cos\theta}{2}\right)-\psi\left(\frac{1+i\cos\theta}{2}\right)\,d\theta.$$ These real parts of digamma functions are extremely well-behaved on $[0,\pi]$, hence any numerical integration algorithm is able to find $S(0)\approx -0.271597$ with arbitrary accuracy. The same approach can be applied to $S(1),S(2),\ldots$ and the sequence $\{S(n)\}_{n\geq 0}$ roughly converges to zero like $\frac{1}{(1+\sqrt{2})^{2n}}$, so it is sufficient to compute $S(n)$ up to a small $n$ with good accuracy, then invoke interpolation to approximate $$\sum_{n\geq 1}\sin(\pi\sqrt{n^2+1})=\sum_{m\geq 0}\frac{\pi^{2m+1}(-1)^m}{(2m+1)!}S(m)\approx \color{red}{-0.566582}.$$