Let $\mu(n)$ be the Möbius function and $S(x)$ be the number of positive integers $n \le x$ such that
$$ \sum_{r = 1}^{n-1} \mu(r)\gcd(n,r) = 0 $$
My experimental data for $n \le 6 \times 10^5 $seems to suggest that the number of solutions $\le x$ is growing at about
$$ S(x) \sim (\zeta(2)-1)\sqrt{x} $$ Is there any explanation for this? The square may come from the growth rate of the Mertens function $M(n) = \sum_{r \le n}\mu(r)$ while the appearance $\zeta(2)$ may be due to the fact that it appears in many sums involving the Pillai function $P(n) = \sum_{r \le n} \gcd(n,r)$.
I also observed that $200$ out of the $230$ zeroes for $n \le 1.3 \times 10^5$, were primes which indicate that the zeroes might be dominated by primes. For a prime $p$, $\sum_{r = 1}^{p-1} \mu(r)\gcd(p,r) = M(p-1)$ so I guess it is more likely that a sequence of $\pm 1$ adds to to $0$ than a sequence of integers of higher absolute value.
Update: Increased graph for the number of zeroes from $1.5 \times 10^5$ to $6 \times 10^5$
Here's an extended comment with some observations.
Quick summary is, I don't think the hypothesis is correct as $S(x)/\sqrt{x}$ starts to drop off as $x$ increases beyond $1,000,000$.
One way of thinking of this is in terms of $\mu(r)$ being "randomly" $\pm1$ for square-free $r$, so the cumulative sum of $\mu(r)$ would be like a random walk ... except the $\mu(r)$ will tend to cancel more frequently as, eg, $\mu(r)$ and $\mu(2r)$ will cancel for odd $r$. This would explain that $S(x)$ increases at a rate proportional to $\sqrt{x}$: the variance of the sum $\sum_{r=1}^n \mu(r)$ would have variance proportional to $n$ and thus a likelihood proportional to $1/\sqrt{n}$ of being zero. However, using this in a naive manner does not seem to work as I'll explain later, however it may still give some indications as to what is going on.
As is noted in the question, a large portion of the zeroes are for $n$ prime. This seems natural as all the terms in the sum then has $\gcd(r,n)=1$ and thus adds less variance. If $n$ has large proper divisors $d|n$, values of $r$ which are multiples of $d$ add more variance as they are multiplied by $d$.
However, if we look at the zeroes which are not prime, many of them are multiples of $3$: often $n=3p$ for some prime $p$. This is also natural, as the terms for multiples of $p$, $r=p$ and $r=2p$, cancel out.
There are more zeroes than just for $n=p$ and $n=3p$, but they tend to have in common that the partial sums for each divisor (ie for each different value of $\gcd(r,n)$) tends to be zero or close to zero.
It is not clear to me how the relative density of the different kinds of zeroes---primes, 3 times primes, and other composites---changes as $n$ increases. It seems to be fairly constant, perhaps with the "other composites" class falling slightly, although I find the numbers hard to judge. If that is the case, the number of zeroes for $n$ prime should be an indication of how fast $S(x)$ increases. However, if the likelihood of a prime $n$ being a zero is proportional to $1/\sqrt{n}$ and the density of primes is $~1/\ln n$, this should make $S(x)$ increase proportional to $\sqrt{x}/\ln x$.
I've run computations up to $n=5,000,000$ using the same formula as provided by Greg Martin at MathOverflow which allows partial sums per divisor to be cached for efficient computation: went reasonably fast even in Sage/Python.
What seems clear is that the ratio $S(x)/\sqrt{x}$ starts to drop as $x$ increases past $1,000,000$: for $x>2,500,000$, the ratio seems to stay below $0.5$, and I suspect it will keep dropping.