Dirichlet Convolution of the Mobius Function with Itself

As you say, this is multiplicative. This means that once you know how to calculate it for numbers of the form $p^a$ you can calculate it for arbitrary $n$ by writing $n$ as a product of powers of different primes $p^aq^b\cdots$, and then multiplying the corresponding values $(\mu*\mu)(n)=((\mu*\mu)(p^a))((\mu*\mu)(q^b))\cdots$.

If $a=1$, $\sum_{d\mid p}\mu(d)\mu(p/d)=-2$ (the two factors each contribute $-1$). If $a=2$, $\sum_{d\mid p^2}\mu(d)\mu(p^2/d)=1$, since the only way for $\mu(d)\mu(p^2/d)$ to be non-zero is if $d\leq p$ and $p^2/d\leq p$, which requires $d=p$. If $a>2$ then $\sum_{d\mid p^a}\mu(d)\mu(p^2/d)=0$, since for each term in the sum either $d$ or $p^a/d$ is divisible by $p^2$.

So your function is $0$ if $n$ is divisible by the cube of any prime. Otherwise it is $(-2)^k$, where $k$ is the number of primes that divide $n$ exactly once (i.e. their squares do not divide $n$).

Note on the Previous Solution: Let $\operatorname{rad}(n)$ denote the radix, or squarefree part, of $n$, i.e., so that if $n = p_1^{\alpha_1} \cdots p_k^{\alpha_k}$ is the factorization of $n$ into powers of distinct primes then $\operatorname{rad}(n) = p_1 p_2 \cdots p_k$. Then we see that $\mu(\operatorname{rad}(n))$ is always non-zero and that $\mu^2(n / \operatorname{rad}(n)$ is the characteristic function of the cube-free integers. Also, if $\omega(n)$ counts the number of distinct prime factors of $n$ (the prime omega function) then the number of primes that divide $n$ exactly once is given by $$\#\{p \text{ prime } : p|n \text{ and } p^2 \nmid n\} = \omega(n) - \omega(n / \operatorname{rad}(n)).$$ This implies that the solution given by Essentially Lime in the previous response has the following exact closed-form expression for any $n \geq 1$: $$(\mu \ast \mu)(n) = (-2)^{\omega(n) - \omega(n / \operatorname{rad}(n))} \cdot \mu^2(\omega(n) - \omega(n / \operatorname{rad}(n)).$$ Curiously enough, parts of this formula for the solution can be re-expressed in terms of 1) the identity that $\sum_{d|n} \mu^2(d) = 2^{\omega(n)}$, or in other words the convolution $\mu^2 \ast 1 = 2^{\omega}$; and 2) the fact that $(-1)^{\omega(n)} = \mu(n)$ whenever $n$ is squarefree.

Alternate Solution Using Dirichlet Inverses: Another, perhaps shorter and more direct way to find this solution is to consider the Dirichlet inverse of the divisor function $d(n) := \sum_{d|n} 1 = (1 \ast 1)(n)$. Recall that the Dirichlet inverse $f^{-1}$ of any arithmetic function $f$ is the unique function (if one exists) satisfying $(f \ast f^{-1})(n) = \varepsilon = \delta_{n,1}$ where $\varepsilon$ is the multiplicative identity of Dirichlet convolutions. This inverse exists precisely when $f(1) \neq 0$.

Now consider the following phrasing of the problem: write $f := mu \ast \mu$. Then since $\mu \ast 1 = \varepsilon$ (an elementary fact which is seen here, for example, or proved by Dirichlet generating functions), we can invert the right-hand-side to obtain that $\varepsilon = f \ast 1 \ast 1 = f \ast d$. We can easily verify that $d(1) = 1 \neq 0$, so that the divisor function has a Dirichlet inverse. This implies that $f = d^{-1}$. The values of this inverse function are defined by $d(1) = 1$ and for $n > 1$ by the recurrence relation: $$d^{-1}(n) = -\sum_{\substack{d|n \\ d>1}} d(n) d^{-1}(n/d).$$ By computation, the first few values of this sum correspond to the sequence $\{d^{-1}(n)\}_{n \geq 1} = \{1,-2,-2,1,-2,4,-2,0,1,4,\ldots\}$ (A007427).