**1. Introduction**

The properties of many families of special numbers have been profitably studied by the use of methods tracing back to the umbral calculus [1]. Within this context, the definition of the associated polynomials naturally emerges as umbral Newton binomial convolutions (see "The Bernoulli Polynomials §4.2.2" in Reference [1]). The formalism is extremely powerful, and has allowed for the extension of the method to generalized forms of special numbers ([2,3]). The use of umbral techniques has been recently employed in the study of harmonic numbers, whose relationship to Bernoulli numbers has been pointed out in Reference [4]. In this paper, we will extend the use of umbral methods to the case of higher-order harmonic numbers.

In a number of previous papers ([5–7]), different problems concerning harmonic numbers and the relevant generating functions have been touched. The already mentioned use of the umbral-like formalism has allowed for the framing of the theory of harmonic numbers within an algebraic context. Some of the points raised in ([5–7]) have been reconsidered, made rigorous, and generalized by means of different technical frameworks in successive research ([8–15]).

The present investigation concerns the application of the method foreseen in ([5–7]) to generalized forms of harmonic numbers, such as (We use the notation *mhn* instead of *H*(*m*) *n* recommended in Reference [4] for continuity with previous papers, where it has been adopted to avoid confusion with higher-order Hermite polynomials):

$$\begin{aligned} \,\_m h\_n &= \sum\_{r=1}^n \frac{1}{r^m}, \quad n \ge 1, \\\,\_m h\_0 &= 0, \end{aligned} \tag{1}$$

namely, "higher-order harmonic numbers" satisfying the property:

$$h\_m h\_{n+1} = \,\_m h\_n + \frac{1}{(n+1)^m} \,\_ {m} \tag{2}$$

whose associated series is provided by the limit lim*n*→∞ *<sup>m</sup>hn*, *m* > 1 is, unlike the ordinary harmonic numbers (*m* = 1), not diverging.

It can be argued that for negative *m* values, the Harmonic numbers reduce to a finite sum of integers, expressible in terms of Bernoulli numbers, as discussed in the concluding part of the paper (Remark 3). In the following, we will derive a number of apparently new properties and the relevant consequences.

As an introductory example, we provide the following:

**Example 1.** *We consider the second-order harmonic numbers (m* = 2*) and write*

$$\,\_2h\_{\mathbb{N}} = \int\_0^1 \frac{1 - \mathbf{x}^n}{\mathbf{x} - 1} \ln(\mathbf{x}) \, d\mathbf{x}, \quad \forall n \in \mathbb{N}, \tag{3}$$

*which is obtained after setting*

$$\frac{1}{r^2} = \int\_0^\infty e^{-sr}s \, ds\tag{4}$$

*by noting that*

$$\varepsilon\_2 \hbar\_n = \int\_0^\infty \frac{e^{-s(n+1)} - e^{-s}}{e^{-s} - 1} \text{ s } ds,\tag{5}$$

*and then by changing the variable of integration.*

*It is worth stressing that the integral representation allows for the extension of harmonic numbers to non-integer values of the index. The second-order harmonic numbers interpolates between the integer and real values of the index, as shown in the plot given in Figure 1, where it is pointed out that the asymptotic limit of the second-order harmonic numbers is π*26*.*

The relevant extension to negative real indices will be considered later in the article.

Let us first consider the generating function associated with the second-order harmonic numbers, which can be cast in the form of an umbral exponential series, as follows.

**Definition 1.** *We introduce*

$$\mathbb{E}\_{2^{\tilde{H}}}e(t) := 1 + \sum\_{n=1}^{\infty} \frac{t^n}{n!} \left( {}\_2\mathbb{I}h\_n \right) = e^{\mathbb{A}t} \mathbb{I}\_{0\prime} \tag{6}$$

*where* 2 ˆ *h is an umbral-like operator acting on the vaccum ξ*<sup>0</sup>*, such that (see [10] for a complete treatment of the umbral method):*

$$\begin{aligned} \,\_2\hat{h}^n\mathbb{J}\_0 &:= \mathbb{J}\_n = {}\_2h\_n & n > 0, \\ \,\_2\hat{h}^0\mathbb{J}\_0 = 1 \neq {}\_2h\_0 & 0 \end{aligned} \tag{7}$$

*and*

$$
\varepsilon\_2 \hat{h}^{\nu} \hat{h}^{\mu} \mathbb{Q}\_0 = \varepsilon \hat{h}^{\nu + \mu} \mathbb{q}\_{0\prime} \qquad \forall \nu, \mu \in \mathbb{R}. \tag{8}
$$

(The action of the operator 2 ˆ *h* should be defined as explained in Reference [10]—namely, as the action of a shift operator on its vacuum, such as here *ξ*0, in the following more rigorous way):

$$\begin{aligned} \,\_2\hat{h} &:= \, ^{\partial\_z} , \\ \,\_2\hat{h}^\mu \,\_5\mathbb{X}\_0 &= \left. \epsilon^{\mu \partial\_z} \, \_z \mathbb{X}\_z \right|\_{z=0} = \, ^\mu\_{z+\mu} \vert\_{z=0} = \int\_0^1 \frac{1 - \mathbf{x}^{z+\mu}}{\mathbf{x} - 1} \ln(\mathbf{x}) \, \mathrm{d}x \bigg|\_{z=0} = \int\_0^1 \frac{1 - \mathbf{x}^{\mu}}{\mathbf{x} - 1} \ln(\mathbf{x}) \, \mathrm{d}x . \end{aligned} $$

**Proposition 1.** ∀*m* ∈ N*,*

$$\delta\_{2^{h}}e(t,m) := \partial\_{t}^{m}{}\_{2^{h}}e(t) = {}\_{2}h\_{m} + \sum\_{n=1}^{\infty} \frac{t^{n}}{n!} \left({}\_{2}h\_{n+m}\right) \,. \tag{9}$$

**Proof.** From Equations (6) and (7) it follows that, ∀*m* ∈ N,

$$\partial\_{2}h e(t,m) = \partial\_{t}^{\mathfrak{m}}{}\_{2}h e(t) = \partial\_{t}^{\mathfrak{m}}\mathfrak{e}^{2\mathfrak{j}t}\zeta\_{0} = {}\_{2}\hat{h}^{\mathfrak{m}}\mathfrak{e}^{2\mathfrak{j}t}\zeta\_{0} = {}\_{2}h\_{\mathfrak{m}} + \sum\_{n=1}^{\infty} \frac{t^{n}}{n!} \left({}\_{2}h\_{n+m}\right) \dots$$

**Corollary 1.** *Limiting ourselves to the first derivative only, it appears evident that the generating function* (6) *satisfies the identity*

$$\begin{cases} \partial\_{t\_{2}h}e(t) = \,\_{2h}e(t) + f\_{2}(t), \\ \,\_{2h}e(0) = 1 \end{cases} , \tag{10}$$

*where*

$$\begin{split} f\_2(t) &= \sum\_{n=1}^{\infty} \frac{t^n}{(n+1)(n+1)!} = \frac{1}{t} \int\_0^t \frac{e^s - s - 1}{s} \, ds = -\frac{Ein(-t) + t}{t}, \\ \operatorname{Ein}(z) &= \int\_0^z \frac{1 - e^{-\overline{\zeta}}}{\overline{\zeta}} \, d\zeta. \end{split} \tag{11}$$

**Observation 1.** *The problem of specifying the generating function of second-order harmonic numbers is reduced to the solution of the first-order differential Equation* (10)*. The solution writes*

$$\sigma\_{2h}e(t) = \varepsilon^{t} \left( 1 + \sum\_{n=1}^{\infty} \frac{1}{(n+1)^2} \left( 1 - \varepsilon^{-t} \varepsilon\_n(t) \right) \right), \tag{12}$$

*where*

$$c\_{\rm nl}(\mathbf{x}) = \sum\_{r=0}^{n} \frac{\mathbf{x}^r}{r!} \tag{13}$$

*are the truncated exponential polynomials [16]. They belong to the family of Appél type polynomials [13] and are defined through the operational identity [17,18]:*

$$e\_n(\mathbf{x}) = \frac{1}{1 - \partial\_{\mathbf{x}}} \frac{\mathbf{x}^n}{n!}. \tag{14}$$

**Corollary 2.** *We can further elaborate on the previous identities (Equations* (13) *and* (14)*), and set*

$$\begin{aligned} \sum\_{n=1}^{\infty} \frac{e\_n(t)}{(n+1)^2} &= Q\_2(t), \\ Q\_2(t) &= \frac{1}{1 - \partial\_t} f\_2(t). \end{aligned} \tag{15}$$

*Furthermore, since*

$$\sum\_{n=1}^{\infty} \frac{1}{(n+1)^2} = \frac{\pi^2}{6} - 1,\tag{16}$$

*we end up with*

$$\begin{aligned} \partial\_{2h}\varepsilon(t) &= \varepsilon^t \, \Sigma\_2(t), \\ \Sigma\_2(t) &= \frac{\pi^2}{6} - Q\_2(t) \, \, e^{-t}. \end{aligned} \tag{17}$$

*This new result can be viewed as an extension of the generating function for the first-order harmonic numbers derived by Gosper (see below) [19].*

It is furthermore evident that the formalism allows the straightforward derivation of other identities, such as

**Lemma 1.** ∀*t* ∈ R

$$\sum\_{n=1}^{\infty} \frac{t^n}{n!} \left( {}\_2\mathfrak{h}\_{n+m} \right) = e^t \sum\_{s=0}^m \binom{m}{s} \Sigma\_2^{(s)}(t) - {}\_2\mathfrak{h}\_{m\prime} \tag{18}$$

*where the upper index (s) denotes a s-order derivative and is a direct consequence of the identity in Equation* (9)*.*

The extension to higher-order harmonic numbers with *m* > 2 follows the same logical steps—namely, the derivation of the associated Cauchy problem.

**Corollary 3.** *Let* ∀*t* ∈ R, ∀*p* ∈ N : *p* > 1*,*

$$f\_p(t) = \sum\_{n=1}^{\infty} \frac{t^n}{(n+1)^{p-1} \ (n+1)!} \tag{19}$$

*and*

$$\begin{cases} \partial\_t \left( \,\_p\mathfrak{k}\mathfrak{c}(t) \right) = \,\_p\mathfrak{k}\mathfrak{c}(t) + f\_p(t) \\\\ \,\_p\mathfrak{k}\mathfrak{c}(0) = 1 \end{cases} , \tag{20}$$

*a Cauchy problem. Then, we can write the solution as:*

$$\varepsilon\_{p, \text{fl}} e(t) = \varepsilon^t \left( 1 + \sum\_{n=1}^{\infty} \frac{1}{(n+1)^p} \left( 1 - e^{-t} \varepsilon\_n(t) \right) \right) \tag{21}$$

*or*

$$
\varepsilon\_{p,h}e(t) = \varepsilon^t \,\,\Sigma\_p(t),
\tag{22}
$$

*with*

$$\begin{aligned} \Sigma\_p(t) &= \zeta(p) - Q\_{\overline{p}}(t)e^{-t}, \\ \zeta'(p) &= \sum\_{n=1}^{\infty} \frac{1}{n^p}, \\ Q\_p(t) &= \sum\_{n=1}^{\infty} \frac{1}{(n+1)^p} e\_n(t) = \frac{1}{1 - \partial\_t} f\_p(t). \end{aligned} \tag{23}$$

*The case p* = 1 *should be treated separately, because the sum on the right-hand side of Equation* (21) *apparently diverges.*

**Observation 2.** *It is accordingly worth noting that, since*

$$f\_1(t) = \sum\_{n=1}^{\infty} \frac{t^n}{(n+1)!} = \frac{1}{t}(e^t - t - 1),\tag{24}$$

*we find*

$$\begin{aligned} \varepsilon\_{1h}\varepsilon(t) &= \varepsilon^t \left( 1 + \int\_0^t \frac{1 - (\tau + 1)\,\varepsilon^{-\tau}}{\tau} d\tau \right) = \varepsilon^t \,\Sigma\_1(t), \\ \Sigma\_1(t) &= \varepsilon^{-t} + E\widehat{m}(t), \end{aligned} \tag{25}$$

*which is a restatement of the Gosper derivation of the generating function of first-order harmonic numbers.*

Further comments on the role played by the functions <sup>Σ</sup>*p*(*t*) will be provided in the final section of the paper.

We conclude this introductory section with the inclusion of a further identity.

**Definition 2.** *Here, we introduce higher-order harmonic number umbral polynomials (for m* = 1*, see Definition* 5 *in Reference [7]):*

$$\begin{aligned} \mathbf{u}\_m \boldsymbol{h}\_\mathbf{n}(\mathbf{x}) &:= \left(\mathbf{x} + \, \_m\hat{\mathbf{h}}\right) \, ^\mathbf{n}\mathbf{\hat{g}}\_0 = \mathbf{x}^n + \sum\_{s=1}^n \binom{n}{s} \, \mathbf{x}^{n-s} \, \_m \mathbf{h}\_{s\boldsymbol{\omega}} \\ \mathbf{u}\_m \boldsymbol{h}\_0(\mathbf{x}) &= 1. \end{aligned} \tag{26}$$

The *mhn*(*x*) are introduced in umbral form in complete analogy with those associated with the Bernoulli polynomials. They belong to the Sheffer family, and the relevant properties can be studied by means of the techniques discussed in Reference [17].

It has already been noted that (see Equation (31) in Reference [7])

$$h\_1 h\_n(-1) = (-1)^n \left(1 - \frac{1}{n}\right),\tag{27}$$

which coincides with an analogous identity derived with Mathematica in Reference [20] and, within the present framework from the recurrence

$$h\_1 h\_{n+1}(\mathbf{x}) = (\mathbf{x} + \mathbf{1}) \,\_1h\_{\mathbb{H}}(\mathbf{x}) + \sum\_{s=1}^{n} \binom{n}{s} \frac{1}{(s+1)} \mathbf{x}^{n-s} = (\mathbf{x} + \mathbf{1}) \,\_1h\_{\mathbb{H}}(\mathbf{x}) + n \int\_0^1 dy \int\_0^y (\mathbf{x} + z)^{n-1} dz. \tag{28}$$

Analogous results can be derived for the higher-order case, as in

$$h\_2 h\_n(-1) = (-1)^n \left( 1 - \frac{{}\_2 h\_n}{n} \right),\tag{29}$$

and can be further generalized.

#### **2. Harmonic Numbers and Integral Transforms**

The identities we have dealt with in the previous section can be further generalized if the umbral procedure is merged with other techniques, involving things such as methods of an operational nature.

**Proposition 2.** *We note that an extension of identities of the type reported in Equation* (9) *is provided by the sum* 

$$\sum\_{n=0}^{\infty} \frac{t^n}{n!} \, ^m{p}{}\_p h\_n = e^t \sum\_{k=0}^m \left\{ \begin{array}{c} m \\ k \end{array} \right\} t^k \sum\_{s=0}^k \binom{k}{s} \Sigma\_p^{(k)}(t) \,. \tag{30}$$

*where nk are Stirling numbers of the second kind, namely*

$$\left\{ \begin{array}{c} n \\ k \end{array} \right\} = \frac{1}{k!} \sum\_{j=0}^{k} (-1)^{k-j} j^n \binom{k}{j}. \tag{31}$$

**Proof.** The result in Equation (30) was obtained by merging the umbral formalism with identities of an operational nature. By noting that (see, e.g., Equation (18) in Reference [21] and Equation (1) in [22], and for a more general use of these numbers, see [23–25])

$$(\mathbf{x}\partial\_{\mathbf{x}})^{\
u} = \sum\_{k=0}^{n} \left\{ \begin{array}{c} n \\ k \end{array} \right\} \mathbf{x}^{k} \hat{\partial}\_{\mathbf{x}\prime}^{k} \tag{32}$$

and that if the following sum

$$\sum\_{n=0}^{\infty} \left( t^n a\_n \right) = \sigma(t) \tag{33}$$

does exist, then the following relation holds true

$$\sum\_{n=0}^{\infty} n^m \left(t^n a\_{\mathbb{H}}\right) = \left(t \partial\_t\right)^m \sum\_{n=0}^{\infty} \left(t^n a\_{\mathbb{H}}\right) = \left(t \partial\_t\right)^m \sigma(t) = \sum\_{k=0}^m \left\{\begin{array}{c} m\\k \end{array}\right\} t^k \sigma^{(k)}(t). \tag{34}$$

The result contained in Equation (30) is, therefore, a consequence of Equations (9) and (34).

Furthermore, we noticed that by combining the umbral, Laplace transform, and integral representation methods, we could make further progress.

**Example 2.** *Let us therefore note that (see Corollary* 1 *in [7])*

$$\frac{1}{1 - \,\_1\hat{h}\,\,t}^{\chi}\zeta\_0 = \sum\_{r=1}^{\infty} \,\_1h\_r\,\,t^r + 1,\qquad \vert \,t \vert < 1. \tag{35}$$

*The use of the Laplace transform allowed us to write the left-hand side of Equation* (35) *in the form of*

$$\frac{1}{1 - \,\_1\hat{\boldsymbol{\mu}} \, t} \mathbf{\check{\boldsymbol{\zeta}}}\_0 = \int\_0^\infty e^{-s} e^{s\_1 \hat{\boldsymbol{\mu}} \, t} \, ds \, \boldsymbol{\zeta}\_{0\prime} \tag{36}$$

*which, on account of Equation* (25)*, allows the conclusion controlla se l'eq e' la 25*

$$\begin{aligned} 1 + \sum\_{n=1}^{\infty} t^n \left( {}\_1 h\_{\mathbb{H}} \right) &= 1 - \frac{\ln(1 - t)}{1 - t} = 1 + \frac{L i\_1(t)}{1 - t}, \quad |t| < 1, \\ L i\_m(x) &= \sum\_{r=1}^{\infty} \frac{x^r}{r^m} \equiv \text{Polyllogarithm.} \end{aligned} \tag{37}$$

*The same procedure applied to higher-order harmonic numbers yields the generating functions:*

$$\sum\_{r=1}^{\infty} \, \_m h\_r t^r + 1 = 1 + \frac{L i\_m(t)}{1 - t}, \qquad | \, t \mid < 1,\tag{38}$$

*already known for the case m* = 2*.*

The following example further underscores the versatility of the procedure we have proposed. We will indeed show that the use of the Gaussian identity

$$\epsilon^{b^2} = \frac{1}{\sqrt{\pi t}} \int\_{-\infty}^{\infty} \epsilon^{-\zeta^2 + 2b\xi} d\xi \tag{39}$$

can be exploited to infer further identities on the properties of harmonic numbers.

**Example 3.** *Now, we consider the following generating function*

$$\mathfrak{a}\_{2^h}\mathfrak{e}\_2(t) := 1 + \sum\_{n=1}^{\infty} \frac{t^n}{n!} \left( {}\_2\mathfrak{h}\_{2^n} \right). \tag{40}$$

*According to our formalism, the corresponding r.h.s. can be written as*

$$\mathfrak{a}\_{2h}\mathfrak{e}\_2(t) = \mathfrak{e}^{\left(\frac{1}{2\bar{h}^2}\right)t}\mathfrak{f}\_{0\prime} \tag{41}$$

*which, on account of the identity* (39)*, can be written as*

$$
\epsilon\_{2h}\epsilon\_2(t) = \frac{1}{\sqrt{\pi t}} \int\_{-\infty}^{+\infty} e^{-s^2 + 2s\sqrt{t}\sqrt{t}} s^{\frac{\alpha}{2}} \delta s \,\,\xi\_0. \tag{42}
$$

*The derivation of the sum in Equation* (40) *is therefore reduced to the evaluation of the following integral*

$$\mathcal{L}\_2 \mathfrak{e}\_2(t) = \frac{1}{\sqrt{\pi t}} \int\_{-\infty}^{+\infty} e^{-s^2} \left( \frac{\pi^2}{6} e^{2s\sqrt{t}} - Q\_2(2s\sqrt{t}) \right) ds = \frac{\pi^2}{6} e^t - q\_2(t),\tag{43}$$

*with*

$$q\_2(t) = \frac{1}{\sqrt{\pi t}} \sum\_{n=1}^{\infty} \frac{1}{(n+1)^2} \int\_{-\infty}^{+\infty} e^{-s^2} e\_n \left(2 \, s \, \sqrt{t} \right) ds. \tag{44}$$

#### **3. Final Comments**

Before concluding this paper, we show that the umbral formalism we have employed can be pushed even further to infer new properties of the harmonic numbers. To emphasise this point, we start from certain identities established in Reference [20] by the use of the Mathematica code, SIGMA. Among the examples discussed in Equation (18) in Reference [20], we pick out the following two:

$$\begin{aligned} 1) \quad &\sum\_{r=1}^{n} \binom{n}{r} \frac{(-1)^r}{r} = -{}\_{1}h\_{n}, \\ 2) \quad &\sum\_{r=1}^{n} \binom{n}{r} \frac{(-1)^r}{r} \,\_1h\_r = -{}\_{2}h\_{n} \end{aligned} \tag{45}$$

We can transform the left-hand side of Equations (45) in a Newton binomial by an appropriate definition of umbra.
