Stirling type formula for Sum on $\ln(n)^2$ Is there a similar formula like the Stirling one on the sum over $\ln(n)$ (take logarithms on its factorial representation):
$$\sum_{n=1}^N \ln(n) = N \ln(N)-N+\ln(N)/2+\ln(2\pi)/2+\mathcal{O}(\ln(N)/N)$$ 
but on the sum over its squares?
$$\sum_{n=1}^N (\ln(n))^2$$ 
I already advanced on getting good approximation on asymptotics integrating $\ln(n)^2$ and arrive to correct terms till $\mathcal{O}(N)$ order. But further advance is becoming hard for me in $\mathcal{O}(\ln(N))$ terms. 
I am specially interested in $\mathcal{O}(1)$ term.

We can use the Euler-Maclaurin formula to obtain the asymptotic expansion
$$ f_N (a) \equiv \sum \limits_{n=1}^N n^a \sim \zeta(-a) + \sum \limits_{k=0}^\infty \frac{B_k^*}{a+1} {a+1 \choose k} N^{a+1-k} $$
for $N \in \mathbb{N}$ and $a \in \mathbb{R} \setminus \{-1\}$ . Here $(B_k^*)_{k\in \mathbb{N}_0}$ are the Bernoulli numbers with $B_1^* = \frac{1}{2}$. We want to find an asymptotic expansion for 
$$f_N''(0) = \sum \limits_{n=1}^N \ln^2 (n) \, .$$ 
Obviously, the contribution of the terms with $k \geq 2$ vanishes as $N \to \infty$ . More precisely, their leading term is $\mathcal{O} (\ln^2 (N) /N)$ . Therefore we have
\begin{align}
\sum \limits_{n=1}^N \ln^2 (n) &\sim \frac{\mathrm{d}^2}{\mathrm{d}a^2} \left[\zeta(-a) + \frac{N^{a+1}}{a+1} + \frac{1}{2} N^{a} \right] \Bigg\rvert_{a=0} + \mathcal{O} \left(\frac{\ln^2 (N)}{N}\right) \\
&\sim N \ln^2 (N) - 2 N \ln(N) + 2N + \frac{1}{2} \ln^2 (N) + \zeta''(0) + \mathcal{O} \left(\frac{\ln^2 (N)}{N}\right)
\end{align}
as $N \to \infty$ .
$\zeta''(0)$ can be found using the series expansions
\begin{align}
(2\pi)^s &= 1 + \ln(2 \pi) s + \frac{1}{2} \ln^2 (2 \pi) s^2 + \mathcal{O}(s^3) \, , \\
\sin \left(\frac{\pi s}{2}\right) &= \frac{\pi}{2} s - \frac{\pi^3}{48} s^3 + \mathcal{O}(s^5) \, , \\
\Gamma(1-s) &= 1 - \gamma s + \frac{6 \gamma^2 + \pi^2}{12} s^2 + \mathcal{O}(s^3) \, , \\
\zeta(1-s) &= - \frac{1}{s} + \gamma + \gamma_1 s + \mathcal{O}(s^2)
\end{align}
near $s = 0$ in the functional equation
\begin{align}
\zeta(s) &= \frac{1}{\pi} (2\pi)^s \sin \left(\frac{\pi s}{2}\right) \Gamma(1-s) \zeta(1-s)\\
&= - \frac{1}{2} - \frac{1}{2} \ln (2 \pi) s - \frac{1}{2} \left[\frac{\pi^2}{24} + \frac{1}{2} \ln^2(2\pi) - \frac{1}{2} \gamma^2 - \gamma_1 \right] s^2 + \mathcal{O}(s^3)\, .
\end{align}
We obtain $\zeta(0) = -\frac{1}{2}$ , $\zeta'(0) = -\frac{1}{2} \ln(2 \pi)$ and
$$ \zeta''(0) = - \left[\frac{\pi^2}{24} + \frac{1}{2} \ln^2(2\pi) - \frac{1}{2} \gamma^2 - \gamma_1 \right] \, . $$
