I can give a partial answer to this. My apologies if I am telling you what you already know.
First, let's review how the explicit formula for the Merten's function is derived. We have that
$$\frac{1}{\zeta(s)} = \sum_{n=1}^\infty \mu(n) n^{-s}.$$
We have
$$M(x) = \frac{1}{2 \pi i}\int_{2-i\infty}^{2+i\infty} \frac{1}{\zeta(s)} \frac{x^s}{s} dx,$$
and the residue theorem gives you the explicit formula you mentioned above. Note that the zeros of $\zeta$ become poles of $\zeta^{-1}$, so the two summations in your formula correspond to the nontrivial zeros and the trivial zeros of Riemann zeta function respectively. Also the $1/s$ term gives us a pole at $s=0$, giving us a residue of $1/\zeta(0) = -2$.
More generally, we can define a function $d_k(n)$ by the Dirichlet series given by $\zeta^k$:
$$\zeta(s)^k = \sum_{n=1}^\infty d_k(n) n^{-s}.$$
Then, for $k\geq 1$, $d_k(n)$ gives the number of representations of $n$ as a product of $k$ positive integers, given explicitly by your product formula above. In particular, $d_2(n)$ gives the number of positive divisors of $n$. Also, as you mentioned, $d_{-1}(n) = \mu(n)$, the Möbius function.
(For $k<-1$, I am have not checked whether my definition using the Dirichlet series will coincide with your product formula.)
Now you can attempt to find an explicit formula by carrying out the contour integral
$$D_k(x) = \sum_{n<x} d_k(n) = \frac{1}{2 \pi i}\int_{2-i\infty}^{2+i\infty} \zeta(s)^k \frac{x^s}{s} dx.$$
Now we can see that the explicit formula you will get will differ greatly depending on whether $k$ is positive or negative. If $k$ is positive, the main term will come from the residue of the pole at $s=1$, and this main term will be of the form $xP(\log x)$, where $P$ is a polynomial of degree $k-1$. If $k$ is negative, you'll pick up residues at all the zeros of $\zeta$, as in the case of $k=-1$.
Best Answer
The number of divisors of a square is the divisor function convolved with the square of the Möbius function (see $g(n)$ here)
$$e(n)=d(n^2)=(d\star\mu^2)(n)$$
by combining these three identities from the linked table
$$e=1\star1\star Q_2,\space d=1\star1,\space Q_2=\mu^2$$
and since
$$\mu^2(n)=\sum_{d^2|n}\mu(d)$$
therefore
$$e(n)=\sum_{a \left| n \right.} d \left( \frac{n}{a} \right) \sum_{b^2 \left| a \right.} \mu \left( b \right)$$
which can be simplified and rewritten as
$$e(n)=\sum_{b^2 \left| n \right.} \mu \left( b \right) d_3 \left( \frac{n}{b^2} \right)$$
where $d_3(n)$ is the number of ways that a given number can be written as a product of three integers. This identity can be verified by noting that $e(n)$ is multiplicative and checking at prime powers which yields $e(p^a)={2a+1}$ and can be compared with $d(p^a)={a+1}$. In particular note that $d_3(p^a)={\binom{a+2}{2}}$ (see $d_k$ here).
Then the summation of the number of divisors of the square numbers can be computed as:
$$E(x)=\sum_{n\le x} d(n^2) =\sum_{n \leq x} \sum_{b^2 \left| n \right.} \mu \left( b \right) d_3 \left( \frac{n}{b^2} \right)$$
which can be reorganized as:
$$E(x)=\sum_{b \leq \sqrt{x}} \mu \left( b \right) \sum_{n \leq x / b^2} d_3 \left( n \right)$$ $$E(x)=\sum_{a \leq \sqrt{x}} \mu \left( a \right) D_3 \left( \frac{x}{a^2} \right)$$
where $D_3$ is the summatory function for $d_3$. Since $D_3(x)$ can be computed in $O(x^{2/3})$ time complexity using the three dimensional analogue of the hyperbola method, $E(x)$ can also therefore be computed in
$$O(\int_{a=1}^\sqrt{x}{(\frac{x}{a^2})}^{2/3}da)=O(x^{2/3})$$
which is better than $O(x)$ as desired.
By taking an $O(x^{1/3})$ algorithm to compute $D(x)$ and using it for $D_3(x)$, this bound can be reduced to $O(x^{5/9})$. You can find such an algorithm and one formula for $D_3(x)$ in my article here which uses the notation $T(n)$ and $T_3(n)$.