The reader's familiarity is assumed with the Euler-Mascheroni constant and the values at $s=\tfrac12$ of the Gamma, digamma and trigamma functions. Warning: you should double-check all my arithmetic, but this code suggests I'm right.
Let $I_k:=\int_{\Bbb R}\frac{1}{\sqrt{2\pi}}\exp(-\tfrac12x^2)\ln^k|x|dx$, so $R$ has mean $nI_1$ and $n(I_2-I_1^2)$. But$$I_k=\sqrt{\frac{2}{\pi}}\int_0^\infty\exp(-\tfrac12x^2)\ln^kxdx=\frac{1}{\sqrt{\pi}}\int_0^\infty y^{-1/2}\exp(-y)\ln^k(\sqrt{2}y^{1/2})dy.$$Next define$$J_l:=\int_0^\infty y^{-1/2}\exp(-y)\ln^lydy=\Gamma^{(l)}(\tfrac12),$$so special values of the gamma, digamma and trigamma functions give$$J_0=\sqrt{\pi},\,\frac{J_1}{J_0}=-\gamma-\ln4,\,\frac{J_2}{J_0}=(\gamma+\ln4)^2+\frac{\pi^2}{2}.$$Hence$$I_1=\frac{1}{\sqrt{\pi}}\left(\ln\sqrt{2}\Gamma(\tfrac12)+\tfrac12\Gamma^\prime(\tfrac12)\right)=\frac{-\gamma-\ln 2}{2}$$and$$I_2=\frac{1}{\sqrt{\pi}}(\ln^2\sqrt{2}\Gamma(\tfrac12)+\ln\sqrt{2}\Gamma^\prime(\tfrac12)+\frac14\Gamma^{\prime\prime}(\tfrac12))=\frac{\gamma\ln 2}{2}+\frac{\ln^22}{4}+\frac{\gamma^2}{4}+\frac{\pi^2}{8}.$$In particular,$$I_2-I_1^2=\frac{\pi^2}{8}.$$