Asymptotic behaviour of sum over the inverse japanese symbol

If we exploit:

$$\mathcal{L}^{-1}\left(\frac{1}{\sqrt{x^2+\omega}}-\frac{1}{x}\right)=-1+J_0(s\sqrt{\omega})=\sum_{n\geq 1}\frac{s^{2n}\omega^n(-1)^n}{4^nn!^2} \tag{1}$$ we have:

$$ c(\omega)-\gamma = \int_{0}^{+\infty}\frac{J_0(s\sqrt{\omega})-1}{e^s-1}\,ds. \tag{2} $$

The Euler-Maclaurin Sum Formula applied to the power series in $\frac1m$ yields

$$ \begin{align} \sum_{m=1}^M\frac1{\sqrt{m^2+\omega}} &=\log(M)+c(\omega)+\frac1{2M}+\frac{3\omega-1}{12M^2}-\frac\omega{4M^3}\\ &+\frac{4+60\omega-45\omega^2}{480M^4}+\frac{3\omega^2}{M^5}+O\!\left(\frac1{M^6}\right)\tag{1} \end{align} $$

and $$ \sum_{m=1}^M\frac1m =\log(M)+\gamma+\frac1{2M}-\frac1{12M^2}+\frac1{120M^4}+O\!\left(\frac1{M^6}\right)\tag{2} $$ Note that as $\omega\to0$, $(1)\to(2)$, as long as $c(\omega)\to\gamma$.

Subtracting $(2)$ from $(1)$ and letting $M\to\infty$, we get $$ c(\omega)-\gamma =\sum_{m=1}^\infty\left(\frac1{\sqrt{m^2+\omega}}-\frac1m\right)\tag{3} $$ Therefore, $$ \begin{align} c(\omega) &=\gamma+\sum_{m=1}^\infty\left(\frac1{\sqrt{m^2+\omega}}-\frac1m\right)\\ &=\gamma+\frac1{\sqrt{1+\omega}}-1+\sum_{m=2}^\infty\left(\frac1{\sqrt{m^2+\omega}}-\frac1m\right)\\ &=\gamma+\frac1{\sqrt{1+\omega}}-1+\sum_{m=2}^\infty\frac1m\sum_{k=1}^\infty\binom{2k}{k}\left(-\frac\omega{4m^2}\right)^k\\ &=\gamma+\frac1{\sqrt{1+\omega}}-1+\sum_{k=1}^\infty\left(-\frac\omega4\right)^k\binom{2k}{k}\underbrace{(\zeta(2k+1)-1)}_{\sim2^{-2k-1}}\tag{4} \end{align} $$ That is,

$$ c(\omega)=\gamma+\frac1{\sqrt{1+\omega}}-1+\sum_{k=1}^\infty\left(-\frac\omega4\right)^k\binom{2k}{k}(\zeta(2k+1)-1)\tag{5} $$

The absolute value of the coefficient of $\omega^k$ in the summation is $\,\sim\!\frac1{4^k\sqrt{4\pi k}}$, so the sum converges for $\left|\omega\right|\lt4$. Therefore, $\lim\limits_{\omega\to0}c(\omega)=\gamma$, as required above.

---

Note that if we incorporate $\frac1{\sqrt{1+\omega}}-1$ into the summation in $(4)$ and move $\gamma$ to the right side, we get $$ c(\omega)-\gamma=\sum_{k=1}^\infty\left(-\frac\omega4\right)^k\binom{2k}{k}\zeta(2k+1)\tag{6} $$ which matches formula $(1)$ from Jack D'Aurizio's answer applied to the integral in $(2)$ from Jack D'Aurizio's answer. However, $(6)$ only converges for $-1\lt\omega\le1$.

---

Here is a plot of $c(\omega)$ for $0\le\omega\le1$:

If you use the Riemann integral trick, one gets: $$\sum_{m=1}^M\frac{1}{\sqrt{m^2+\omega}} \sim \ln\frac{2M}{1+\sqrt{\omega+1}}+\frac{1}{4}\frac{\omega}{M^2}+\mathcal{O}(\frac{\omega^2}{M^4})$$

To see how well the RHS function (which let us call it $fa(M)$, as the asymptotic form of the sum) approximates the sum, here is the ratio of the sum to its asymptotic (?) formula (ie, the RHS of our equation):

($N$ in the figure is our $M$). The red line is a fit to the ratio ($\omega=0.43$) .

Hope this helps you a bit.

Added note: if in the Riemann sum we sum up the trapezes (with areas slightly bigger than the area under the curve) instead of rectangles, the new formula (with absolute error almost an order of magnitude less) reads now: $$\sum_{m=1}^M\frac{1}{\sqrt{m^2+\omega}} \sim \ln\frac{2M}{1+\sqrt{\omega+1}}+\frac{1}{2\sqrt{\omega+1}} + \frac{1}{2M} + \frac{1}{4}\frac{\omega}{M^2}+\mathcal{O}(\frac{\omega}{M^3})$$