Power Series (Cambridge Tripos 1900)
If $\alpha,\beta$ are roots of $az^2+bz+c$ then we have $(a+bz+cz^2)/a=(1-\alpha z)(1-\beta z)$, so $$ p_n=\alpha^n+\alpha^{n-1}\beta+\dots+\alpha\beta^{n-1}+\beta^n=\frac{\alpha^{n+1}-\beta^{n+1}}{\alpha-\beta} $$ Thus $$ p_n^2=\frac{1}{(\alpha-\beta)^2}\left[\alpha^{2n+2}-2\alpha^{n+1}\beta^{n+1}+\beta^{2n+2}\right] $$ so $$ \begin{align*} \color{red}{(\alpha-\beta)^2}\sum_{n=0}^\infty p_n^2z^n &=\sum_{n=0}^\infty \alpha^{2n+2}z^n-2\sum_{n=0}^\infty (\alpha\beta)^{n+1}z^n+\sum_{n=0}^\infty\beta^{2n+2}z^n\\ &=\frac{\alpha^2}{1-\alpha^2z}-\frac{2\alpha\beta}{1-\alpha\beta z}+\frac{\beta^2}{1-\beta^2z}\\ &=\frac{\alpha^2+\beta^2-2\alpha^2\beta^2z}{(1-\alpha^2z)(1-\beta^2z)}-\frac{2\alpha\beta}{1-\alpha\beta z}\\ &=\frac{\color{red}{(\alpha-\beta)^2}(1+\alpha\beta z)}{(1-\alpha^2z)(1-\beta^2z)(1-\alpha\beta z)} \end{align*} $$ So $$ \sum_{n=0}^\infty p_n^2z^n=\frac{1+\alpha\beta z}{(1-(\alpha^2+\beta^2)z+\alpha^2\beta^2z^2)(1-\alpha\beta z)} $$ and this works even for the limiting case $\alpha=\beta$. Substituting $\alpha\beta=\dfrac ca$ and $\alpha^2+\beta^2=\dfrac{b^2-2ac}{a^2}$ gives $$ \begin{align*} \sum_{n=0}^\infty p_n^2z^n&=\frac{1+\dfrac{c}{a}z}{\left(1-\left(\dfrac{b^2-2ac}{a^2}\right)z+\dfrac{c^2}{a^2}z^2\right)\left(1-\dfrac{c}{a} z\right)}\\ &=\frac{a+cz}{a-cz}\frac{a^2}{a^2-(b^2-2ac)z+c^2z^2}\\ \end{align*} $$
If we start with the original function, we can reformulate is conveniently as $$ \frac{a}{a+bz+cz^2} = \frac{1}{1+(\frac{b}{a})z+(\frac{c}{a})z^2} = \frac{1}{(1-\frac{z}{z_+})(1-\frac{z}{z_-})} = \frac{z_+}{(z_+-z_-)} \frac{1}{(1-\frac{z}{z_-})} - \frac{z_-}{(z_+-z_-)} \frac{1}{(1-\frac{z}{z_+})} $$ with $z_\pm$ the roots of the denominator, i.e., $$ z_\pm = \frac{- b \pm \sqrt{b^2 - 4 a c}}{2 c} $$ If we expand the two expressions on the left as a series in $z$, we get $$ \frac{z_+}{(z_+-z_-)} \frac{1}{(1-\frac{z}{z_-})} = \sum_{n\geq0} l_n z^n \quad \text{ with } \quad l_n = \frac{z_+}{(z_+-z_-)} \frac{1}{z_-^n} $$ $$ \frac{z_-}{(z_+-z_-)} \frac{1}{(1-\frac{z}{z_+})} = \sum_{n\geq0} r_n z^n \quad \text{ with } \quad r_n = \frac{z_-}{(z_+-z_-)} \frac{1}{z_+^n} $$ and hence we identify $p_n=l_n-r_n$ for $n\geq0$. Now let us construct a series with coefficients $p_n^2 = (l_n-r_n)^2$. It is relatively easy to recognise that we can make the following identifications $$ \sum_{n\geq0} l_n^2 z^n = \frac{z_+^2}{(z_+-z_-)^2} \frac{1}{(1-\frac{z}{z_-^2})} $$ $$ \sum_{n\geq0} r_n^2 z^n = \frac{z_-^2}{(z_+-z_-)^2} \frac{1}{(1-\frac{z}{z_+^2})} $$ $$ \sum_{n\geq0} l_n r_n z^n = \frac{z_+ z_-}{(z_+-z_-)^2} \frac{1}{(1-\frac{z}{z_+ z_-})} $$ and hence that the required series $$ \sum_{n\geq0} p_n^2 z^n = \frac{1}{(z_+-z_-)^2} \left(\frac{z_+^2}{(1-\frac{z}{z_-^2})} + \frac{z_-^2}{(1-\frac{z}{z_+^2})} - \frac{2 z_+ z_-}{(1-\frac{z}{z_+ z_-})}\right) = \frac{z_+^2 z_-^2 (z_+ z_- + z)}{(z_+ z_- -z)(z_-^2 - z)(z_+^2 - z)} $$ Using that $z_+ z_- = a/c$ and $z_+^2+z_-^2=(b^2-2 ac)/c^2$ it then follows that $$ \sum_{n\geq0} p_n^2 z^n = \frac{a+c z}{a - c z} \frac{a^2}{a^2 - (b^2-2ac)z + c^2 z^2} $$