I've revisited an old post of mine--Dirac's Delta Functions and Riemann's Jump Function J(x) for the Primes--dealing with Riemann's "jump" or "staircase" function (aka, Π(x)) that has unit steps for each prime along the horizontal axis and smaller steps of size $1/n$ for an $n$-th power of a prime. This function is derived as an integral of the inverse Mellin transform of $\log(\zeta(1-s))$. The inverse Mellin transform can also be realized as a differential operator acting on a delta function:

$$\frac{J(x)}{dx} =\log[\zeta(1+xD_x)] \delta(x-1)=\sum_p \sum_{n>0}\frac{1}{n} \delta(x-p^n),$$

where the sum is over the primes $p$ and $D_x=d/dx.$

Another instance in which values of the Riemann zeta appear in a differential operator is presented in the MO-Q "Riemann zeta function at positive integers and an Appell sequence of polynomials related to fractional calculus"

$$R_x = -\log(x)-\Psi(1+xD_x) = -\log(x) +\gamma + \sum_{n=1}^{\infty } (-1)^n\zeta (n+1)(xD_x)^n,$$

where $\gamma$ is the Euler-Mascheroni constant and $\Psi$ the digamma function. $R_x$ is an infinitesimal generator for a class of fractional calculus operators; i.e.,

$$ e^{\beta R_x} \frac{x^{\alpha}}{\alpha!}= D^\beta \frac{x^{\alpha}}{\alpha!}= \frac{x^{\alpha-\beta}}{(\alpha-\beta)!}.$$

In addition, with a the change of variable $x=e^z$, it becomes the raising operator for an Appell sequence of polynomials $p_n(z)$ related to gamma classes as shown in the MO-Q, which also has the generator

$$\frac{1}{D_z!}z^n =\exp\left [-\gamma D_z -\sum_{k=2}^{\infty } \frac{\zeta (k)D_z^k}{k} \right ]z^n = p_n(z).$$

Another example of the occurrence of the zeta function, disguised as the Bernoulli numbers $b_n$ is a generator of the Bernoulli polynomials $B_n(x)$, which can be related to differentiation: umbrally,

$$T \; x^n = \frac{D_x}{e^{D_x} -1} x^n = \exp(b.D_x)x^n= (b.+x)^n = B_n(x).$$

[Edit Jan 11, 2021: This is the Todd operator, whose relation to combinatorics, volumes of polytopes, summation/trace formulas, and the calculus in general is indicated with references in answers and comments to this MO-Q and this one. Also used by Hirzebruch to construct the Todd (characteristic) class.]

Values for the zeta function at the negative integers, the Bernoullis again, also pop up in raising operators $R$ for several Appell Sheffer polynomial sequences $p_n(x)$ defined by

$$R p_n(x) = p_{n+1}(x).$$

They also satisfy an evolution equation

$$\frac{d}{dt}\exp[tp.(x)]= R \exp[tp.(x)],$$

where $$\exp[tp.(x)]=e^{tp.(0)}e^{xt}$$

is the umbral representation of the e.g.f. of the Appell sequence. Umbrally, $(p.(x))^n=p_n(x).$

1. The raising op for the reversed face polynomials of the simplices normalized by an integer,

$$p_n(x)=\frac{(x+1)^{n+1}-x^{n+1}}{n+1}$$

(OEIS A074909, a truncated Pascal triangle),

$R = x - D_t \; \frac{t}{e^t-1} \; |_{t=D_x} = x - \exp[\omega. D_x] $

[1/16/21, corrected notation with $(\omega.)^n = \omega_n =-\tfrac{b_{n+1}}{n+1} = \zeta(-n) $]

2. The raising op for the Bernoulli polynomials

$R = x + D_t \; \frac{t}{e^t-1} \; |_{t=D_x} = x + \exp[\omega. D_x].$

(These polynomials and the ones above are an inverse pair under umbral composition, which means their lower triangular coefficient matrices are also an inverse pair. This property holds for Appell sequences whose raising ops differ only by the one sign as above.)

3. The raising op for the integral, normalized Euler polynomials A081733

$$R = x - \frac{2}{e^{-2D}+1}$$

with $$\frac{2}{e^{2t}+1} = 2 \sum_{n g.e. 0} \eta(-n) (-2t)^n/n!,$$

where $\eta(s)$ is the Dirichlet eta function, and

$$2(-2)^n \eta(-n) = (-1)^n [2^{n+1}-4^{n+1}] \zeta(-n) = [2^{n+1}-4^{n+1}] \tfrac{b_{n+1}}{n+1}.$$

4. Its umbral compositional inverse is A119468 with

$$R = x + \frac{2}{e^{-2D}+1}.$$

This generates a set of polynomials which when multiplied by 2 gives essentially the face polynomials of the hypercubes A038207, which is the square of the Pascal triangle.

5. The raising op for the Swiss-knife polynomials A119879, from which the Bernoulli, Gennochi, Euler, tangent, and Springer numbers can be computed.

$$R = x + \tanh(D).$$

$\tanh(x)$ is the e.g.f. of the zag numbers A000182,

$$zag(n) = 2^{2n} (2^{2n} - 1) \frac{|b_{2n}|}{2n}=2^{2n} (2^{2} - 1)|\zeta(-2n+1)|.$$

Note that the matrix of coefficients is a signed, masked Pascal triangle.

6. Its umbral inverse is A119467, the same masked Pascal triangle unsigned, with

$$R = x - \tanh(D).$$

7. (added 1/11/21) Dragovich presents operators of the form

$$\zeta(a\square+b),$$

where $\square= \partial^2_t + \triangledown^2$ is the D-dimensional d'Alembert operator, in "From p-Adic to Zeta Strings" and "p-Adic Mathematical Physics: the First 50 Years" along with Khrenikov, Kozyrev, Volovich, and Zelenov.

Question: In what other differential operators does the Riemann zeta function play an important role?