I was exploring some raising and lowering operators related to an infinitesimal generator for fractional integro-derivatives and found an Appell sequence of polynomials, i.e., an infinite sequence of polynomials for which $\frac{d}{dx}p_n(x)=np_{n-1}(x)$, that is defined by the following recursion relation:
$p_{0}(x)=1$, $p_{1}(x)=x+\gamma$, and for $n>0$ $$p_{n+1}(x)=(x+\gamma)p_{n}(x)+\sum_{j=1}^{n}(-1)^j\binom{n}{j}j!\zeta (j+1)p_{n-j}(x)$$
where $\gamma=-\frac{\mathrm{d} }{\mathrm{d} \beta }\beta !\mid_{\beta =0 }$, the Euler-Mascheroni constant, and $\zeta(s)$ is the Riemann zeta function.
They satisfy $$p_{n}(x)=\frac{\mathrm{d^n} }{\mathrm{d} \beta^n }\frac{\exp(\beta x)}{\beta !} \mid_{\beta =0 }.$$
Explicitly,
$$p_2(x)=(x+\gamma)^2-\zeta(2)$$ $$p_3(x)=(x+\gamma)^3-3\zeta(2)(x+\gamma)+2\zeta(3)$$ $$p_4(x)=(x+\gamma)^4-6\zeta(2)(x+\gamma)^2+8\zeta(3)(x+\gamma)+3[\zeta^2(2)-2\zeta(4)]$$ $$p_5=p_1^5-10\zeta(2)p_1^3+20\zeta(3)p_1^2+15[\zeta^2(2)-2\zeta(4)]p_1+4[-5\zeta(2)\zeta(3)+6\zeta(5)]$$
Can anyone provide a reference for these polynomials or point out an interesting combinatorial interpretation?
Background: Rich associations with fractional calculus, digamma function, ladder operators
The fractional integro-derivative can be represented as an exponentiated convolutional infinitesimal generator (cf. MSE-Q125343):
$\displaystyle\frac{d^{-\beta}}{dx^{-\beta}}\frac{x^{\alpha}}{\alpha!}= \displaystyle\frac{x^{\alpha+\beta}}{(\alpha+\beta)!} = exp(-\beta R_x) \frac{x^{\alpha}}{\alpha!}$
where
$$R_xf(x)=\frac{1}{2\pi i}\displaystyle\oint_{|z-x|=|x|}\frac{-ln(z-x)+\lambda}{z-x}f(z)dz$$
$$=(-ln(x)+\lambda)f(x)+\displaystyle\int_{0}^{x}\frac{f\left ( x\right )-f(u)}{x-u}du.$$
with $\lambda=d\beta!/d\beta|_{\beta=0}$. (Note the integrand is related to the q (Jackson) derivative, and the Pincherle derivative / commutator is $[R_x,x]=D_x^{-1}$.)
Then $$exp(-\beta R_x) 1 =\displaystyle\frac{x^\beta}{\beta!} = exp(-\beta\psi_{.}(x)), $$
with $(\psi_{.}(x))^n=\psi_n(x)$, which implies
$$\psi_{n}(x)=(-1)^n \frac{d^n}{d\beta^n}\frac{x^\beta}{\beta!}|_{\beta=0},$$ $$L_x\psi_{n}(x)=n\psi_{n-1}(x)=-x\displaystyle\frac{d}{dx}\psi_{n}(x),$$ $$R_x\psi_{n}(x)=\psi_{n+1}(x).$$
Let $x=e^z$ and $p_n(z)=(-1)^n \psi_{n}(e^z)$. Then
$$exp(-\beta R_z) 1 =\displaystyle\frac{exp(\beta z)}{\beta!} = exp(\beta p_{.}(z)), $$
$$L_z p_{n}(z)=n p_{n-1}(z)=\displaystyle\frac{d}{dz} p_{n}(z),$$ $$R_z p_{n}(z)= p_{n+1}(z)= (z+\gamma)p_n(z)-\displaystyle\int_{-\infty}^{z}\frac{p_n\left ( z\right )-p_n(u)}{e^z-e^u} e^u du$$
with $\gamma=-d\beta!/d\beta|_{\beta=0}$, the Euler-Mascheroni constant.
Since $p_n(z)$ is an Appell sequence and, consequently, $p_n(x+y)=(p.(x)+y)^n$, umbrally, a change of integration variables $\omega=z-u$ gives
$$R_z p_{n}(z)= p_{n+1}(z)= (z+\gamma)p_n(z)-\displaystyle\int_{0}^{\infty}[p_n(z)-(p_{.}(z)-\omega)^n] \frac{1}{e^{\omega}-1}d\omega$$
from which the recursion formula follows.
In addition, using the operator formalism for Sheffer sequences, of which the Appell is a special case,
$$R_z=z-\frac{\mathrm{d} }{\mathrm{d} \beta}ln[\beta!]\mid _{\beta=\frac{\mathrm{d} }{\mathrm{d} z}=D_z}=z-\Psi(1+D_z)$$
where $\Psi(x)$ is the digamma or Psi function.
UPDATE (Nov. 16, 2012): Just found this exact sequence in the thesis "Regularized Equivariant Euler Classes and Gamma Functions" by R. Lu with a discussion of the relationships to Chern and Pontrjagin classes.