How to calculate the asymptotic expansion of $\sum \sqrt{k}$? Denote $u_n:=\sum_{k=1}^n \sqrt{k}$. We can easily see that
$$ k^{1/2} = \frac{2}{3} (k^{3/2} - (k-1)^{3/2}) + O(k^{-1/2}),$$
hence $\sum_1^n \sqrt{k} = \frac{2}{3}n^{3/2} + O(n^{1/2})$, because $\sum_1^n O(k^{-1/2}) =O(n^{1/2})$.
With some more calculations, we get
$$ k^{1/2} = \frac{2}{3} (k^{3/2} - (k-1)^{3/2}) + \frac{1}{2} (k^{1/2}-(k-1)^{-3/2}) +  O(k^{-1/2}),$$
hence $\sum_1^n \sqrt{k} = \frac{2}{3}n^{3/2} + \frac{1}{2} n^{1/2} + C + O(n^{1/2})$ for some constant $C$, because $\sum_n^\infty O(k^{-3/2}) = O(n^{-1/2})$.
Now let's go further. I have made the following calculation
$$k^{1/2} = \frac{3}{2} \Delta_{3/2}(k) + \frac{1}{2} \Delta_{1/2}(k) + \frac{1}{24} \Delta_{-1/2}(k) + O(k^{-5/2}),$$
where $\Delta_\alpha(k) = k^\alpha-(k-1)^{\alpha}$. Hence :
$$\sum_{k=1}^n \sqrt{k} = \frac{2}{3} n^{3/2} + \frac{1}{2} n^{1/2} + C + \frac{1}{24} n^{-1/2} + O(n^{-3/2}).$$
And one can continue ad vitam aeternam, but the only term I don't know how to compute is the constant term.

How do we find $C$ ?


Let us substitute into the sum
$$\sqrt k=\frac{1}{\sqrt \pi }\int_0^{\infty}\frac{k e^{-kx}dx}{\sqrt x}. $$
Exchanging the order of summation and integration and summing the derivative of geometric series, we get
\begin{align*}
\mathcal S_N:=
\sum_{k=1}^{N}\sqrt k&=\frac{1}{\sqrt \pi }\int_0^{\infty}\frac{\left(e^x-e^{-(N-1)x}\right)-N\left(e^{-(N-1)x}-e^{-Nx}\right)}{\left(e^x-1\right)^2}\frac{dx}{\sqrt x}=\\&=\frac{1}{2\sqrt\pi}\int_0^{\infty}
\left(N-\frac{1-e^{-Nx}}{e^x-1}\right)\frac{dx}{x\sqrt x}=\\
&=\frac{1}{2\sqrt\pi}\int_0^{\infty}
\left(N-\frac{1-e^{-Nx}}{e^x-1}\right)\frac{dx}{x\sqrt x}.
\end{align*}
To extract the asymptotics of the above integral it suffices to slightly elaborate the method used to answer this question. Namely
\begin{align*}
\mathcal S_N&=\frac{1}{2\sqrt\pi}\int_0^{\infty}
\left(N-\frac{1-e^{-Nx}}{e^x-1}+\left(1-e^{-Nx}\right)\left(\frac1x-\frac12\right)-\left(1-e^{-Nx}\right)\left(\frac1x-\frac12\right)\right)\frac{dx}{x\sqrt x}=\\
&={\color{red}{\frac{1}{2\sqrt\pi}\int_0^{\infty}\left(1-e^{-Nx}\right)\left(\frac1x-\frac12-\frac{1}{e^x-1}\right)\frac{dx}{x\sqrt x}}}+\\&+
{\color{blue}{\frac{1}{2\sqrt\pi}\int_0^{\infty}
\left(N-\left(1-e^{-Nx}\right)\left(\frac1x-\frac12\right)\right)\frac{dx}{x\sqrt x}}}.
\end{align*}
The reason to decompose $\mathcal S_N$ in this way is that


*

*the red integral has an easily computable finite limit: since $\frac1x-\frac12-\frac{1}{e^x-1}=O(x)$ as $x\to 0$, we can simply neglect the exponential $e^{-Nx}$.

*the blue integral can be computed exactly.
Therefore, as $N\to \infty$, we have
$$\mathcal S_N={\color{blue}{\frac{\left(4n+3\right)\sqrt n}{6}}}+
{\color{red}{\frac{1}{2\sqrt\pi}\int_0^{\infty}\left(\frac1x-\frac12-\frac{1}{e^x-1}\right)\frac{dx}{x\sqrt x}+o(1)}},$$
and the finite part you are looking for is given by
$$C=\frac{1}{2\sqrt\pi}\int_0^{\infty}\left(\frac1x-\frac12-\frac{1}{e^x-1}\right)\frac{dx}{x\sqrt x}=\zeta\left(-\frac12\right).$$
