The Legendre function of the second kind, $Q_n(z)$, along with the usual Legendre polynomial $P_n(z)$, are the two linearly independent solutions of the Legendre differential equation.

$Q_n(z)$ can be expressed in the following form:

$$Q_n(z)=P_n(z)\mathrm{artanh}\,z-W_{n-1}(z)$$

where $W_{n-1}(z)$ can be expressed either as

$$W_{n-1}(z)=\sum_{k=1}^n \frac{P_{k-1}(z) P_{n-k}(z)}{k}$$

or as

$$W_{n-1}(z)=\sum_{k=0}^{n-1} \frac{(H_n-H_k)(n+k)!}{2^k (n-k)! (k!)^2} (z-1)^k$$

where $H_k$ is the $k$-th harmonic number, $H_k=\sum\limits_{j=1}^k \frac1{j}$.

My questions:

1. Mathematica returns a rather complicated expression for $W_{n-1}(z)$ involving the unknown solution of a certain recurrence (i.e. DifferenceRoot[]). Is there possibly a simpler form for this polynomial?

2. Might there be a (hopefully simple) $n$-term recurrence that generates these polynomials?

Addendum:

After staring long and hard at Pietro's answer, I feel now that my second question was sorta kinda dumb; I already knew that both Legendre functions satisfied the same difference equation, so it stands to reason that a linear combination of them should also be a solution to that recurrence.

I now would like to expand my first question a bit: is it possible to express $W_n(z)$ as a single hypergeometric function (e.g. ${}_p F_q$ or some of the fancy multivariate ones), perhaps with one of the parameters being a negative integer? For instance, $P_n(z)$ is expressible as a Gaussian hypergeometric function ${}_2 F_1$ with one of the numerator parameters being a negative integer. Might there be something similar for the $W_n$?

I would also like to consider an additional question: are the $W_n$ orthogonal polynomials with respect to some weight function $\omega(x)$ and an associated support interval $(a,b)$? That is, if

$$\int_a^b\omega(t)W_j(t)W_k(t)\mathrm dt=0,\qquad j\neq k$$

for some $\omega(x)$ and some interval $(a,b)$, what is this weight function and its support interval?