Find the limiting value of $S=a^{\sqrt{1}}+a^{\sqrt{2}}+a^{\sqrt{3}}+a^{\sqrt{4}}+…$ for $0 \leq a < 1$

Question

How can I find the sum of the series
$$S=a^{\sqrt{1}}+a^{\sqrt{2} }+a^{\sqrt{3}}+a^{\sqrt{4}}+…$$ under the condition $0 \leq a < 1$

Short version

I think the sum of the series should be about $$S=\frac{2} {(\ln a)^2}+c$$
where c is a correction term. So the question is how should I find the value of $c$?

Long version

My approach so far:

The sum can be restated in terms of $y=f(x)$ where
$$y_1=\sum_{n=1}^{\left \lfloor{x}\right \rfloor}a^{\sqrt{n}}$$
so the problem now is to find the limiting value of $y_1$ as $x \rightarrow \infty$. Instead of finding the limit of $y_1(x)$ as defined now, I tried to find another function with a similar growth rate:
$$\frac{\mathrm{d} y_2}{\mathrm{d} x}=a^{\sqrt{x}}$$
$${\textrm{i.e. }} y_2=\int a^{\sqrt{x}}{\mathrm{d} x} $$
which on solving yields
$$y_2=\frac{2 \sqrt{x} a^{\sqrt{x}}} {\ln a} – \frac{2a^{\sqrt{x}}} {(\ln a)^2}+\frac{2} {(\ln a)^2}$$
where the constant of integration is so chosen to make the curve pass through the origin.

Now plotting both these functions on Desmos screenshot
shows that $y_1$ and $y_2$ resemble each other closely. However $y_2$ has faster initial growth, and therefore it’s limiting value is slightly larger than the limiting value of $y_1$ which is what I’m after. The limiting value of $y_2$ is
$$S_2=\lim_{x \rightarrow \infty}{y_2(x)} =\frac{2} {(\ln a)^2}$$
Therefore, I think the actual sum of the original series should be:
$$S=\frac{2} {(\ln a)^2}+c$$
where $c$ is a correction term, possibly depending on $a$.

I’m stuck on how to find the correction term $c$ and what it’s form should be. Is my approach so far correct? Any help would be appreciated.

Solutions Collecting From Web of "Find the limiting value of $S=a^{\sqrt{1}}+a^{\sqrt{2}}+a^{\sqrt{3}}+a^{\sqrt{4}}+…$ for $0 \leq a < 1$"

We can keep going on the line of the comments for something that converges as close as necessary, starting from
$$
S = \frac{2}{(\log a)^2} – \frac12 + \frac1\pi \int_0^\infty \frac{e^{-\beta\sqrt{t}}\sin(\beta\sqrt{t})}{e^t-1}dt \quad\text{ where }\quad \beta = \frac{-\log a}{2\sqrt{\pi}}
$$
as given by @achille hui in the comments. We can try and get a form for the integral,
$$
I(\beta) = \int_0^\infty \frac{e^{-\beta\sqrt{t}}\sin(\beta\sqrt{t})}{e^t-1}dt
$$
by taking a Mellin transform of the integrand with respect to $\beta$ we get
$$
\int_0^\infty \beta^{s-1}I(\beta)\; d\beta = \Gamma(s)\sin\left(\frac{\pi s}{4}\right)\int_0^\infty \frac{(2t)^{-s/2}}{e^t-1} \; dt
$$
this is then
$$
\int_0^\infty \beta^{s-1}I(\beta)\; d\beta = 2^{-s/2}\Gamma(s)\sin\left(\frac{\pi s}{4}\right)\Gamma\left(1-\frac{s}{2}\right)\text{Li}_{1-\frac{s}{2}}(1)
$$
we want to take the inverse Mellin transform of both sides to get $I(\beta)$ but it’s not trivial to do this for the RHS, so invoke the Ramanujan Master Theorem, that the Mellin transform of $I(\beta)$ is written
$$
\int_0^\infty \beta^{s-1}I(\beta)\; d\beta = \Gamma(s)\phi(-s)
$$
such that
$$
I(\beta)=\sum_{n=0}^\infty \frac{(-1)^n}{n!}\phi(n)\beta^n
$$
then we have
$$
\phi(s)= \sin\left(\frac{-\pi s}{4}\right)2^{s/2}\Gamma\left(1+\frac{s}{2}\right)\text{Li}_{1+\frac{s}{2}}(1)
$$
and using $$
\text{Li}_{1+\frac{s}{2}}(1) = \zeta\left(1+\frac{s}{2}\right), \;\; s>0
$$
we can write
$$
I(\beta)=\sum_{n=0}^\infty \frac{2^{n/2}}{n!}\sin\left(\frac{-\pi n}{4}\right)\Gamma\left(1+\frac{n}{2}\right)\zeta\left(1+\frac{n}{2}\right)\left(\frac{\log{a}}{2\sqrt{\pi}}\right)^n
$$
finally
$$
\sum_{k=1}^\infty a^{\sqrt{n}}=\frac{2}{(\log a)^2} – \frac12 -\frac{\log(a)\zeta\left(\frac{3}{2}\right)}{4\pi} -\frac{\log(a)^2\zeta\left(2\right)}{4\pi^2}-\frac{\log(a)^3\zeta\left(\frac{5}{2}\right)}{32\pi^2}+\frac{\log(a)^5\zeta\left(\frac{7}{2}\right)}{512\pi^3}+\cdots
$$