Here is a generating function approach. Note that the RHS is almost a convolution $\sum_{k=0}^n a_{k}b_{n-k}$, which suggests the use of generating functions, after some bricolage on the the recursion aimed to get a convolution form (whence a Cauchy product of power series).
To start with, note that the first term on the RHS may be naturally included in the sum for $k=1$. The binomial has already a suitable form; and then note the exponent of $(1-p)$: it can be written $k(n-k)= {n^2\over2}-{ (n-k)^2\over 2} -{k^2\over2} $, which suggests to multiply the relation by $ (1-p)^{-{n^2\over2}}{x^{n-1}\over(n-1)!}$ : $$(1-p)^{-{n^2\over2}}P_{n}{x^{n-1}\over (n-1)!}= p^{k(k-1)\over 2 }(1-p)^{-{k^2\over 2}}{x^{k-1}\over(k-1)!}\cdot P_{n-k}{x^{n-k}\over (n-k)!} $$
If we introduce the generating functions
$$g(x):=\sum_{n=1}^\infty p^{n(n-1)\over 2 }(1-p)^{-{n^2\over 2}}{x^{n}\over n!}$$ $$f(x):=\sum_{n=0}^\infty (1-p)^{-{n^2\over2}}P_{n}{x^{n}\over n!} $$ the recurrence relations write $f'=g'f$, whence $$f(x)=\exp\Big( \sum_{n=1}^\infty p^{n(n-1)\over 2 }(1-p)^{-{n^2\over 2}}{x^{n}\over n!}\Big).$$ This is possibly still far from your needs, but could hopefully be of use. Further steps would require some reference (hopefully somebody here does have it!) e.g. on the function $$G(x,z):=\sum_{n=1}^\infty {z^{n^2}x^{n}\over n!}$$ since $g(x)=G\big(x p^{-{1\over 2}},\big({p\over 1-p}\big)^{1\over2}\ \big)$.