Skip to main content


Front. Phys., 01 October 2013
Sec. Statistical and Computational Physics
Volume 1 - 2013 |

Umbral Vade Mecum

Thomas L. Curtright1 Cosmas K. Zachos2*
  • 1Department of Physics, University of Miami, Coral Gables, FL, USA
  • 2High Energy Physics Division, Argonne National Laboratory, Argonne, IL, USA

In recent years the umbral calculus has emerged from the shadows to provide an elegant correspondence framework that automatically gives systematic solutions of ubiquitous difference equations—discretized versions of the differential cornerstones appearing in most areas of physics and engineering—as maps of well-known continuous functions. This correspondence deftly sidesteps the use of more traditional methods to solve these difference equations. The umbral framework is discussed and illustrated here, with special attention given to umbral counterparts of the Airy, Kummer, and Whittaker equations, and to umbral maps of solitons for the Sine-Gordon, Korteweg–de Vries, and Toda systems.

1. Introduction

Robust theoretical arguments have established an anticipation of a fundamental minimum measurable length in Nature, of order LPlanckGN/c3=1.6162×1035m, the corresponding mass and time being MPlanckc/GN=2.1765×108kg and LPlanck/c=5.3911×1044s. The essence of such arguments is the following (in relativistic quantum geometrical units, wherein ħ, c, and MPlanck are all unity).

In a system or process characterized by energy E, no lengths smaller than L can be measured, where L is the larger of either the Schwarzschild horizon radius of the system (~ E) or, for energies smaller than the Planck mass, the Compton wavelength of the aggregate process (~1/E). Since the minimum of max(E, 1/E) lies at the Planck mass (E = 1), the smallest measurable distance is widely believed to be of order LPlanck. Thus, continuum laws in Nature are expected to be deformed, in principle, by modifications at that minimum length scale.

Remarkably, however, if a fundamental spacetime lattice of spacing a = O(LPLanck) is the structure that underlies conventional continuum physics, then it turns out that continuous symmetries, such as Galilei or Lorentz invariance, can actually survive unbroken under such a deformation into discreteness, in a non-local, umbral realization (10, 11, 18).

Umbral calculus, pioneered by Rota and associates in a combinatorial context (4, 16), specifies, in principle, how functions of discrete variables in infinite domains provide systematic “shadows” of their familiar continuum limit properties. By preserving Leibniz's chain rule, and by providing a discrete counterpart of the Heisenberg algebra, observables built from difference operators shadow the Lie algebras of the standard differential operators of continuum physics. [For a review relevant to physics, see (13).] Nevertheless, while the continuous symmetries and Lie algebras of umbrally deformed systems might remain identical to their continuum limit, the functions of observables themselves are modified, in general, and often drastically so.

Traditionally, the controlling continuum differential equations of physics are first discretized (2, 5, 18), and then those difference equations are solved to yield umbral deformations of the continuum solutions. But quite often, routine methods to solve such discrete equations become unwieldy, if not intractable. On the other hand, some technical difficulties may be bypassed by directly discretizing the continuum solutions. That is, through appropriate umbral deformation of the continuum solutions, the corresponding discrete difference equations may be automatically solved. However, as illustrated below for the simplest cases of oscillations and wave propagation, the resulting umbral modifications may present some subtleties when it comes to extracting the underlying physics.

In (21) the linearity of the umbral deformation functional was exploited, together with the fact that the umbral image of an exponential is also an exponential, albeit with interesting modifications, to discretize well-behaved functions occurring in solutions of physical differential equations through their Fourier expansion. This discrete shadowing of the Fourier representation functional should thus be of utility in inferring wave disturbance propagation in discrete spacetime lattices. We continue to pursue this idea here with some explicit examples. We do this in conjunction with the umbral deformation of power series, especially those for hypergeometric functions. We compare both Fourier and power series methods in some detail to gain further insight into the umbral framework.

Overall, we utilize essentially all aspects of the elegant umbral calculus to provide systematic solutions of discretized cornerstone differential equations that are ubiquitous in most areas of physics and engineering. We pay particular attention to the umbral counterparts of the Airy, Kummer, and Whittaker equations, and their solutions, and to the umbral maps of solitons for the Sine-Gordon, Korteweg–de Vries, and Toda systems.

2. Overview of the Umbral Correspondence

For simplicity, consider discrete time, t = 0, a, 2a, …, na, …. Without loss of generality, broadly following the summary review of (13), consider an umbral deformation defined by the forward difference discretization of ∂t,


and whence of the elementary oscillation Equation, (t) = −x(t), namely,


Now consider the solutions of this second-order difference equation. Of course, (2) can be easily solved directly by the textbook Fourier-component Ansatz x(t) ∝ rt, (2), to yield (1 ± ia)t/a. However, to illustrate instead the powerful systematics of umbral calculus (13, 18), we produce and study the solution in that framework.

The umbral framework considers associative chains of operators, generalizing ordinary continuum functions by ultimately acting on a translationally-invariant “vacuum” 1, after manipulations to move shift operators to the right and have them absorbed by that vacuum, which we indicate by T · 1 = 1. Using the standard Lagrange-Boole shift generator

Teat, so that Tf(t)1=f(t+a)T1=f(t+a),(3)

the umbral deformation is then

  ttT1,     (5)

so that [t]0 = 1, and, for n > 0, [0]n = 0. The [t]n are called “basic polynomials”1 for positive n (5, 13, 16), and they are eigenfunctions of tT−1Δ.

A linear combination of monomials (a power series representation of a function) will thus transform umbrally to the same linear combination of basic polynomials, with the same series coefficients, f(t) ↦ f(tT−1). All observables in the discretized world are thus such deformation maps of the continuum observables, and evaluation of their direct functional form is in order. Below, we will be concluding the correspondence by casually eliminating translation operators at the very end, first through operating on the vacuum and then leaving it implicit, so that F(t) ≡ f(tT−1) · 1.

The umbral deformation relies on the respective umbral entities obeying operator combinatorics identical to their continuum limit (a → 0), by virtue of obeying the same Heisenberg commutation relation (18),


Thus, e.g., by shift invariance, TΔT−1 = Δ,

[t,tn]=ntn1  [Δ,[t]nTn]=n[t]n1T1n,(8)

so that, ultimately, Δ[t]n = n[t]n − 1. For commutators of associative operators, the umbrally deformed Leibniz rule holds (10, 11),


ultimately to be dotted onto 1. Formally, the umbral deformation reflects (unitary) equivalences of the unitary irreducible representation of the Heisenberg-Weyl group, provided for by the Stone-von Neumann theorem. Here, these equivalences reflect the alternate consistent realizations of all continuum physics structures through systematic maps such as the one we have chosen. It is worth stressing that the representations of this algebraic relation on the real or complex number fields can only be infinite dimensional, that is, the lattices covered must be infinite.

Now note that, in this case the basic polynomials [t]n are just scaled falling factorials, for n ≥ 0, i.e., generalized Pochhammer symbols, which may be expressed in various ways:

[t]n(tT1)n·1=t(ta)(t(n1)a)=an(t/a)!(t/an)!      =anΓ(ta+1)Γ(tan+1)=(a)nΓ(nta)Γ(ta) .(10)

Thus [−t]n = (−)n[t + a(n − 1)]n. Furthermore, [an]n = ann!; [t]m[tam]nm = [t]n for 0 ≤ mn; and for integers 0 ≤ m < n, [am]n = 0. Thus, Δm[t]n = [an]m[t]nm/am.

Negative umbral powers, by contrast, are the inverse of rising factorials, instead:

[ 1t ]n=(T1t)n1=1(t+a)(t+2a)(t+na)       =an(t/a)!(t/a+n)!       =anΓ(ta+1)Γ(ta+n+1)=(a)nΓ(tan)Γ(ta) .(11)

These correspond to the negative eigenvalues of tT−1Δ.

The standard umbral exponential is then natural to define as (6, 11, 16)2

E(λt,λa)eλ[t]eλtT11=n=0λnn![t]n             =n=0(λa)n(t/an)=(1+λa)t/a,(12)

the compound interest formula, with the proper continuum limit (a → 0). N.B. There is always a 0 at λ = −1/a.

Evidently, since Δ · 1 = 0,


and, as already indicated, one could have solved this equation directly3 to produce the above Et, λa).

Serviceably, the umbral exponential E happens to be an ordinary exponential,


and it actually serves as the generating function of the umbral basic polynomials,


Conversely, then, this construction may be reversed, by first solving directly for the umbral eigenfunction of Δ, and effectively defining the umbral basic polynomials through the above parametric derivatives, in situations where these might be more involved, as in the next section.

As a consequence of linearity, the umbral deformation of a power series representation of a function is given formally by


This may not always be easy to evaluate, but, in fact, the same argument may be applied to linear combinations of exponentials, and hence the entire Fourier representation functional, to obtain

F(t)=dτ f(τ)dω2πeiωτ(1+iωa)t/a      =(1+aτ)t/af(τ)|τ=0 .(17)

The rightmost equation follows by converting iω into ∂τ derivatives and integrating by parts away from the resulting delta function. Naturally, it identifies with Equation (16) by the (Fourier) identity f(∂x)g(x)|x = 0 = g(∂x)f(x)|x = 0. It is up to individual ingenuity to utilize the form best suited to the particular application at hand.

It is also straightforward to check that this umbral transform functional yields

tfΔF ,(18)

and to evaluate the umbral transform of the Dirac delta function, which amounts to a cardinal sine or sampling function,


or to evaluate umbral transforms of rational functions, such as


to obtain an incomplete Gamma function (1), and so on. Note how the last of these is distinctly, if subtly, different from the umbral transform of negative powers, as given in (11).

In practical applications, evaluation of umbral transforms of arbitrary functions of observables may be more direct, at the level of solutions, through this deforming functional, Equation (17). For example, one may evaluate in this way the umbral correspondents of trigonometric functions,

Sin[t]ei[t]ei[t]2i, Cos[t]ei[t]+ei[t]2,(21)

so that

ΔSin[t]=Cos[t], ΔCos[t]=Sin[t].(22)

As an illustration, consider phase-space rotations of the oscillator. The umbral deformation of phase-space rotations,

x˙=p, p˙=x  ΔX(t)=P(t), ΔP(t)=X(t),(23)

readily yields, by directly deforming continuum solutions, the oscillatory solutions,


In view of (14), and also


the umbral sines and cosines in (24) are seen to amount to discrete phase-space spirals,


with a frequency decreased from the continuum value (i.e., 1) to


So the frequency has become, effectively, the inverse of the cardinal tangent function.4 Note that the umbrally conserved quantity is,


such that Δ=0, with the proper energy as the continuum limit.

3. Reduction from Second-Order Differences to Single Term Recursions

In this section and the following, to conform to prevalent conventions, the umbral variable will be denoted by x, instead of t. In this case there is a natural way to think of the umbral correspondence that draws on familiar quantum mechanics language (16): The discrete difference equations begin as operator statements, for operator xs and Ts, but are then reduced to equations involving classical-valued functions just by taking the matrix element 〈x|…|vac〉 where |vac〉 is translationally invariant. The overall x-independent non-zero constant 〈x|vac〉 is then ignored.

To be specific, consider Whittaker's equation (1) for μ = 1/2,


This umbrally maps to the operator statement


Considering either y(xT−1)·1 ≡ Y(x), or else 〈x|y(xT−1)|vac〉 = Y(x) 〈x|vac〉, this operator statement reduces to a classical difference equation,


Before using umbral mapping to convert continuous solutions of (29) into discrete solutions (14, 15) of (31), here we note a simplification of the latter equation upon choosing a = 2, which amounts to setting the scale of x. With this choice (31) collapses to a mere one-term recursion. Shifting xx − 2 this is


Despite being a first-order difference equation, however, the solutions of this equation still involve two independent “constants of summation” even for x restricted to only integer values, because the choice a = 2 has decoupled adjacent vertical strips of unit width on the complex x plane. To be explicit, for integer x > 0, forward iteration gives (2)

Y(2k+1)=2k(j=1kj2kj)Y(1)   andY(2k+2)=2k(j=1kjkj)Y(2)    for integer k0,(33)

with Y(1) and Y(2) the two independent constants that determine values of Y for all larger odd and even integer points, respectively.

Or, if generic x is contemplated, the Equation (32) has elementary solutions, for arbitrary complex constants C1 and C2, given by

Y(x)=2x/2Γ(x2κ)Γ(x2)C1+(2)x/2Γ(x2)Γ(1x2+κ) C2(34)

In the second expression, we have used Γ(z)Γ(1 − z) = π/sin πz. Note the C2 part of this elementary solution differs from the C1 part just through multiplication by a particular complex function with period 2. This is typical of solutions to difference equations since any such periodic factors are transparent to Δ, as mentioned in an earlier footnote (12).

As expected, even for generic x the constants C1 and C2 may be determined given Y(x) at two judiciously chosen points, not necessarily differing by an integer. For example, if 0 < κ < 1,

C1=Γ(1+κ)21+κY(2+2κ) ,   C2=πsinπκC112Γ(κ)Y(2).(36)

Moreover, poles and zeros of the solution are manifest either from the Γ functions in (34), or else from continued product representations such as (33). For the latter, either forward or backward iterations of the first-order difference Equation (32) may be used. Schematically,


or alternatively,


Although both terms in (34) have zeroes, the C1 term also has poles while the C2 term has none—it is an entire function of x—and it is complex for any nonzero choice of C2. Of course, since the Equation (32) is linear, real and imaginary parts may be taken as separate real solutions. All this is evident in the following plots for various selected integer κ.

2x/2Γ(12xκ)Γ(12x) for κ = 1, 2, and 3 in red, blue, and green.

2x/2cosπ(12x)Γ(12x)Γ(112x+κ) for κ = 1, 2, and 3 in red, blue, and green.

2x/2sinπ(12x)Γ(12x)Γ(112x+κ) for κ = 1, 2, and 3 in red, blue, and green.

Collapse to a mere one-term recursion also occurs for an inverse-square potential,


For μa2 = 1, which amounts to setting the scale of the energy of the solution, the umbral version of this equation reduces to

Y(x)=12(1+κa2x(x+a))Y(x+a)      =12(1+aκxaκa+x)Y(x+a).(40)

That is to say,

Y(x+a)=2(1+xa)xa(xa+1+14κ2)(xa+114κ2)Y(x) .(41)

Elementary solutions for generic x, for arbitrary complex constants C1 and C2, are given by


Again, the C2 part of this elementary solution differs from the C1 part just through multiplication by a particular complex function with period a. And again, poles and zeros of these and other solutions are manifest either from those of the Γ functions, or else from a continued product form, e.g.


It is not surprising that (29) and (39) share the privilege to become only first-order difference equations for specific choices of a, as in (32) and (41), because they are both special cases of Whittaker's differential equation, as discussed in the next section. No other linear second-order ODEs lead to umbral equations with this property.

4. Discretization Through Hypergeometric Recursion

In this section we discuss several examples using umbral transform methods to convert solutions of continuum differential equations directly into solutions of the corresponding discretized equations. We use both Fourier and power series umbral transforms.

As an explicit illustration of the umbral transform functional (17), inserting the Fourier representation of the Airy function (1) yields

AiryAi(x)  UmAiryAi(x,a)Re(1π0+e13ik3(1+ika)xadk).(45)

This integral is expressed in terms of hypergeometric functions and evaluated numerically in Appendix A.

Likewise, gaussians also map to hypergeometric functions, as may be obtained by formal series manipulations:

ex2G(x,a)         n=0()n[x]2nn!         =n=01n!(1)na2nΓ(xa+1)Γ(xa2n+1)(46)

where the reflection and duplication formulas were used to write


While the series (47) actually has zero radius of convergence, it is Borel summable, and the resulting regularized hypergeometric function is well-defined. See Appendix B for some related numerics.

For another example drawn from the familiar repertoire of continuum physics, consider the confluent hypergeometric equation of Kummer (A&S 13.1.1):

x y+(βx)yαy=0,(50)

whose regular solution at x = 0, expressed in various dialects, is


with series and integral representations

1F1(α;β;x)=n=0Γ(α+n)Γ(α)Γ(β)Γ(β+n)xnn!                 =Γ(β)Γ(α)Γ(βα)01exssα1(1s)βα1ds                 =1+αβx+12α(α+1)β(β+1)x2+16α(α+1)(α+2)β(β+1)(β+2)x3+O(x4).(52)

The second, independent solution of (50), with branch point at x = 0, is given by Tricomi's confluent hypergeometric function (1), sometimes known as HypergeometricU:

U(α,β,x)=πsinπβ (M(α,β,x)Γ(1+αβ)Γ(β)x1βM(1+αβ,2β,x)Γ(α)Γ(2β)).(53)

Invoking the umbral calculus for x, either of these confluent hypergeometric functions can be mapped onto their umbral counterparts using

1F1(α;β;x)  2F1(α,xa;β;a),(54)

where 2F1 is the well-known Gauss hypergeometric function (1). This map from 1F1 to 2F1 follows from the basic monomial umbral map,


and from the series (52). When combined, these give the well-known series representation of 2F1.

Next, reconsider the one-dimensional Coulomb problem defined by Whittaker's equation for general μ (1):


Since κ and μ are both arbitrary, this also encompasses the inverse-square potential, (39). Exact solutions of this differential equation are

whittakerW(κ,μ,x)=xμ+1/2ex/2            (Γ(2μ)Γ(μκ+12) 1F1(μκ+12;2μ+1;x)+Γ(2μ)Γ(μκ+12)x2μ1F1(μκ+12;2μ+1;x)).(59)

Umbral versions of these solutions are complicated by the exponential and overall power factors in the classical relations between the 1F1's and the Whittaker functions, but this complication is manageable. (In part this is because in the umbral calculus there are no ordering ambiguities (20).)

To obtain the umbral version of the Whittaker functions, we begin by evaluating

e12xT11F1(α;β;xT1)·1=m=0n=0(12)mΓ(α+n)Γ(α)Γ(β+n)Γ(β)[ x ]m+nm!n!                                       =(1a2)xa 2F1(α,xa;β;2aa2),(60)

where we have performed the sum over m first, to obtain


The sum over n then gives the Gauss hypergeometric function in (60).

Next, to deal with the umbral deformations of the Whittaker functions, we need to use the continuation of (10) and (11) to an arbitrary power of xT−1, namely,


This continuation leads to the following:

(xT1)γe12xT11F1(α;β;xT1)·1    =aγΓ(xa+1)Γ(xaγ+1)Tγe12xT11F1(α;β;xT1)·1    =aγΓ(xa+1)Γ(xaγ+1)e12(xγa)T11F1(α;β;(xγa)T1)·1.(63)

Thus we obtain the umbral map

xγe12x1F1(α;β;x)Γ(xa+1)Γ(xaγ+1) aγ(1a2)xaγ 2F1(α,γxa;β;2aa2).(64)

Finally then, specializing to the relevant α, β, and γ, we find the umbral Whittaker functions. In particular,

whittakerM(κ,μ,x)Γ(xa+1)Γ(xaμ+12)aμ+1/2(1a2)xaμ12 2F1(μ+12κ,μ+12xa;2μ+1;2aa2).(65)

This result for general a exhibits what is special about the choice a = 2, as exploited in the previous section. To realize that choice from (65) requires taking a limit a ↗ 2, hence it requires the asymptotic behavior of the Gauss hypergeometric function (1):


Now with sufficient care, a = 2 solutions can be coaxed from the umbral version of whittakerM in (65), and/or the corresponding umbral counterpart of whittakerW, upon taking lima↗2 and making use of (66). Moreover, in principle the umbral correspondents of both Whittaker functions could be used to obtain from this limit a solution with two arbitrary constants.

On the other hand, for a = 2, the umbral equation corresponding to (56) again reduces to a one-term recursion, namely,


For generic x, solutions for arbitrary complex constants C1 and C2 are then given by

      =2x/2Γ(1+x2)Γ(x2κ)Γ(x2+12+μ)Γ(x2+12μ)(C1+1π2 C2sin(πx2)sinπ(x2κ)),(69)

which agrees with (34) when μ = 1/2, of course. As in that previous special case, the C2 part of (68) differs from the C1 part just through multiplication by a particular complex function with period 2 (12).

We graph some examples to show the differences between the Whittaker functions and their umbral counterparts, for a = 1.

whittakerM(κ, 1/2, x) for κ = 1, 2, and 3 in red, blue, and green.

Umbral whittakerM(κ, 1/2, x) for a = 1, and for κ = 1, 2, and 3 in red, blue, and green.

The examples above are specific illustrations of combinatorics that may be summarized in a few umbral hypergeometric mapping lemmata, the simplest being

Lemma 1:

pFq(α1,,αp;β1,,βq;x)p+1Fq(α1,,αp,x/a;β1,,βq;a) ,(70)

where the series representation of the generalized hypergeometric function pFq is5


A proof of (70) follows from formal manipulations of these series.

The umbral version of a more general class of functions is obtained by replacing xxT−1 in functions of xk for some fixed positive integer k. Thus, again for hypergeometric functions, we have

Lemma 2:

pFq(α1,,αp;β1,,βq;xk)p+kFq( α1,,αp,1k(xa),1k(1xa),, 1k(k1xa);β1,,βq;(ak)k ).(72)

And again, a proof follows from formal series expansions.

Multiplication by exponentials produces only minor modifications of these general results, as was discussed above in the context of Whittaker functions, namely,

Lemma 3:

eλxpFq(α1,,αp;β1,,βq;xk)(1+aλ)xap+kFq( α1,,αp,1k(xa),1k(1xa),, 1k(k1xa);β1,,βq;(ak1+aλ)k ).(73)

In addition, multiplication by an overall power of x gives

Lemma 4:

xγeλxpFq(α1,,αp;β1,,βq;xk)Γ(xa+1)aγ(1+aλ)xaγΓ(xaγ+1) p+kFq( α1,,αp,γxak,1+γxak,,k1+γxak; β1,,βq;(ak1+aλ)k ).(74)

5. Wave Propagation

Given the umbral features of discrete time and space equations discussed above, separately, it is natural to combine the two.

For example, the umbral version of simple plane waves in 1+1 spacetime would obey an equation of the type (6, 11, 12),

(Δx2Δt2) F=0 ,(75)

on a time-lattice with spacing a and a space-lattice with spacing b, not necessarily such that b = a in all spacetime regions. For generic frequency, wavenumber and velocity, the basic solutions are


For right-moving waves, say, these have phase velocity

v(ω,k)=ωk aarcsin(b)barcsin(a).(77)

Thus, the effective index of refraction in the discrete medium is (b arcsin(a))/(a arcsin(b)), i.e., modified from 1. Small inhomogeneities of a and b in the fabric of spacetime over large regions could therefore yield interesting effects.

Technically, a more challenging application of umbral methods involves nonlinear, solitonic phenomena (21), such as the one-soliton solution of the continuum Sine-Gordon equation,

(x2t2)f(x,t)=sin(f(x,t)) ,  fSG(x,t)=4arctan(mexvt1v2).(78)

The corresponding umbral deformation of the PDE itself would now also involve a deformed potential sin(f(xT−1x, tT−1t))·1. But rather than tackling this difficult nonlinear difference equation, one may instead use the umbral transform (17) to infer that fSG(x, t) maps to


The continuum Korteweg–de Vries soliton is likewise mapped:


Closed-form evaluations of these Fourier integrals are not available, but the physical effects of the discretization could be investigated numerically, and compared to the Lax pair integrability machinery of (13), or to the results on a variety of discrete KdVs in (17), or to other studies (8, 9).

However, a more accessible example of umbral effects on solitons may be found in the original Toda lattice model (19). For this model the spatial variable is already discrete, usually with spacing b = 1 so x = n is an integer, while the time t is continuous. The equations of motion in that case are


for integer n. Though x = n is discrete, nevertheless there are exact multi-soliton solutions valid for all continuous t, as is well-known.

Specific one-soliton Toda solutions are given for constant α, β, γ, and q0 by


provided that


So the soliton's velocity is just v=±2βsinh(β2).

While obtained only for discrete x = n, for plotting purposes q(n, t) may be interpolated for any x (see graph below). To carry out the complete umbral deformation of this system, it is then only necessary to discretize t in the equations of motion (81). Consider what effects this approach to discrete time has on the specified one-soliton solutions.

To that end, expand the exact solutions in (82) as series,


Upon umbralizing t, the one-soliton solutions then map as


and these are guaranteed to give solutions to the umbral operator equations of motion,

Δp(n,tT1)1a(T1)p(n,tT1)                   =(e(q(n+1,tT1)q(n,tT1))e(q(n,tT1)q(n1,tT1))),(88)

upon projecting onto a translationally invariant “vacuum” (i.e., Q(n, t) ≡ q(n, tT−1) · 1).

Now, for integer time steps, t/a = m, consider the series at hand:

S(m,c,z)=k=1zkk(ekβ1)(1+ck)m             =ln(1z1zeβ)+j=1mcj(mj)R(j,z),(89)

where c = γa, z = −αe−βn, and where for j > 0,

R(j,z)=k=0(ekβ1)zkkj1         Φ(eβz,1j,0)Φ(z,1j,0).(90)

Fortunately, for positive integer t/a, we only need the Lerch transcendent function,


for those cases where the sums are expressible as elementary functions. For example,


The ln(…) term on the RHS of (89) then reproduces the specified classical one-soliton solutions at t = 0, while the remaining terms give umbral modifications for t ≠ 0.

Altogether then, we have


These umbral results are compared to some time-continuum soliton profiles for t/a = 0, 1, 2, 3, and 4 in the following Figure (with q0 = 0, α = 1 = β, and γ = 2 sinh(1/2) = 1.042).

Toda soliton profiles q interpolated for all x ∈ [−5, 5] at integer time slices superimposed with their time umbral maps Q (thicker curves) for a = 1.

Thus, the umbral-mapped solutions no longer evolve just by translating the profile shape. Rather, they develop oscillations about the classical fronts that dramatically increase with time, that evince not only dispersion but also generation of harmonics, and that, strictly speaking, disqualify use of the term soliton for their description. Be that as it may, this model is referred to in some studies as integrable (8, 9).

These umbral effects on wave propagation evoke scattering and diffraction by crystals. But here the “crystal” is spacetime itself. It is tempting to speculate based on this analogy. In particular, were a well-formed wave packet to pass through a localized region of crystalline spacetime, with sufficiently large lattice spacings, the packet could undergo dramatic deformations in shape, wavelength, and frequency—far greater than and very different from what would be expected just from the dispersion of a free packet propagating through continuous space and time.

6. Concluding Remarks

We have emphasized how the umbral calculus has visibly emerged to provide an elegant correspondence framework that automatically gives solutions of ubiquitous difference equations as maps of well-known continuous functions. This correspondence systematically sidesteps the use of more traditional methods to solve these difference equations.

We have used the umbral calculus framework to provide solutions to discretized versions of several differential equations that are widespread building-blocks in many if not all areas of physics and engineering, thereby avoiding the rather unwieldy frontal assaults often engaged to solve such discrete equations directly.

We have paid special attention to the Airy, Kummer, and Whittaker equations, and illustrated several basic principles that transform their continuum solutions to umbral versions through the use of hypergeometric function maps. The continuum limits thereof are then manifest.

Finally, we have applied the solution-mapping technique to single solitons of the Sine-Gordon, Korteweg–de Vries, and Toda systems, and we have noted how their umbral counterparts—particular solutions of corresponding discretized equations—evince dispersion and other non-solitonic behavior, in general. Such corrections to the continuum result may end up revealing discrete spacetime structure in astrophysical wave propagation settings.

We expect to witness several applications of the framework discussed and illustrated here.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


This work was supported in part by NSF Award PHY-1214521; and in part, the submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory. Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. Thomas L. Curtright was also supported in part by a University of Miami Cooper Fellowship.


1. ^We stress that the notation [t]n is shorthand for the product t(ta)…(t − (n − 1)a). It is not just the nth power of [t] = t.

2. ^Again we stress that eλ[t] is a short-hand notation, and not just the usual exponential of λ[t] = λt.

3. ^N.B. There is an infinity of “non-umbral” extensions of the Et, λa) solution (12): Multiplying the umbral exponential by an arbitrary periodic function g(t + a) = g(t) will pass undetected through Δ, and thus will also yield an eigenfunction of Δ. Often, such extra solutions have either a vanishing continuum limit, or else an ill-defined one.

4. ^That is, for Θ ≡ arctan(a), the spacing of the zeros, period, etc, are scaled up by a factor of tanc(Θ)tan(Θ)Θ1. For complete periodicity on the time lattice, one further needs return to the origin in an integral number of N steps, thus a solution of N = 2πn/arctan a. Example: For a = 1, the solutions' radius spirals out as 2t/2, while ω = π/4, and the period is τ = 8.

5. ^Recall results from using the ratio test to determine the radius of convergence for the pFq1, …, αp1,…,βq;x) series:

If p < q + 1 then the ratio of coefficients tends to zero. This implies that the series converges for any finite value of x.

If p = q + 1 then the ratio of coefficients tends to one, hence the series converges for |x| < 1 and diverges for |x| > 1.

If p > q + 1 then the ratio of coefficients grows without bound. The series is then divergent or asymptotic, and is a symbolic shorthand for the solution to a differential equation.


1. Abramowitz M, Stegun I. Handbook of Mathematical Functions., National Bureau of Standards. AMS 55, (1964).

2. Bender C, Orszag S. Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill (1978).

3. Cholewinski F, Reneke J. Electron J Diff Eq. (2003) 2003:1–64.

4. Di Bucchianico A, Loeb D. A selected survey of umbral calculus. Electron J Combin. (1995) DS3.

5. Dimakis A, Müller-Hoissen F, Striker T. Umbral calculus, discretization, and quantum mechanics on a lattice. J Phys. (1996) A 29:6861–76.

6. Floreanini R, Vinet L. Lie symmetries of finite-difference equations. J Math Phys. (1995) 36:7024–42

7. Levi D, Negro J, del Olmo M. Discrete q-derivatives and symmetries of q-difference equations. J Phys. (2004) 37:3459–73.

8. Grammaticos B, Kosmann-Schwarzbach Y, Tamizhmani T. (eds.). Discrete Integrable Systems. In: Lect Notes Phys 644:Springer (2004). doi: 10.1007/b94662.

CrossRef Full Text

9. Grammaticos B, Ramani A, Willox R. A sine-Gordon cellular automaton and its exotic solitons. J Phys A Math Theor. (2013) 46:145204.

10. Levi D, Negro J, del Olmo M. Discrete derivatives and symmetries of difference equations. J. Phys. (2001) 34:2023–2030.

11. Levi D, Tempesta P, Winternitz P. Lorentz and galilei invariance on lattices. Phys Rev. (2004a) D69:105011.

12. Levi D, Tempesta P, Winternitz P. Umbral calculus, difference equations and the discrete Schrödinger equation. J Math Phys. (2004b) 45: 4077–4105.

13. Levi D, Winternitz P. Continuous symmetries of difference equations Schiff: Loop groups and discrete KdV equations. J Phys. (2006) A39:R1–R63.

14. López-Sendino J, Negro J, Del Olmo M, Salgado E. “Quantum mechanics and umbral calculus” J. Phys. (2008) 128:012056. doi: 10.1088/1742-6596/128/1/012056

CrossRef Full Text

15. López-Sendino J, Negro J, Del Olmo M. “Discrete coulomb potential” Phys Atom. Nuclei (2010) 73.2:384–90.

16. Rota G-C. Finite Operator Calculus, Academic Press (1975).

17. Schiff J. Loop groups and discrete KdV equations. Nonlinearity (2003) 16:257–75.

18. Smirnov Y, Turbiner A. Lie algebraic discretization of differential equations. Mod Phys Lett. (1995) A10:1795 [Erratum-ibid. A10:3139 (1995)].

19. Toda M. Theory of Nonlinear Lattices (2nd Edn.) Springer (1989).

20. Ueno K. Umbral calculus and special functions; Hypergeometric series formulas through operator calculus. Adv Math. (1988) 67:174–229; Funkcialaj Ekvacioj (1990) 33:493–518.

21. Zachos CK. Umbral deformations on discrete space-time. Int J Mod Phys. (2008) A23:2005.

Appendix A: Umbral Airy Functions

Formally, these can be obtained by expressing the Airy functions in terms of hypergeometric functions and then umbral mapping the series. The continuum problem is given by

yxy=0, y(x)=C1AiryAi(x)+C2AiryBi(x),(94)


AiryAi(x)=132/3Γ(2/3) 0F1(;23;19x3)131/3Γ(1/3) 0F1(;43;19x3),(95)
AiryBi(x)=131/6Γ(2/3) 0F1(;23;19x3)+31/6zΓ(1/3) 0F1(;43;19x3).(96)

The yY umbral images of these, solving the umbral discrete difference equation (3, 12)


are then given by (72) for k = 3. In particular,

UmAiryAi(x,a)=132/3Γ(2/3) 3F1(13xa,13(1xa),13(2xa);23;3a3)                       131/3Γ(1/3) 3F1(13xa,13(1xa),13(2xa);43;3a3).(98)

Since the number of “numerator parameters” in the hypergeometric function 3F1 exceeds the number of “denominator parameters” by 2, the series expansion is at best asymptotic. However, the series is Borel summable. In this respect, the situation is the same as for the umbral gaussian (see Appendix B).

Alternatively, as previously mentioned in the text, using the familiar integral representation of AiryAi(x), the umbral map devolves to that of an exponential. That is to say,

  AiryAi(xT1)=12π+exp(13is3+isxT1)ds                        (99)

Just as AiryAi(x) is a real function for real x, UmAiryAi(x, a) is a real function for real x and a,


After some hand-crafting, the final result may be expressed in terms of just three 2F2 generalized hypergeometric functions. To wit,


where the hypergeometric functions 2F2 (a, b; c, d; z) appear in the expression as

H1(w,z)=Γ(13z)Γ(13+13z) 2F2(13z,13+13z;13,23;13w3),(103)
H2(w,z)=Γ(13+13z)Γ(23+13z) 2F2(13+13z,23+13z;23,43;13w3),(104)
H3(w,z)=Γ(23+13z)Γ(1+13z) 2F2(23+13z,1+13z;43,53;13w3),(105)

and where the coefficients in (102) are

C3(w,z)=336cos(12πz12πzsignum(w))636cos(16πz+12πzsignum(w))                +636cos(56πz+12πzsignum(w))+3×323sin(12πz12πzsignum(w))                336cos(12πz12πzsignum(w))+2×323sin(16πz+12zπsignum(w))                +323sin(12πz+12πzsignum(w))+2×323sin(56πz+12πzsignum(w)).(109)

While the coefficient functions C0−3 are not pretty, they are comprised of elementary functions, and they are nonsingular functions of z. On the other hand, the hypergeometric functions do have singularities and discontinuities for negative z. However, the net result for UmAiryAi is reasonably well-behaved.

We plot UmAiryAi(x, a) for a=0, ±14, ±12, and ±1.

UmAiryAi(x, a) for a = ±1, ±1/2, and ±1/4 (red, blue, & green dashed/solid curves, resp.) compared to AiryAi(x) = UmAiryAi(x, 0) (black curve).

Appendix B: Umbral Gaussians

As discussed in the text, straightforward discretization of the series yields the umbral gaussian map:


(NB G(x, a) ≠ G(−x, a).) Now, it is clear that term by term the series (110) reduces back to the continuum gaussian as a → 0. Nonetheless, since the series is asymptotic and not convergent for |a| > 0, it is interesting to see how this limit is obtained from other representations of the hypergeometric function in (111), in particular from using readily available numerical routines to evaluate 2F0 for specific small values of a. Some examples are shown here.

G(x, 1/2n) vs. x ∈ [−3,2], for n = 1, 2, and 3, in red, blue, and green, respectively, compared to G(x, 0) = exp(−x2), in black.

Mathematicaő code is available online to produce similar graphs, for those interested. It is amusing that Mathematica manipulates the Borel regularized sum to render the 2F0 in question in terms of Tricomi's confluent hypergeometric function U, as discussed above in the context of Kummer's Equation, cf. (53). Thus G can also be expressed in terms of 1F1s. The relevant identities are:


Keywords: umbral correspondence, discretization, difference equations, umbral transform, hypergeometric functions

Citation: Curtright TL and Zachos CK (2013) Umbral Vade Mecum. Front. Physics 1:15. doi: 10.3389/fphy.2013.00015

Received: 27 June 2013; Paper pending published: 04 September 2013;
Accepted: 10 September 2013; Published online: 01 October 2013.

Edited by:

Manuel Asorey, Universidad de Zaragoza, Spain

Reviewed by:

Apostolos Vourdas, University of Bradford, UK
An Huang, Harvard University, USA
Mariano A. Del Olmo, Universidad de Valladolid, Spain

Copyright © 2013 Curtright and Zachos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Cosmas K. Zachos, High Energy Physics Division 362, Argonne National Laboratory, Argonne, IL 60439-4815, USA e-mail: