Bernoulli umbra is some object $B$, an element of a commutative ring, such that there is an “index lowering” linear operator $\operatorname{eval}$ which applied to $B^n$ will give $B_n$, the $n$-th Bernoulli number. At first sight, Bernoulli umbra can be represented as formal power series:
$B= \varepsilon^{-1} -\frac{1}{2}-\frac{\varepsilon }{24}+\frac{3 \varepsilon ^3}{640}-\frac{1525 \varepsilon ^5}{580608}+\dotsb$, with lowering index operator corresponding to taking the coefficient of $1=\varepsilon^0$ of the power series. The numerators of the terms are given in OEIS A118050 and the denominators are in OEIS A118051.
This indeed will work well with positive integer powers of $B$.
But if we generalize the Bernoulli numbers with (Hurwitz) Zeta function, this definition will fail.
For instance, the expected value of $\operatorname{eval} (B+1)^{-1}$ is $\pi^2/6$ and in general, $\operatorname{eval}f(B+x)=\frac{D}{e^D-1}[f](x)=D\Delta^{-1}[f](x)$. As an example, it it desired that $\operatorname{eval}\ln(B+x)=\psi(x)$ (digamma function).
So, my question is: is there an accurate representation of Bernoulli umbra via (infinite) matrices such that the generalization of Bernoulli numbers via Zeta function will be accurately reflected?
