Finding asymptotic for $S(X)$ amounts to estimation of $\sum_{(a,q)=1} aa^*$ for any given $q$.
From the rearrangement inequality, it follows that $$\frac{1}{2}q^2\varphi(q) - T = \sum_{(a,q)=1} a(q-a) \le \sum_{(a,q)=1} aa^* \le \sum_{(a,q)=1} a^2 = T,$$ where $$T = \sum_{(a,q)=1} a^2 = \sum_{d\mid q} \mu(d) d^2 \frac{q/d(q/d+1)(2q/d+1)}{6} = \frac{1}3 q^2\varphi(q) + O(q^2).$$ Hence, $$\frac{1}6 q^2\varphi(q) \lesssim \sum_{(a,q)=1} aa^* \lesssim\frac{1}3 q^2\varphi(q),$$ and the rest should follow as under$$\frac{1}6 q^2\varphi(q) \lesssim \sum_{(a,q)=1} aa^* \lesssim\frac{1}3 q^2\varphi(q).$$
UPDATE. Notice that the assumption of independent $a$ and $a^*$ (it essentially suggests that $\sum_{(a,q)=1} aa^* \sim\frac{1}4 q^2\varphi(q)$, which deviates from the actual asymptotic bounds above by constant factors). However, these asymptotic bounds are not tight enough to obtain the anticipated asymptotic formula for $S(X)$.