what is the sum of square of all elements in $U(n)$? I know that $\sum\limits_{a\in U(n)} a=\frac{n\varphi(n)}{2}$ where $U(n):=\{1\leq r\leq n: (r, n)=1\}$ is a multiplicative group. And I know how to prove this result.
What I was willing to know was this $\sum\limits_{a\in U(n)} a^2$. is it possible to find in closed form?
what I tried is the following:
Let $S=\sum\limits_{a\in U(n)} a^2$. Now $(n, a)=1$ shows that $(n, n-a)=1$ which again under the fact  $(a, b)=1, (a, c)=1\Rightarrow (a, bc)=1$, shows that $(n, (n-a)^2)=1$. Hence $\{(n-a_1)^2, \cdots, (n-a_{\varphi(n)})^2\}$ is nothing but a permutation of the original set $\{a_1^2, \cdots, a_{\varphi(n)}^2\}$ in some order. Hence $S=\sum\limits_{a\in U(n)} (n-a)^2$. 
In other words we must have:
\begin{align*}
S=&\sum\limits_{a\in U(n)} (n^2-2an+a^2)\\
 =&n^2 \sum\limits_{a\in U(n)} 1-2n\sum\limits_{a\in U(n)} a+S\\
 =&n^2\varphi(n) -2n\times \frac{n\varphi(n)}{2}+S\\
 =&S
\end{align*}
and no result is obtained. 
What to do ? Please help me. 
Thanks
EDIT:
After the link has been provided below by Robert Israel
the formula reads $ n = p_1^{a_1}p_2^{a_2}\cdots p_r^{a_r}$ then $S = n^2\frac{\varphi(n)}{3}+(-1)^r p_1p_2\cdots p_r\frac{\varphi(n)}{6}$ but how to establish this ?

Lemma 1: For any function $f$ defined on rational values in $[0,1]$, denoting, $\displaystyle F(n) = \sum\limits_{k=1}^{n}f\left(\frac{k}{n}\right)$ and $\displaystyle F^{*}(n) = \sum\limits_{1\le k\le n;(k,n)=1} f\left(\frac{k}{n}\right)$.
We have, $\displaystyle 1*F^{*} = F$
Proof: First we observe that each integer $m \in \{1,2,\cdots,n\}$ has a unique representation in the form $m = k.\frac{n}{d}$, where, $d|n$, $1 \le k\le d$ and $(k,d) = 1$. 
Since, $\displaystyle \frac{m}{n} = \frac{k}{d}$, be the reduced fraction i.e., $(k,d) = 1$, each integer $m \le n$ has a representation in the required form, and representation is unique. 
Now, $\displaystyle \begin{align} (1*F^{*})(n) = \sum\limits_{d|n}F^{*}(d) = \sum\limits_{d|n}\sum\limits_{1\le k\le n;(k,d)=1} f\left(\frac{k}{d}\right) &= \sum\limits_{d|n}\sum\limits_{1\le k\le n;(k,d)=1} f\left(\frac{k\frac{n}{d}}{n}\right) \\ &= \sum\limits_{k=1}^{n}f\left(\frac{k}{n}\right) = F(n)\end{align}$
Now, take $\displaystyle f(x) = x^2$, we have $\displaystyle F(n) = \frac{\sum\limits_{k=1}^{n} k^2}{n^2} = \sum\limits_{d|n}\sum\limits_{1\le k\le n;(k,d)=1} \left(\frac{k}{d}\right)^2 = \sum\limits_{d|n} \frac{\phi_2(d)}{d^2}$
Where, $\displaystyle \phi_2(n) = \sum\limits_{k \le n;(k,n) = 1} k^2$
That is, $\displaystyle \frac{1}{6}n(n+1)(2n+1) = (\phi_2*f)(n) = g(n)$, (say)
Applying Dirichlet Inversion to the above identity, we have,
$\displaystyle \begin{align} \phi_2(n) = \sum\limits_{d|n} g(d)f^{-1}\left(\frac{n}{d}\right) \end{align}$
Since, $f$ is an absolutely multiplicative function the inverse is given by $f^{-1}(n) = \mu(n)f(n)$
(a proof of this fact can be found in T. Apostol, Introduction to Analytic Number Theory, page 36)
Thus, 
$\displaystyle \begin{align} \phi_2(n) &= \sum\limits_{d|n} g(d)\left(\frac{d}{n}\right)^2\mu\left(\frac{n}{d}\right) \\ &= \frac{n^2}{3}\sum\limits_{d|n} d\mu\left(\frac{n}{d}\right) + \frac{n^2}{2}\sum\limits_{d|n} \mu\left(\frac{n}{d}\right) + \frac{n}{6}\sum\limits_{d|n} \left(\frac{n}{d}\right)\mu\left(\frac{n}{d}\right) \tag{1} \\ &= \frac{n^2}{3}\phi(n) + \frac{n}{6}\sum\limits_{d|n}d\mu(d) \end{align}$
Where, in $(1)$, we used $\sum\limits_{d|n} d\mu\left(\frac{n}{d}\right) = \phi(n)$ and $\displaystyle \sum\limits_{d|n} \mu(d) = 0$, for $n > 1$.
Now, $\displaystyle \sum\limits_{d|n}d\mu(d) = \prod\limits_{p^\alpha || n} \left(\sum\limits_{e|p^{\alpha}}e\mu(e)\right) = \prod\limits_{p | n}(1-p)$
Thus, $\displaystyle \phi_2(n) = \frac{n^2}{3}\phi(n) + \frac{n}{6}\prod\limits_{p | n}(1-p)$ as required.
