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.