Boson operator ordering identities from generalized Stirling and Eulerian numbers
Robert S. Maier

Depts. of Mathematics and Physics, University of Arizona, Tucson,
AZ 85721, USA

Ordering identities in the Weyl–Heisenberg algebra generated by
single-mode boson operators are investigated.  A boson string composed
of creation and annihilation operators can be expanded as a linear
combination of other such strings, the simplest example being a normal
ordering.  The case when each string contains only one annihilation
operator is already combinatorially nontrivial.  Two kinds of
expansion are derived: (i) that of a power of a string $\Omega$ in
lower powers of another string $\Omega'$, and (ii) that of a power
of $\Omega$ in twisted versions of the same power of $\Omega'$.  The
expansion coefficients are shown to be, respectively, generalized
Stirling numbers of Hsu and Shiue, and certain generalized Eulerian
numbers.  Many examples are given.  These combinatorial numbers are
binomial transforms of each other, and their theory is developed,
emphasizing schemes for computing them: summation formulas,
Graham–Knuth–Patashnik (GKP) triangular recurrences, terminating
hypergeometric series, and closed-form expressions.  The results on
the first type of expansion subsume a number of previous results on
the normal ordering of boson strings.

Boson operator Weyl–Heisenberg algebra ordering identity generalized Stirling numbers generalized Eulerian numbers triangular recurrence binomial transform hypergeometric series
[2020] 11B73 11B37 05A10 05A15 81R05

§ INTRODUCTION

Single-mode boson creation and annihilation operators $a^\dag,a$ do
not commute, but satisfy the commutation relation $[a,a^\dag]=1$.
Over the complex numbers, they generate the Weyl–Heisenberg
algebra $W_{\mathbb{C}}$.  Similarly noncommuting pairs of operators
appear in many places in quantum mechanics and quantum field theory,
and many examples could be adduced [Fleury and Turbiner, 1994].  For instance, the
position operator $x$ and the operator $D=D_x={\rm d}/{\rm d}x$
satisfy the equivalent relation $[D,x]=1$.  They too generate a
$\mathbb{C}$-algebra, but negative integer powers of $x$
and indeed arbitrary powers of $x$ will be allowed here, so the
algebra is a larger one, $\tilde W_{\mathbb{C}}$.  The algebra
$W_{\mathbb{C}}$ has a representation on Fock space, and $\tilde
W_{\mathbb{C}}$ can be represented on any convenient space of
functions of $x$, but this paper will deal with their algebraic
properties rather than their representations.

Reordering the operators in any algebra element built from $a^\dag,a$
(or $x,D$), to satisfy constraints, may be combinatorially nontrivial.
A standard problem is that of normal ordering: moving all
operators $a$ (or $D$) to the right, to facilitate the computing of
quantum expectation values.  It has long been known [Katriel, 1974]
that the normal ordering of $(a^\dag a)^n$ and $(a \adag)^n$ involves
Stirling numbers of the second kind [Graham et al., 1994], which are denoted
by $S(n,k)$ or $\stirsub{n}k$.  In present notation,

\begin{equation}
\label{eq:katriel}
(\adag a)^n = \sum_{k=0}^n S_{n,k}(0,1;0)\, \adag^k (a)^k,
\qquad
(a \adag)^n = \sum_{k=0}^n S_{n,k}(0,1;1)\, \adag^k (a)^k,
\end{equation}

where the coefficients $S_{n,k}(0,1;0)$ and $S_{n,k}(0,1;1)$
respectively equal $\stirsub{n}k$ and $\stirsub{n+1}{k+1}$.  These are
ordering identities, and they have been greatly generalized
[Mansour and Schork, 2016, Schork, 2021].  For example, if the word raised to the
$n$'th power is an arbitrary word $\Omega=\Omega(a^\dag,a)$ in the
noncommuting $a^\dag,a$, or a linear combination, it may be possible
to compute the coefficients in the normal ordering of $\Omega^n$
iteratively or recursively, or view them as counting combinatorial
structures [Błasiak and Flajolet, 2012].

This paper derives ordering identities in the Weyl–Heisenberg algebra
of two new kinds, not previously considered in detail.  The derivation
is largely algebraic rather than combinatorial.  In the first, one
expands the $n$'th power of a word $\Omega=\Omega(a^\dag, a)$ in
powers $\Omega{}^\prime{}^k$, $k=0,\dots,n$, of another specified
word $\Omega'$.  For simplicity, each of $\Omega,\Omega'$ will be
required to contain exactly one annihilation operator $a$, so $\Omega
= \adag^L a\adag^R$ and $\Omega' = \adag^{L'} a\adag^{R'}$, or
correspondingly, $\Omega = x^L D x^R$ and $\Omega'= x^{L'} D x^{R'}$.
The first theorem dealing with this situation contains the following
parametric operator identity in $\tilde W_{\mathbb{C}}$ (see
Theorem <ref>):

\begin{equation}
\begin{aligned}
&\left(x^{1-L-R}\right)^n \left(x^L D x^R\right)^n\\
&\qquad {}= \sum_{k=0}^n S_{n,k}\left(1-L-R,1-L'-R';\,R-R'\right)
\left(x^{1-L'-R'}\right)^k \left(x^{L'}Dx^{R'}\right)^k.
\label{eq:firstmain}
\end{aligned}
\end{equation}

Here, $L,R,L',R'\in\mathbb{C}$ are arbitrary, and the coefficients
$S_{n,k}$, $0\le k\le n$, with three complex parameters, are
“unified” generalized Stirling numbers of Hsu and
Shiue [Hsu and Shiue, 1998].  They are one among many generalizations of the
classical Stirling numbers; for others, see
[Mansour and Schork, 2016, Sándor and Crstici, 2004].

One can formally substitute $\adag,a$ for $x,D$
in (<ref>), but for the result to be an identity in
$W_{\mathbb{C}}$, it must contain only nonnegative integer powers of
$\adag$ and $a$.  So one must require $L,R,L',R'\in\mathbb{N}$, and to
remove negative powers, left-multiply both sides
of (<ref>) by $(x^E)^n$, where $E\in\mathbb{N}$ is
sufficiently large.  With $x,D$ replaced by $a^\dag,a$, this yields a
parametric identity in $W_{\mathbb{C}}$:

\begin{equation}
\begin{aligned}
\left({a^\dag}^{1-L-R+E}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n
=& \sum_{k=0}^n S_{n,k}\left(1-L-R,1-L'-R';\,R-R'\right)\\
\left({a^\dag}^E\right)^{n-k}   \left({a^\dag}^{1-L'-R'+E}\right)^k \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k,
\label{eq:firstmainmod}
\end{aligned}
\end{equation}

where taking $E=\max\left(L+R-1,L'+R'-1,0\right)$ suffices (see
Theorem <ref>).

The normal orderings in (<ref>) are the
$(L,R;L',R')=(1,0;0,0)$ and $(0,1;0,0)$ cases
of (<ref>), with $E=0$, and several other known
identities [Kim and Kim, 2022, Lang, 2000, Mikhailov, 1985, Mohammad-Noori, 2011] are
also subsumed by (<ref>).  This includes the normal
ordering of any monomial $({a^\dag}^L a{a^\dag}^R)^n$, which comes by
taking $(L',R')=(0,0)$ and $\Omega'=\nobreak a$.  It is known that
normally ordering $\Omega^n$ when $\Omega$ is a “single annihilator”
word yields generalizations of the classical Stirling
numbers [Duchamp et al., 2004, Duchamp et al., 2010], and the case when
$\Omega^n$ is replaced by a product $\adag^{r_m}a^{s_m}\dots
\adag^{r_1}a^{s_1}$, which includes $({a^\dag}^L a {a^\dag}^R)^n$, has
been treated [Méndez et al., 2005]. But the
expansion (<ref>) and the extension to words
$\Omega'=\adag^{L'} a \adag^{R'}$ other than $a$ are new.

It is explained at length below how the Hsu–Shiue Stirling numbers
$S_{n,k}$ appearing as coefficients can be computed: by summation
formulas, triangular recurrences, or identities relating them at
different indices $n,k$ and different values of their parameters.
Closed-form expressions at certain parameter values are also derived.
A sample application of eq. (<ref>) is

\begin{equation}
\label{eq:sampleappl}
\begin{aligned}
\left( {a^\dag}^3 a\right)^n  &=  {a^\dag}^n\,\sum_{k=0}^n
{a^\dag}^{(n-k)} \left({a^\dag}^2 a\right)^k\\
&=  {a^\dag}^n \, \sum_{k=0}^n 2^{-(n-k)}k\,\binom{n}{k} \frac{(2n-k-1)!}{n!}\:  {a^\dag}^{(n-k)} \left({a^\dag}^2 a\right)^k, \qquad n\ge1,
\end{aligned}
\end{equation}

which does not normally-order $\Omega^n=({a^\dag}^3a)^n$, but rather
expands it in powers of $\Omega'={a^\dag}^2a$.  Here,
$(L,R;L',R')=(3,0;2,0)$ and $E=2$.  This formula incorporates a
closed-form expression for the numbers $S_{n,k}(-2,-1;0)$ with ${0\le
k\le n}$, which applies when $n\ge1$.  They are also called the
unsigned Bessel numbers of the first kind, and are denoted below
by $\hat b_{n,k}$.

There has been some work [Kim et al., 2023, Mansour et al., 2011, Mansour et al., 2012] on
the combinatorics of the algebra $W_{\mathbb{C}}(s)$ over $\mathbb{C}$
defined by $UV-\nobreak VU=V^s$.  When $s=0$, this
is $W_{\mathbb{C}}$, and when $s\in\mathbb{N}$ it can be identified
with a subalgebra of $W_{\mathbb{C}}$, its embedding
in $W_{\mathbb{C}}$ being given by $V=\adag$, $U=\adag^s a$.
In particular, $W_{\mathbb{C}}(1)$ is the well-known shift algebra,
and $W_{\mathbb{C}}(2)$ is the Jordan plane, also called the
meromorphic Weyl algebra <cit.>.  Many
ordering identities in $W_{\mathbb{C}}(s)$ follow from the parametric
identity (<ref>), and the appearance of Hsu–Shiue
generalized Stirling numbers in the context of $W_{\mathbb{C}}(s)$ has
already been noted [Mansour et al., 2012].

The second kind of ordering identity derived in this paper is based on
a generalization of the classical Eulerian numbers, which are
traditionally denoted by $A_{n,k}$ and have more recently been denoted
by $\eulerian{n}k = A_{n,k+1}$.  These numbers are not well
known to physicists, their combinatorial interpretation having been
discovered as late as the 1950's [Comtet, 1974, Graham et al., 1994, Riordan, 1958].
Identities in the Weyl–Heisenberg algebra $W_{\mathbb{C}}$ based on
the classical numbers are (in present notation)

\begin{align}
n!\,(\adag a)^n &= \sum_{k=0}^n E_{n,k}(0,1;\,0)\, \adag^k (a)^n \adag^{(n-k)},\label{eq:euleriank1}\\
n!\,(a \adag)^n &= \sum_{k=0}^n E_{n,k}(0,1;\,1)\, \adag^k (a)^n \adag^{(n-k)},\label{eq:euleriank2}
\end{align}

where the coefficients $E_{n,k}(0,1;0)$ and $E_{n,k}(0,1;1)$
respectively equal $A_{n,k}$ and $\eulerian{n}k=A_{n,k+1}$.

Generalized Eulerian numbers $E_{n,k}$, $0\le k\le n$, with three
complex parameters, which are of the type needed for
extending (<ref>), have been considered in
[Carlitz, 1979, Charalambides, 1982, Hsu and Shiue, 1999] and by the
author [Maier, 2023]; for other generalizations,
see [Sándor and Crstici, 2004].  These numbers appear as coefficients in
expansions that are generalized versions of (<ref>)
and (<ref>).  Given any single-annihilator word
$\Omega=\adag^L a \adag^R $, one can expand its $n$'th power in
“twistings” of the $n$'th power of another such word,
$\Omega'=\adag^{L'} a \adag^{R'}$.  More precisely, for sufficiently
large $p,q\in\mathbb{N}$, one can expand
$\left({a^\dag}^p\right)^n\Omega{}^n\left({a^\dag}^q\right)^n$ as a
combination of terms of the form $\adag^{p'} {\Omega'}{}^n
\adag^{q'}$, with each term having the same value of $p'+\nobreak q'$,
which is a certain multiple of $n$.  In each identity
in (<ref>), $p=q=0$, with $\Omega'=a$ and $p'+q'=n$.

A general expansion of this type (the first of four in
Theorem <ref>) is

\begin{equation}
\label{eq:introE2}
\begin{aligned}
&n!\, (-e')^n \left({a^\dag}^{(E_L-e)}\right)^n \left({a^\dag}^L a
{a^\dag}^R\right)^n \left(\adag^{E_R}\right)^n
\\ &\qquad{}=\sum_{k=0}^n E_{n,k}\left(-e,-e';\,R-R'\right)
\\ &\qquad\qquad{}\times\left({a^\dag}^{E_L}\right)^{n-k}
\left({a^\dag}^{({E_L}-e')}\right)^k
\left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^n
\left(\adag^{(E_R-e')}\right)^{n-k} \left(\adag^{E_R}\right)^k,
\end{aligned}
\end{equation}

where $e=L+R-1$, $e'=L'+R'-1$, with $E_L=\max(e,e',0)$ and
$E_R=\max(e',0)$.  In (<ref>), one sees that
$p=E_L-\nobreak e$, $q=E_R$, $p'+\nobreak q' = (E_L+\nobreak
E_R-\nobreak e')n$.  Similarly to (<ref>), this
expansion is obtained from a simpler one in $\tilde W_{\mathbb{C}}$
(see Theorem <ref>) by left-multiplying
by $(x^{E_L})^n$ and right-multiplying by $(x^{E_R})^n$, with
$E_L,E_R$ sufficiently large that each negative power of $x$ on either
side is removed or made positive.  The identities
(<ref>),(<ref>) come
from (<ref>) by setting $(L,R;\allowbreak
L',R')=(1,0;\allowbreak 0,0)$ and $(0,1;0,0)$ respectively, so that
(in each case) $e=\nobreak0$, $e'=\nobreak-1$, and $E_L=\nobreak

As with the $S_{n,k}$, the generalized Eulerian numbers $E_{n,k}$
appearing as coefficients can be computed by various means: summation
formulas, triangular recurrences, or identities relating them at
different indices $n,k$ and different values of their parameters.
Closed-form expressions at certain parameter values are also derived
below.  Another sample application of (<ref>) is the

\begin{equation}
\label{eq:sampleeulerian}
2^n\,(\adag)^n (\adag^2 a)^n (\adag^2)^n = \sum_{k=0}^n
\binom{n+1}{2k+1}\, (\adag^2)^k (\adag^3 a)^n (\adag^2)^{n-k},
\end{equation}

in which $\Omega^n = (\adag^2 a)^n $ and $\Omega'{}^n = (\adag^3
a)^n$.  This comes by setting $(L,R;L',R')=(2,0;3,0)$, so that $e=1$,
$e'=2$, $E_L=E_R=2$, and replacing the summation index $k$
by $n-\nobreak k$.  The binomial coefficients on the right-hand side
come from the closed-form expression $(-1)^n\, n!\binom{n+1}{2n-2k+1}$
for $E_{n,k}(-1,-2;0)$, and the only nonvanishing terms in the sum are
those with $0\le k\le\lfloor\frac{n}2\rfloor$.

This paper is heavily influenced by the early work of Toscano, who in
a series of papers published in the 1930's (mostly cited
in [Mansour and Schork, 2016]) obtained many singly parametrized differential
operator expansions, the coefficients in which satisfy triangular
recurrences.  Among his tools was a generalized Boole identity
(Lemma <ref>), which will play a role below.  An entry point
to his work is provided by [Toscano, 1939].  Other than being
$1$-indexed rather than $0$-indexed, most of
his expansion coefficients are one-parameter specializations of the
three-parameter $S_{n,k}$ and $E_{n,k}$ used here.  For instance, his
$B_{n+1,k+1}^{(u)}$ equals $E_{n,k}(1-\nobreak u, 1; u)$.

The plan of the paper is as follows.  Section <ref> explains
notation and mathematical conventions.  Section <ref>
contains some basic facts about the algebras $W_{\mathbb{C}}$
and $\tilde W_{\mathbb{C}}$, including the generalized Boole identity.
In Section <ref> the parametric numbers $S_{n,k}$
and $E_{n,k}$ are introduced and investigated in depth, with an
emphasis on recurrences, summation formulas, and other computationally
useful results, such as the fact that they are binomial transforms of
each other.  Section <ref> derives some
closed-form expressions for the $S_{n,k}$ and comments on
combinatorial interpretations.  Ordering identities of the generalized
Stirling kind are derived in Section <ref>.
Section <ref> derives some closed-form
expressions for the $E_{n,k}$ and briefly discusses combinatorial
interpretations.  Finally, ordering identities of the generalized
Eulerian kind are derived in Section <ref>.

§ NOTATION AND CONVENTIONS

Knuth's overline and underline notations [Knuth, 1992] are used for
rising and falling factorial powers.  For all $r\in\mathbb{C}$
and $m\in\mathbb{N}$, $(r)^{\overline m}$ equals
$(r)(r+\nobreak1)\dots\allowbreak(r+\nobreak m-\nobreak1)$, with
$(r)^{\overline 0}=1$.  Likewise, $(r)^{\underline m}$ equals
$(r)(r-\nobreak1)\dots\allowbreak(r-\nobreak m+\nobreak1)$, with
$(r)^{\underline 0}=1$.

These notations are extended here to include a non-unit “stride.”
For all $\alpha\in\mathbb{C}$, $(r)^{\overline m,\alpha}$ equals $(r)(r+\nobreak
\alpha)\dots\allowbreak[r+\nobreak (m-\nobreak 1)\alpha]$, and likewise,
$(r)^{\underline m,\alpha}$ equals $(r)(r-\nobreak
\alpha)\dots\allowbreak[r-\nobreak (m+\nobreak 1)\alpha]$.  Note also that
$\binom{n}k$ is defined for all $n,k\in\mathbb{N}$, with
$\binom{n}k=0$ if $k>n$, and that $0^0=1$.

Knuth's notations $\stirsub{n}k$, $\stircyc{n}k$, and $\eulerian{n}k$,
where $0\le k\le n$, are used for the classical Stirling numbers of
the second kind, the unsigned ones of the first kind, and the Eulerian
numbers. The last differs from the classical definition in having its
lower index shifted by unity.  That is, $\eulerian{n}k=A_{n,k+1}$.

The three-parameter generalized Eulerian numbers $E_{n,k}(\alpha,\beta;r)$
employed here are specializations of four-parameter ones of the
author.  In previously used notation and terminology [Maier, 2023],
the $E_{n,k}(\alpha,\beta;r)$ would be written as $E_{n,k}(\alpha,\beta;\allowbreak
r,\beta-\nobreak r)$ and called “single progression” generalized
Eulerian numbers.

The triangular recurrences derived below, satisfied by $E_{n,k}$
and $S_{n,k}$, are of the GKP type, $c_{n+1,k+1}=[\alpha n+\beta (k+1)
+ \gamma]\,c_{n,k+1} + [\alpha' n+\beta'k + \gamma']\,c_{n,k}$.  The
importance of such recurrences was noted by Graham, Knuth, and
Patashnik <cit.>, after whom they are

Several products $\mathcal{A}\mathcal{B}$ of infinite matrices appear,
the rows and columns of each matrix being indexed by $\mathbb{N}$.
In all cases $\mathcal{A}$ is row-finite (each of its rows has finite
support), or $\mathcal{B}$ is column-finite.  Hence each element
of $\mathcal{A}\mathcal{B}$ can be computed as a finite sum.

In a standard definition, terminating hypergeometric series are those
in which an “upper” parameter is a nonpositive integer, such as the
Gauss hypergeometric series

\begin{equation}
\biggl[
\begin{array}{cc}
-N, & b
\\
& c
\end{array}
\biggm|
\biggr]
\defeq
\sum_{k=0}^\infty \frac{(-N)^{\overline k}(b)^{\overline k}}{(c)^{\overline k}}\, \frac{z^k}{k!}\,.
\end{equation}

If $N\in\mathbb{N}$, one expects the series to terminate with the
$k=N$ term.  Usually the lower parameter $c$ is required not to be a
nonpositive integer, as the series is not defined due to division(s)
by zero if $c$ is an integer in the range ${-N\le c\le 0}$.  To cover
such cases, a desingularized function ${}_2\hat F_1$, defined
informally as $(c)^{\overline N}\,{}_2F_1$, is used here.  That is,

\begin{equation}
{}_2\hat F_1
\biggl[
\begin{array}{cc}
-N, & b
\\
& c
\end{array}
\biggm|
\biggr]
\defeq
\sum_{k=0}^N (-N)^{\overline k}(b)^{\overline k} (c+k)^{\overline {N-k}} \: \frac{z^k}{k!}\,,
\end{equation}

because $(c)^{\overline N}/(c)^{\overline k}=(c+k)^{\overline{N-k}}$
if there is no division by zero.  ${}_2\hat F_1(-N,b;c;z)$ when
$N\in\mathbb{N}$ is defined for all $c\in\mathbb{C}$, and is
continuous in $c$ as well as in $b$.  This desingularization of the
parametric function ${}_2F_1$ differs from that of
Olver <cit.>, in which it is multiplied by
$1/\Gamma(c)$ rather than $(c)^{\overline N}$.

§ THE WEYL–HEISENBERG ALGEBRA

The (one-dimensional) Weyl–Heisenberg algebra $W_{\mathbb{C}}$ is a
unital associative algebra over the field of scalars $\mathbb{C}$
generated by noncommuting indeterminates $\adag,a$.  Any element
of $W_{\mathbb{C}}$ is a combination over $\mathbb{C}$ of finite
products of $\adag,a$, and the unit (multiplicative identity) of the
algebra, conventionally written as $1$, is the empty product
multiplied by the scalar $1$.

The representation of elements of $W_{\mathbb{C}}$ is not unique, due
to the existence of the commutation relation $[a,\adag]=1$.  One can

\begin{equation}
W_{\mathbb{C}} = \mathbb{C}\langle \adag,a\rangle \,/\, \mathbb{C}(a\adag - \adag
a - 1),
\end{equation}

where $\mathbb{C}\langle \adag,a\rangle$ is the free associative
unital algebra over $\mathbb{C}$ generated by $\adag,a$, and the ideal
quotiented out is the two-sided one generated over $\mathbb{C}$ by the
element $a\adag - \nobreak\adag a - \nobreak1$, where $1$ is the unit
of the free algebra.  By definition, an ordering identity
in $W_{\mathbb{C}}$ is an element of this two-sided ideal.  For
example, one can prove by induction that for all $m,n\in\mathbb{N}$,

\begin{equation}
\label{eq:difflr}
a^m \adag^n= \sum_{\ell=0}^{\min(m,n)} \ell\,!\, \binom{m}{\ell}\binom{n}{\ell}\,
\adag^{(n-\ell)} a^{m-\ell}
\end{equation}

holds as an identity in $W_\mathbb{C}$.  This expresses an
anti-normally ordered monomial in $\adag,a$ in terms of normally
ordered ones.  Viewed as an element of $\mathbb{C}\langle
\adag,a\rangle$, the difference between the left-hand and right-hand
sides of (<ref>) is an element of the ideal
$\mathbb{C}( a\adag -\nobreak \adag a -\nobreak 1)$.

The polynomial algebra over $\mathbb{C}$ generated by $x,D=D_x={\rm
d}/{\rm d}x$ with $[D,x]=1$ is isomorphic to $W_{\mathbb{C}}$, but
including negative integer powers of $x$ as elements extends it to an
algebra that has $x^{-1}$ as an additional generator, subject to the
relation $[D,x^{-1}]=-(x^{-1})^2$.  The extension algebra $\tilde
W_{\mathbb{C}}$ used here is a larger, infinitely generated one: it
has each power $x^\alpha$, $\alpha\in\mathbb{C}$, as an element,
subject to $[D,x^\alpha]=\alpha(x^{\alpha-1})$.  It is also useful to
consider the deformation $W_{\mathbb{C}}(s)$ of $W_{\mathbb{C}}$
defined by $UV-\nobreak VU=V^s$, where $s\in\mathbb{C}$.  As noted,
when $s\in\mathbb{N}$ it is isomorphic to a subalgebra
of $W_{\mathbb{C}}$, with $V,U$ and $\adag,\adag^s a$ identified.
When $s\in\mathbb{C}\setminus{\mathbb{N}}$ it is a subalgebra
of $\tilde W_{\mathbb{C}}$, with $V,U$ and $x,x^s D$ identified.
Owing to the embedding $W_{\mathbb{C}}(s)\hookrightarrow W_{\mathbb{C}}$
(resp. $\tilde W_{\mathbb{C}}$), any ordering identity
in $W_{\mathbb{C}}(s)$ can be viewed as an ordering identity
in $W_{\mathbb{C}}$ (resp. $\tilde W_{\mathbb{C}}$).

The algebras $W_{\mathbb{C}}$ and $\tilde W_{\mathbb{C}}$ are
respectively $\mathbb{Z}$-graded and
$\mathbb{C}$-graded.  For any $\omega\in{W}_{\mathbb{C}}$
that is a nonzero scalar multiple of a product of $\adag$'s and $a$'s
(i.e., a monomial), let $|\omega|_\adag,\allowbreak |\omega|_a$ denote
the number of $\adag$'s and $a$'s.  Similarly, if an element $\tilde
\omega\in\tilde{W}_{\mathbb{C}}$ is a nonzero multiple of a product of
powers of $x$ and $D$, let $|\tilde \omega|_x,\allowbreak |\tilde
\omega|_D$ denote the sum of those powers: for instance,
$|x^{1/2}Dx^{-1/2}|_x=0$ and $|x^{1/2}Dx^{-1/2}|_D=1$.  The
excess (or degree) of $\omega$ is
$e(\omega)\defeq|\omega|_\adag-\nobreak |\omega|_a\in\mathbb{Z}$, and
that of $\tilde \omega$ is $e(\tilde \omega)\defeq|\tilde
\omega|_x-\nobreak |\tilde \omega|_D\in\mathbb{C}$.  These algebras
are graded by excess: any element can be written as a sum of groups of
terms, each group containing terms with a distinct excess.
For $W_{\mathbb{C}}$, the quotienting out of the ideal
$\mathbb{C}(a\adag -\nobreak \adag a -\nobreak 1)$
respects the grading, as each term in $a\adag -\nobreak \adag a
-\nobreak 1$ has zero excess; and similarly for $\tilde

As any ordering identity (i.e., element of the ideal) can be grouped
into terms with distinct excesses, it is reasonable to focus on ones
in which only a single excess is represented.  For instance, each term
in (<ref>) has excess $n-\nobreak m\in\mathbb{Z}$, and in
the identities (<ref>), (<ref>),
(<ref>), (<ref>) of the introduction,
the common excess is respectively $0,0,\allowbreak En,2n$.  The
algebra $W_{\mathbb{C}}$ has a
subalgebra $W_{\mathbb{C}}^{\textrm{bal}}$ (and $\tilde
W_{\mathbb{C}}$ has a subalgebra $\tilde
W_{\mathbb{C}}^{\textrm{bal}}$) that comprises all linear combinations
of “balanced” or “scale invariant” monomials: ones with excess $0$.
It can be shown that $W_{\mathbb{C}}^{\textrm{bal}}$ is singly
generated: it is $\mathbb{C}[\adag a]$, the algebra of
complex-coefficient polynomials in the boson number operator
$\mathcal{N}\defeq\adag a$, though it could also be written as
$\mathbb{C}[a \adag]$ or $\mathbb{C}[\adag a+ a\adag]$, etc.

Identities relevant to single-annihilator $\Omega(\adag,a)$ include
what is properly called the Tait–Toscano–Viskov
identity [Mansour and Schork, 2016, Toscano, 1955], which is $(\adag a \adag)^n =
\adag^na^n\adag^n$ in $W_{\mathbb{C}}$ and $(xDx)^n=x^nD^nx^n$
in $\tilde W_{\mathbb{C}}$; it holds for all $n\in\mathbb{N}$.
Another pair of identities is

\begin{align}
& \adag^L a \adag^R + \adag^R a \adag^L = 2(\adag)^{(L+R)/2} a (\adag)^{(L+R)/2},&&\text{in $W_{\mathbb{C}}$},\label{eq:otherpaira}\\
&x^L D x^R + x^R D x^L = 2x^{(L+R)/2} D x^{(L+R)/2},&&\text{in $\tilde W_{\mathbb{C}}$}.\label{eq:otherpairb}
\end{align}

These hold respectively for all $L,R\in\mathbb{N}$ with $L+\nobreak R$
even, and for all $L,R\in\mathbb{C}$.  The following will play a major
role in deriving operator ordering identities from polynomial

In $\tilde W_{\mathbb{C}}^{\rm{bal}}$, the identities

\begin{align}
\label{eq:majorrole1a}
(xD)^{\underline n,\alpha} &= (x^\alpha)^n (x^{1-\alpha}D)^n,\\
&= (xDx^\alpha)^n (x^{-\alpha})^n
\label{eq:majorrole1b}
\end{align}

hold for all $n\in\mathbb{N}$ and $\alpha\in\mathbb{C}$.  More generally,

\begin{align}
\label{eq:majorrole2a}
(xD+r)^{\underline n,\alpha} &= x^{-r}\left[(x^\alpha)^n (x^{1-\alpha}D)^n\right]x^r
\\
&= x^{-r}\left[ (xDx^\alpha)^n (x^{-\alpha})^n \right] x^r
= (x^{1-r}Dx^{\alpha+r})^n (x^{-\alpha})^n
\label{eq:majorrole2b}
\end{align}

hold for all $n\in\mathbb{N}$ and $\alpha,r\in\mathbb{C}$, with

\begin{equation}
\label{eq:powerfullemma}
(x^{1-L-R})^n(x^LDx^R)^n = (x^{1-R}Dx^{1-L})^n (x^{L+R-1})^n
\end{equation}

being a reparametrized statement of the equality between
(<ref>) and (<ref>), which
holds for all $n\in\mathbb{N}$ and $L,R\in\mathbb{C}$.

When $\alpha=1$, (<ref>) specializes to $(xD)^{\underline
n}=x^nD^n$, which is the well-known identity of
Boole [Mansour and Schork, 2016]: in the $W_{\mathbb{C}}^{\text{bal}}$ form
$(\adag a)^{\underline n}=\adag^n a^n$, it normally-orders any falling
factorial power of the number operator.  But
(<ref>),(<ref>) are not so well known,
though (<ref>) appears as Proposition 5
of [Mohammad-Noori, 2011].  The versions
(<ref>),(<ref>) come from them by
similarity transformations.  All these identities can be proved by
induction, or (in $\tilde W_{\mathbb{C}}$) operationally: showing that
the left and right-hand sides act identically on monomials or
polynomials in $x$ [Carlitz, 1976].

The algebras $W_{\mathbb{C}}$ and $\tilde W_{\mathbb{C}}$ have
nontrivial automorphisms, but the only one that will appear here is
the formal adjoint operation $*$ on $\tilde W_{\mathbb{C}}$ (strictly,
an anti-automorphism), which treats its elements as differential
operators.  For instance, $\left[ (x^L D x^R)^n\right]^*$ equals
$(-1)^{n} (x^R D x^L)^n$.  Up to parameter renamings, the adjoint
operation $*$ transforms (<ref>)
to (<ref>), and vice versa, and
(<ref>) is self-adjoint: it is transformed to itself.

§ GENERALIZED STIRLING AND EULERIAN NUMBERS

The unified generalized Stirling numbers
$S_{n,k}=S_{n,k}(\alpha,\beta;r)$ of Hsu and Shiue [Hsu and Shiue, 1998], and
related generalized Eulerian numbers
$E_{n,k}=E_{n,k}(\alpha,\beta;r)$, both of which are indexed by $n,k$
with $0\le k\le n$, can be defined in several ways.  They are given by
finite sums (see Theorem <ref> below), and they satisfy
triangular recurrences (see Theorem <ref>), which can
be used to compute them when $n$ and $k$ are small.  But they can also
be defined abstractly, as coefficients that connect alternative bases
of the linear space of polynomials of degree $\le n$ in an
indeterminate $x$.  This leads to many identities in polynomial

These numbers will be viewed as elements of infinite, lower triangular
matrices $\mathcal{S}=[S_{n,k}]$ and $\mathcal{E}=[E_{n,k}]$, indexed
by $n,k\in\mathbb{N}$.  The triangles $\mathcal{S}$ and $\mathcal{E}$
are parametrized by $\alpha,\beta,r\in\mathbb{C}$, which specify the
bases.  As will be seen, the modified triangle
$\hat{\mathcal{S}}=[\hat S_{n,k}]\defeq[\beta^kk!\,S_{n,k}]$ and the
triangle $\mathcal{E}$ are in a sense dual to each other.  To reveal
the symmetry, many identities in this section are written in terms of
$\hat S_{n,k}$ rather than $S_{n,k}$.  Though a (weak) case can be
made that the generalized Stirling numbers should be redefined to
include the factor $\beta^kk!$, the original
normalization [Hsu and Shiue, 1998] will be maintained here.

In the linear-algebraic approach, the numbers $\hat S_{n,k}$
and $E_{n,k}$, $0\le k\le n$, are defined to satisfy

\begin{align}
\label{eq:Sdef2}
(\beta x+r)^{\underline n,\alpha} &= \sum_{k=0}^n \hat
\label{eq:Edef2}
(\beta x+r)^{\underline n,\alpha} &= \sum_{k=0}^n
\end{align}

The sets of polynomials $\left[\binom{x}k\right]_{k=0}^n$ and
$\left[\binom{x+n-k}n\right]_{k=0}^n$ are bases of the space of
polynomials of degree $\le n$ in $x$ (the latter statement requiring a
proof [Carlitz, 1978]), so $\hat S_{n,k}$ and $E_{n,k}$ are uniquely
defined.  They are continuous functions of $\alpha,\beta,r$.  If
$\beta\neq0$, replacing $x$ by $(x-\nobreak r)/\beta$ and defining
$S_{n,k}\defeq \hat S_{n,k}/\beta^kk!$ produces the less symmetrical pair
of definitions

\begin{align}
\label{eq:Sdef1}
(x)^{\underline n,\alpha} &= \sum_{k=0}^n S_{n,k}(\alpha,\beta;\,r)\,(x-r)^{\underline k,\beta},
\\
\label{eq:Edef1}
\beta^nn!\,(x)^{\underline n,\alpha} &= \sum_{k=0}^n E_{n,k}(\alpha,\beta;\,r)\,\left[x-r+\beta(n-k)\right]^{\underline n,\beta}.
\end{align}

The definition (<ref>) is that of Hsu and Shiue, and it
defines the $S_{n,k}$ uniquely even when $\beta=0$, as continuous
functions of $\alpha,\beta,r$.  (Note that if $\beta=0$, $\hat
S_{n,k}$ vanishes unless $k=0$, but the same is not true
of $S_{n,k}$.)  A consequence of (<ref>) is that the
parametric triangle $\mathcal{S}=[S_{n,k}]$ satisfies

\begin{equation}
\label{eq:Smultrule}
\mathcal{S}(\alpha,\gamma;\,r_1+r_2)= \mathcal{S}(\alpha,\beta;\,r_1)\,\mathcal{S}(\beta,\gamma;\,r_2),
\end{equation}

matrix multiplication being meant.  In accordance
with (<ref>), $\mathcal{S}(\alpha,\alpha;0)$ for
any $\alpha$ equals the infinite identity matrix $\mathcal{I}$, and
for any $\alpha,\beta,r$, the inverse matrix
$\mathcal{S}(\alpha,\beta;r)^{-1}$ exists and equals

For the $E_{n,k}$, the definition (<ref>) has been used by
the author [Maier, 2023].  It defines the $E_{n,k}$ uniquely when
$\beta\neq0$, though (<ref>) must be used when $\beta=0$.
When $(\alpha,\beta;r)=(0,1;0)$, it reduces to Worpitzky's
identity [Graham et al., 1994], which is satisfied by the classical Eulerian
numbers.  The case when $\alpha=0$ but $\beta,r$ are arbitrary, which
is intermediate between the classical case and the present one, has
recently been treated [Mansour et al., 2019, Mező and Ramírez, 2016, Ramírez et al., 2018].

Eq. (<ref>) has a matrix-algebraic interpretation.  Let
$\mathcal{V}(\alpha,\beta;r)=\allowbreak\left[(\beta x+\nobreak r)^{\underline
n,\alpha}\right]$ be an infinite matrix of generalized Vandermonde type
indexed by $n,x\in\mathbb{N}$, and let
$\mathcal{P}=\left[{\binom{x}{k}}\right]$ be an infinite Pascal matrix
indexed by $x,k\in\mathbb{N}$, which is lower triangular (nonzero
elements having $0\le k\le x$).  Also, let ${\mathcal{P}}^t$ be its
upper-triangular transpose.  Then (<ref>) says that

\begin{equation}
\label{eq:LU}
\mathcal{V}(\alpha,\beta;\,r) = \hat{\mathcal{S}}(\alpha,\beta;\,r)\,{\mathcal{P}}^t,
\end{equation}

which defines $\hat{\mathcal{S}}(\alpha,\beta;r)$ as the
$\text{L}$-factor in a certain $\text{LU}$ decomposition
of $\mathcal{V}(\alpha,\beta;r)$.  For a related interpretation
of (<ref>), see the remark after Theorem <ref> below.

It follows by inspection that irrespective of the choice of parameters
$\alpha,\beta,r$, the apex elements $S_{0,0}$, $\hat S_{0,0}$, and
$E_{0,0}$ of the triangles $\mathcal{S}$,
$\hat{\mathcal{S}}$, $\mathcal{E}$ equal unity.  The homogeneity

\begin{equation}
\begin{gathered}
S_{n,k}(\lambda \alpha,\lambda \beta;\,\lambda r) = \lambda^{n-k} S_{n,k}(\alpha,\beta;\,r)
\\
\hat S_{n,k}(\lambda \alpha,\lambda \beta;\,\lambda r) = \lambda^n \hat S_{n,k}(\alpha,\beta;\,r),\qquad
E_{n,k}(\lambda \alpha,\lambda \beta;\,\lambda r) = \lambda^n E_{n,k}(\alpha,\beta;\,r)
\label{eq:homogeneity}
\end{gathered}
\end{equation}

also follow by inspection.  Each $S_{n,k}$ and $E_{n,k}$ turns out to
be a polynomial in $\alpha,\beta,r$ with integral coefficients (see
below); so if $\alpha,\beta,r$ are integers, each $S_{n,k}$ and
$E_{n,k}$ will be an integer.

For some choices of $\alpha,\beta,r\in\mathbb{Z}$, all $S_{n,k}$
and/or $E_{n,k}$ may be non-negative integers, in which case they may
have combinatorial interpretations.  This is the case for many
$\alpha,\beta,r$, as can be confirmed by comparing the recurrences
satisfied by the $S_{n,k}$ and $E_{n,k}$ (see
Theorem <ref> below) with ones that have appeared in
the literature [Comtet, 1974, OEIS Foundation, Inc., 2023].  Examples of such
triangularly indexed combinatorial numbers include the Stirling subset
numbers $\stirsub{n}k=S_{n,k}(0,1;0)$, which count the partitions of
an $n$-set into $k$ non-empty blocks.  They are often
called the Stirling numbers of the second kind and are also denoted by
$S(n,k)$.  (When $(\alpha,\beta;r)=(0,1;0)$, the Vandermonde
factorization (<ref>) is well known [Cheon and Kim, 2001].) Another
example is provided by the Eulerian numbers $\eulerian{n}k =
E_{n,k}(0,1;1)$, which count the permutations of an
$n$-set that have exactly $k$ descents
[Graham et al., 1994, Riordan, 1958]. The traditionally indexed Eulerian numbers
$A_{n,k}$ are horizontally shifted versions: $A_{n,k+1}$
equals $\eulerian{n}k$, with $A_{n,0}=0$ unless $n=0$; and

Known $S_{n,k}$ with $\beta=0$ include $\binom{n}k=S_{n,k}(0,0;1)$ and
the Stirling cycle numbers $\stircyc{n}k=S_{n,k}(-1,0;0)$, which count
the permutations of an $n$-set, the cycle decompositions
of which have exactly $k$ cycles.  They are often called the unsigned
Stirling numbers of the first kind and are also denoted by
$(-1)^{n-k}s(n,k)$, the numbers $s(n,k)=S_{n,k}(1,0;0)$ being the
signed versions.  It follows from (<ref>) that
$\mathcal{S}(\alpha,\beta;\,r) = \mathcal{S}(\alpha,0;\,0)\,
\mathcal{S}(0,0;\,r)\, \mathcal{S}(0,\beta;\,0)$, i.e.,

\begin{equation}
\label{eq:elegant}
S_{n,k}(\alpha,\beta;\,r) = \sum_{j=k}^n\sum_{p=k}^j
(-\alpha)^{n-j}\stircyc{n}{j} r^{j-p}\binom{j}{p}
\beta^{p-k}\stirsub{p}k.
\end{equation}

The existence of this elegant decomposition
formula [Spivey, 2011, Pan, 2012] justifies keeping the Hsu–Shiue
definition of $S_{n,k}$, as $\stircyc{n}j$ can be expressed in terms
of $\hat S_{n,j}$ only by taking a $\beta\to0$ limit.

Relatively simple formulas for the $\hat S_{n,k}$ and $E_{n,k}$ can be
obtained by extracting coefficients from (<ref>)
and (<ref>).  The following can be used to compute $\hat
S_{n,k}$ (and $S_{n,k}$ if $\beta\neq0$), and $E_{n,k}$, for all $k$ with
$0\le k\le n$, though the effort needed grows rapidly with $n$
and $k$.

The triangles $\hat{\mathcal{S}}=[\hat S_{n,k}]$ and
$\mathcal{E}=[{E}_{n,k}]$ are alternatively defined by

\begin{align}
\label{eq:Sform}
\begin{aligned}
\hat S_{n,k}(\alpha,\beta;\,r) &= \Delta_x^k\left[(\beta x+r)^{\underline n,\alpha}\right]\bigm|_{x=0}
\\
&= \sum_{x=0}^k (-1)^{k-x} \binom{k}{x} (\beta x+r)^{\underline n,\alpha},
\end{aligned}
\\
\shortintertext{and}
\label{eq:Eform}
\begin{aligned}
E_{n,k}(\alpha,\beta;\,r) &= \nabla_x^{n+1}\left[(\beta x+r)^{\underline n,\alpha}\right]\bigm|_{x=k}
\\
&= \sum_{x=0}^k (-1)^{k-x} \binom{n+1}{k-x} (\beta x+r)^{\underline n,\alpha},
\end{aligned}
\end{align}

with each right-hand side equaling zero unless $0\le k\le n$.

Acting on functions of $x=0,1,2,\dots$, the forward and backward
difference operators, $\Delta_x$ and $\nabla_x$, are defined by
$(\Delta_x f)(x) =\allowbreak f(x+\nobreak 1)-\nobreak f(x)$ and
$(\nabla_x f)(x) =\allowbreak f(x) -\nobreak f(x-\nobreak 1)$, with
$f(-1)$ taken to equal zero in the latter.  From
$\Delta_x^j\left[(x)^{\underline k}\right] = k^{\underline
j}(x)^{\underline{k-j}}$, it follows readily that

\begin{align}
&\frac1{j!}\Delta_x^j\left[(x)^{\underline k}\right]\bigm|_{x=0} = \delta_{j,k},\label{eq:extract1}\\
&\frac1{n!}\nabla_x^{n+1}\Bigl\{\left[ x+(n-k)\right]^{\underline n}\Bigr\}\bigm|_{x=j} = \delta_{j,k},
\label{eq:extract2}
\end{align}

where $\delta_{j,k}$ is the Kronecker delta.  These statements hold
for all $j,k,n\in\mathbb{N}$, with $k\le n$ assumed in the latter.
Applying (<ref>),(<ref>) to
(<ref>),(<ref>) yields

Formulas (<ref>) and (<ref>) are generalizations of
known ones for the Stirling subset and classical Eulerian
numbers [Graham et al., 1994]. For each value of the row
index $n\in\mathbb{N}$, they express $\hat S_{n,k} =
\beta^kk!\,S_{n,k}$ and $E_{n,k}$, viewed as functions of the column
index $k=0,\dots,n$, as binomial transforms of $(\beta x+\nobreak
r)^{\underline n,\alpha}$, viewed as a function of $x=0,\dots,n$.  In
particular, (<ref>) says that
Having the factor $(-1)^{k-x}$ in their summands,
(<ref>) and (<ref>) are properly called
inverse binomial transforms.  According to the general theory
of binomial transforms [Boyadzhiev, 2018], they are inverse to the
forward transforms (<ref>) and (<ref>), which
include no such factor, and are their immediate corollaries.

In fact, corresponding rows of the triangles $\hat{\mathcal{S}}=[\hat
S_{n,k}]$ and $\mathcal{E}=[E_{n,k}]$ are binomial transforms of
each other, so each triangle can be computed from the other.  The
following theorem, the $\alpha=0$ case of which has previously
appeared [Mező and Ramírez, 2016, Ramírez et al., 2018], generalizes the known pair of
relations [Graham et al., 1994] between the Stirling subset numbers
$\stirsub{n}k=S(n,k)=S_{n,k}(0,1;0)$ and the traditionally indexed
Eulerian numbers $A_{n,k}=\eulerian{n}{k-1}=E_{n,k}(0,1,0)$.  The
involutive “rev” transformation
$\mathcal{M}\mapsto{\mathcal{M}}^{\rm rev}$ is defined to produce from
any infinite lower-triangular matrix a like matrix in which each row
is reversed.  If $\mathcal{M}=[M_{n,k}]$, then ${\mathcal{M}}^{\rm
rev}=[M^{\rm rev}_{n,k}]$ has $M^{\rm rev}_{n,k} = M_{n,n-k}$,

For all $n\ge0$ and $\alpha,\beta,r$, one has

\begin{align}
\hat S_{n,k} &= \sum_{j=0}^k \binom{n-j}{n-k} E_{n,j},\label{eq:EtoS}\\
E_{n,k} &= \sum_{j=0}^k (-1)^{k-j} \binom{n-j}{n-k}\, \hat S_{n,j},\label{eq:StoE}
\end{align}

which are a binomial transform and its inverse.  That is,
$\hat{\mathcal{S}}^{\rm{rev}} = {\mathcal{E}}^{\rm{rev}}\,\mathcal{P}$
and ${\mathcal{E}}^{\rm{rev}} =
\hat{\mathcal{S}}^{\rm{rev}}\,\mathcal{P}^{-1}$.

The basis functions (of $x$) that appear in the original definitions
(<ref>) and (<ref>) are related by the identities

\begin{align}
&\binom{x+n-k}{n} = \sum_{j=k}^n \binom{n-k}{n-j}\binom{x}{j},\label{eq:elident1}\\
&\binom{x}{k} = \sum_{j=k}^n (-1)^{j-k}\binom{n-k}{n-j}\binom{x+n-j}{n}.\label{eq:elident2}
\end{align}

The first is a version of Vandermonde's identity, which has an
elementary combinatorial proof; one can also view it as a binomial
transform, in which case the second is its inverse.  (It will be
recalled that if $\mathcal{P} = \left[{\binom{x}k}\right]$, indexed
by $x,k$, then $\mathcal{P}^{-1} =
\left[(-1)^{k-x}{\binom{k}x}\right]$, indexed by $k,x$.)  To verify
(<ref>), substitute it into (<ref>) and use
(<ref>) to reduce the double sum to a single sum; one
obtains (<ref>).  To verify (<ref>), substitute it
into (<ref>) and use (<ref>) similarly; one
obtains (<ref>).

As $\hat{\mathcal{S}}^{\text{rev}} =
{\mathcal{E}}^{\text{rev}}\,\mathcal{P}$, the Vandermonde
decomposition (<ref>) can be written as

\begin{equation}
\mathcal{V}(a,b;\,r) = \left({\mathcal E}^\text{rev}\,\mathcal{P}\right)^{\text{rev}}
\mathcal{P}^t.
\end{equation}

This is the promised matrix-algebraic restatement of (<ref>),
which defines the matrix $\mathcal{E}=\mathcal{E}(a,b;r)$.
Equivalently, ${\mathcal{E}}^{\rm rev} = \left[\mathcal{V}
\left(\mathcal{P}^{-1}\right)^t\right]^{\rm rev}\mathcal{P}^{-1}$,
which is a restatement of the summation formula (<ref>).

One can work out generating functions for the numbers $\hat S_{n,k}$
and $E_{n,k}$, and triangular recurrences which they satisfy.  Define
the parametric row polynomials $\hat S_n(t)$ and $E_n(t)$,
$n\in\mathbb{N}$, by

\begin{equation}
\hat S_n(t) = \sum_{k=0}^n \hat S_{n,k}\,t^k,\qquad
E_n(t) = \sum_{k=0}^n E_{n,k}\,t^k,
\end{equation}

and the parametric bivariate exponential generating functions (EGF's)
$\hat S(t,z)$ and $E(t,z)$ by

\begin{equation}
\hat S(t,z) = \sum_{n=0}^\infty\sum_{k=0}^n \hat S_{n,k}\,t^k\frac{z^n}{n!},\qquad
E(t,z) = \sum_{n=0}^\infty\sum_{k=0}^n E_{n,k}\,t^k\frac{z^n}{n!}.
\end{equation}

(Each depends on $\alpha,\beta,r$.)  It is not difficult to see that the pair

\begin{equation}
\hat S_n(t) = (1+t)^nE_n\left(t/(1+t)\right),\qquad
E_n(t) = (1-t)^n\hat S_n\left(t/(1-t)\right)
\end{equation}

and the pair

\begin{equation}
\label{eq:relationpair}
\hat S(t,z) = E\left( t/(1+t), (1+t)z\right),\qquad
E(t,z) = \hat S\left( t/(1-t), (1-t)z\right)
\end{equation}

are each equivalent to (<ref>),(<ref>).

When $\alpha\neq0$ and $\beta\neq0$, one has the explicit EGF's

\begin{align}
\hat S(t,z) &=\frac{(1+\alpha z)^{r/\alpha}}{1-t\left[(1+\alpha z)^{\beta/\alpha}-1\right]}, \label{eq:SEGF} \\
E(t,z) &= \frac{(1-t)(1+\alpha z-\alpha tz)^{r/\alpha}}{1-t(1+\alpha z-\alpha tz)^{\beta/\alpha}}.
\label{eq:EEGF}
\end{align}

(Taking the $\alpha\to0$ limit is straightforward; the $\beta\to0$ limit is not
discussed here.)

Define the forward difference operator $\Delta_{r,\beta}$, acting on
functions of $r$, by $(\Delta_{r,\beta} f)(r) =\allowbreak f(r+\nobreak
\beta)-\nobreak f(r)$.  Then by (<ref>),

\begin{equation}
\hat S_{n,k}(\alpha,\beta;\,r) = \Delta_x^k\left[(\beta x+r)^{\underline n,\alpha}\right]\bigm|_{x=0}
= \Delta_{r,\beta}^k\left[(r)^{\underline n,\alpha}\right]
\end{equation}

which vanishes if $k>n$, so one has the “vertical” univariate EGF

\begin{equation}
\begin{aligned}
\sum_{n=k}^\infty \hat S_{n,k}(\alpha,\beta;\,r)\frac{z^n}{n!}
\sum_{n=0}^\infty
(r)^{\underline n,\alpha}
\frac{z^n}{n!}
\right]
\Delta^k_{r,\beta}\left[
(1+\alpha z)^{r/\alpha}
\right]
\\
\qquad=
\sum_{j=0}^k(-1)^{k-j}\binom{k}j (1+\alpha z)^{(r+j\beta)/\alpha}
=(1+\alpha z)^{r/\alpha}\left[(1+\alpha z)^{\beta/\alpha}-1\right]^k.
\end{aligned}
\end{equation}

Multiplying this by $t^k$ and summing from $k=0$ to $\infty$
yields the bivariate EGF (<ref>), and (<ref>) comes by
applying (<ref>).

By direct calculation, $\hat S(t,z)$ and $E(t,z)$ given by
(<ref>) and (<ref>) satisfy the partial differential
equations (PDE's)

\begin{align}
\frac{\partial\hat S}{\partial z}
&= -\alpha z \frac{\partial\hat S}{\partial z} + \beta(1+t)t\frac{\partial \hat S}{\partial t} +(r+\beta t)\hat S,\\
\frac{\partial E}{\partial z} &=
\left[-\alpha+(\alpha+\beta)t\right]z\frac{\partial E}{\partial z} + \beta(1-t)t\frac{\partial E}{\partial t} + \left[r+(\beta-r)t\right]E.
\end{align}

The numbers $\hat S_{n,k}$ and $E_{n,k}$ satisfy

\begin{align}
\hat S_{n+1,k+1} &= \left[-\alpha n+\beta (k+1)+r\right] \hat S_{n,k+1} + \beta(k+1)\hat S_{n,k},\label{eq:SrecGKP}\\
E_{n+1,k+1} &= \left[-\alpha n+\beta(k+1)+r\right] E_{n,k+1} + \left[(\alpha+\beta)n-\beta k+(\beta-r)\right] E_{n,k},\label{eq:ErecGKP}
\end{align}

-and the unmodified numbers $S_{n,k}$ satisfy

\begin{equation}
\label{eq:SGKPrec}
S_{n+1,k+1} = \left[-\alpha n+\beta (k+1)+r\right]  S_{n,k+1} +  S_{n,k},
\end{equation}

all three parametric recurrences being of the triangular, GKP type.

As $\hat S_{n,k}$, $E_{n,k}$ are the coefficients of $t^kz^n/n!$ in
$\hat S(t,z)$, $E(t,z)$, (<ref>) and (<ref>)
follow from (<ref>) by equating coefficients.  If $\beta\neq0$,
(<ref>) follows from (<ref>) by $\hat S_{n,k} =
\beta^k k!\, S_{n,k}$; and if $\beta=0$, by continuity.

The recurrences (<ref>) and (<ref>) can be used
as alternative definitions of the $\hat S_{n,k}$ and $E_{n,k}$.  By
starting with the apex elements ${\hat S_{0,0}=1}$ and ${E_{0,0}=1}$,
and iterating, one can compute any desired $\hat S_{n,k}$
and $E_{n,k}$, though the effort needed grows rapidly with $n$
and $k$.  Each $\hat S_{n,k}$ and $E_{n,k}$ will be a homogeneous
polynomial in $\alpha,\beta,r$, of joint degree $n$ and with integral
coefficients.  The useful reversal (or “reflection”) identity
$\mathcal{E}^{\text{rev}}(\alpha,\beta;r) = \mathcal{E}(-\alpha,\beta;\beta-\nobreak r)$,

\begin{equation}
\label{eq:reflectionprop}
E_{n,n-k}(\alpha,\beta;\,r) = E_{n,k}(-\alpha,\beta;\,\beta-r),
\end{equation}

follows from (<ref>) by substituting $n-\nobreak k$
for $k$.  When $(\alpha,\beta;r)=(0,1;1)$, it reduces to the known reversal
identity $\eulerian{n}{n-k}=A_{n,k}$, i.e., the property
$\eulerian{n}{n-k}=\eulerian{n}{k-1}$ of the classical Eulerian
numbers [Graham et al., 1994].

$S_{n,k}$ is similarly a homogeneous polynomial in $\alpha,\beta,r$
with integral coefficients, but it is of joint degree $n-\nobreak k$,
as follows from the decomposition formula (<ref>) and also
from (<ref>).  By applying (<ref>) and
(<ref>),(<ref>) repeatedly, one finds that the
elements on the left edges of these triangles are $S_{n,0} = \hat
S_{n,0} = E_{n,0} = (r)^{\underline n,\alpha}$, and those on the right
edges are $S_{n,n}=1$, $\hat S_{n,n}=\beta^nn!$, and
$E_{n,n}=(\beta-\nobreak r)^{\overline n,\alpha}$.

It is worth noting here that the
generalized Vandermonde factorization (<ref>) can be written
in terms of the unmodified matrix $\mathcal{S}$ as

\begin{equation}
\label{eq:LDU}
\mathcal{V}(\alpha,\beta;\,r) = \mathcal{S}(\alpha,\beta;\,r)
\mathcal{D}(\beta){\mathcal{P}}^t,
\end{equation}

where $\mathcal{D}(\beta)$ is the infinite diagonal matrix with
$(k,k)$ element equal to $\beta^k/k!$.  This has the form of a
standard $\text{LDU}$ decomposition, because each diagonal element of
the $\text{L}$-factor $\mathcal{S}$ equals unity, as is
the case with the $\text{U}$-factor $\mathcal{P}^t$.  When
$(\alpha,\beta;r)$ equals $(0,1;0)$ and
$\mathcal{V}(\alpha,\beta;r)=\allowbreak \left[(\beta x+\nobreak
r)^{\underline n,\alpha}\right]$ becomes the infinite Vandermonde
matrix $\left[x^n\right]$, indexed by $n,x\in\mathbb{N}$, it
specializes to the well-known Vandermonde factorization mentioned
above, in which $\mathcal{S}$ is the matrix of Stirling numbers of the
second kind [Cheon and Kim, 2001].

Additional identities can be deduced from the bivariate EGF's
in (<ref>).  Again by direct calculation, they satisfy
difference equations on $r$, viz.,

\begin{gather}
\label{eq:Sdiffrec}
\begin{aligned}
&\Delta_{r,\alpha}\hat S(t,z) = \alpha z\hat S(t,z),\qquad  t\Delta_{r,\beta}\hat S(t,z) = \hat S(t,z)-\hat S(0,z),
\end{aligned}
\\
\shortintertext{and}
\label{eq:Ediffrec}
\begin{aligned}
\begin{aligned}
&\Delta_{r,\alpha} E(t,z) = \alpha(1-t)z\,E(t,z),\\
t\Delta_{r,\beta}E(t,z) &= (1-t)\left[E(t,z) -E\left(0, (1-t)z\right)\right]
\\
&= (1-t)\left[E(t,z)-(1+\alpha z-\alpha tz)^{r/\alpha}\right],
\end{aligned}
\end{aligned}
\end{aligned}
\end{gather}

with (<ref>) and (<ref>) related to each other
by (<ref>).  Equating coefficients yields the

The parametric triangles $\hat{\mathcal S}=[\hat S_{n,k}]$ and
${\mathcal E}=[E_{n,k}]$ satisfy the difference equations

\begin{gather}
\begin{aligned}
& \Delta_{r,\alpha} \hat S_{n+1,k}(\alpha,\beta;\,r) = \alpha(n+1)\hat S_{n,k}(\alpha,\beta;\,r),
\label{eq:Sdiffonr}
\\
& \Delta_{r,\beta} \hat S_{n,k}(\alpha,\beta;\,r) = \hat S_{n,k+1}(\alpha,\beta;\,r),
\end{aligned}
\\
\shortintertext{and}
\begin{aligned}
& \Delta_{r,\alpha} E_{n+1,k+1}(\alpha,\beta;\, r) = \alpha(n+1)\left[E_{n,k+1}-E_{n,k}\right](\alpha,\beta;\,r),
\\
\begin{aligned}
\Delta_{r,\beta} E_{n,k}(\alpha,\beta;\,r) &= \left[E_{n,k+1}-E_{n,k}\right](\alpha,\beta;\,r) +(-1)^k \binom{n+1}{k+1} E_{n,0}(\alpha,\beta;\,r)   \\
& = \left[E_{n,k+1}-E_{n,k}\right](\alpha,\beta;\,r) +(-1)^k \binom{n+1}{k+1}(r)^{\underline n,\alpha}
\end{aligned}
\end{aligned}
\label{eq:gaa}
\end{gather}

on the third parameter $r$.

These difference equations can be iterated, so as to increment or
decrement the parameter $r$ of any $\hat{S}_{n,k}(\alpha,\beta;r)$ or
$E_{n,k}(\alpha,\beta;r)$ by any desired positive integer multiple of
$\alpha$ or $\beta$.

One can also construct Newton series, which are finite-difference
analogues of Maclaurin or Taylor series.  For a function $f(r)$, the
Newton series at $r=0$ is

\begin{equation}
\label{eq:diffMacl}
f(r) = \sum_{m=0}^\infty \frac{\Delta_{r,c}^mf(0)}{c^mm!} (r)^{\underline m,c},
\end{equation}

where $c\neq0$ is arbitrary.  It follows from (<ref>)

\begin{equation}
\begin{aligned}
& \Delta_{r,\alpha}^m \hat S_{n,k}(\alpha,\beta;\,0) = \alpha^m(n)^{\underline m}\,\hat S_{n-m,k}(\alpha,\beta;\,0),
\\
& \Delta_{r,\beta}^m \hat S_{n,k}(\alpha,\beta;\,0) = \hat S_{n,k+m}(\alpha,\beta;\,0),
\end{aligned}
\end{equation}

and (<ref>) therefore yields the formulas

\begin{equation}
\begin{aligned}
\hat{S}_{n,k}(\alpha,\beta;\,r) = \sum_{m=0}^{n-k} \binom{n}{m} \hat S_{n-m,k}(\alpha,\beta;\,0)\,(r)^{\underline m,\alpha},
\\
\hat{S}_{n,k}(\alpha,\beta;\,r) = \sum_{m=0}^{n-k} \frac{1}{\beta^mm!} \hat S_{n,k+m}(\alpha,\beta;\,0)\,(r)^{\underline m,\beta},
\end{aligned}
\label{eq:Shatsumformulas}
\end{equation}

it being assumed in the latter than $\beta\neq0$.  Thus from any row
or column of the matrix $\hat{\mathcal{S}}(\alpha,\beta;0)$, any
element of the corresponding row or column
of $\hat{\mathcal{S}}(\alpha,\beta;r)$ can be computed as a finite
sum.  In terms of the unmodified numbers $S_{n,k}(\alpha,\beta;r)$,
the formulas (<ref>) become

\begin{equation}
\begin{aligned}
{S}_{n,k}(\alpha,\beta;\,r) = \sum_{m=0}^{n-k} \binom{n}{m} S_{n-m,k}(\alpha,\beta;\,0)\,(r)^{\underline m,\alpha},
\\
{S}_{n,k}(\alpha,\beta;\,r) = \sum_{m=0}^{n-k} \binom{k+m}m   S_{n,k+m}(\alpha,\beta;\,0)\,(r)^{\underline m,\beta},
\end{aligned}
\end{equation}

with no restriction on $\beta$.  As $S_{n,k}(\beta,\beta;\,r) =
\binom{n}k (r)^{\underline{n-k},\beta}$ (see Theorem <ref>
below), these can be seen as specializations of the matrix product
formula (<ref>).

§ STIRLING-NUMBER FORMULAS AND COMBINATORICS

When deriving and interpreting operator ordering identities with
parametric matrix elements $S_{n,k}(\alpha,\beta;r)$ as coefficients,
it is desirable to be able to express them in closed form, and/or
interpret them combinatorially.  For certain choices of $\alpha,\beta$
and $r$, both are possible.

The classical Stirling numbers $S_{n,k}(0,1;0)\eqdef\stirsub{n}k$ and
$S_{n,k}(-1,0;0)\eqdef\stircyc{n}k$ have already been mentioned, and
from the recurrence (<ref>) it can be proved by induction
that $S_{n,k}(0,1;1)=\stirsub{n+1}{k+1}$ and
$S_{n,k}(-1,0;1)=\stircyc{n+1}{k+1}$.  In fact when $r\in\mathbb{N}$,
the so-called $r$-Stirling numbers
$S_{n,k}(0,1;r)\eqdef{\stirsub{n}k}_r$ and
$S_{n,k}(-1,0;r)\eqdef{\stircyc{n}k}_r$ have combinatorial
interpretations [Broder, 1984, Maier, 2023].  For instance,
${\stirsub{n}{k}}_r$ counts the number of partitions of an
$(n+\nobreak r)$-set into $k+\nobreak r$ blocks, such that
$r$ distinguished elements are placed in distinct blocks.  But the
focus in this section is on triples $\alpha,\beta,r$ for which a
closed-form expression for $S_{n,k}(\alpha,\beta;r)$ exists, obviating
the need for an iterative or recursive calculation.  For deriving
additional expressions, homogeneity should be recalled:
$S_{n,k}(\lambda \alpha,\lambda \beta;\lambda r)$ equals

One has the formulas

(i) $S_{n,k}(1,1;\,0) = \delta_{n,k}$,

(ii)  $S_{n,k}(-1,1;\,0) = \binom{n}k (k)^{\overline{n-k}}$,

(iii) $S_{n,k}(1,2;\,0) = 2^{-(n-k)}\,\frac{n!}{k!}\binom{k}{n-k}$,

(iv) $S_{n,k}(-2,-1;\,0) =2^{-(n-k)}\, \frac{(n)^{\overline{n-k}}\,(k)^{\overline{n-k}}}{(n-k)!}$,

which can be extended to $r\in\mathbb{C}$:

(i $^\prime$) $S_{n,k}(\beta,\beta;\,r) = \binom{n}k (r)^{\underline{n-k},\beta}$,

(ii $^\prime$)  $S_{n,k}(-\beta,\beta;\,r) = \binom{n}k (\beta k+r)^{\overline{n-k},\beta}$,

(iii $^\prime$) $S_{n,k}(1,2;\,r) = 2^{-(n-k)}\,\binom{n}{k} \, {}_2{\hat F}_1\Bigl[
\begin{smallarray}{c}-(n-k),\:-r \\ -n+2k+1\end{smallarray}
\Bigm|2\,
\Bigr]$,

(iv $^\prime$) $S_{n,k}(-2,-1;\,r) = (-2)^{-(n-k)}\,\binom{n}k \, {}_2{\hat
\begin{smallarray}{c}-(n-k),\:r-1 \\ -2n+k\end{smallarray}
\Bigm|2\, \Bigr]$.

Formulas (i)–(iv), (i$^\prime$)–(ii$^\prime$) are verified by
confirming that they satisfy the recurrence (<ref>).
Formulas equivalent to (iii$^\prime$),(iv$^\prime$) were derived
in [Cheon et al., 2013] and restated in ${}_2F_1$ terms by the
author [Maier, 2023]; they are stated here in terms of ${}_2\hat
F_1$.  For an alternative derivation of
(iii$^\prime$),(iv$^\prime$), see the proof of
Theorem <ref> below.

Combinatorial interpretations of these $S_{n,k}(a,b;r)$ can be deduced
from the triangular recurrence (<ref>), and a wealth of
information on these and other number triangles can be found in the
OEIS [OEIS Foundation, Inc., 2023].

The $S_{n,k}(-1,1;0) \eqdef \hat L_{n,k}$ in (ii) are the unsigned Lah
numbers.  The $S_{n,k}(1,2;0)\eqdef B_{n,k}$ in (iii) are the Bessel
numbers of the second kind, and the $S_{n,k}(-2,-1;0)\eqdef \hat
b_{n,k}$ in (iv) are the unsigned Bessel numbers of the first kind.
They appear respectively in the OEIS as A271703,
A122848 (reversed as A100861),
and A132062 (reversed as A001498).  The
formulas (ii),(iii),(iv) are equivalent to those given by
Comtet <cit.>.  When $n\ge1$, the formula (iv) can be
rewritten in terms of factorials in the form that appears in
eq. (<ref>) of the introduction.

These three triangles are the $r=0$ cases of the more general ones in
(ii$^\prime$),(iii$^\prime$),(iv$^\prime$), which have combinatorial
interpretations when $r\in\mathbb{N}$.

* The $S_{n,k}(-1,1;r)\eqdef \hat L_{n,k}^{(r)}$ in (ii$^\prime$)
are what will be called the unsigned $r$-Lah numbers,
which when $r$ is even, count the partitions of a set with
$(n+\nobreak r/2)$ elements into $k+\nobreak r/2$ ordered blocks or
lists, such that $r/2$ distinguished elements of the set are placed
in distinct lists [Nyul and Rácz, 2015]. (Unsignedness is indicated by a
hat; the signed versions are $S_{n,k}(1,-1;-r)\eqdef
L_{n,k}^{(r)}$.)  Notably, $\hat L_{n,k}^{(2)}= \hat
L_{n+1,k+1}^{(0)}=\hat L_{n+1,k+1}$.  The triangles $\hat
L_{n,k}^{(r)}$ with $r=1,2,3,4$ appear separately as
A021009 (signed), A105278,
A062139 (signed), A143497.

* The $S_{n,k}(1,2;r) \eqdef B_{n,k}^{(r)}$ in (iii$^\prime$) are
the $r$-Bessel numbers of the second kind, which count
the partitions of an $(n+\nobreak r)$-set into
$k+\nobreak r$ blocks of size $1$ or $2$, such that
$r$ distinguished elements of the set are placed in distinct
blocks [Cheon et al., 2013, Jung et al., 2018].  (This implies that when
$r\in\mathbb{N}$, $S_{n,k}(1,2;r)$ is nonzero if and only if $0\le
n-\nobreak k \le \left\lfloor\frac{n+r}2\right\rfloor$.)  Notably,
$B_{n,k}^{(1)}= B_{n+1,k+1}^{(0)}=B_{n+1,k+1}$.  The triangles
$B_{n,k}^{(r)}$ with $r=1,2$ appear separately as
A049403, A330209.

* The $S_{n,k}(-2,-1;r) \eqdef \hat b_{n,k}^{(r)}$
in (iv$^\prime$) are the unsigned $r$-Bessel numbers of
the first kind, which count certain so-called
G-permutations of a multiset depending on $n,k$
and $r$ [Cheon et al., 2013].  More simply, $\hat b_{n,k}^{(0)}= \hat
b_{n,k}$ counts $n$-vertex rooted plane forests composed
of $k$ rooted plane trees, the vertices being labeled by
$1,\dots,n$, such that each of the $k$ trees is “increasing” in
that the labels on its vertices increase as one moves away from its
root [Bergeron et al., 1992].  Notably, $\hat b_{n,k}^{(1)}= \hat
b_{n+1,k+1}^{(0)}=\hat b_{n+1,k+1}$.  The triangle $\hat
b_{n,k}^{(r)}$ with $r=1$ appears separately as

Additional expressions are given in the following theorem and

For all $r\in\mathbb{C}$,

(v) $S_{n,k}(-1,2;\,r) = 2^{-(n-k)}\, \binom{n}{k}\,
{}_2{\hat F}_1\Bigl[
\begin{smallarray}{c}-(n-k),\:-n-r+1 \\ -n+2k+1\end{smallarray}
\Bigm|2\,
\Bigr]$,

(vi) $S_{n,k}(-2,1;\,r) = 2^{-(n-k)}\, \binom{n}{k}\,
{}_2{\hat F}_1\Bigl[
\begin{smallarray}{c}-(n-k),\:-2n-r+1 \\ -2n+k\end{smallarray}
\Bigm|2\,
\Bigr]$.

Each of (v),(vi), and in fact (iii$^\prime$),(iv$^\prime$), can be
proved by showing that it yields the correct vertical univariate
EGF, which was found in the proof of Theorem <ref>
to be

\begin{equation}
\sum_{n=k}^\infty \hat S_{n,k}(\alpha,\beta;\,r)
\frac{z^n}{n!}
=(1+\alpha z)^{r/\alpha} \left[ (1+\alpha z)^{\beta/\alpha} - 1 \right]^k,
\end{equation}

where $\hat S_{n,k} = \beta^kk!\,S_{n,k}$.  In doing this, a useful
lemma is the Srivastava–Singhal generating function for Jacobi
polynomials [Srivastava and Singhal, 1973],

\begin{equation}
\sum_{m=0}^\infty
P_m^{(\alpha-\lambda m,\beta-\mu m)}(x) \,z^m =
\frac{(1+\xi)^{\alpha+1}\,(1+\eta)^{\beta+1}}{1+\lambda\xi + \mu\eta - (1-\lambda-\mu)\xi\eta}\,,
\end{equation}

in which the degree-$m$ Jacobi polynomial $P_m^{(\alpha,\beta)}(x)$
is defined as usual by

\begin{equation}
{}_2{\hat F}_1
\Bigl[
\begin{smallarray}{c}-m,\:m+\alpha+\beta+1 \\ \alpha+1\end{smallarray}
\Bigm|\frac{1-x}2\,
\Bigr],
\end{equation}

and $\xi,\eta$ are defined implicitly as functions of (or strictly,
as formal power series in) $z$ by

\begin{align}
\xi &= \tfrac12 (x+1)z\, (1+\xi)^{1-\lambda}(1+\eta)^{1-\mu},\\
\eta &= \tfrac12 (x-1)z\, (1+\xi)^{1-\lambda}(1+\eta)^{1-\mu}.
\end{align}

One sets $m=n-k$ and $x=-3$.  To prove (iii$^\prime$),(iv$^\prime$)
and (v),(vi), one lets the tuple
$(\lambda,\mu;\allowbreak\alpha,\beta)$ equal $(1,0;\allowbreak
k,-k-\nobreak r-\nobreak1)$, $(2,-1; \allowbreak
-k-\nobreak1,\allowbreak k+\nobreak r-\nobreak1)$ and
$(1,1;\allowbreak k,-2k-\nobreak r)$, $(2,1; -k-\nobreak
1,-k-\nobreak r+\nobreak1)$, respectively.

The following formula for $\hat S_{n,k}(-1,2;r)=2^{k}k!\,
S_{n,k}(-1,2;r)$ was found heuristically and may merit further study.

For all $r\in\mathbb{Z}$,

n-jn-k n!n+12j+r-1.

This superficially resembles the binomial transform
formula (<ref>) of Theorem <ref>, but the range of the
summation index $j$ is not $0\le j\le k$, and the upper argument
$n-\nobreak j$ of the first binomial coefficient may be negative.  So
it is not easy to extract the generalized Eulerian numbers
$E_{n,k}(-1,2;r)$ corresponding to the $\hat S_{n,k}(-1,2;r)$ from
this conjectured formula.

The generalized Stirling numbers $S_{n,k}(-1,2;0)$ and
$S_{n,k}(-2,1;0)$ obtained by setting $r=0$ in
Theorem <ref> could be called numbers of the second
kind, and unsigned ones of the first kind, belonging to an as yet
unnamed combinatorial family.  They appear as A046089
and A035342 in the OEIS, where they are interpreted as
counting $n$-vertex rooted forests composed of $k$ rooted
trees, suitably restricted and labeled.  The more general versions
$S_{n,k}(-1,2;r)$ and $S_{n,k}(-2,1;r)$ seem not to have distinctive
interpretations.  $S_{n,k}(-1,2;1)$ and $S_{n,k}(-2,1;1)$ appear
separately as A321966 and A265649, and
notably, one has $S_{n,k}(-1,2;3)=S_{n+1,k+1}(-1,2;0)$ and

It has been observed (see A321966 in the OEIS) that for
all $n\ge0$,

\begin{align}
\label{eq:Hermiteidentity}
& z^n\Hermiteprob_n(z) = \sum_{k=0}^n S_{n,k}(-1,2;\,1)\Hermiteprob_{2k}(z),\\
& \Hermiteprob_{2n}(z) = \sum_{k=0}^n (-1)^{n-k}\,S_{n,k}(-2,1;\,1)\,z^k \Hermiteprob_{k}(z),
\label{eq:Hermiteidentityinverse}
\end{align}

the two identities being inverse to each other, as the inverse of the
matrix $\mathcal{S}(-1,2;1)$ is $\mathcal{S}(2,-1;-1)$, and
$S_{n,k}(2,-1;-1)$ equals $(-1)^{n-k} S_{n,k}(-2,1;1)$.  Here
$\Hermiteprob_n$, $n\ge0$, are the probabilists' version of the
physicists' Hermite polynomials $H_n$, i.e.,
$\Hermiteprob_n(z)=2^{-n/2}\,H_n(z/\sqrt2) =(-1)^n{\rm e}^{z^2/2}
D_z^n {\rm e}^{-z^2/2}$.

There may be simpler expressions for the connection coefficients in
these identities than are provided by formulas (v) and (vi) of
Theorem <ref>.  The $S_{n,k}(-1,2;1)$
in (<ref>) could be computed instead by

\begin{equation}
\label{eq:lesscomb}
S_{n,k}(-1,2;\,1) = 2^{-k}\sum_{j=0}^k (-1)^{k-j}\,
\frac{(2j+1)^{\overline n}}{j!\,(k-j)!},
\end{equation}

which is an application of the summation formula (<ref>) of
Theorem <ref>.  Also, one should notice that
$\mathcal{S}(-1,2;1) = \mathcal{S}(-1,1;1)\,\mathcal{S}(1,2;0)$, which
is an instance of the matrix product formula (<ref>).  It
implies that

\begin{equation}
S_{n,k}(-1,2;\,1) = \sum_{j=k}^n \hat L_{n,j}^{(1)}\,B_{j,k}
\frac{(n!)^2}{k!} \sum_{j=k}^n \frac{2^{-(j-k)}}{j!\, (n-j)!}\,\binom{k}{j-k}.
\label{eq:seemssimpler}
\end{equation}

Equation (<ref>) is more combinatorial
than (<ref>) in the sense that none of its summands is
negative, and it generalizes to

\begin{equation}
\label{eq:callrewrite}
S_{n,k}(-1,2;\,r) = \sum_{j=k}^n \hat L_{n,j}^{(r)}\,B_{j,k}
\frac{n!\,(n+r-1)!}{k!} \sum_{j=k}^n \frac{2^{-(j-k)}}{(j+r-1)!\, (n-j)!}\,\binom{k}{j-k},
\end{equation}

where $r$ is any positive integer.  But
rewriting (<ref>) in hypergeometric terms and extending
it to all $r\in\mathbb{C}$ would yield an expression comparable to the
${}_2{\hat F}_1$-based one in part (v) of
Theorem <ref>.  The numbers $S_{n,k}(-2,1;1)$ in the
inverse identity (<ref>) will reappear at
the end of the next section.

As regards combinatorial interpretations, it must be said that when
$\alpha,\beta$ and $r$ are integers, many quite broadly applicable
interpretations of the numbers $S_{n,k}(\alpha,\beta;r)$ are known,
such as ones involving the placement of balls in cells and
compartments, or the drawing of balls from an urn.  (Some are reviewed
in <cit.>.)  The ones mentioned in this
section are specific to the choices of $\alpha,\beta,r$ for which a
closed-form expression for $S_{n,k}(\alpha,\beta;r)$ is available,
terminating hypergeometric series being regarded as closed-form

§ ORDERING IDENTITIES (GENERALIZED STIRLING KIND)

Operator ordering identities involving the generalized Stirling
numbers of Hsu and Shiue [Hsu and Shiue, 1998] will now be derived.  They are of
two types: (1) the normal ordering of a power $\Omega^n$ of any
single-annihilator word in $\adag,a$, and (2) the expansion
of $\Omega^n$ in powers $\Omega'{}^k$, $k=0,\dots,n$, of any other
such word.  Normal orderings are subsumed by expansions, which are
more general.  The respective main theorems are
Theorem <ref>, which can be viewed as a special case of a
normal ordering result of Méndez et al. [Méndez et al., 2005], and
Theorem <ref>, which is new.

The basic technique used here is a long known and used
one [Kim and Kim, 2022, Mikhailov, 1985].  It is the transformation of an
identity in (commutative) polynomial algebra over $\mathbb{C}$ to one
in operator algebra by replacing an indeterminate $x$ by $xD$.  What
results is an identity in $\tilde W_{\mathbb{C}}$, but if the only
exponents appearing in it are elements of $\mathbb{N}$, it can also be
viewed as an identity in $W_{\mathbb{C}}$.

For all $L,R\in\mathbb{C}$, one has a normal ordering identity in $\tilde

\begin{equation}
\label{eq:normord}
\left(x^{-e}\right)^n \left(x^L D x^R\right)^n
= \sum_{k=0}^n S_{n,k}\left(-e,1;\,R\right)
\end{equation}

where $e=L+R-1$.

By the definition (<ref>) of the Hsu–Shiue numbers
$S_{n,k}(\alpha,\beta;r)$, with $x$ replaced by $xD+r$ and $\beta$
set equal to $1$,

\begin{equation}
(xD+r)^{\underline n,\alpha} = \sum_{k=0}^n
S_{n,k}(\alpha,1;\,r)\,(xD)^{\underline k}.
\end{equation}

Applying Boole's identity $(xD)^{\underline k}=x^kD^k$ and the
generalization (<ref>) given in Lemma <ref>
converts this to

\begin{equation}
= \sum_{k=0}^n
\end{equation}

and letting $\alpha=1-L-R$, $r=R$ yields (<ref>).

For all $L,R\in\mathbb{N}$ with $L+R\ge1$, one has a normal ordering
identity in $W_{\mathbb{C}}$:

\begin{equation}
\begin{aligned}
\label{eq:untilnow}
\left(\adag^L a \adag^R\right)^n &= \left(\adag^{e}\right)^n \,\sum_{k=0}^n
S_{n,k}\left(-e,1;\,R\right) \adag^k a^k \\
&= \left(\adag^e\right)^n
\,\sum_{k=0}^n
\left[
\sum_{j=0}^k (-1)^{k-j} \,\frac{(R+j)^{\overline n,e}}{j!(k-j)!}
\right] \adag^k a^k,
\end{aligned}
\end{equation}

where $e=L+R-1\ge0$ is the excess of the word $\Omega={a^\dag}^L a

Left-multiply both sides of (<ref>) by
$\left(x^e\right)^n=\left(x^{L+R-1}\right)^n$.  If all powers of $x$
and $D$ are nonnegative integers, which is the case if
$L,R\in\mathbb{N}$ with $L+\nobreak R\ge\nobreak 1$, substituting
$\adag,a$ for $x,D$ yields an identity in $W_\mathbb{C}$.  The formula
for $S_{n,k}(-e,1;r)$ is taken from Theorem <ref>, it
being recalled that $\hat S_{n,k}(\alpha,\beta;r) = \beta^k k!
\,S_{n,k}(\alpha,\beta;r)$.

This is a special case of the abovementioned result of Méndez
et al. [Méndez et al., 2005], who normally-ordered an arbitrary “boson
string” of the form $\adag^{r_m}a^{s_m}\dots \adag^{r_1}a^{s_1}$.
(Also see <cit.>.)  Typical applications
include the following.  Setting $(L,R)=(1,0)$ and $(0,1)$
in (<ref>) yields the ordering identities shown in
eq. (<ref>) of the introduction, which readily generalize
with the aid of Lemma <ref> to

\begin{equation}
\label{eq:katrielplus}
(\adag a)^{\underline n,\alpha} = \sum_{k=0}^n S_{n,k}(\alpha,1;0)\, \adag^k (a)^k,
\qquad
(a \adag)^{\underline n,\alpha} = \sum_{k=0}^n S_{n,k}(\alpha,1;1)\, \adag^k (a)^k,
\end{equation}

for all $\alpha\in\mathbb{C}$.  Another application involves the
unsigned $r$-Lah numbers $\hat
L_{n,k}^{(r)}=S_{n,k}(-1,1;r)$, which by formula (ii${}^\prime$) of
Theorem <ref> equal $\binom{n}k (k+r)^{\overline{n-k}}$.
Setting $(L,R)$ successively equal to $(0,2)$, $(1,1)$, and $(2,0)$,
for each of which $e=\nobreak 1$, yields the triple of normal

\begin{align}
\left(\adag^0 a \adag^2\right)^n &= \adag^n\:
\sum_{k=0}^n \binom{n}k (k+2)^{\overline{n-k}}\,\adag^k a^k,\\
\left(\adag^1 a \adag^1\right)^n &= \adag^n\:
\sum_{k=0}^n \binom{n}k (k+1)^{\overline{n-k}}\,\adag^k a^k,\\
\left(\adag^2 a \adag^0\right)^n &= \adag^n\:
\sum_{k=0}^n \binom{n}k (k+0)^{\overline{n-k}}\,\adag^k a^k.
\end{align}

An additional application is this.

For all $L,R\in\mathbb{N}$ with $L+R\ge1$, one has a special
normal-ordering identity in $W_{\mathbb{C}}$:

\begin{equation}
\begin{split}
&\left(\adag^L a \adag^{R} + \adag^R a \adag^L\right)^n\\
&\qquad =\left(\adag^{(L+R-1)}\right)^n \,\sum_{k=0}^n \,2^k S_{n,k}\left(2-2L-2R,2;\,L+R\right)
\adag^k a^k.
\end{split}
\end{equation}

If $L+R$ is even, the left-hand side equals $2(\adag)^{(L+R)/2} a
(\adag)^{(L+R)/2}$ by the identity (<ref>)
in $W_{\mathbb{C}}$, and the claim follows from
Theorem <ref>, if one also uses the homogeneity property
$S_{n,k}(\alpha,\beta;r) = 2^{-(n-k)} S_{n,k}(2\alpha,2\beta;2r)$.
(See eqs. (<ref>).)  If $L+R$ is odd, one must go
back to Theorem <ref>, and apply the counterpart
identity (<ref>) in $\tilde W_{\mathbb{C}}$.

The following subsumes Theorem <ref>, because its first
identity, eq. (<ref>), reduces to
Theorem <ref> when $(L',R')=(0,0)$.

For all $L,R,L',R'\in\mathbb{C}$, one has four Stirling-type expansion identities
in $\tilde W_{\mathbb{C}}^{\rm{bal}}$:

\begin{align}
\left(x^{-e}\right)^n \left(x^L D x^R\right)^n
&= \sum_{k=0}^n S_{n,k}\left(-e,-e';\,R-R'\right)\left(x^{-e'}\right)^k \left(x^{L'}Dx^{R'}\right)^k   \label{eq:firstmain2a}\\[\jot]
&= \sum_{k=0}^n S_{n,k}\left(-e,e';\,R+L'-1\right)\left(x^{L'}Dx^{R'}\right)^k\left(x^{-e'}\right)^k,\label{eq:firstmain2b}
\\
\shortintertext{and}
\stepcounter{parentequation}\gdef\theparentequation{\arabic{parentequation}}\setcounter{equation}{0}
\left(x^L D x^R\right)^n\left(x^{-e}\right)^n
&= \sum_{k=0}^n S_{n,k}\left(e,-e';\,1-L-R'\right)\left(x^{-e'}\right)^k \left(x^{L'}Dx^{R'}\right)^k   \label{eq:firstmain2c}\\[\jot]
&= \sum_{k=0}^n S_{n,k}\left(e,e';\,L'-L\right)\left(x^{L'}Dx^{R'}\right)^k\left(x^{-e'}\right)^k,\label{eq:firstmain2d}
\end{align}

where $e=L+R-1$ and $e'=L'+R'-1$.

Up to parameter renamings, the identities
(<ref>),(<ref>) are related by the
formal adjoint operation $*$ on $\tilde W_{\mathbb{C}}$, which
preserves $\tilde W_{\mathbb{C}}^{\rm bal}$; and so are

The identity (<ref>) of Theorem <ref> is
based on the matrix $\mathcal{S}(1-\nobreak L-\nobreak R,1;R)$, the
inverse of which is $\mathcal{S}\bigl(1,1-\nobreak L-\nobreak
R;-R\bigr)$.  So the inverse identity

\begin{equation}
\label{eq:inversenormord}
x^nD^n = \sum_{k=0}^n S_{n,k}\left(1,1-{L'}-{R'};\,\,{-R'}\right)
\left(x^{1-{L'}-{R'}}\right)^k \left(x^{L'} D x^{R'}\right)^k
\end{equation}

also holds.  Composing (<ref>)
with (<ref>) yields (<ref>),
because by the matrix product formula (<ref>), the
product of $\mathcal{S}(1-\nobreak L-\nobreak R,1;R)$ and
$\mathcal{S}\bigl(1,\allowbreak 1-\nobreak L'-\nobreak R';-R'\bigr)$
is $\mathcal{S}(1-\nobreak L-\nobreak R,\allowbreak 1-\nobreak
L'-\nobreak R';\allowbreak R-\nobreak R')$.  To obtain
(<ref>) from (<ref>), apply
(<ref>) to the right-hand side and perform the
renamings $L'\leftarrow 1-\nobreak R'$, $R'\leftarrow 1-\nobreak
L'$.  To obtain (<ref>) from (<ref>),
apply (<ref>) to the left-hand side and perform the
renamings $L\leftarrow 1-\nobreak R$, $R\leftarrow 1-\nobreak L$.
To obtain (<ref>) from (<ref>), do

The following four-identity theorem subsumes Theorem <ref>,
as the $e=\allowbreak L+\nobreak R-\nobreak 1\ge0$ case of its first
identity reduces to Theorem <ref> if one sets
$(L',R')=\allowbreak (0,0)$, so that $\Omega'=\nobreak a$ and
$e'=\nobreak -1$.  That is, it produces a normal ordering of its
left-hand side.  The fourth identity similarly produces an anti-normal
ordering of its left-hand side.

For all $L,R,L',R'\in\mathbb{N}$, one has four Stirling-type expansion identities in $W_{\mathbb{C}}$:

\begin{align}
&\left({a^\dag}^{(E_L-e)}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n \label{eq:main1a}\\
&{}=\sum_{k=0}^n S_{n,k}\left(-e,-e';\,R-R'\right)
\left({a^\dag}^{E_L}\right)^{n-k}  \left({a^\dag}^{({E_L}-e')}\right)^k  \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k   ,\nonumber\\[\jot]
&\left({a^\dag}^{(E_L-e)}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n   \left(\adag^{E_R}\right)^n\label{eq:main1b}\\
&{}=\sum_{k=0}^n S_{n,k}\left(-e,e';\,R+L'-1\right)
\left({a^\dag}^{E_L}\right)^{n}   \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k    \left({a^\dag}^{(E_R-e')}\right)^k  \left(\adag^{E_R}\right)^{n-k}  ,\nonumber
\\
\shortintertext{and}
\stepcounter{parentequation}\gdef\theparentequation{\arabic{parentequation}}\setcounter{equation}{0}
&\left({a^\dag}^{E_L}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n \left({a^\dag}^{(E_R-e)}\right)^n\label{eq:main2a}\\
&{}=\sum_{k=0}^n S_{n,k} \left(e,-e';\,1-L-R'\right)
\left({a^\dag}^{E_L}\right)^{n-k}  \left({a^\dag}^{(E_L-e')}\right)^k   \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k   \left(\adag^{E_R}\right)^n ,\nonumber\\[\jot]
&\left({a^\dag}^L a {a^\dag}^R\right)^n \left({a^\dag}^{(E_R-e)}\right)^n \label{eq:main2b}\\
&{}=\sum_{k=0}^n S_{n,k}\left(e,e';\,L'-L\right)
\left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k  \left(\adag^{(E_R-e')}\right)^k   \left(\adag^{E_R}\right)^{n-k} ,\nonumber
\end{align}

where $e=L+R-1$ and $e'=L'+R'-1$ are the excesses of the words
$\Omega={a^\dag}^L a {\adag}^R$ and $\Omega'={a^\dag}^{L'} a
{\adag}^{R'}$, and the quantities $E_L=E_L(e,e')$, $E_R=E_R(e,e')$
are as follows:

Identity     $\hfil E_L$     $\hfil E_R$

(<ref>)     $\max(e,e',0)$     $\hfil \text{---}\hfil$

(<ref>)     $\max(e,0)$     $\max(e',0)$

(<ref>)     $\max(e',0)$     $\max(e,0)$

(<ref>)     $\hfil \text{---}\hfil$     $\max(e,e',0)$.

Restrict each of (<ref>), (<ref>),
(<ref>), (<ref>) by requiring
$L,R,\allowbreak L',R'\in\nobreak\mathbb{N}$, and also left-multiply
by $(x^{E_L})^n$ and right-multiply by $(x^{E_R})^n$, where
$E_L,E_R$ are chosen to make all powers of $x$ non-negative
integers; and then substitute $\adag, a$ for $x,D$.

The expansion identities of Theorem <ref>, being
parametric, yield many individual operator ordering identities
in $W_{\mathbb{C}}$.  Their coefficients $S_{n,k}$ can be computed
from the summation formula (<ref>) in
Theorem <ref> or the triangular
recurrence (<ref>) in Theorem <ref>.  But
by using the closed-form expressions for certain $S_{n,k}$ derived in
Section <ref>, one can make these identities more
explicit, perhaps facilitating new combinatorial interpretations.  The
following discussion focuses on the first and second identities.

It should be noted that the first identity, eq. (<ref>),
splits into two or three cases.  Suppose that $e,e'\ge0$ (which
incidentally rules out the normal ordering case, when $e'=-1$).  Then
if $e\ge e'$, so that $E_L=e$, (<ref>) says that

\begin{equation}
\left({a^\dag}^L a {a^\dag}^R\right)^n
\sum_{k=0}^n S_{n,k}\left(-e,-e';\,R-R'\right)\,
\adag^{(en-e'k)} \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k,
\end{equation}

and if $e\le e'$, so that $E_L=e'$, it says that

\begin{equation}
{a^\dag}^{(e'-e)n} \left({a^\dag}^L a {a^\dag}^R\right)^n
\sum_{k=0}^n S_{n,k}\left(-e,-e';\,R-R'\right)
{a^\dag}^{e'(n-k)}   \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k.
\end{equation}

In both cases,

\begin{equation}
S_{n,k}(-e,-e';\,R-R') =
\sum_{j=0}^k (-1)^{k-j}\, \frac{(R-R'-e'j)^{\overline n,e}}{j!(k-j)!},
\end{equation}

which comes from Theorem <ref>.  If $e=e'\ge0$, i.e.,
$R+L=R'+L'\ge1$, (<ref>) becomes an expansion identity with
quite simple coefficients,

\begin{equation}
\left({a^\dag}^L a {a^\dag}^R\right)^n
\sum_{k=0}^n \binom{n}k (R-R')^{\overline{n-k},e}\:
\adag^{e(n-k)} \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^k,
\end{equation}

where $e=R+L-1=R'+L'-1$.  The coefficients come from
formula (i$^\prime$) of Theorem <ref>.

Explicit examples of the expansion of a power of a word $\Omega$ in
powers of another word $\Omega'$ appeared in eqs. (<ref>)
and (<ref>) of the introduction.  The former is the case
$(L,R;\allowbreak L',R') = (1,0;\allowbreak0,0)$ of the
identity (<ref>), and the latter is the case
$(L,R;\allowbreak L',R') = (3,0;\allowbreak2,0)$.  The closed-form
expression for $\hat b_{n,k}=S_{n,k}(-2,-1;0)$ appearing
in (<ref>) comes not from Theorem <ref> but
from formula (iv) of Theorem <ref>, having been
rewritten in terms of ordinary factorials.  The identity inverse
to (<ref>) is the case $(L,R;\allowbreak
L',R')=(2,0;\allowbreak 3,0)$ of (<ref>), which is

\begin{equation}
\begin{aligned}
\adag^n \left( {a^\dag}^2 a\right)^n
&=  \sum_{k=0}^n  S_{n,k}(-1,-2;0)\: {a^\dag}^{2(n-k)} \left({a^\dag}^3 a\right)^k\\
&=  \sum_{k=0}^n  (-2)^{-(n-k)}\, \frac{n!}{k!}\,\binom{k}{n-k}  \: {a^\dag}^{2(n-k)} \left({a^\dag}^3 a\right)^k.
\end{aligned}
\label{eq:viewedas}
\end{equation}

Here the coefficient $S_{n,k}(-1,-2;0)$ equals
$(-1)^{n-k}S_{n,k}(1,2;0)$ by homogeneity, and the expression for
$B_{n,k}=S_{n,k}(1,2;0)$ is taken from formula (iii) of
Theorem <ref>.  The only nonvanishing terms in the sum
are those with $0\le n-\nobreak k\le
\left\lfloor{\frac{n}2}\right\rfloor$.

Equation (<ref>) has a companion identity, which is the
case $(L,R;\allowbreak L',R')=\allowbreak (2,1;2,0)$ of
eq. (<ref>).  It is

\begin{equation}
\begin{aligned}
\left( {a^\dag}^2 a\adag\right)^n &= \adag^n\: \sum_{k=0}^n S_{n,k}(-2,-1;1)\:
{a^\dag}^{(n-k)} \left({a^\dag}^2 a\right)^k\\ &= \adag^n\:\sum_{k=0}^n
2^{-(n-k)} \frac{(2n-k)!}{k!(n-k)!} \: {a^\dag}^{(n-k)}
\left({a^\dag}^2 a\right)^k.
\end{aligned}
\end{equation}

The expression used here for the numbers $\hat
b_{n,k}^{(1)}=S_{n,k}(-2,-1;1)$, which are the unsigned
$1$-Bessel numbers of the first kind, comes from
formula (iv$^\prime$) of Theorem <ref>.  When $r=1$, the
desingularized hypergeometric series ${}_2\hat F_1(2)$ on which
formula (iv$^\prime$) is based has zero as an upper parameter.  It
therefore terminates with its first term, which is $(-2n+\nobreak
k)^{\overline{n-k}}$; whence the factor $(2n-\nobreak k)!$.

An interesting example of normal ordering is the following.

One has the normal-ordering identity

(^3 a)^n =
^2n  ∑_k=0^n
{2^-(n-k)  nk
c-(n-k), -2n+1
^k a^k

in $W_{\mathbb{C}}$, with desingularized hypergeometric series as

This is the case $(L,R;L',R')=(3,0;0,0)$ of eq. (<ref>),
or equivalently the case $(L,R)=(3,0)$ of Theorem <ref>,
with the bracketed expression for $S_{n,k}(-2,1;0)$ coming from
formula (vi) of Theorem <ref>.

A seemingly simpler expression for the bracketed coefficient
$S_{n,k}(-2,1;0)$ in this proposition, based on an explicit summation,
is provided by Theorem <ref>: it is

\begin{equation}
\label{eq:lesscomb2}
S_{n,k}(-2,1;\,0) = \sum_{j=0}^k(-1)^{k-j} \frac{(j)^{\overline n,2}}{j!(k-j)!}.
\end{equation}

A third summation formula for $S_{n,k}(-2,1;0)$, which is more
combinatorial than (<ref>), follows from

$\mathcal{S}(-2,1;0) = \mathcal{S}(-2,-1;0) \,
\mathcal{S}(-1,1;0)$.  It is

\begin{equation}
\begin{aligned}
\sum_{j=k}^n \hat b_{n,j}\hat L_{j,k}
\sum_{j=k}^n 2^{-(n-j)}\, \frac{(n)^{\overline{n-j}}(j)^{\overline{n-j}}}{(n-j)!}
\cdot
\binom{j}k (k)^{\overline{j-k}}\\
&= (k)^{\overline{n-k}}\,
\sum_{j=k}^n
2^{-(n-j)}\, \frac{(n)^{\overline{n-j}}}{(n-j)!} \binom{j}k,
\end{aligned}
\end{equation}

because $S_{n,k}(-2,-1;0)=\hat b_{n,k}$ and $S_{n,k}(-1,1;0)=\hat
L_{n,k}$.  Which formula is the most efficient for computing
$S_{n,k}(-2,1;0)$ depends on the (relative) sizes of $k$, $n-\nobreak
k$, and $n$.

A sample application of the second expansion identity in
Theorem <ref>, eq. (<ref>), is the following.
As was noted in Section <ref>, the numbers
$S_{n,k}(-2,1;1)$, which are A265649 in the
OEIS [OEIS Foundation, Inc., 2023], appear as coefficients in an identity satisfied
by the Hermite polynomials $\Hermiteprob_n(z)$.  Setting
$(L,R;\allowbreak L',R')$ in eq. (<ref>) successively equal
to $(1,2;\allowbreak 0,2)$, $(2,1;\allowbreak1,1)$, and
$(3,0;\allowbreak 2,0)$, for each of which $E_L=e=\nobreak2$ and
$E_R=e'=\nobreak1$, yields the triple of expansions

\begin{align}
%    \left(\adag^L a \adag^R\right)^n \left(\adag\right)^n &= \sum_{k=0}^n S_{n,k}(-2,1;\,1)\, \left(\adag^2\right)^n \left(\adag^{L'} a \adag^{R'}\right)^k \left(\adag\right)^{n-k}\\
\left(\adag^1 a \adag^2\right)^n \adag^n &= \adag^{2n} \:\sum_{k=0}^n S_{n,k}(-2,1;\,1)\, \left(\adag^{0} a \adag^{2}\right)^k \adag^{(n-k)},\\
\left(\adag^2 a \adag^1\right)^n \adag^n &= \adag^{2n} \:\sum_{k=0}^n S_{n,k}(-2,1;\,1)\, \left(\adag^{1} a \adag^{1}\right)^k \adag^{(n-k)},\\
\left(\adag^3 a \adag^0\right)^n \adag^n &= \adag^{2n} \:\sum_{k=0}^n S_{n,k}(-2,1;\,1)\, \left(\adag^{2} a \adag^{0}\right)^k \adag^{(n-k)},
\end{align}

where the $S_{n,k}(-2,1;1)$ appear as coefficients.  The connection to
the Hermite polynomial identity (<ref>), if
any, is unclear.

§ EULERIAN-NUMBER FORMULAS AND COMBINATORICS

As with the numbers $S_{n,k}(\alpha,\beta;r)$, it is desirable to be
able to compute the generalized Eulerian numbers
$E_{n,k}(\alpha,\beta;r)$ efficiently, and to interpret them
combinatorially if possible.  They can be computed from the summation
formula (<ref>) of Theorem <ref>, or the
triangular recurrence (<ref>) of
Theorem <ref>.  They can also be computed with the aid
of Theorem <ref>, according to which each row of the matrix
${\mathcal{E}}(\alpha,\beta;r)=[E_{n,k}(\alpha,\beta;r)]$ is the
(inverse) binomial transform of the corresponding row of
$\hat{\mathcal{S}}(\alpha,\beta;r)=[\hat S_{n,k}(\alpha,\beta;r)]$.

Similarly to the $S_{n,k}(\alpha,\beta;r)$, closed-form expressions
are available for certain choices of $\alpha,\beta$ and $r$.  For
deriving additional closed-form expressions one can use homogeneity:
$E_{n,k}(\lambda \alpha,\lambda \beta;\lambda r)$ equals
$\lambda^{n}E_{n,k}(\alpha,\beta;r)$.  The reflection
property (<ref>), according to which
$E_{n,k}(\alpha,\beta;r)$ equals
$E_{n,n-k}(-\alpha,\beta;\beta-\nobreak r)$, can also be used.

For all $r\in\mathbb{C}$,

(i) $E_{n,k}(\beta,\beta;\,r) = \binom{n}k
(r)^{\underline{n-k},\beta}(\beta n-r)^{\underline{k},\beta}$,

(ii)  $E_{n,k}(-\beta,\beta;\,r) = \binom{n}k (\beta k+r)^{\overline{n-k},\beta}(\beta-r)^{\underline{k},\beta}$.

Additionally, one has the formulas

(iii) $E_{n,k}(-1,2;\,0) =
\begin{cases}
\delta_{n,0}, & $k=0$,\\
n!\,\binom{n+1}{2k-1}, & $k=1,\dots,n$,
\end{cases}$

(iv) $E_{n,k}(-1,2;\,1) = n!\,\binom{n+1}{2k}$,

(v) $E_{n,k}(-1,2;\,2) = n!\,\binom{n+1}{2k+1}$.

Moreover, $E_{n,k}(-1,2;r)$ for any $r\in\mathbb{N}$ can be computed from the following.

For all $p\in\mathbb{N}$ and $\zeta\in\{0,1\}$,

\begin{align*}
& {}=n!\,\binom{n+1}{2k+2p+1-\zeta}
- (-1)^{k+p} \sum_{\ell=0}^{p-1}\, (-1)^{\ell} (2-\zeta+2\ell)^{\overline n}
\binom{n+1}{k+p-\ell}.
\end{align*}

Formulas (i)–(v) are verified by confirming that they satisfy the
recurrence (<ref>).  Also, by the reflection property
(i),(ii) are equivalent.

The cases $\zeta=1,0$ of formula (vi) are proved separately, by
induction on $p$.  Their respective base subcases (with $p=0$) are
formulas (iv) and (v), and the inductive step in each proof
(an incrementing of the third parameter by $\beta=2$) uses the second
difference equation in eq. (<ref>), which is a formula for

For computing not individual numbers $E_{n,k}$ but entire row
polynomials, the following simple formula may be of use.  The $n$'th
row polynomial is

\begin{equation}
E_n(t)= \sum_{k=0}^n E_{n,k}(\alpha,\beta;\,r)\,t^k =(1-t)^{n+1}
\sum_{j=0}^\infty (\beta j+r)^{\underline n,\alpha} \,t^j.
\end{equation}

(It is not immediately obvious that the last expression, viewed as a
formal power series in $t$, is in fact a polynomial in $t$.)  This
follows from the summation formula (<ref>), or alternatively
the EGF (<ref>).  When $(\alpha,\beta;r)=(0,1;0)$, it
specializes to a known formula for the classical Eulerian
polynomials <cit.>, and a version with $\alpha=0$ but
$\beta,r$ arbitrary was obtained previously [Mező and Ramírez, 2016].

For some choices of the parameters $\alpha,\beta,r$, combinatorial
interpretations of the numbers $E_{n,k}(\alpha,\beta;r)$ are known.
The standard interpretation (when $n>0$) of the classical numbers
$E_{n,k}(0,1;0)$ and $E_{n,k}(0,1;1)$ (i.e., $A_{n,k}$ and
$\eulerian{n}k = A_{n,k+1}$, identical except for a horizontal shift
by unity) has already been mentioned.  The former counts the number of
permutations of the ordered set $\{1,\dots,n\}$ that have exactly
$k$ descents, and the latter is an incremented version.

The numbers $E_{n,k}(0,2;1) \eqdef {\eulerian{n}{k}}_B$ are less well
known.  They are the MacMahon or type-$B$ Eulerian
numbers, and count signed permutations of the ordered set
$\{-n,\dots,-1\}\cup\{1,\dots,n\}$ with exactly $k$ “signed
descents.”  (See [Bagno et al., 2022, Brenti, 1994] and A060187
in the OEIS.)  They are related by the case
$(\alpha,\beta;r)=(0,2;1)$ of the generalized Worpitzky
identity (<ref>) to the type-$B$ Stirling numbers
$S_{n,k}(0,2;1)\eqdef {\stirsub{n}k}_B$, which can also be interpreted
combinatorially [Bagno et al., 2019].

The number triangles $E_{n,k}(0,3;2)$ and $E_{n,k}(0,4;3)$ can be
found in the OEIS as well, as A225117 and
A225118.  Like the classical and type-$B$
Eulerian numbers, they belong to the family of triangles of the form
$E_{n,k}(0,r+\nobreak 1;\nobreak r)$, $r\in\mathbb{N}$, each of which
has an enumerative interpretation [Wang and Zhu, 2023].  In fact, so does
any triangle $E_{n,k}(0,\beta;r)$ in which $\beta,r\in\mathbb{N}$ and
$r\le\nobreak\beta$ (see [Ramírez et al., 2018]).

An interesting interpretation of $E_{n,k}(-1,\beta;r)/n!$, when
$\beta$ is a positive integer and $r\in\mathbb{N}$ with $r\le\beta$,
should also be mentioned.  It counts the number of words
$s_1s_2,\dots,\nobreak s_n$ with
$s_i\in\{0,\dots,\nobreak\beta-\nobreak 1\}$ that have exactly $k$
descents, provided each word is prefixed with a hidden, “zeroth”
letter $s_0=\beta\nobreak-r$ (which affects the count of the number of
descents).  Closed-form expressions for $E_{n,k}(-1,2;r)/n!$ with
$r=0,1,2$ are supplied by parts (iii),(iv),(v) of
Theorem <ref> above.  The latter two triangles appear in
the OEIS as A228955 and A131980, and the
triangles $E_{n,k}(-1,3;r)/n!$ with $r=1$ and $r=3$ as
A178618 and A120906.

Additionally worth mentioning is the recent
interpretation [Herscovici, 2020] of the “degenerate” Eulerian
numbers of Carlitz [Carlitz, 1979], which in present notation would
be denoted by $E_{n,k}(\alpha,1;r)$, in the case when $r=1$.  But
in general, combinatorial interpretations of the
$E_{n,k}(\alpha,\beta;r)$ have yet to be developed.

§ ORDERING IDENTITIES (GENERALIZED EULERIAN KIND)

Operator ordering identities in $W_{\mathbb{C}}$ based on the
generalized Eulerian numbers $E_{n,k}$ will now be derived.  Each is
the expansion of the $n$'th power $\Omega^n$ of a single-annihilator
word $\Omega$ in twistings of the $n$'th power $\Omega{{}^\prime{}^n}$
of another such word.  More precisely, each provides an expansion of
$\left(\adag^p\right)^n\Omega^n\left(\adag^q\right)^n$, where $p,q$
are sufficiently large, as a combination of terms of the form
$\adag^{p'}\Omega'{}^n\adag^{q'}$, with each term having the same
value of $p'+\nobreak q'$, which is a certain multiple of $n$.

The main theorem is Theorem <ref>, and like its Stirling
counterpart Theorem <ref>, it arises from a simpler
version in $\tilde W^{\rm bal}_{\mathbb{C}}$ (here,
Theorem <ref>).  The version in $\tilde W^{\rm
bal}_{\mathbb{C}}$ comes in turn from an identity in (commutative)
polynomial algebra, by replacing an indeterminate $x$ by $xD$.  Here,
this is the generalized Worpitzky identity (<ref>), which in
the linear-algebraic approach of this paper serves as the definition
of the numbers $E_{n,k}$.

For all $L,R,L',R'\in\mathbb{C}$, one has four Eulerian-type expansion identities
in $\tilde W_{\mathbb{C}}^{\rm bal}$:

\begin{align}
\begin{aligned}
\left(x^{-e}\right)^n \left(x^L D x^R\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(-e,-e';\,R-R'\right)\left(x^{-e'}\right)^k \left(x^{L'}Dx^{R'}\right)^n \left(x^{-e'}\right)^{n-k},
\end{aligned}
\label{eq:secondmain2a}
\\[\jot]
\begin{aligned}
\left(x^{-e}\right)^n \left(x^L D x^R\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(-e,e';\,R+L'-1\right) \left(x^{-e'}\right)^{n-k}   \left(x^{L'}Dx^{R'}\right)^n \left(x^{-e'}\right)^k ,
\end{aligned}
\label{eq:secondmain2b}
\\
\shortintertext{and}
\stepcounter{parentequation}\gdef\theparentequation{\arabic{parentequation}}\setcounter{equation}{0}
\begin{aligned}
\left(x^L D x^R\right)^n\left(x^{-e}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(e,-e';\,1-L-R'\right)\left(x^{-e'}\right)^k \left(x^{L'}Dx^{R'}\right)^n  \left(x^{-e'}\right)^{n-k} ,
\end{aligned}
\label{eq:secondmain2c}
\\[\jot]
\begin{aligned}
\left(x^L D x^R\right)^n\left(x^{-e}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}
\left(e,e';\,L'-L\right) \left(x^{-e'}\right)^{n-k}   \left(x^{L'}Dx^{R'}\right)^n  \left(x^{-e'}\right)^{k} ,
\end{aligned}
\label{eq:secondmain2d}
\end{align}

where $e=L+R-1$ and $e'=L'+R'-1$.

Substitute $xD+r_1$ for $x$ in the identity (<ref>), with
$r=r_1-\nobreak r_2$.  This yields

\begin{equation}
n!\,\beta^n(xD+r_1)^{\underline n,\alpha} = \sum_{k=0}^n E_{n,k}(\alpha,\beta;\,r)\,\left[xD+\beta(n-k)+r_2\right]^{\underline n,\beta}.
\end{equation}

Applying (<ref>) to both sides converts this to

\begin{equation}
\begin{aligned}
&\qquad{}= \sum_{k=0}^n E_{n,k}(\alpha,\beta;\,r)
x^{-r_2}x^{-\beta(n-k)} \left[(x^\beta)^n (x^{1-\beta}D)^n\right]
x^{\beta(n-k)} x^{r_2} \\
&\qquad{}= \sum_{k=0}^n E_{n,k}(\alpha,\beta;\,r)
\left(x^{\beta}\right)^k\left[x^{1-\beta-r_2}Dx^{r_2}\right]^n (x^\beta)^{n-k},
\end{aligned}
\end{equation}

and (<ref>) follows by letting $L=1-\nobreak
\alpha-\nobreak r_1$, $R=r_1$, $L'=1-\nobreak \beta-\nobreak r_2$,
$R'=r_2$.  The expansions
(<ref>),(<ref>) are equivalent by
homogeneity and the reflection property (<ref>),
as one sees by replacing the summation index $k$ by $n-\nobreak k$.
The expansions (<ref>),(<ref>) are
too, and up to parameter renamings they are the respective adjoints
of (<ref>),(<ref>).  They also
follow from (<ref>),(<ref>) by
applying the self-adjoint identity (<ref>).

For all $L,R,L',R'\in\mathbb{N}$, one has four Eulerian-type expansion
identities in $W_{\mathbb{C}}$:

\begin{align}
\begin{aligned}
&n!\, (-e')^n \left({a^\dag}^{(E_L-e)}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n    \left(\adag^{E_R}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(-e,-e';\,R-R'\right) \\
&\qquad\qquad{}\times\left({a^\dag}^{E_L}\right)^{n-k} \left({a^\dag}^{({E_L}-e')}\right)^k   \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^n \left(\adag^{(E_R-e')}\right)^{n-k} \left(\adag^{E_R}\right)^k,
\end{aligned}
\label{eq:2main1a}
\\[\jot]
\begin{aligned}
&n!\, (e')^n \left({a^\dag}^{(E_L-e)}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n   \left(\adag^{E_R}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(-e,e';\,R+L'-1\right)\\
&\qquad\qquad{}\times\left({a^\dag}^{E_L}\right)^{k} \left(\adag^{(E_L-e')}\right)^{n-k}   \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^n     \left({a^\dag}^{(E_R-e')}\right)^k  \left(\adag^{E_R}\right)^{n-k},
\end{aligned}
\label{eq:2main1b}
\\
\shortintertext{and}
\stepcounter{parentequation}\gdef\theparentequation{\arabic{parentequation}}\setcounter{equation}{0}
\begin{aligned}
&n!\, (-e')^n \left({a^\dag}^{E_L}\right)^n \left({a^\dag}^L a {a^\dag}^R\right)^n \left({a^\dag}^{(E_R-e)}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(e,-e';\,1-L-R'\right) \\
&\qquad\qquad{}\times\left({a^\dag}^{E_L}\right)^{n-k} \left({a^\dag}^{(E_L-e')}\right)^k     \left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^n   \left(\adag^{(E_R-e')}\right)^{n-k} \left(\adag^{E_R}\right)^{k},
\end{aligned}
\label{eq:2main2a}
\\[\jot]
\begin{aligned}
&n!\, (e')^n \left(\adag^{E_L}\right)^n  \left({a^\dag}^L a {a^\dag}^R\right)^n \left({a^\dag}^{(E_R-e)}\right)^n \\
&\qquad{}=\sum_{k=0}^n E_{n,k}\left(e,e';\,L'-L\right)\\
&\qquad\qquad{}\times\left(\adag^{E_L}\right)^k  \left(\adag^{(E_L-e')}\right)^{n-k}
\left({a^\dag}^{L'}a{a^\dag}^{R'}\right)^n    \left(\adag^{(E_R-e')}\right)^k \left(\adag^{E_R}\right)^{n-k},
\end{aligned}
\label{eq:2main2b}
\end{align}

where $e=L+R-1$ and $e'=L'+R'-1$ are the excesses of the words
$\Omega={a^\dag}^L a {\adag}^R$ and $\Omega'={a^\dag}^{L'} a
{\adag}^{R'}$, and the quantities $E_L=E_L(e,e')$, $E_R=E_R(e,e')$
are as follows:

Identity     $\hfil E_L$     $\hfil E_R$

(<ref>),(<ref>)       $\max(e,e',0)$      $\max(e',0)$

(<ref>),(<ref>)       $\max(e',0)$      $\max(e,e',0)$.

By homogeneity of the $E_{n,k}$ and the reflection
property (<ref>), the expansions
(<ref>),(<ref>) are equivalent, as one sees by
replacing the summation index $k$ by $n-\nobreak k$; and so are

Restrict each of (<ref>), (<ref>),
(<ref>), (<ref>) by requiring
$L,R,\allowbreak L',R'\in\mathbb{N}$, and also left- and
right-multiply by $(x^{E_L})^n$, $(x^{E_R})^n$, where $E_L,E_R$ are
chosen to make all powers of $x$ non-negative integers; and
substitute $\adag, a$ for $x,D$.

The identities of Theorem <ref>, being parametric, yield
many individual operator ordering identities in $W_{\mathbb{C}}$.
Their coefficients $E_{n,k}$ can be computed from the summation
formula (<ref>) in Theorem <ref> or the
triangular recurrence (<ref>) in
Theorem <ref>.  But by using the closed-form
expressions for certain $E_{n,k}$ in
Section <ref>, one can make these identities more
explicit, perhaps facilitating combinatorial interpretations.  The
following brief discussion focuses on the first
identity (<ref>).

Typical examples of an expansion of a word $\Omega^n$ in twistings of
a word $\Omega'{}^n$ include the following.  Setting $(L,\nobreak
R;\allowbreak L',\nobreak R')$ in (<ref>) equal to
$(1,\nobreak 0;\allowbreak 0,\nobreak 0)$ and $(0,\nobreak
1;\allowbreak 0,\nobreak 0)$ respectively yields the identities
(<ref>),(<ref>) that appeared in the
introduction.  More generally, when $(L,R)\neq(0,0)$, setting
$(L',R')=(0,0)$, so that $e'=-1$ and $E_L=E_R=0$, yields

\begin{equation}
n!\,\left(\adag^La\adag^R\right)^n = {\left(\adag^e\right)}^n
\sum_{k=0}^n E_{n,k}(1-L-R,1;\,R)\, \adag^k (a)^n \adag^{(n-k)},
\end{equation}

where $e=L+R-1$.  This representation of any
$\Omega^n=\left(\adag^La\adag^R\right)^n$ with $(L,R)\neq(0,0)$ as a
combination of twistings of ${\Omega'}^n = a^n$ subsumes

Another example, with explicit coefficients, is

\begin{equation}
\label{eq:last}
(\adag a \adag)^n \adag^n = \sum_{k=0}^n (-1)^{n-k} \binom{n+1}k
\adag^{(n-k)} (\adag^2 a)^n \adag^k.
\end{equation}

Here, $\Omega^n= (\adag a \adag)^n$ and $\Omega'{}^n= (\adag^2 a)^n$.
This is the case $(L,R;L',R')=(1,1;2,0)$ of (<ref>), with
$e=e'=1$ and $E_L=E_R=1$, and the expression used for
$E_{n,k}(-1,-1;\allowbreak1)$ comes from formula (i) of
Theorem <ref>.

An additional example of the expansion of a word $\Omega^n$ in
twistings of a word $\Omega'{}^n$ appeared in
eq. (<ref>) of the introduction.  There, $\Omega^n=
(\adag^2 a)^n$ and $\Omega'{}^n = (\adag^3 a)^n$.
Equation (<ref>) is the case $(L,R;\allowbreak
L',R') = (2,0;\allowbreak3,0)$ of (<ref>), with $e=1$,
$e'=2$, and $E_L=E_R=2$.  The expression for
$E_{n,k}(-1,-2;0)=(-1)^{n} E_{n,k}(1,2;0)$ incorporated
in (<ref>) comes from formula (v) of
Theorem <ref>, with the aid of the reflection
identity (<ref>).

A final observation is that if a word $\Omega^n=({a^\dag}^L a
{\adag}^R)^n$ or $\Omega'{}^n=({a^\dag}^{L'} a {\adag}^{R'})^n$ in the
Eulerian ordering identities of Theorem <ref> is of the
form $({a^\dag} a {\adag})^n$, it can be rewritten by the
Tait–Toscano–Viskov identity as $\adag^n a^n \adag^n$.  An instance
occurs in (<ref>), and in consequence, its left-hand side
equals $\adag^n a^n \adag^{2n}$.  A similar observation applies to the
words $\Omega^n=({a^\dag}^L a {\adag}^R)^n$ and
$\Omega'{}^k=({a^\dag}^{L'}a{\adag}^{R'})^k$ in the Stirling-kind
ordering identities of Theorem <ref>.

§ ACKNOWLEDGEMENTS / AUTHOR DECLARATIONS

This research did not receive any
specific grant from funding agencies in the public, commercial, or
not-for-profit sectors.  Declarations of interest: none.

[Bagno et al., 2019]
authorE. Bagno, authorR. Biagioli,
authorD. Garber, titleSome identities involving second
kind Stirling numbers of types $B$ and $D$,
journalElectron. J. Combin.
volume26 (number3) (year2019)
pagesarticle 3.9 [20 pp.], noteavailable as
arXiv:1901.07830 [math.CO].

[Bagno et al., 2022]
authorE. Bagno, authorD. Garber,
authorM. Novick, titleThe Worpitzky identity for the
groups of signed and even-signed permutations,
journalJ. Algebraic Combin.
volume55 (number2) (year2022)
pages413–428, noteavailable as arXiv:2004.03681

[Bergeron et al., 1992]
authorF. Bergeron, authorP. Flajolet,
authorB. Salvy, titleVarieties of increasing trees,
in: booktitleCAAP '92. Proceedings of the Seventeenth Colloquium
on Trees in Algebra and Programming, no. number581 in
seriesLecture Notes in Computer Science,
publisherSpringer-Verlag, addressNew York/Berlin,
pages24–48, year1992.

[Błasiak and Flajolet, 2012]
authorP. Błasiak, authorP. Flajolet,
titleCombinatorial models of creation–annihilation,
journalSém. Lothar. Combin. volume65
(year2010/12/2012) pagesarticle B65c [78 pp.],
noteavailable as arXiv:1010.0354.

[Boyadzhiev, 2018]
authorK. N. Boyadzhiev, titleNotes on the Binomial
Transform. Theory and Table with Appendix on Stirling Transform,
publisherWorld Scientific, addressHackensack, NJ,

[Brenti, 1994]
authorF. Brenti, title$q$-Eulerian polynomials
arising from Coxeter groups, journalEuropean J. Combin.
volume15 (number5) (year1994)

[Broder, 1984]
authorA. Z. Broder, titleThe $r$-Stirling numbers,
journalDiscrete Math. volume49 (number3)
(year1984) pages241–259.

[Carlitz, 1976]
authorL. Carlitz, titleA theorem on differential
operators, journalAmer. Math. Monthly
volume83 (number5) (year1976)

[Carlitz, 1978]
authorL. Carlitz, titlegeneralized Stirling and
related numbers, journalRiv. Mat. Univ. Parma (4)
volume4 (year1978) pages79–99.

[Carlitz, 1979]
authorL. Carlitz, titleDegenerate Stirling,
Bernoulli and Eulerian numbers, journalUtilitas Math.
volume15 (year1979) pages51–88.

[Charalambides, 1982]
authorC. A. Charalambides, titleOn the Enumeration of
Certain Compositions and Related Sequences of Numbers,
journalFibonacci Quart.
volume20 (number2) (year1982)

[Cheon et al., 2013]
authorG.-S. Cheon, authorJ.-H. Jung,
authorL. W. Shapiro, titleGeneralized Bessel numbers
and some combinatorial settings, journalDiscrete Math.
volume313 (number20) (year2013)

[Cheon and Kim, 2001]
authorG.-S. Cheon, authorJ.-S. Kim,
titleStirling matrix via Pascal matrix,
journalLinear Algebra Appl. volume329
(year2001) pages49–59.

[Comtet, 1974]
authorL. Comtet, titleAdvanced Combinatorics. The Art of
Finite and Infinite Expansions, publisherReidel,
addressBoston/Dordrecht, year1974.

[Duchamp et al., 2004]
authorG. H. E. Duchamp, authorK. A. Penson,
authorA. I. Solomon, authorA. Horzela,
authorP. Błasiak, titleOne-parameter groups and
combinatorial physics, in: editorJ. Govaerts,
editorM. N. Hounkonnou, editorA. Z. Msezane (Eds.),
booktitleContemporary Problems in Mathematical Physics,
publisherWorld Scientific, pages436–449,

[Duchamp et al., 2010]
authorG. H. E. Duchamp, authorL. Poinsot,
authorA. I. Solomon, authorK. A. Penson,
authorP. Błasiak, authorA. Horzela,
titleLadder operators and endomorphisms in combinatorial physics,
journalDiscrete Math. Theor. Comput. Sci,
volume12 (number2) (year2010)
pages23–45, noteavailable as arXiv:0908.2332.

[Fleury and Turbiner, 1994]
authorN. Fleury, authorA. Turbiner,
titlePolynomial relations in the Heisenberg algebra,
journalJ. Math. Phys. volume35 (number11)
(year1994) pages6144–6149.

[Graham et al., 1994]
authorR. L. Graham, authorD. E. Knuth,
authorO. Patashnik, titleConcrete Mathematics,
publisherAddison–Wesley, addressBoston,
edition2nd edn., year1994.

[Herscovici, 2020]
authorO. Herscovici, titleGeneralized permutations
related to the degenerate Eulerian numbers, notepreprint, year2020,
available as https://doi.org/10.48550/arXiv.2007.13205 and arXiv:2007.13205 [math.CO].

[Hsu and Shiue, 1998]
authorL. C. Hsu, authorP. J.-S. Shiue,
titleA unified approach to generalized Stirling numbers,
journalAdv. in Appl. Math.
volume20 (number3) (year1998)

[Hsu and Shiue, 1999]
authorL. C. Hsu, authorP. J.-S. Shiue,
titleOn certain summation problems and generalizations of
Eulerian polynomials and numbers, journalDiscrete Math.
volume204 (number1–3) (year1999)

[Jung et al., 2018]
authorJ.-H. Jung, authorI. Mező,
authorJ. L. Ramírez, titleThe $r$-Bessel
and restricted $r$-Bell numbers, journalAustralas.
J. Combin. volume70 (year2018)

[Katriel, 1974]
authorJ. Katriel, titleCombinatorial aspects of boson
algebra, journalLett. Nuovo Cimento volume10
(year1974) pages565–567.

[Kim and Kim, 2022]
authorT. Kim, authorD. S. Kim,
titleDegenerate $r$-Whitney numbers and degenerate
$r$-Dowling polynomials via boson operators, journalAdv. in
Appl. Math. volume140 (year2022)
pagesarticle 102394 [21 pp.], noteavailable as

[Kim et al., 2023]
authorT. Kim, authorD. S. Kim, authorH. K.
Kim, titleNormal Ordering Associated with $\lambda$-Stirling
numbers in $\lambda$-shift algebra, journalDemonstr. Math.
volume56 (number1) (year2023)
pagesarticle 20220250 [10 pp.], noteavailable as

[Knuth, 1992]
authorD. E. Knuth, titleTwo notes on notation,
journalAmer. Math. Monthly
volume99 (number5) (year1992)

[Lang, 2000]
authorW. Lang, titleOn generalizations of the Stirling
number triangles, journalJ. Integer Seq.
volume3 (number2) (year2000)
pagesarticle 00.2.4 [19 pp.].

[Maier, 2023]
authorR. S. Maier, titleTriangular recurrences,
generalized Eulerian numbers, and related number triangles,
journalAdv. in Appl. Math. volume146
(year2023) pagesarticle 102485 [62 pp.],
noteavailable as arXiv:2207.10224.

[Mansour et al., 2019]
authorT. Mansour, authorJ. L. Ramírez,
authorM. Shattuck, authorS. N. Villamarín,
titlesome combinatorial identities of the
$r$-Whitney-Eulerian numbers, journalAppl. Anal. Discrete
Math. volume13 (number2) (year2019)

[Mansour and Schork, 2016]
authorT. Mansour, authorM. Schork,
titleCommutation Relations, Normal Ordering, and Stirling
Numbers, publisherCRC Press, addressBoca Raton, FL,

[Mansour et al., 2011]
authorT. Mansour, authorM. Schork,
authorM. Shattuck, titleOn a new family of generalized
Stirling and Bell numbers, journalElectron. J. Combin.
volume18 (number1) (year2011)
pagesarticle 77 [33 pp.].

[Mansour et al., 2012]
authorT. Mansour, authorM. Schork,
authorM. Shattuck, titleThe generalized Stirling and
Bell numbers revisited, journalJ. Integer Seq.
volume15 (number8)
(year2012a) pagesarticle 12.8.3
[47 pp.].

[Mansour et al., 2012]
authorT. Mansour, authorM. Schork,
authorM. Shattuck, titleOn the Stirling numbers
associated with the meromorphic Weyl algebra, journalAppl.
Math. Lett. volume25 (number11)
(year2012b) pages1767–1771.

[Méndez et al., 2005]
authorM. A. Méndez, authorP. Błasiak,
authorK. A. Penson, titleCombinatorial approach to
generalized Bell and Stirling numbers and boson normal ordering problem,
journalJ. Math. Phys. volume46 (number8)
(year2005) pagesarticle 083511 [8 pp.].

[Mező and Ramírez, 2016]
authorI. Mező, authorJ. L. Ramírez,
titleSome identities of the $r$-Whitney numbers,
journalAequationes Math.
volume90 (number2) (year2016)

[Mikhailov, 1985]
authorV. V. Mikhailov, titleNormal ordering and
generalised Stirling numbers, journalJ. Phys. A
volume18 (year1985) pages231–235.

[Mohammad-Noori, 2011]
authorM. Mohammad-Noori, titleSome remarks about the
derivation operator and generalized Stirling numbers,
journalArs Combin. volume100 (year2011)

[Nyul and Rácz, 2015]
authorG. Nyul, authorG. Rácz, titleThe
$r$-Lah numbers, journalDiscrete Math.
volume338 (number10) (year2015)

[OEIS Foundation, Inc., 2023]
authorOEIS Foundation, Inc., titleThe On-Line
Encyclopedia of Integer Sequences, howpublishedPublished
electronically at https://oeis.org, year2023.

[Olver et al., 2010]
editorF. W. J. Olver, editorD. W. Lozier,
editorR. F. Boisvert, editorC. W. Clark (Eds.),
titleNIST Handbook of Mathematical Functions,
publisherU.S. Department of Commerce, National Institute of
Standards and Technology, addressWashington, DC,

[Pan, 2012]
authorJ. Pan, titleMatrix decomposition of the unified
generalized Stirling numbers and inversion of the generalized factorial
matrices, journalJ. Integer Seq.
volume15 (number6) (year2012)
pagesarticle 12.6.6 [9 pp.].

[Ramírez et al., 2018]
authorJ. L. Ramírez, authorS. N.
Villamarín, authorD. Villamizar, titleEulerian
numbers associated with arithmetical progressions,
journalElectron. J. Combin.
volume25 (number1) (year2018)
pagesarticle P1.48 [12 pp.].

[Riordan, 1958]
authorJ. Riordan, titleAn Introduction to Combinatorial
Analysis, publisherWiley, addressNew York,

[Sándor and Crstici, 2004]
authorJ. Sándor, authorB. Crstici,
titleHandbook of Number Theory. II, publisherKluwer,
addressDordrecht, year2004.

[Schork, 2021]
authorM. Schork, titleRecent developments in
combinatorial aspects of normal ordering, journalEnumer. Comb.
Appl. volume1 (number2) (year2021)
pagesarticle S2S2 [40 pp.].

[Spivey, 2011]
authorM. Z. Spivey, titleOn solutions to a general
combinatorial recurrence, journalJ. Integer Seq.
volume14(9) (year2011) pagesarticle 11.9.7
[19 pp.], noteavailable as arXiv:1307.2010 [math.CO].

[Srivastava and Singhal, 1973]
authorH. M. Srivastava, authorJ. P. Singhal,
titleNew generating functions for Jacobi and related
polynomials, journalJ. Math. Anal. Appl.
volume41 (number3) (year1973)

[Toscano, 1939]
authorL. Toscano, titleNumeri di Stirling
generalizzati, operatori differenziali e polinomi, journalPont.
Acad. Sci. Comment. volume3 (year1939)

[Toscano, 1955]
authorL. Toscano, titleNuove formule sugli operatori
permutabili di secondo ordine, journalMatematiche (Catania)
volume10 (year1955) pages37–43.

[Wang and Zhu, 2023]
authorZ. Wang, authorZ.-Y. Zhu, titleThe
spiral property of $q$-Eulerian numbers of type $B$,
journalAustralas. J. Combin.
volume87 (number1) (year2023)