Riemann zeta function at positive integers and an Appell sequence of polynomials related to fractional calculus 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.
 A: Let $P_i$ be the power sum symmetric function. In your $p_n$, Replace $x+\gamma$ by $P_1$ and $\zeta(i)$ by $P_i$. Then divide the result by $n!$. What you get looks like a well-known symmetric function, which corresponds to the sign representation of the symmetric group $S_n$.
A: Follow-up on Rupinski's and Chapoton's observations:
To nail down the identification of the $p_n(x)$ with the cycle index polynomials for $S_n$ (or the partition polynomials of the refined Stirling numbers of the first kind A036039), look at the Taylor series rep of the digamma operator for the raising / creation operator for the $p_n(z)$ basis 
$$R_z = z-\Psi(1+D_z) = z+\gamma + \sum_{n=1}^{\infty } (-1)^n\zeta (n+1)D_z^n.$$ 
This is precisely the raising operator for the cycle index polynomials as presented on page 23 of Lagrange à la Lah Part I with $c_1=z+\gamma=p_1(x)$ and $c_{n+1}=(-1)^n\zeta(n+1)$ for $n>0$
$$D^{-1}_{c_1}= :\frac{c_{.}}{1-c_{.}D_{c_1}}: = c_1+\sum_{n=1}^{\infty } c_{n+1}D_{c_1}^n.$$
Alternatively, the Newton identities extrapolated to an entire function as an infinite order polynomial using the Weierstrass factorization maneuver can be applied to see the connections to the power and elementary symmetric polynomial formalism: 
$$\exp\left (-\beta p_{.}(z)\right )=\frac{\exp\left (-\beta z  \right )}{\left (-\beta  \right )!}=\exp\left (-\beta(z+\gamma)  \right )\prod_{k=1}^{\infty }\left ( 1-\frac{\beta}{k} \right )\exp\left (\frac{\beta}{k}  \right )$$
$$=\exp\left [-(z+\gamma)\beta -\sum_{k=2}^{\infty } \frac{\zeta (k)\beta ^k}{k}  \right ]=\exp\left [ :ln(1-a\beta ) :\right ]$$ where $a^1=a_{1}=(z+\gamma)$ and $a^k=a_k=\zeta(k)$ for $k>1$, but this is precisely the umbral form of the e.g.f. for the cycle index polynomials (mod signs).
(Also there are connections to rational zeta series.)
Update (Nov. 16, 2012): The generating series appears on pg. 58 in "Hodge theoretic aspects of mirror symmetry" by L. Katzarkov, M. Kontsevich, and T. Pantev (following Lu's references).
