A recent question on math overflow on sums of the primes squared was answered/put on hold by pointing to an old paper by T. Šalát and S. Znám, On the sums of the prime powers. Šalát and Znám's (SZ) asymptotic formula for $ S(n, k) = 2^k + 3^k + 5^k + \ldots p_n^k$ is $n^{k+1}(\ln(n))^{k}/(k+1)$

Coding this result up, I found that this result was always smaller by a factor that grows to around 8% (for $n < 78498$) than the real sum. So I decided to try a different approach inspired by Dusart who has a better approximation for the $n$-th prime as $n (\ln n + \ln \ln n - 1)$. So I (heuristically) replaced $\ln n$ in the SZ formula by $(\ln n + \ln \ln n - 1)$. But then I found that the approximate sums now were too large (by sometimes as much as the values were too small for the SZ formula). So I decided to balance between these two extremes using a factor $c$ as in $\ln n + c ( \ln \ln n - 1)$. I found for the prime sums that $c = 0.6$ seemed to work very well, but for the squares, I found $c = 0.7$ seemed to work well and for the cubes $c = 0.8$ worked well. So I wonder if the "best" constant (i.e. value of $c$) was approaching $1$ for larger powers. The results for prime sums, for example, had errors less than 2% for sums over 100 primes (and less than 1% for $1959 < n < 78498$) whereas either the SZ or Dusart methods had errors of over 10% at times. Of course, the idea of a best constant $c$ is vague: the errors in my data go up and down as one advances through the primes and oscillates between overestimating and underestimating the real sum (for example, at $n = 15885$ for $k = 1$ and $c = 0.6$).

So just as there is a great deal of error analysis for the N-th prime (e.g., Rosser and Schoenfield, Approximate Formulas for some functions of prime numbers, Dusart etc.) is there some reason why:

$S(n, k)$ is approximated well by $n^{k+1}(\ln(n) + c_k \ln (\ln (n) - c_k)^{k}/(k+1)$ for values of $c_1 \approx 0.6$ and $c_2 \approx 0.7$ and where $c_k$ seems to approach $1$ as $k$ gets larger,

(where "approximated well" needs a definition :=)) Note: The sum of primes is A007504 in the OEIS; there is a more refined error analysis for $k=1$ by Shevelev of $S(n) \le n(n+1)(\ln(n) + \ln(\ln(n))+ 1)/2$ but even this upper bound is closer to the Dusart-based approximation. There is a similar error analysis by Shevelev for the case of $k = 2$ in A024450 but again this does not seem to answer the question above.