Asymptotic expansion of $f(x)=\sum _{n=1}^{\infty } \frac{\sin \left(\sqrt{n}x\right)}{n}$ at the origin

If, instead of the sum, we compute the integral (play later with inequalities) $$\int \frac{\sin \left(\sqrt{n} t\right)}{n} \,dn=2 \,\text{Si}\left(\sqrt{n} t\right)$$ $$\int_1^\infty \frac{\sin \left(\sqrt{n} t\right)}{n} \,dn=\pi\, \text{sgn}(t)-2 \,\text{Si}(t)$$

Close to $t=0$ $$\text{Si}(t)=t-\frac{t^3}{18}+O\left(t^5\right)$$

Initial Problem ...before considering the more general series $\;\displaystyle\sum _{n=1}^{\infty } \frac{\sin \left(n^a x\right)}{n^b}$.

First an approximate picture of $\,\displaystyle f(x):=\sum_{n=1}^{\infty}\frac{\sin \left(\sqrt{n}\,x\right)}n\,$ for $\displaystyle\,x\in\left(\frac 1{10},10\right)$ :

(the oscillations increase quickly at the right) $$-$$ We could continue by observing that $\;\displaystyle f''(x)=-\sum _{n=1}^{\infty }\sin \left(\sqrt{n}x\right)\;$ and see what was done for $x=1$ or follow this neat thread but I'll rather try the Euler–Maclaurin formula proposed by Mark Viola to get an actual asymptotic expansion at zero.

Euler–Maclaurin formula applied to the function $\;\displaystyle g_x(n):=\frac{\sin \left(\sqrt{n}\,x\right)}{n}\;$ using $\,p$ even Bernoulli terms is : $$\sum_{i=m}^n g_x(i) = \int^n_m g_x(t)\,dt + \frac{g_x(n) + g_x(m)}{2} + \sum_{k=1}^{p} \frac{B_{2k}}{(2k)!} (g_x^{(2k - 1)}(n) - g_x^{(2k - 1)}(m)) + R_{2p}(x)$$

From the definition of the sine integral we get $\;\displaystyle\lim_{n\to\infty}\int_m^n \frac{\sin \left(\sqrt{t}\,x\right)}t \,dt=\pi-2 \,\text{Si}(\sqrt{m}\,x)$
(the special case $m=1$ was indicated by Claude) $\;\displaystyle g_x'(n):=\frac{\sqrt{n}\,x\cos \left(\sqrt{n}\,x\right)-2\sin \left(\sqrt{n}\,x\right)}{2\,n^2}\;$ and so on so that at the limit $\,n\to +\infty\,$ we obtain : $$f_m(x)=\sum _{n=1}^{m-1} \frac{\sin \left(\sqrt{n}x\right)}{n}+\pi-2 \operatorname{Si}\left(\sqrt{m}\,x\right)+\frac {\sin \left(\sqrt{m}\,x\right)}{2\,m}-\frac 1{12}g_x'(n)+\frac 1{720}g_x'''(n)+\cdots$$

Notice that $\,g_x(n)\,$ as well as its derivatives and $\operatorname{Si}\left(\sqrt{m}\,x\right)$ may all be expanded as odd powers of $x$ so that the constant term $\pi$ will be kept unchanged.
I didn't fix $\,m=1\,$ to allow higher precision computation of the numerical factors of the $x^k$ appearing (the precision will decrease by adding too many $p$ binomial terms but will increase with $m$).

I won't provide these numerical factors since they all posses this (unproved) closed form that should please you : $$\boxed{\sum_{n=1}^{\infty}\frac{\sin \left(\sqrt{n}\,x\right)}n=\pi+\sum_{k=0}^{\infty}\zeta\left(\frac 12-k\right)\frac{(-1)^k\,x^{2k+1}}{(2k+1)!}}$$

...because this is nearly 'dual' to my other answer from your link (could be Deep!) : $$\boxed{\sum_{n=1}^\infty \frac{\sin(nx)}{\sqrt{n}}=\sqrt{\frac{\pi}{2\,x}}+\sum_{k=0}^\infty\zeta\left(-\frac 12-2k\right)\frac{(-1)^k\,x^{2k+1}}{(2k+1)!}}$$

---

GENERALIZATION and conjecture

"could be Deep" or perhaps not... but still...

Let's search the asymptotic expansion of $\;\displaystyle f_{a,b}(x):=\sum _{n=1}^{\infty } \frac{\sin \left(n^a x\right)}{n^b},\quad a,b\in\mathbb{R^+}$
(whose derivative should return the corresponding expansion for the $\,\cos\,$ function)

The expansion of the sine function is well known and has an infinite radius of convergence so let's somewhat carelessly change the order of summation : \begin{align} f_{a,b}(x)&=\sum _{k=1}^{\infty } \frac{\sin \left(k^a x\right)}{k^b}\\ &=\sum _{k=1}^{\infty } \frac 1{k^b}\sum _{n=0}^{\infty } \frac{(-1)^n}{(2n+1)!}\left(k^a x\right)^{2n+1}\\ &=\sum _{n=0}^{\infty } \frac{(-1)^nx^{2n+1}}{(2n+1)!}\sum _{k=1}^{\infty } \frac{\left(k^a\right)^{2n+1}}{k^b}\\ &\overset{?}=\sum _{n=0}^{\infty } \frac{(-1)^nx^{2n+1}}{(2n+1)!}\zeta\left(b-a(2n+1)\right)\\ \end{align} of course the previous equality is wrong since the last series at the right will converge to $\,\zeta(b-a(2n+1))\,$ only if $\;b-a(2n+1)>1\;$ which can't be satisfied for $a>0\;$ for all values of $n$ but the point seems to be that for different values of $a$ and $b$ we will simply have to use the previous formula but with a "corrective (regularized)" additional term :

• $\displaystyle\pi\quad$ for $\;\displaystyle a=\frac 12,\;b=1$
• $\displaystyle\sqrt{\frac{\pi}{2\,x}}\quad$ for $\;a=1,\;b=\displaystyle\frac 12$
• a more general expression for the general case

After some work let's expose my actual conjecture (with numerical verification for small fractions using Euler Maclaurin but no formal proof. The "regularized term" finally appearing as the only term of the Laurent series of $\;\displaystyle \int_1^{\infty} \frac{\sin \left(t^a x\right)}{t^b} \,dt\;$ not of the $\;a_n\,x^{2n+1},\;n\ge 0\;$ form) :

General expression for $\;a,\,b\in \mathbb{Q}^+$ and $\;\displaystyle f:=\frac {b-1}a\;$ :

• For $\;f= 2n\;$ with $\;n\ge 0\;$ integer : $$\tag{1}\boxed{\sum _{k=1}^{\infty } \frac{\sin \left(k^a x\right)}{k^b}=\frac{\pi}{2\,a}\frac{(-x^2)^n}{(2n)!}+\sum _{k=0}^{\infty } \zeta\left(b-a(2k+1)\right)\frac{(-1)^kx^{2k+1}}{(2k+1)!}}$$

• For $\;f= 2n+1\;$ with $\;n\ge 0\;$ integer (we have to take care of the singularity of $\zeta$ at $1$) : $$\tag{2}\boxed{\sum _{k=1}^{\infty } \frac{\sin \left(k^a x\right)}{k^b}=\sum _{k=0}^{\infty } \frac{(-1)^kx^{2k+1}}{(2k+1)!}\begin{cases} \zeta\left(b-a(2k+1)\right)&\text{if}\ \;\small b-a(2k+1)\neq 1\\ \gamma-\large\frac{\log(x)+\gamma\,-H_{2k+1}}a&\normalsize \text{else} \\ \end{cases} }$$

• else $$\tag{3}\boxed{\sum _{k=1}^{\infty } \frac{\sin \left(k^a x\right)}{k^b}=-\;\sin\left(\frac {\pi}2\,f\right)\,\Gamma(-f)\;\frac {x^{\,f}}a+\sum _{k=0}^{\infty } \zeta\left(b-a(2k+1)\right)\frac{(-1)^kx^{2k+1}}{(2k+1)!}}$$

where $\Gamma$ is the gamma function, $\gamma$ the Euler constant and $H_n$ the $n$-th harmonic number.

The first case is the limit of the general expression $(3)\,$ as $\,f\to 2n$.

In the second case we obtain $\;\displaystyle\gamma-\frac{\log(x)+\gamma\,-H_{2n+1}}a\;$ as the limit of $\;\displaystyle -\frac{(2n+1)!}{(-1)^nx^{2n+1}}\sin\left(\frac {\pi}2\,f\right)\,\Gamma(-f)\;\frac {x^{\,f}}a+\zeta(b-a(2n+1))\;$ as $\;f\to 2n+1\;$

(replace $\,b\,$ by $\,b+\epsilon\,$ and study the collision of the $\,x^{\,f}$ and $\,x^{2n+1}\,$ coefficients in $(3)\,$ as $\,\epsilon\to 0$).