I'm trying to find the ordinary generating function $f(n)$ such that $$f(n) = \sum_{k=0}^\infty\mu^2(n)\, x^k$$
I have a nasty looking answer that involves Hadamard products, and I was hoping there was a cleaner, canonical answer.
The Hadamard solution, for those interested:
\begin{align*} e_{n}(x) & = \sum_{k=0}^\infty \mu(x) \, \mathbf{1}_{k \leq n} x^k \\ f_n(x) & = \sum_{k=0}^\infty \mu(x) \, \mathbf{1}_{k \leq n} x^k \\ g_n(x) & = \sum_{k=0}^\infty \mu^2(x) \,\mathbf{1}_{k \leq n} x^k \\ & = \sum_{k=0}^\infty\bigg(\mu(x) \, \mathbf{1}_{k \leq n}\bigg)\bigg( \mu(x)\, \mathbf{1}_{k \leq n} \bigg)x^k \end{align*}
where we recognize the last line as the Hadamard product of the power series $e_n, \, f_n$ (i.e. the term-wise product).
According to the Wikipedia page above, if closed forms are known for $e_n, f_n$ then
$$ g_n(x) =\frac1{2\pi} \int_{0}^{2\pi} e_n\big(\sqrt{x} e^{xt}\big) \, f_n\big(\sqrt{x} e^{-xt}\big) dt$$ where
$$e_n(x) = f_n(x) = \sum \limits_{n=1}^{\infty} \mu(n)x^n = x - \frac{x^{2}}{1-x} + \sum \limits_{a=2}^{\infty} \frac{x^{2a}}{1-x^{a}} - \sum \limits_{b=2}^{\infty} \sum \limits_{a=2}^{\infty} \frac{x^{2ab}}{1-x^{ab}} + \sum \limits_{c=2}^{\infty} \sum \limits_{b=2}^{\infty} \sum \limits_{a=2}^{\infty} \frac{x^{2abc}}{1-x^{abc}} - \sum \limits_{d=2}^{\infty} \sum \limits_{c=2}^{\infty} \sum \limits_{b=2}^{\infty} \sum \limits_{a=2}^{\infty} \frac{x^{2abcd}}{1-x^{abcd}} + ...,$$ a known result.
Here are some ideas which I got, starting with Lambert series.
This first part is in https://en.wikipedia.org/wiki/Lambert_series
A zeta series of a sequence is $Z_f(s) =\sum_{n=1}^{\infty} \dfrac{f_n}{n^s} $.
A Lambert series is of the form $S_a(q) =\sum_{n=1}^{\infty} a_n\dfrac{q^n}{1-q^n} =\sum_{m=1}^{\infty} b_mq^m $ where $b_m =\sum_{n|m} a_n $.
Examples are for the Mobius function, $S_{\mu}(q) =\sum_{n=1}^{\infty} \mu(n)\dfrac{q^n}{1-q^n} =q $ and for Liouville's $\lambda$ function, $S_{\lambda}(q) =\sum_{n=1}^{\infty} \lambda(n)\dfrac{q^n}{1-q^n} =\sum_{m=1}^{\infty} q^{m^2} $ (a theta function).
The reason I mentioned the last series is because of this which is in https://en.wikipedia.org/wiki/Liouville_function:
"The Liouville function's Dirichlet inverse is the absolute value of the Möbius function."
If you use another result in the article above, $\dfrac{\zeta(2s)}{\zeta(s)} =\sum_{n=1}^{\infty} \dfrac{\lambda(n)}{n^s} $, this means that, if $Z_{\mu^2}(s) =\sum_{n=1}^{\infty} \dfrac{\mu^2(n)}{n^s} $, then $Z_{\mu^2}(s)\dfrac{\zeta(2s)}{\zeta(s)} =1 $ so $Z_{\mu^2}(s) =\dfrac{\zeta(s)}{\zeta(2s)} $.
Since $\mu^2(n) = |\mu(n)|$, this is the zeta series for your function.
Combining all these, you might be able to get what you want.
I'll stop here.