Babusci and Dattoli, On the logarithm of the derivative operator, arXiv:1105.5978, gives some great results: \begin{align*} (\ln D) 1 & {}= -\ln x -\gamma \\ (\ln D) x^n & {}= x^n (\psi (n+1)-\ln x) \\ (\ln D) \ln x & {}= -\zeta(2) -(\gamma+\ln x)\ln x. \end{align*} I wonder, what is its matrix, or otherwise, is there a method of applying it to a function?
What is its intuitive role in various fields of math?
[Edit Oct. 13, 2022: Corrected $1$ to $z^n$ in item 3).]
The interpretation of a $\ln(D)$ depends on the interpolation that one chooses of the usual derivative operator and its positive integer powers to a fractional integro-derivative operator (FID), i.e., an interpretation of $D$ exponentiated by any real (or complex number via analytic continuation), which in turn, depends on the functions the FID is to act upon. The extension described below produces B & Ds three identities and is consistent with the properties that Pincherle imposed on any legitimate family of FIDs (see this MO-Q on a 1/2 derivative and this MO-Q on fractional calculus). It can be defined by the action on a 'basis set' of entire functions in the complex variable $\omega$ as
$$D_x^\alpha \; H(x) \; \frac{x^\omega}{\omega!} = H(x) \frac{x^{\omega-\alpha}}{(\omega-\alpha)!} ,$$
where $H(x)$ is the Heaviside step function, and $\alpha$ and $\omega$ may be any complex numbers with the usual identification in the theory of generalized functions and distributions of
$$(-1)^n \delta^{(n)}(x) = H(x) \frac{x^{-n-1}}{(-n-1)!},$$
with $n=0,1,2,3,\ldots$.
Note this has little to do with a Fourier transform over the real line or any pseudo-diff op/symbol associated with such. In particular, $D^\alpha$ here is NOT associated with multiplication by $(i 2 \pi f)^\alpha$ in frequency space. Elsewhere I show various equivalent convolutional reps of this FID as 1) a FT over a circle via a transformation of a regularized Cauchy complex contour integral, 2) the analytic continuation of the integral rep of the Euler beta function either through a blow-up into the complex plane of the integral along the real line segment or regularization via the Hadamard finite part or via the Pochhammer contour, 3) the Mellin interpolation of the standard derivative operator via the action of the generating function $e^{tD_x}$, an operator application of Ramanujan's master formula, or 4) a sinc function/cardinal series interpolation of the generalized binomial coefficients.
Let's see how viable the above definition of the FID is; its connection to an infinitesimal generator (infinigen) of the FID and the three B & D identities; a connection to the formalism of Appell Sheffer polynomial sequences and, therefore, symmetric polynomial/function theory; and matrix reps of the infinigen and FID.
If we assume that an infinitesimal generator $IG$ exists such that
$$ e^{\alpha \; IG} \; H(x) \; \frac{x^\omega}{\omega!} = D_x^{\alpha} \; H(x) \; \frac{x^\omega}{\omega!} = H(x) \frac{x^{\omega-\alpha}}{(\omega-\alpha)!} = e^{-\alpha D_\omega} \; H(x) \; \frac{x^\omega}{\omega!},$$
then formally
$$D_\alpha \; e^{\alpha IG} \; H(x) \; \left. \frac{x^\omega}{\omega!} \right|_{\alpha =0} = IG \; H(x) \; \frac{x^\omega}{\omega!} = \ln(D_x) \; H(x) \; \frac{x^\omega}{\omega!}$$
$$ = D_{\alpha} \; H(x) \; \left. \frac{x^{\omega-\alpha}}{(\omega-\alpha)!} \right|_{\alpha =0} = -D_{\omega} \;\frac{x^{\omega}}{\omega!}$$
$$ = [\; -\ln(x) + \psi(1+\omega) \;] H(x) \; \frac{x^{\omega}}{\omega!} $$
$$ = [ \; -\ln(x) + \psi(1+xD_x) \;] \; H(x) \; \frac{x^\omega}{\omega!}, $$
and the infinigen is
$$ \ln(D_x) := IG = -\ln(x) + \psi(1+xD_x),$$
where $\psi(x)$ is the digamma function, which can be defined over the complex plane as a meromorphic function and is intimately related to the values of the Riemann zeta function at $s = 2,3,4,\ldots$.
Some reps (that give the same identities as in B & D) are
$$IG \; f(x)=\frac{1}{2\pi i} \oint_{|z-x|=|x|}\frac{-\ln(z-x)+\lambda}{z-x}f(z) \; dz$$
$$=(-\ln(x)+\lambda) \; f(x)+ \int_0^x \frac{f(x)-f(u)}{x-u}\, du$$
$$ = \left[\; -\ln(x)+ \left. \frac{\mathrm{d} }{\mathrm{d} \beta} \ln[\beta!]\right| _{\beta =xD} \; \right] \; f(x)= \left[ \; -\ln(x)+\Psi(1+xD) \;\right] \; f(x)$$
$$ = \left[ \; -\ln(x)+\lambda - \sum_{n=1}^\infty (-1)^n\zeta (n+1) \; (xD)^n \;\right] \; f(x)$$
where $\lambda$ is related to the Euler–Mascheroni constant via $\lambda=D_\beta \; \beta! \;|_{\beta=0}$.
Other reps and other ways of arriving at the reps above are given in the refs below.
Let's look at a way via the formalism of Appell Sheffer polynomial sequences, which settles any issues of convergence upon exponentiation of the explicit diff op formula for the infinigen and allows connections to the theory of symmetric polynomials/functions.
The relevant Appell sequence of polynomials $p_n(z) = (p.(z))^n$ has the exponential generating function, entire in the complex variable $t$, i.e., with its Taylor series globally convergent,
$$\frac{1}{t!} \; e^{zt} = e^{a.t} \; e^{zt} = e^{(a.+z)t} = e^{p.(z)t} = \sum_{n\geq 0} p_n(z) \frac{t^n}{n!}$$
with the reciprocal polynomial sequence defined in four consistent ways $\hat{p}(z)$
1) $t! \;e^{zt} = e^{\hat{a}.t} \; e^{zt} = e^{(\hat{a}.+z)t} = e^{\hat{p}.(z)t} $, an e.g.f.,
2) $M_p \cdot M_{\hat{p}} = I $, in terms of the lower triangular coefficient matrices of the two sequences in the monomial power basis $z^n$ with unit diagonal,
3) $p_n(\hat{p}.(z)) = \hat{p}_n(p.(z)) = (a. + \hat{a.}+z)^n = z^n$, an umbral convolutional inversion,
4) $D_z! \; z^n = e^{\hat{a.}D_z} \; z^n = (\hat{a.}+z)^n = \hat{p}_n(z)$, an operational generator.
It follows that the raising op of the Appell polynomials $p_n(z)$ defined by
$$R_z \; p_n(z) = p_{n+1}(z)$$
is given by
$$ R_z \; p_n(z) = \frac{1}{D_z!} \; z \; D_z! \; p_n(z) = \frac{1}{D_z!} \; z \; p_n(\hat{p}.(z))$$
$$ = \frac{1}{D_z!} \; z \; z^n = \frac{1}{D_z!} \; z^{n+1} = p_{n+1}(z),$$
an operator conjugation, or 'gauge transformation', of the raising operator $z$ for the power monomials.
In addition, with the operator commutator $[A,B] = AB - BA$,
$$R_z = \frac{1}{D_z!} \; z \; D_z! = z + \left[\frac{1}{D_z!},z \right] \; D_z! .$$
Now re-enter Pincherle and the eponymous operator derivative, which Rota touted for the finite operator calculus. The Graves–Pincherle derivative derives its power from the Graves-Lie-Heisenberg-Weyl commutator $[D_z,z] = 1$ from which, by normal re-ordering, implies for any function expressed as a power series in $D_z$
$$[f(D_z),z] = f'(D_z) = D_t \; f(t) \; \Big|_{t = D_z}.$$
This is an avatar of the Pincherle derivative (PD) that follows from the action $$[D^n,z] \; \frac{z^\omega}{\omega!} = \left[\;\frac{\omega+1}{(\omega+1-n)!} - \frac{1}{(\omega-n)!}\;\right] \; z^{\omega+1-n} = n \; D_z^{n-1} \; \frac{z^\omega}{\omega!},$$
but the PD is valid for more general lowering and raising (ladder) ops that satisfy $[L,R]= 1$.
Then
$$R_z = \frac{1}{D_z!} \; z \; D_z! = z + \left[\frac{1}{D_z!},z \right] \; D_z! = z + D_{t = D_z}\; \ln\left[\frac{1}{t!}\right] $$
$$ = z - \psi(1+D_z).$$
With the substitution $ z = \ln(x)$
$$R_z = R_x = \ln(x) - \psi(1+ x D_x) = -IG = -\ln(D_x).$$
The raising op is defined such that
$$ e^{t \; R_z} \; 1 = \sum_{n \geq 0} \frac{t^n}{n!} R_z^n \; 1 = e^{tp.(z)} = \frac{1}{t!} \; e^{zt},$$
an entire function for $t$ complex; therefore,
$$e^{-t \; IG} \;1 = e^{t \;R_x} \; 1 = e^{t \; p.(\ln(x))} = \frac{x^t}{t!},$$
so
$$e^{-(\alpha+\beta) \; IG} \;1 = e^{(\alpha+\beta) \; R_x} \; 1 = e^{(\alpha+\beta) \; p.(\ln(x))} = \frac{x^{\alpha+\beta}}{(\alpha+\beta)!}, $$
$$ = e^{-\alpha \; IG} e^{-\beta \; IG} \;1 = e^{-\alpha \; IG} \; \frac{x^\beta}{\beta!} , $$
and we can identify that indeed
$$e^{-\alpha \; IG} = D_x^{-\alpha}$$
and
$$IG = \ln(D_x).$$
Now apply the PD to $\ln(D)$, as a check of the formalism and an avenue to a matrix rep, giving formally
$$ [\ln(D),x] = [\ln(1-(1-D)),x] = \frac{1}{1-(1-D)} = \frac{1}{D} = D^{-1}.$$
This is given an explicit meaning by evaluating the commutator for a general function $g(x)$ analytic at the origin (which generalizes to our 'basis' set) using the integral rep for $R_x = -\ln(D_x)$, giving
$$[\ln(D_x),x] \; g(x) = [-R_x,x] \; g(x) = (-\ln(x)+\lambda) \; [x,g(x)]$$
$$ + \int_{0}^{x}\frac{xg(x)-ug(u)}{x-u} \; du - x \int_{0}^{x}\frac{g(x)-g(u)}{x-u} \; du$$
$$ = \int_{0}^{x} \; g(u) \; du = D_x^{-1} g(x).$$
So, we have
$$[\ln(D_x),x] = [-R_x,x] = D_x^{-1} = [-\ln([-R_x,x]),x]$$
and
$$-R_x = \ln(D_x) = -\ln(D_x^{-1}) = -\ln([-R_x,x]),$$
implying
$$e^{R_x} =\exp[\ln([-R_x,x])] = [-R_x,x] = D_x^{-1}.$$
In addition, with
$$\bigtriangledown^{s}_{n} \; c_n=\sum_{n=0}^{\infty}(-1)^n \binom{s}{n}c_n,$$
then
$$R_x = -\ln(D_x) = \ln(D_x^{-1}) = \ln[1-(1-D_x^{-1})]$$
$$ = - \sum_{n \geq 1} \frac{1}{n} \; \bigtriangledown^{n}_{k} D_x^{-k}, $$
where
$$D_x^{-1} \frac{x^\omega}{\omega!} = \frac{x^{\omega+1}}{(\omega+1)!}.$$
The finite difference op series is embedded in the derivative $D_{\alpha =0}$ of the Newton interpolator
$$ \frac{x^{\alpha+\omega}}{(\alpha+\omega)!} = \bigtriangledown^{\alpha}_n \bigtriangledown^n_k \frac{x^{\omega+k}}{(\omega+k)!}$$
$$ = \bigtriangledown^\alpha_n \bigtriangledown^n_k D_x^{-k} \;\frac{x^\omega}{\omega!}$$
$$ = [1-(1-D_x^{-1})]^\alpha \; \;\frac{x^\omega}{\omega!} = D_x^{-\alpha}\;\frac{x^\omega}{\omega!}. $$
For $\alpha = -m$ with $m = 1,2,\ldots$ and $\omega = 0$, this Newton interpolator gives
$$D^m_x \; H(x) = \delta^{(m-1)}(x) = H(x) \; \frac{x^{-m}}{(-m)!} = \bigtriangledown^{-m}_{n}\bigtriangledown^{n}_{k} D_x^{-k} \; H(x)$$
$$ = \sum_{n \geq 0} (-1)^n \binom{-m}{n} \bigtriangledown^{n}_{k} \; H(x) \frac{x^k}{k!} = H(x) \; \sum_{n \geq 0} (-1)^n \binom{-m}{n} \; L_n(x)$$
$$ = H(x) \; \sum_{n \geq 0} \binom{m-1+n}{n} \; L_n(x), $$
which agrees in a distributional sense with the Laguerre polynomial resolutions of $f(x) = \delta^{(m-1)}(x)$ in the formulas of this MO-Q since, with $c_n = f_n$ in the notation there,
$$ f(x) = \sum_{n \geq 0} c_n \; L_n(x)$$
with
$$\sum_{n \geq 0} t^n \; c_n = \frac{1}{1-c.t} = \int_0^{\infty} e^{-x} \sum_{n \geq 0} t^n \; L_n(x) f(x) \; dx$$
$$ = \int_0^{\infty} e^{-x} \frac{e^{-\frac{t}{1-t}x}}{1-t} f(x) \; dx = \int_0^{\infty} \frac{e^{-\frac{1}{1-t}x}}{1-t} f(x) \; dx,$$
so, for the $m$-th derivative of the Heaviside function,
$$\frac{1}{1-c_{m,.}t}= \int_0^{\infty} e^{-x} \frac{e^{-\frac{t}{1-t}x}}{1-t} f(x) \; dx = \int_0^{\infty} \frac{e^{-\frac{1}{1-t}x}}{1-t} \delta^{(m-1)}(x) \; dx = \frac{1}{(1-t)^{m}},$$
and, therefore, the coefficients of the Laguerre series resolution of the $m$-th derivative of the Heaviside function are
$$c_{m,n} =(-1)^n \binom{-m}{n} = \binom{m-1+n}{n},$$
in agreement with the Newton interpolator.
Applying $D_x^{-1}$ iteratively to both sides of this identity establishes convergent interpolations for $\omega = 1,2,3,...$, and acting on the power basis within the binomial expansion of $\frac{x^{\omega}}{\omega!} = \frac{(1-(1-x))^{\omega}}{\omega!}$ should give convergent expressions as well.
Similarly for $\omega=0$, we have the Laplace transform (or more accurately, the modified Mellin transform central to Ramanujan's master formula via which the FIDs may be cast as Mellin interpolations of the standard derivatives),
$$\frac{1}{1-c.t} = \int_0^\infty \frac{e^{-\frac{1}{1-t}x}}{1-t} \frac{x^\alpha}{\alpha!} \; dx = (1-t)^\alpha,$$
for $\operatorname{Re}(\alpha) > -1$, giving
$$c_n = (-1)^n \binom{\alpha}{n}.$$
This Laplace transform and, therefore, the Newton interpolator can be analytically continued in several standard ways (e.g., blow-up from the real line to the complex plane via a Hankel contour, Hadamard finite part) to the full complex plane for $\alpha$. For the negative integer exponents, the Hankel contour contracts to the usual Cauchy contour rep for differentiation. The Hadamard-finite-part approach allows the Newton interpolator to be appropriately modified strip by strip to give the intended results.
Returning to the finite difference rep for $\ln(D_x)$, action of the infinigen on 1 then gives, for $x > 0$,
$$\ln(D_x) 1 = \sum_{n \geq 1} \frac{1}{n} \; \bigtriangledown^{n}_{k} D_x^{-k} 1$$
$$ = \sum_{n \geq 1} \frac{1}{n} \; \bigtriangledown^{n}_{k} \frac{x^k}{k!}$$
$$ = \sum_{n \geq 1} \frac{1}{n} \; L_n(x) = -\ln(x)-.57721... , $$
where $L_n(x)$ are the Laguerre polynomials, in agreement with the first equation of B & D in the question.
Plots of the results of evaluation of the operator series truncated at $n=80$, or so, acting on $x^2$ and $x^3$ match the analytic results as well.
The matrix rep $M$ of the action of this integration op $D_x^{-1}$ on $x^n$ is simple enough in the power basis--a matrix with all zeros except for the first subdiagonal, or superdiagonal, depending on left or right matrix multiplication, with elements $(1,1/2,1/3,...)$.
The matrix rep for $R_x$ is then
$$ R_M = \ln[I-(I-M)] = - \sum_{n \geq 1} \frac{1}{n} \; \bigtriangledown^{n}_{k} M^k. $$
Exponentiating,
$$D_x^{-\beta} = \exp(-\beta R_x)= (1-(1-D_x^{-1} ) )^{\beta} = \bigtriangledown^{\beta}_{n} \bigtriangledown^{n}_{k} (D_x^{-1})^k.$$
The associated matrix rep is
$$ \exp(-\beta R_M)= \bigtriangledown^{\beta}_{n} \bigtriangledown^{n}_{k} M^k.$$
(I haven't checked these matrix computations numerically as I normally would since my MathCad disc is in storage in another state.)
To act on non-integer powers of $x$, you must represent them as superpositions of the integer power basis as in the binomial expansion
$$x^{\alpha} = [1 - (1-x)]^{\alpha} = \bigtriangledown^{\alpha}_{n} \bigtriangledown^{n}_{k} x^k .$$
Alternatively, return to the $z$ rep and write down the matrix rep of the raising op $R_z$. This is a simple transformation of the infinite lower triangular Pascal matrix augmented with a first superdiagonal of all ones. OEIS A039683 has an example of the matrix equivalent of a raising op in the monomial power basis, also known as a production matrix in another approach (Riordan?) to polynomial sequences. Better in this case to switch to the divided power basis $z^n/n!$. Then the augmented Pascal matrix becomes the simple summation matrix of all ones. Multiply along the n-th diagonal by $c_n$ where $(c_0,c_1,..) = (1-\lambda,-\zeta(2),...,(-1)^k \; \zeta(k+1),...)$ to generate the matrix rep for the raising op, but since, e.g., $x^2=e^{2z}$, this quickly becomes a messy algorithm to apply compared to the finite difference rep.
Further references (not exhaustive):
Upon Fourier transformation $x\mapsto k$, this becomes a diagonal operator with matrix elements $\langle k|\ln D|k'\rangle=2\pi \delta(k-k')\ln k$. So to find the matrix elements in the $x$-representation we would need to invert the Fourier transform of the logarithm $\ln k$. From this MSE answer for the Fourier transform of $\ln |k|$ (with absolute value signs) I would conclude that $$\langle x|\ln D|x'\rangle=\left(\frac{i \pi}{2}-\gamma\right) \delta (x-x')+\text{P.V.}\left(\frac{1}{2 (x-x')}-\frac{1}{2 | x-x'| }\right).$$
This notation means that $\ln D$ acting on a function $f(x)$ produces a new function $g(x)$ given by $$g(x)=\int_{-\infty}^\infty \left[\left(\frac{i \pi}{2}-\gamma\right) \delta (x-x')+\text{P.V.}\left(\frac{1}{2 (x-x')}-\frac{1}{2 | x-x'| }\right)\right]f(x')\,dx'$$ $$=\left(\frac{i \pi}{2}-\gamma\right) f(x)+\frac{1}{2}\,\text{P.V.}\int_{-\infty}^\infty \left(\frac{1}{x-x'}-\frac{1}{| x-x'| }\right)\,f(x')\,dx'.$$
