What is the mean value of a pair of Ramanujan Sums when summed over squares? Does anyone know of the mean value of two Ramanujan Sums when summed over the square of integers?
In my research on the Landau problem regarding nearly square primes, I have run into the mean value of a product of Ramanujan Sums:
\begin{equation*}
\lim\limits_{N \to \infty}
\frac{1}{N}
\sum\limits_{n =1}^{N} 
c_r\left( n^2 \right)
c_s\left( n^2 \right)
\end{equation*}
Does anyone know of a listing of properties of Ramanujan Sums similar to the this one, but with summations taken over the squares of positive numbers instead summation over the integers?  
The closest I can find is the function, $E_G$ found in Sums of products of Ramanujan sums by Toth. I am interested in the particular case where $g_1(n) = n^2$ though and not general case case represented by $E_G$ for and arbitrary set of functions $g_i(n)$
 A: Here is a partial answer (Added: completed below). Assuming $(r,s)=1$ and denoting $q=rs$, we have
$$ \sum_{n=1}^N c_r\left( n^2 \right) c_s\left( n^2 \right) 
= \sum_{n=1}^N c_q\left( n^2 \right) = \sum_{n=1}^N \sum_{d\mid(q,n^2)}\mu\left(\frac{q}{d}\right)d = \sum_{d\mid q}\mu\left(\frac{q}{d}\right)d\sum_{\substack{{1\leq n\leq N}\\{d\mid n^2}}} 1. $$
Now let $f(d)$ be the number of residue classes modulo $d$ whose square is zero modulo $d$. Then 
$$ \sum_{n=1}^N c_r\left( n^2 \right) c_s\left( n^2 \right) 
= \sum_{d\mid q}\mu\left(\frac{q}{d}\right)d f(d)\left(\frac{N}{d}+O(1)\right)
= N \sum_{d\mid q}\mu\left(\frac{q}{d}\right) f(d) + O_q(1),$$
whence
$$\lim_{N\to\infty}\frac{1}{N}\sum_{n=1}^N c_r\left( n^2 \right) c_s\left( n^2 \right) 
= \sum_{d\mid q}\mu\left(\frac{q}{d}\right) f(d).$$
The right hand side is the multiplicative convolution $g:=\mu\ast f$ evaluated at $q$, in particular it is multiplicative in $q$. Hence it suffices to determine $g$ at prime powers $p^k$, which is straightforward:
$$ g(p^k) = f(p^k)-f(p^{k-1}) = p^{\lfloor\frac{k}{2}\rfloor}-p^{\lfloor\frac{k-1}{2}\rfloor}.$$
We infer that the sought mean value equals
$$ g(q) = \begin{cases}\phi(\sqrt{q}),&q=\square\\0,&q\neq\square\end{cases} \tag{$\ast$}$$
Note that under our initial assumptions, $q$ is a square if and only if both $r$ and $s$ are squares.
Added. In the general case, i.e. without the coprimality condition on $r$ and $s$, we can proceed as follows. The function $n\mapsto c_r(n^2)c_s(n^2)$ is periodic mod $[r,s]$, hence the mean value exists and equals
$$ h(r,s):=\frac{1}{[r,s]}\sum_{n=1}^{[r,s]}c_r(n^2)c_s(n^2). $$
This function satisfies the multiplicativity relation
$$ h(rr',ss') = h(r,s)h(r',s'),\qquad (rs,r's')=1, $$
as follows from the Chinese Remainder Theorem and the invariance property
$$ c_q(m)=c_q(km),\qquad (k,q)=1. $$
Therefore,
$$ h(r,s) = \prod_{p\mid rs}h(p^{v_p(r)},p^{v_p(s)}), $$
so that it suffices to evaluate $h(p^k,p^l)$ for any prime $p$ and any exponents $k,l\geq 0$. This is straightforward, given the explicit formulae for $c_{p^k}(m)$ here. In particular, for $l>k\geq 0$ we obtain
$$ h(p^k,p^l)=\phi(p^k)h(1,p^l)
=\begin{cases}\phi(p^k)\phi(p^{l/2}),&\text{$l$ even}\\0,&\text{$l$ odd}\end{cases}$$
which generalizes the explicit formula ($\ast$) above.
