Let $$ P_0(x) = x - \lfloor x\rfloor\ - \frac{1}{2} = \left\{x\right\} - \frac{1}{2} $$
You can use Euler-Maclaurin summation:
$$ \int_1^{n} x^{1/x} dx = \sum_{k=1}^{n-1} \int_k^{k+1} x^{1/x} dx = \sum_{k=1}^{n-1} \bigg|_k^{k+1} P_0(x) x^{1/x} - \int_k^{k+1} P_0(x) \ x^{1/x} \ \frac{1- \ln(x)}{x^2} \ dx $$
$$ = \sum_{k=1}^{n-1} \left( \frac{1}{2} (k+1)^{1/(k+1)} + \frac{1}{2} k^{1/k} \right) + \int_1^{n} P_0(x) \ x^{1/x} \ \frac{\ln(x)-1}{x^2} \ dx $$
$$ = \sum_{k=1}^n k^{1/k} - \frac{1}{2} 1^{1/1} - \frac{1}{2} n^{1/n} + \int_1^{n} P_0(x) \ x^{1/x} \ \frac{\ln(x)-1}{x^2} \ dx $$
and so $$ \sum_{k=1}^n k^{1/k} = \int_1^n x^{1/x} dx + \frac{1}{2} + \frac{1}{2} n^{1/n} - \int_1^{n} P_0(x) \ x^{1/x} \ \frac{\ln(x)-1}{x^2} \ dx $$ Now, $P_0(x)x^{1/x}$ is bounded, and $\frac{\ln(x)-1}{x^2}$ behaves like $\frac{\ln(x)}{x^2}$ for large $x$, whose integral converges, so the right hand integral converges, and $n^{1/n} \to 1$ as $n \to \infty$, so:
$$ \sum_{k=1}^n k^{1/k} = \int_1^n x^{1/x} dx + \left( 1 - \int_1^{\infty} P_0(x) \ x^{1/x} \ \frac{\ln(x)-1}{x^2} \ dx \right) + o(1) $$
(note the term in the brackets is just a constant), so all that is left is to approximate the first integral:
$$ \int_1^n x^{1/x} dx = \int_1^n dx + \int_1^n \frac{\ln(x)}{x} dx + \int_1^n \left( x^{1/x} - 1 - \frac{\ln(x)}{x} \right) dx $$
$$ = n - 1 + \frac{\ln(n)^2}{2} + \int_1^n \left( x^{1/x} - 1 - \frac{\ln(x)}{x} \right) dx $$
A Taylor expansion shows that $$ \left( x^{1/x} - 1 - \frac{\ln(x)}{x} \right) \sim \frac{\ln(x)^2}{2x^2} $$ as $x \to \infty$, and so the final integral converges, so
$$ \int_1^n x^{1/x} dx = n + \frac{\ln(n)^2}{2} -1 + \int_1^{\infty} \left( x^{1/x} - 1 - \frac{\ln(x)}{x} \right) dx + o(1)$$
Putting it all together:
$$ \sum_{k=1}^n k^{1/k} = n + \frac{\ln(n)^2}{2} + \int_1^{\infty} \left( x^{1/x} - 1 - \frac{\ln(x)}{x} \right) dx - \int_1^{\infty} P_0(x) \ x^{1/x} \ \frac{\ln(x)-1}{x^2} \ dx +o(1) $$
i.e. $$\sum_{k=1}^n k^{1/k} = n + \frac{\ln(n)^2}{2} + c + o(1) $$ for a constant $c$. If you want even $\textit{further}$ terms in the asymptotic expansion, you can use this same method, and simply approximate each of the integrals to further accuracy, for the first integral by using more terms of the taylor expansion, for the second integral by integrating by parts using the Bernoulli polynomials, you can look them up in the Wikipedia page on Euler Maclaurin Summation.