Determining the Lambert series for $xq+x^2q^4+x^3q^9+...+x^nq^{n^2}+...$ I am trying to determine the polynomials $P_n(x)$ from
$$
xq+x^2q^4+x^3q^9+...+\ x^nq^{n^2}+...=\sum_{n\geqslant1}\frac{P_n(x)q^n}{1-xq^n};
$$
that is,
$$
\sum_{d|n}x^{\frac nd-1}P_d(x)=\begin{cases}x^k,&n=k^2,\\0&\text{otherwise.}\end{cases}
$$
Or if you prefer, the Dirichlet generating function is
$$
\sum_{n\geqslant1}\frac{P_n(x)}{n^s}=\frac{x\operatorname{Li}_{2s}(x)}{\operatorname{Li}_s(x)}.
$$
In particular ($p$, $p_1$, $p_2$, $p_3$ are primes),
\begin{align*}
P_1(x)&=x\\
P_p(x)&=-x^p\\
P_{p^2}(x)&=-x^{p^2}+x^{2p-1}+x^p\\
P_{p_1p_2}(x)&=-x^{p_1p_2}+2x^{p_1+p_2-1}\\
P_{p^3}(x)&=-x^{p^3}+2x^{p^2+p-1}-x^{3p-2}-x^{2p-1}\\
P_{p_1p_2^2}(x)&=-x^{p_1p_2^2}+2x^{p_1+p_2^2-1}+2x^{p_1p_2+p_2-1}-3x^{p_1+2p_2-2}-x^{p_1+p_2-1}\\
P_{p_1p_2p_3}(x)&=-x^{p_1p_2p_3}+2x^{p_1p_2+p_3-1}+2x^{p_1p_3+p_2-1}+2x^{p_2p_3+p_1-1}-6x^{p_1+p_2+p_3-2}
\end{align*}
but how to handle the whole thing?
In fact ultimately I only need an explicit (well...) expression for $P_n'(1)$. Note that $P_n(1)$ is the Liouville function.
 A: This problem crucially involves the analogue of the Mobius function $\mu(n)$ for polylogarithms. However, the "polylogarithmic" analogues $\mu_n(x)$, say with $\mu(n)=\mu_n(1)$, are generally not multiplicative, so one cannot use the Fundamental Theorem of Arithmetic to determine them (this is one of the properties Lambert series work nicely with). I'm pretty sure no non-recursive (see below) formula is known for $\mu_n(x)$ and, therefore, neither for the polynomials $P_n(x)$. 
Let $\mu_n(x)$, $n\in\mathbb{N}$, be the unique sequence of polynomials in $x$ of degree $n$ with the property that 
$$\sum_{d|n}\mu_{n/d}(x)x^{d-1}=\delta(n)$$
so $\mu_n(x)$ is the Dirichlet inverse of $x^{n-1}$, in other words,
$$x{Li}^{-1}_s(x)=\sum_1^{\infty}\frac{\mu_n(x)}{n^s}.$$
Forming the Dirichlet product of the Dirichlet product in the question (your second displayed equation) with $\mu_{n}$ gives an expression for $P_n$ explicitly in terms of $\mu_n$, that is
$$P_n(x)=\sum_{d|n: d=k^2}\mu_{n/d}(x)x^{\sqrt{d}}.$$
Now let 
$$f*g(n)=\sum_{d|n}f\left(\frac{n}{d}\right)g(d)$$
 denote the Dirichlet convolution of $f$ and $g$, and 
$$f^{* j}(n)=\underbrace{f* f* \cdots *f}_{j\text{ copies of }f}(n)$$
denote Dirichlet self-convolution to order $j$, that is, the Dirichlet convolution of $f$ with itself $j-1$ times. The reason I claim that no non-recursive solution to your problem is known is simply that, when 
$f(n)={n\choose m}$, I believe a closed form for the Dirichlet self convolution to order $j$ is not known, i.e, one must actually carry out $j$ nested summations, and this is related to the fact that the binomial coefficient is not multiplicative (of course this could be incorrect so if anyone knows otherwise please comment on the solution here). 
Anyway, if 
$$\frac{\mu_n(x)}{x}=\sum_0^{n}c_k(n)(x-1)^k,$$
then 
$$c_k(n)=\sum_{j_1+2j_2+\cdots=k}\frac{(-1)^{j_1+j_2+\cdots}}{j_1!j_2!\cdots}\prod_{p|.}{j_1+j_2+\cdots\choose \alpha_p}*{.\choose 1}^{* j_1}*{.\choose 2}^{* j_2}*\cdots(n)$$
where $\alpha_p$ is the exponent of $p$ in the factorisation of $n$ and the summation extends over all partitions of $k$. This is derived using Faa di Bruno's formula for derivatives of composite functions applied to the reciprocal of $Li_s(x)$. From here you can easily get a formula for $P_n(x)$. 
A: Decided to add an expression for $P_n'(1)$ as I finally got it simplified to certain extent. I decided still to leave no answers approved as (1) I would be grateful for any possible further alternative expressions for $P_n'(1)$ and (2) the question itself is broader and seems to be interesting to others in this generality. Besides, this expression for $P_n'(1)$ somehow gives hopes for (and maybe hints at) further simplifying $P_n(x)$ along the lines of the answer given by Kevin Smith: after all, knowing $n$ values like $P_n(k)$, $P_n'(l)$, $P_n''(m)$, ... one could reconstruct $P_n(x)$, right?
So what I obtained is
$$
P_n'(1)=\lambda(n)\surd(n)-\sum_{1<d|n}\lambda\left(\frac nd\right)\varphi(d)
$$
where $\lambda$ is the Liouville function, $\varphi$ is the Euler totient function and $\surd(n)$ is the square root of the largest square dividing $n$.
This may be seen from
\begin{align*}
\sum_{n\geqslant1}\frac{P_n'(1)}{n^s}=\left.\frac d{dx}\left(\frac{x\text{Li}_{2s}(x)}{\text{Li}_s(x)}\right)\right|_{1\leftarrow x}&=\left.\frac{\text{Li}_{2s-1}(x)+\text{Li}_{2s}(x)}{\text{Li}_s(x)}
-\text{Li}_{2s}(x)\frac{\text{Li}_{s-1}(x)}{\text{Li}_s(x)^2}\right|_{1\leftarrow x}\\
&=\frac{\zeta(2s-1)}{\zeta(s)}+\frac{\zeta(2s)}{\zeta(s)}-\frac{\zeta(2s)}{\zeta(s)}\frac{\zeta(s-1)}{\zeta(s)}.
\end{align*}
