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. **P.S.** I am sure the general case can be treated similarly, but I have not attempted this (for lack of time), and I expect the final formula to be less elegant. **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][1]. 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. [1]: https://en.wikipedia.org/wiki/Ramanujan's_sum