Here is a partial answer. 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} $$ 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.