# Orthogonal polynomials with respect to the weighting function $\omega(x)=\frac{x}{e^x-1}$

I'd like to build a family of orthogonal polynomials with respect to the weighting function $$\omega(x)=\frac{x}{e^x-1}$$ on $$[0,+\infty)$$, i.e the inner product

$$=\displaystyle\int_{0}^{\infty}\frac{xP_n(x)P_m(x)}{e^x-1}dx$$

And more specifically, I would like to have an explicit expression of these polynomials, of the form $$P_n(x)=\displaystyle\sum\limits_{k=0}^{n}p_{k,n}x^k$$

Considering that we have the identity $$\displaystyle\int_{0}^{\infty}\frac{x^n}{e^x-1}dx=n!\zeta(n+1)$$, this seems doable.

I tried orthogonalizing $$(1,X,X^2...)$$ with Gram-Schmidt, but I end with a very bad-looking recurrence relation that is impossible to work with. The thing is, it seems to me that Gram-Schmidt always produces those horrible-looking recurrence relations that can't lead any further.

But every known orthogonal polynomials (Legendre, Chebyshev, Laguerre, Hermite...) have nice explicit expressions, and I've even seen casually thrown in a paper "and by the way, the orthogonal polynomials with respect to this particular weighting function on $$(0,4)$$ can be expressed as follows" (about the function $$x:\to\sqrt{\frac{4+x}{4-x}}$$) with an explicit expression.

All of this leads me to believe that their is a generic method to find it. But I could find any references for their derivation, they just seems to be considered as "well-known facts".

I thought that maybe there are derived from their Rodrigues formula, but the Chebyshev polynomials don't have one and still have nice expressions.

Do you have any idea how to find the explicit expression ?

A little too long for a comment : there is no generic method to find a closed form for orthogonal polynomials wrt to a particular inner-product. It is the opposite, to obtain Hermite,Laguerre,Legendre they started from a closed-form two variable analytic function $$f(x,t) = \sum_n g_n(x) t^n$$ satisfying $$\int_a^b f(x,t) f(x,u) w(x) dx = h(ut)$$ for some $$w,h$$, which implies $$\int_a^b g_n(x)g_m(x)w(x)dx = \cases{ h^{(n)}(0) \ \text{ if } n=m\\ 0 \text{ otherwise}}$$ The triplet $$f,w,h$$ exists in closed-form only in a very few special cases.
Given the Laguerre closed-form for $$w(x) = e^{-x}$$ to find the orthogonal polynomials for your $$w_2(x)=\frac1{e^x-1}$$ inner product you'll need to invert the infinite matrix $$M_{nm} = _{w_2}$$