Séminaire Lotharingien de Combinatoire 44 (2000), Article B44d

Mathemagics^{⋆} (A Tribute to L. Euler and R. Feynman)

Pierre Cartier^{⋆⋆}

CNRS, Ecole Normale Supérieure de Paris, 45 rue d’Ulm, 75230 Paris Cedex 05

To the memory of Gian-Carlo Rota,

modern master of mathematical magical tricks

Table of contents

1. Introduction

2. A new look at the exponential 2.1. The power of exponentials 2.2. Taylor’s formula and exponential 2.3. Leibniz’s formula 2.4. Exponential vs. logarithm 2.5. Infinitesimals and exponentials 2.6. Differential equations

3. Operational calculus 3.1. An algebraic digression: umbral calculus 3.2. Binomial sequences of polynomials 3.3. Transformation of polynomials 3.4. Expansion formulas 3.5. Signal transforms 3.6. The inverse problem 3.7. A probabilistic application 3.8. The Bargmann-Segal transform 3.9. The quantum harmonic oscillator

4. The art of manipulating infinite series 4.1. Some divergent series 4.2. Polynomials of infinite degree and summation of series

⋆ Lectures given at a school held in Chapelle des Bois (April 5-10, 1999) on “Noise, oscillators

and algebraic randomness”. ⋆⋆ cartier@ihes.fr

4.3. The Euler-Riemann zeta function 4.4. Sums of powers of numbers 4.5. Variation I: Did Euler really fool himself? 4.6. Variation II: Infinite products

5. Conclusion: From Euler to Feynman References

1 Introduction

The implicit philosophical belief of the working mathematician is today the Hilbert-Bourbaki formalism. Ideally, one works within a closed system: the basic principles are clearly enunciated once for all, including (that is an addition of twentieth century science) the formal rules of logical reasoning clothed in mathematical form. The basic principles include precise definitions of all mathematical objects, and the coherence between the various branches of mathematical sciences is achieved through reduction to basic models in the universe of sets. A very important feature of the system is its non-contradiction ; after Gödel, we have lost the initial hopes to establish this non-contradiction by a formal reasoning, but one can live with a corresponding belief in non-contradiction. The whole structure is certainly very appealing, but the illusion is that it is eternal, that it will function for ever according to the same principles. What history of mathematics teaches us is that the principles of mathematical deduction, and not simply the mathematical theories, have evolved over the centuries. In modern times, theories like General Topology or Lebesgue’s Integration Theory represent an almost perfect model of precision, flexibility and harmony, and their applications, for instance to probability theory, have been very successful.

My thesis is: there is another way of doing mathematics, equally successful, and the two methods should supplement each other and not fight.

This other way bears various names: symbolic method, operational calculus, operator theory … Euler was the first to use such methods in his extensive study of infinite series, convergent as well as divergent. The calculus of differences was developed by G. Boole around 1860 in a symbolic way, then Heaviside created his own symbolic calculus to deal with systems of differential equations in electric circuitry. But the modern master was R. Feynman who used his diagrams, his disentangling of operators, his path integrals … The method consists in stretching the formulas to their extreme consequences, resorting to some internal feeling of coherence and harmony. They are obvious pitfalls in such methods, and only experience can tell you that for the Dirac δ-function an expression like xδ(x) or δ′(x) is lawful, but not δ(x)/x or δ(x)^{2}. Very often, these so-called symbolic methods have been substantiated by later rigorous developments, for instance Schwartz distribution theory gives a rigorous meaning to δ(x), but physicists used sophisticated formulas in “momentum space” long before Schwartz codified the Fourier transformation for distributions. The Feynman “sums over histories” have been immensely successful in many problems, coming from physics as well from mathematics, despite the lack of a comprehensive rigorous theory.

To conclude, I would like to offer some remarks about the word “formal”. For the mathematician, it usually means “according to the standard of formal rigor, of formal logic”. For the physicists, it is more or less synonymous with “heuristic” as opposed to “rigorous”. It is very often a source of misunderstanding between these two groups of scientists.

2 A new look at the exponential

2.1 The power of exponentials The multiplication of numbers started as a shorthand for repeated additions, for instance 7 times 3 (or rather “seven taken three times”) is the sum of three terms equal to 7

7 × 3 = \underbrace{7 + 7 + 7}_{3 times} .

In the same vein 7^{3} (so denoted by Viete and Descartes) means \underbrace{7 × 7 × 7}_{3 factors}.

There is no difficulty to define x^{2} as xx or x^{3} as xxx for any kind of multiplication (numbers, functions, matrices …) and Descartes uses interchangeably xx or x^{2}, xxx or x^{3}.

In the exponential (or power) notation, the exponent plays the role of an operator. A great progress, taking approximately the years from 1630 to 1680 to accomplish, was to generalize a^{b} to new cases where the operational meaning of the exponent b was much less visible. By 1680, a well defined meaning has been assigned to a^{b} for a, b real numbers, a > 0. Rather than to retrace the historical route, we shall use a formal analogy with vector algebra. From the original definition of a^{b} as a × … × a (b factors), we deduce the fundamental rules of operation, namely

(a × a′)^{b} = a^{b} × a′^{b}, a^{b+b′} = a^{b} × a^{b′}, (a^{b})^{b′} = a^{bb′}, a^{1} = a.  (1)

The other rules for manipulating powers are easy consequences of the rules embodied in (1). The fundamental rules for vector algebra are as follows:

(v + v′).λ = v.λ + v′.λ, v.(λ + λ′) = v.λ + v.λ′,
(v.λ).λ′ = v.(λλ′), v.1 = v.  (2)

The analogy is striking provided we compare the product a × a′ of numbers to the sum v + v′ of vectors, and the exponentiation a^{b} to the scaling v.λ of the vector v by the scalar λ.

In modern terminology, to define a^{b} for a, b real, a > 0 means that we want to consider the set R_{+}^{×} of real numbers a > 0 as a vector space over the field of real numbers R. But to vectors, one can assign coordinates: if the coordinates of the vector v(v′) are the v_{i}(v′_{i}), then the coordinates of v + v′ and v.λ are respectively v_{i} + v′_{i} and v_{i}.λ. Since we have only one degree of freedom in R_{+}^{×}, we should need one coordinate, that is a bijective map L from R_{+}^{×} to R such that

L(a × a′) = L(a) + L(a′).  (3)

Once such a logarithm L has been constructed, one defines a^{b} in such a way that L(a^{b}) = L(a).b. It remains the daunting task to construct a logarithm. With hindsight, and using the tools of calculus, here is the simple definition of “natural logarithms”

ln(a) = ∫_{1}^{a} dt/t for a > 0.  (4)

In other words, the logarithm function ln(t) is the primitive of 1/t which vanishes for t = 1. The inverse function exp s (where t = exp s is synonymous to ln(t) = s) is defined for all real s, with positive values, and is the unique solution to the differential equation f′ = f with initial value f(0) = 1. The final definition of powers is then given by

a^{b} = exp(ln(a).b).  (5)

If we denote by e the unique number with logarithm equal to 1 (hence e = 2.71828 …), the exponential is given by exp a = e^{a}.

The main character in the exponential is the exponent, as it should be, in complete reversal from the original view where 2 in x^{2}, or 3 in x^{3} are mere markers.

2.2 Taylor’s formula and exponential

We deal with the expansion of a function f(x) around a fixed value x_{0} of x, in the form

f(x_{0} + h) = c_{0} + c_{1}h + ⋯ + c_{p}h^{p} + · · · .  (6)

This can be an infinite series, or simply a finite order expansion (include then a remainder). If the function f(x) admits sufficiently many derivatives, we can deduce from (6) the chain of relations

f′(x_{0} + h) = c_{1} + 2c_{2}h + ⋯
f′′(x_{0} + h) = 2c_{2} + 6c_{3}h + ⋯
f′′′(x_{0} + h) = 6c_{3} + 24c_{4}h + · · · .

By putting h = 0, deduce

f(x_{0}) = c_{0}, f′(x_{0}) = c_{1}, f′′(x_{0}) = 2c_{2}, …

and in general f^{(p)}(x_{0}) = p!c_{p}. Solving for the c_{p}’s and inserting into (6) we get Taylor’s expansion

f(x_{0} + h) = ∑_{p≥0} \frac{1}{p!}f^{(p)}(x_{0})h^{p}.  (7)

Apply this to the case f(x) = exp x, x_{0} = 0. Since the function f is equal to its own derivative f′, we get f^{(p)} = f for all p’s, hence f^{(p)}(0) = f(0) = e^{0} = 1. The result is

exp h = ∑_{p≥0} \frac{1}{p!}h^{p}.  (8)

This is one of the most important formulas in mathematics. The idea is that this series can now be used to define the exponential of large classes of mathematical objects: complex numbers, matrices, power series, operators. For the modern mathematician, a natural setting is provided by a complete normed algebra A, with norm satisfying ||ab|| ≤ ||a|| · ||b||. For any element a in A, we define exp a as the sum of the series ∑_{p≥0} a^{p}/p!, and the inequality

||a^{p}/p!|| ≤ ||a||^{p}/p!  (9)

shows that the series is absolutely convergent.

But this would not exhaust the power of the exponential. For instance, if we take (after Leibniz) the step to denote by Df the derivative of f, D^{2}f the second derivative, etc… (another instance of the exponential notation!), then Taylor’s formula reads as

f(x + h) = ∑_{p≥0} \frac{1}{p!}h^{p}D^{p}f(x).  (10)

This can be interpreted by saying that the shift operator T_{h} taking a function f(x) into f(x+h) is equal to ∑_{p≥0} \frac{1}{p!}h^{p}D^{p}, that is, to the exponential exp hD (question: who was the first mathematician to cast Taylor’s formula in these terms?). Hence the obvious operator formula T_{h+h′} = T_{h} · T_{h′} reads as

exp(h + h′)D = exp hD · exp h′D.  (11)

Notice that for numbers, the logarithmic rule is

ln(a · a′) = ln(a) + ln(a′)  (12)

according to the historical aim of reducing via logarithms the multiplications to additions. By inversion, the exponential rule is

exp(a + a′) = exp(a) · exp(a′).  (13)

Hence formula (11) is obtained from (13) by substituting hD to a and h′D to a′.

But life is not so easy. If we take two matrices A and B and calculate exp(A + B) and exp A · exp B by expansion we get

exp(A + B) = I + (A + B) + \frac{1}{2}(A + B)^{2} + \frac{1}{6}(A + B)^{3} + ⋯ (14)

exp A. exp B = I + (A + B) + \frac{1}{2}(A^{2} + 2AB + B^{2})

+\frac{1}{6}(A^{3} + 3A^{2}B + 3AB^{2} + B^{3}) + · · · .  (15)

If we compare the terms of degree 2 we get

\frac{1}{2}(A + B)^{2} = \frac{1}{2}(A^{2} + AB + BA + B^{2})  (16)

in (14) and not \frac{1}{2}(A^{2} +2AB +B^{2}). Harmony is restored if A and B commute: indeed AB = BA entails

A^{2} + AB + BA + B^{2} = A^{2} + 2AB + B^{2}  (17)

and more generally the binomial formula

(A + B)^{n} = ∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) A^{i}B^{n−i}  (18)

for any n ≥ 0. By summation one gets

exp(A + B) = exp A · exp B  (19)

if A and B commute, but not in general. The success in (11) comes from the obvious fact that hD commutes to h′D since numbers commute to (linear) operators.

2.3 Leibniz’s formula

Leibniz’s formula for the higher order derivatives of the product of two functions is the following one

D^{n}(fg) = ∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) D^{i}f.D^{n−i}g.  (20)

The analogy with the binomial theorem is striking and was noticed early. Here are possible explanations. For the shift operator, we have

T_{h} = exp hD  (21)

by Taylor’s formula and

T_{h}(fg) = T_{h}f · T_{h}g  (22)

by an obvious calculation. Combining these formulas we get

∑_{n≥0} \frac{1}{n!}h^{n}D^{n}(fg) = ∑_{i≥0} \frac{1}{i!}h^{i}D^{i}f. ∑_{j≥0} \frac{1}{j!}h^{j}D^{j}g;  (23)

equating the terms containing the same power h^{n} of h, one gets

D^{n}(fg) = ∑_{i+j=n} \frac{n!}{i!j!}D^{i}f · D^{j}g  (24)

that is, Leibniz’s formula.

Another explanation starts from the case n = 1, that is

D(fg) = Df · g + f · Dg.  (25)

In a heuristic way it means that D applied to a product fg is the sum of two operators D_{1} acting on f only and D_{2} acting on g only. These actions being independent, D_{1} commutes to D_{2} hence the binomial formula

D^{n} = (D_{1} + D_{2})^{n} = ∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) D_{1}^{i} · D_{2}^{n−i}.  (26)

By acting on the product fg and observing that D_{1}^{i} · D_{2}^{j} transforms fg into D^{i}f ·D^{j}g, one recovers Leibniz’s formula. In more detail, to calculate D^{2}(fg), one applies D to D(fg). Since D(fg) is the sum of two terms Df · g and f · Dg apply D to Df · g to get D(Df)g + Df · Dg and to f · Dg to get Df · Dg + f · D(Dg), hence the sum

D(Df) · g + Df · Dg + Df · Dg + f · D(Dg)
= D^{2}f · g + 2Df · Dg + f · D^{2}g.

This last proof can rightly be called “formal” since we act on the formulas, not on the objects: D_{1} transforms f ·g into Df ·g but this doesn’t mean that from the equality of functions f_{1} · g_{1} = f_{2} · g_{2} one gets Df_{1} · g_{1} = Df_{2} · g_{2} (counterexample: from fg=gf, we cannot infer Df ·g = Dg ·f). The modern explanation is provided by the notion of tensor products: if V and W are two vector spaces (over the real numbers as coefficients, for instance), equal or distinct, there exists a new vector space V ⊗ W whose elements are formal finite sums ∑_{i} λ_{i}(v_{i} ⊗ w_{i}) (with scalars λ_{i} and v_{i} in V, w_{i} in W); we take as basic rules the consequences of the fact that v ⊗ w is bilinear in v, w, but nothing more. Taking V and W to be the space C^{∞}(I) of the functions defined and indefinitely differentiable in an interval I of R, we define the operators D_{1} and D_{2} in C^{∞}(I) ⊗ C^{∞}(I) by

D_{1}(f ⊗ g) = Df ⊗ g, D_{2}(f ⊗ g) = f ⊗ Dg.  (27)

The two operators D_{1}D_{2} and D_{2}D_{1} transform f ⊗ g into Df ⊗ Dg, hence D_{1} and D_{2} commute. Define D̄ as D_{1} + D_{2} hence

D̄(f ⊗ g) = Df ⊗ g + f ⊗ Dg.  (28)

We can now calculate D̄^{n} = (D_{1} + D_{2})^{n} by the binomial formula as in (26) with the conclusion

D̄^{n}(f ⊗ g) = ∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) D^{i}f ⊗ D^{n−i}g.  (29)

The last step is to go from (29) to (20). The rigorous reasoning is as follows. There is a linear operator μ taking f ⊗ g into f · g and mapping C^{∞}(I) ⊗ C^{∞}(I) into C^{∞}(I); this follows from the fact that the product f · g is bilinear in f and g. The formula (25) is expressed by D ∘ μ = μ ∘ D̄ in operator terms, according to the diagram:

C^{∞}(I) ⊗ C^{∞}(I) ⟶^{μ} C^{∞}(I)
D̄ ↓  ↓ D
C^{∞}(I) ⊗ C^{∞}(I) ⟶^{μ} C^{∞}(I).

An easy induction entails D^{n} ∘ μ = μ ∘ D̄^{n}, and from (29) one gets

D^{n}(fg) = D^{n}(μ(f ⊗ g)) = μ(D̄^{n}(f ⊗ g))
= μ(∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) D^{i}f ⊗ D^{n−i}g) = ∑_{i=0}^{n} ( \begin{matrix}n \\ i \end{matrix} ) D^{i}f · D^{n−i}g.  (30)

In words: first replace the ordinary product f · g by the neutral tensor product f ⊗ g, perform all calculations using the fact that D_{1} commutes with D_{2}, then restore the product · in place of ⊗.

When the vector spaces V and W consist of functions of one variable, the tensor product f ⊗ g can be interpreted as the function f(x)g(y) in two variables x, y; moreover D_{1} = ∂/∂x, D_{2} = ∂/∂y and μ takes a function F (x, y) of two variables into the one-variable function F (x, x) hence f(x)g(y) into f(x)g(x) as it should. Formula (25) reads now

\frac{∂}{∂x}(f(x)g(x)) = (\frac{∂}{∂x} + \frac{∂}{∂y} ) f(x)g(y)|_{y=x}.  (31)

The previous “formal” proof is just a rephrasing of a familiar proof using Schwarz’s theorem that \frac{∂}{∂x} and \frac{∂}{∂y} commute.

Starting from the tensor product H_{1} ⊗ H_{2} of two vector spaces, one can iterate and obtain

H_{1} ⊗ H_{2} ⊗ H_{3}, H_{1} ⊗ H_{2} ⊗ H_{3} ⊗ H_{4}, ….

Using once again the exponential notation, H^{⊗n} is the tensor product of n copies of H, with elements of the form ∑ λ · (ψ_{1} ⊗ … ⊗ ψ_{n}). In quantum physics, H represents the state vectors of a particle, and H^{⊗n} represents the state vectors of a system of n independent particles of the same kind. If H is an operator in H representing for instance the energy of a particle, we define n operators H_{i} in H^{⊗n} by

H_{i}(ψ_{1} ⊗ … ⊗ ψ_{n}) = ψ_{1} ⊗ ⋯ ⊗ Hψ_{i} ⊗ ⋯ ⊗ ψ_{n}  (32)

(the energy of the i-th particle). Then H_{1}, …, H_{n} commute pairwise and H_{1} + ⋯ + H_{n} is the total energy if there is no interaction. Usually, there is a pair interaction represented by an operator V in H ⊗ H; then the total energy is given by ∑_{i=1}^{n} H_{i} + ∑_{i<j} V_{ij} with

V_{12}(ψ_{1} ⊗ ψ_{2} ⊗ ⋯ ⊗ ψ_{n}) = V (ψ_{1} ⊗ ψ_{2}) ⊗ ψ_{3} ⊗ ⋯  (33)
V_{23}(ψ_{1} ⊗ ⋯ ⊗ ψ_{n}) = ψ_{1} ⊗ V (ψ_{2} ⊗ ψ_{3}) ⊗ ⋯ ⊗ ψ_{n}  (34)

etc… There are obvious commutation relations like

H_{i}H_{j} = H_{j}H_{i}
H_{i}V_{jk} = V_{jk}H_{i} if i, j, k are distinct.

This is the so-called “locality principle”: if two operators A and B refer to disjoint collections of particles (a) for A and (b) for B, they commute.

Faddeev and his collaborators made an extensive use of this notation in their study of quantum integrable systems. Also, Hirota introduced his so-called bilinear notation for differential operators connected with classical integrable systems (solitons).

2.4 Exponential vs. logarithm

In the case of real numbers, one usually starts from the logarithm and invert it to define the exponential (called antilogarithm not so long ago). Positive numbers have a logarithm; what about the logarithm of −1 for instance?

Things are worse in the complex domain. For a complex number z, define its exponential by the convergent series

exp z = ∑_{n≥0} \frac{1}{n!}z^{n}.  (35)

From the binomial formula, using the commutativity zz′ = z′z one gets

exp(z + z′) = exp z · exp z′  (36)

as before. Separating real and imaginary part of the complex number z = x + iy gives Euler’s formula

exp(x + iy) = e^{x}(cos y + i sin y)  (37)

subsuming trigonometry to complex analysis. The trigonometric lines are the “natural” ones, meaning that the angular unit is the radian (hence sin δ ≃ δ for small δ).

From an intuitive view of trigonometry, it is obvious that the points of a circle of equation x^{2} + y^{2} = R^{2} can be uniquely parametrized in the form

x = R cos θ, y = R sin θ  (38)

with −π < θ ≤ π, but the subtle point is to show that the geometric definition of sin θ and cos θ agree with the analytic one given by (37). Admitting this, every complex number u ≠ 0 can be written as an exponential exp z_{0}, where z_{0} = x_{0} + iy_{0}, x_{0} real and y_{0} in the interval ] − π, π]. The number z_{0} is called the principal determination of the logarithm of u, denoted by Ln(u). But the general solution of the equation exp z = u is given by z = z_{0} + 2πin where n is a rational integer. Hence a nonzero complex number has infinitely many logarithms. The functional property (36) of the exponential cannot be neatly inverted: for the logarithms we can only assert that Ln(u_{1} ⋯ u_{p}) and Ln(u_{1}) + … + Ln(u_{p}) differ by the addition of an integral multiple of 2πi.

The exponential of a (real or complex) square matrix A is defined by the series

exp A = ∑_{n≥0} \frac{1}{n!}A^{n}.  (39)

There are two classes of matrices for which the exponential is easy to compute:

a) Let A be diagonal A = diag(a_{1}, …, a_{n}). Then exp A is diagonal with elements exp a_{1}, …, exp a_{n}. Hence any complex diagonal matrix with non zero diagonal elements is an exponential, hence admits a logarithm, and even infinitely many ones.

b) Suppose that A is a special upper triangular matrix, with zeroes on the diagonal, of the type

A = ( \begin{matrix}0 & a & b & c \\ & 0 & d & e \\ & & 0 & f \\ & & & 0 \end{matrix} ) .

Then A^{d} = 0 if A is of size d × d. Hence exp A is equal to I + B where B is of the form A + \frac{1}{2}A^{2} + \frac{1}{6}A^{3} + ⋯ + \frac{1}{(d−1)!}A^{d−1}. Hence B is again a special upper triangular matrix and A can be recovered by the formula

A = B − \frac{B^{2}}{2} + \frac{B^{3}}{3} − ⋯ + (−1)^{d} \frac{B^{d−1}}{d − 1}.  (40)

This is just the truncated series for ln(I + B) (notice B^{d} = 0). Hence in the case of these special triangular matrices, exponential and logarithm are inverse operations.

In general, A can be put in triangular form A = UT U^{−1} where T is upper triangular. Let λ_{1}, …, λ_{d} be the diagonal elements of T, that is the eigenvalues of A. Then

exp A = U · exp T · U^{−1}  (41)

where exp T is triangular, with the diagonal elements exp λ_{1}, …, exp λ_{d}. Hence

det(exp A) = ∏_{i=1}^{d} exp λ_{i} = exp ∑_{i=1}^{d} λ_{i} = exp(Tr(A)).  (42)

The determinant of exp A is therefore nonzero. Conversely any complex matrix M with a nonzero determinant is an exponential: for the proof, write M in the form UT U^{−1} where T is composed of Jordan blocks of the form

T_{s} = ( \begin{matrix}λ & 1 & . & . & 0 \\ & . & . & . & . \\ & 0 & & . & 1 \\ . & . & . & . & λ \end{matrix}) with λ ≠ 0 .

From the existence of the complex logarithm of λ and the study above of triangular matrices, it follows that T_{s} is an exponential, hence T and M = UT U^{−1} are exponentials.

Let us add a few remarks:

a) A complex matrix with nonzero determinant has infinitely many logarithms; it is possible to normalize things to select one of them, but the conditions are rather artificial.

b) A real matrix with nonzero determinant is not always the exponential

of a real matrix; for example, choose M = ( \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ). This is not surprising since −1 has no real logarithm, but many complex logarithms of the form kπi with k odd.

c) The noncommutativity of the multiplication of matrices implies that in general exp(A + B) is not equal to exp A · exp B . Here the logarithm of a product cannot be the sum of the logarithms, whatever normalization we choose.

2.5 Infinitesimals and exponentials

There are many notations in use for the higher order derivatives of a function f. Newton uses ḟ, f̈, …, the customary notation is f′, f′′, . . .. Once again, the exponential notation can be systematized, f^{(m)} or D^{m}f denoting the m-th derivative of f, for m = 0, 1, . . .. This notation emphasizes that the derivation is a functional operator, hence

(f^{(m)})^{(n)} = f^{(m+n)}, or D^{m}(D^{n}f) = D^{m+n}f.  (43)

In this notation, it is cumbersome to write the chain rule for the derivative of a composite function

D(f ∘ g) = (Df ∘ g) · Dg.  (44)

Leibniz’s notation for the derivative is dy/dx if y = f(x). Leibniz was never able to give a completely rigorous definition of the infinitesimals dx, dy, …^{1}. His explanation of the derivative is as follows: starting from x, increment it by an infinitely small amount dx; then y = f(x) is incremented by dy, see Figure 1.

f(x + dx) = y + dy.  (45)

Then the derivative is f′(x) = dy/dx, hence according to (45),

f(x + dx) = f(x) + f′(x)dx.  (46)

This cannot be literally true, otherwise the function f(x) would be linear. The true formula is

f(x + dx) = f(x) + f′(x)dx + o(dx)  (47)

1 In modern times, Abraham Robinson has vindicated them using the tools of formal logic. There has been many interesting applications of his nonstandard analysis, but one has to admit that it remains too cumbersome to provide a viable alternative to the standard analysis. Maybe in the 21st century!

Fig. 1. Geometrical description: an infinitely small portion of the curve y = f(x), after zooming, becomes infinitely close to a straight line, our function is “smooth”, not fractal-like.

with an error term o(dx) which is infinitesimal, of a higher order than dx, meaning o(dx)/dx is again infinitesimal. In other words, the derivative f′(x), independent of dx, is infinitely close to \frac{f(x+dx)−f(x)}{dx} for all infinitesimals dx. The modern definition, as well as Newton’s point of view of fluents, is a dynamical one: when dx goes to 0, \frac{f(x+dx)−f(x)}{dx} tends to the limit f′(x). Leibniz’s notion is statical: dx is a given, fixed quantity. But there is a hierarchy of infinitesimals: η is of higher order than ϵ if η/ϵ is again infinitesimal. In the formulas, equality is always to be interpreted up to an infinitesimal error of a certain order, not always made explicit.

We use these notions to describe the logarithm and the exponential. By definition, the derivative of ln x is \frac{1}{x}, hence

\frac{d ln x}{dx} = \frac{1}{x}, that is ln(x + dx) = ln(x) + \frac{dx}{x}.

Similarly for the exponential

\frac{d exp x}{dx} = exp x, that is exp(x + dx) = (exp x)(1 + dx).

This is a rule of compound interest. Imagine a fluctuating daily rate of interest, namely ϵ_{1}, ϵ_{2}, …, ϵ_{365} for the days of a given year, every daily rate being of the order of 0.0003. For a fixed investment C, the daily reward is Cϵ_{i} for day i, hence the capital becomes C + Cϵ_{1} + … + Cϵ_{365} = C · (1 + ∑_{i} ϵ_{i}), that is approximately C(1+.11). If we reinvest every day our profit, invested capital changes according to the rule:

C_{i+1} =  C_{i}  + C_{i}ϵ_{i} = C_{i}(1 + ϵ_{i}).
↑  ↑  ↑
capital  capital  profit
at day i + 1 at day i during day i

At the end of the year, our capital is C · ∏_{i}(1 + ϵ_{i}). We can now formulate the “bankers rule”:

if S = ϵ_{1} + … + ϵ_{N}, then exp S = (1 + ϵ_{1}) ⋯ (1 + ϵ_{N}). (B)

Here N is infinitely large, and ϵ_{1}, …, ϵ_{N} are infinitely small; in our example, S = 0.11, hence exp S = 1 + S + \frac{1}{2}S^{2} + … is equal to 1.1163 …: by reinvesting daily, the yearly profit of 11% is increased to 11.63%.

Formula (B) is not true without reservation. It certainly holds if all ϵ_{i} are of the same sign, or more generally if ∑_{i} |ϵ_{i}| is of the same order as ∑ ϵ_{i} = x. For a counter-example, take N = 2p^{2} with half of the ϵ_{i} being equal to +\frac{1}{p}, and the other half to −\frac{1}{p} (hence ∑_{i} ϵ_{i} = 0 while ∏_{i}(1 + ϵ_{i}) is infinitely close to 1/e = exp(−1)).

To connect definition (B) of the exponential to the power series expansion exp S = 1 + S + \frac{1}{2!}S^{2} + ⋯ one can proceed as follows: by algebra we get

∏_{i=1}^{N}(1 + ϵ_{i}) = ∑_{k=0}^{N} S_{k},  (48)

where S_{0} = 1, S_{1} = ϵ_{1} + … + ϵ_{N} = S, and generally

S_{k} = ∑_{i_{1}<...<i_{k}} ϵ_{i_{1}} ⋯ ϵ_{i_{k}}.  (49)

We have to compare S_{k} to \frac{1}{k!}S^{k} = \frac{1}{k!}(ϵ_{1} + ⋯ + ϵ_{N})^{k}. Developing the kth power of S by the multinomial formula, we obtain S_{k} plus error terms each containing at least one of the ϵ_{i}’s to a higher power, ϵ_{i}^{2}, ϵ_{i}^{3}, …, hence infinitesimal compared to the ϵ_{i}’s. The general principle of compensation of errors^{2} is as follows: given a sum of infinitesimals

Σ = η_{1} + ⋯ + η_{M}  (50)

2 This terminology was coined by Lazare Carnot in 1797. Our formulation is more precise than his!

and new summands η′_{j} = η_{j} + o(η_{j}) with an error o(η_{j}) of higher order than η_{j}, we obtain that

Σ′ = η′_{1} + ⋯ + η′_{M}  (51)

is equal to Σ plus an error term o(η_{1}) + ⋯ + o(η_{M}). If the η_{j} are of the same sign, the error is o(Σ), that is negligible compared to Σ.

Fig. 2. Leibniz’ continuum: by zooming, a finite segment of line is made of a large number of atoms of space: a fractal.

The implicit view of the continuum underlying Leibniz’s calculus is as follows: a finite segment of a line is made of an infinitely large number of geometric atoms of space which can be arranged in a succession, each atom x being separated by dx from the next one. Hence in the definition of the logarithm

ln a = ∫_{1}^{a} \frac{dx}{x} (for a > 1),  (52)

we really have ∑_{1≤x≤a} \frac{dx}{x}. Similarly, the bankers rule (B) should be interpreted as

exp a = ∏_{0≤x≤a}(1 + dx) (for a > 0).  (53)

2.6 Differential equations

The previous formulation of the exponential suggests a method to solve a differential equation, for instance y′ = ry. In differential form

dy = r(x)ydx,  (54)

that is

y + dy = (1 + r(x)dx)y.  (55)

The solution is

y(b) = ∏_{a≤x≤b}(1 + r(x)dx) · y(a).  (56)

What is the meaning of this product? Putting ϵ(x) = r(x)dx, an infinitesimal, and expanding the product as in (48), we get

∏_{x}(1 + ϵ(x)) = ∑_{k≥0} ∑_{a≤x_{1}<...<x_{k}≤b} ϵ(x_{1}) ⋯ ϵ(x_{k});  (57)

reinterpreting the multiple sum as a multiple integral, this is

∑_{k≥0} ∫ ⋯ ∫_{Δ_{k}} r(x_{1}) ⋯ r(x_{k})dx_{1} ⋯ dx_{k}.  (58)

The domain of integration Δ_{k} is given by the inequalities

a ≤ x_{1} ≤ x_{2} ≤ … ≤ x_{k} ≤ b.  (59)

The classical solution to the differential equation y′ = ry is given by

y(b) = (exp ∫_{a}^{b} r(x)dx) · y(a).  (60)

Let us see how to go from (58) to (60). Geometrically, consider the hypercube C_{k} given by

a ≤ x_{1} ≤ b, ⋯, a ≤ x_{k} ≤ b  (61)

in the euclidean space R^{k} of dimension k with coordinates x_{1}, …, x_{k}. The group S_{k} of the permutations σ of {1, …, k} acts on R^{k}, by transforming the vector x with coordinates x_{1}, …, x_{k} into the vector σ.x with coordinates x_{σ^{−1}(1)}, …, x_{σ^{−1}(k)}. Then the cube C_{k} is the union of the k! transforms σ(Δ_{k}). Since the function r(x_{1}) … r(x_{k}) to be integrated is symmetrical in the variables x_{1}, …, x_{k} and moreover two distinct domains σ(Δ_{k}) and σ′(Δ_{k}) overlap by a subset of dimension < k, hence of volume 0, we see that the integral of r(x_{1}) ⋯ r(x_{k}) over C_{k} is k! times the integral over Δ_{k}. That is

∫ ⋯ ∫_{Δ_{k}} r(x_{1}) ⋯ r(x_{k})dx_{1} ⋯ dx_{k} =

\frac{1}{k!} ∫_{a}^{b} dx_{1} ⋯ ∫_{a}^{b} dx_{k} r(x_{1}) ⋯ r(x_{k}) = \frac{1}{k!}(∫_{a}^{b} r(x)dx)^{k}.

Summing over k, and using the definition of an exponential by a series, we conclude

∑_{k≥0} ∫ ⋯ ∫_{Δ_{k}} r(x_{1}) ⋯ r(x_{k})dx_{1} ⋯ dx_{k} = exp ∫_{a}^{b} r(x)dx.  (62)

as promised.

The same method applies to the linear systems of differential equations. We cast them in the matrix form

y′ = A · y,  (63)

that is the differential form

dy = A(x)ydx.  (64)

Here A(x) is a matrix depending on the variable x, and y(x) is a vector (or matrix) function of x. From (64) we get

y(x + dx) = (I + A(x)dx)y(x).  (65)

Formally the solution is given by

y(b) = ∏_{a≤x≤b}(I + A(x)dx) · y(a).  (66)

We have to take into account the noncommutativity of the products A(x)A(y)A(z) …. Explicitly, if we have chosen intermediate points

a = x_{0} < x_{1} < … < x_{N} = b,

with infinitely small spacing

dx_{1} = x_{1} − x_{0}, dx_{2} = x_{2} − x_{1}, …, dx_{N} = x_{N} − x_{N−1},

the product in (66) is

(I + A(x_{N})dx_{N})(I + A(x_{N−1})dx_{N−1}) ⋯ (I + A(x_{1})dx_{1}).

We use the notation \overleftarrow{∏}_{1≤i≤N}U_{i} for a reverse product U_{N}U_{N−1} ⋯ U_{1}; hence

the previous product can be written as \overleftarrow{∏}_{1≤i≤N}(I + A(x_{i})dx_{i}) and we should replace ∏ by \overleftarrow{∏} in equation (66). The noncommutative version of equation (48) is

\overleftarrow{∏}_{1≤i≤N}(I + A_{i}) = ∑_{k=0}^{N} ∑_{i_{1}>···>i_{k}} A_{i_{1}} ⋯ A_{i_{k}}.  (67)

Let us define the resolvant (or propagator) as the matrix

U(b, a) = \overleftarrow{∏}_{a≤x≤b}(I + A(x)dx).  (68)

Hence the differential equation dy = A(x)ydx is solved by y(b) = U(b, a)y(a) and from (67) we get

U(b, a) = ∑_{k≥0} ∫ ⋯ ∫_{Δ_{k}} A(x_{k}) ⋯ A(x_{1})dx_{1} ⋯ dx_{k}  (69)

with the factors A(x_{i}) in reverse order

A(x_{k}) ⋯ A(x_{1}) for x_{1} < … < x_{k}.  (70)

One owes to R. Feynman and F. Dyson (1949) the following notational trick. If we have a product of factors U_{1}, ⋯, U_{N}, each attached to a point x_{i} on a line, we denote by T (U_{1} ⋯ U_{N}) (or more precisely by \overleftarrow{T}(U_{1} ⋯ U_{N})) the product U_{i_{1}} ⋯ U_{i_{N}} where the permutation i_{1} … i_{N} of 1 … N is such that x_{i_{1}} > ⋯ > x_{i_{N}}. Hence in the rearranged product the abscisses attached to the factors increase from right to left. We argue now as in the proof of (62) and conclude that

∫ ⋯ ∫_{Δ_{k}} A(x_{k}) ⋯ A(x_{1})dx_{1} ⋯ dx_{k}
= \frac{1}{k!} ∫_{a}^{b} dx_{1} ⋯ ∫_{a}^{b} dx_{k}T (A(x_{1}) ⋯ A(x_{k})).  (71)

We can rewrite the propagator as

U(b, a) = T exp ∫_{a}^{b} A(x)dx,  (72)

with the following interpretation:

a) First use the series exp S = ∑_{k≥0} \frac{1}{k!}S^{k} to expand exp ∫_{a}^{b} A(x)dx. b) Expand S^{k} = (∫_{a}^{b} A(x)dx)^{k} as a multiple integral

∫_{a}^{b} dx_{1} ⋯ ∫_{a}^{b} dx_{k} A(x_{1}) ⋯ A(x_{k}).

c) Treat T as a linear operator commuting with series and integrals, hence

T exp S = ∑_{k≥0} \frac{1}{k!}T (S^{k}) = ∑_{k≥0} \frac{1}{k!}T \{∫_{a}^{b} dx_{1} ⋯ ∫_{a}^{b} dx_{k} A(x_{1}) ⋯ A(x_{k})\}

= ∑_{k≥0} \frac{1}{k!} ∫_{a}^{b} dx_{1} ⋯ ∫_{a}^{b} dx_{k} T (A(x_{1}) ⋯ A(x_{k})).

We give a few properties of the T (or time ordered) exponential:

a) Parallel to the rule

∫_{a}^{c} A(x)dx = ∫_{a}^{b} A(x)dx + ∫_{b}^{c} A(x)dx (for a < b < c)  (73)

we get

T exp ∫_{a}^{c} A(x)dx = T exp ∫_{b}^{c} A(x)dx · T exp ∫_{a}^{b} A(x)dx.  (74)

Notice that, in (73), the two matrices

L = ∫_{a}^{b} A(x)dx, M = ∫_{b}^{c} A(x)dx

do not commute, hence exp(L+M) is in general different from exp L·exp M. Hence formula (74) is not in general valid for the ordinary exponential.

b) The next formula embodies the classical method of “variation of constants” and is known in the modern literature as a “gauge transformation”. It reads as

S(b) · T exp ∫_{a}^{b} A(x)dx · S(a)^{−1} = T exp ∫_{a}^{b} B(x)dx  (75)

with

B(x) = S(x)A(x)S(x)^{−1} + S′(x)S(x)^{−1},  (76)

where S(x) is an invertible matrix depending on the variable x. The general formula (75) can be obtained by “taking a continuous reverse product” \overleftarrow{∏}_{a≤x≤b} over the infinitesimal form

S(x + dx)(I + A(x)dx))S(x)^{−1} = I + B(x)dx  (77)

(for the proof, write S(x + dx) = S(x) + S′(x)dx and neglect the terms proportional to (dx)^{2}). We leave it as an exercise to the reader to prove (75) from the expansion (69) for the propagator.

c) There exists a complicated formula for the T -exponential T exp ∫_{a}^{b} A(x) dx when A(x) is of the form \frac{A_{1}(x)+A_{2}(x)}{2}. Neglecting terms of order (dx)^{2}, we get

I + A(x)dx = (I + A_{2}(x)\frac{dx}{2}) (I + A_{1}(x)\frac{dx}{2})  (78)

and we can then perform the product \overleftarrow{∏}_{a≤x≤b}. This formula is the foundation of the multistep method in numerical analysis: starting from the value y(x) at time x of the solution to the equation y′ = Ay, we split the infinitesimal interval [x, x + dx] into two parts

I_{1} = [x, x + \frac{dx}{2}], I_{2} = [x + \frac{dx}{2}, x + dx];

we move at speed A_{1}(x)y(x) during I_{1} and then at speed A_{2}(x)y(x + \frac{dx}{2}) during I_{2}. Let us just mention one corollary of this method, the so-called Trotter-Kato-Nelson formula:

exp(L + M) = lim_{n→∞}(exp(L/n) exp(M/n))^{n}.  (79)

d) If the matrices A(x) pairwise commute, the T -exponential of ∫_{a}^{b} A(x)dx is equal to the ordinary exponential. In the general case, the following formula holds

T exp ∫_{a}^{b} A(x)dx = exp V (b, a)  (80)

where V (b, a) is explicitly calculated using integration and iterated Lie brackets. Here are the first terms

V (b, a) = ∫_{a}^{b} A(x)dx + \frac{1}{2} ∫ ∫_{Δ_{2}}[A(x_{2}), A(x_{1})]dx_{1}dx_{2}

+\frac{1}{3} ∫ ∫ ∫_{Δ_{3}}[A(x_{3}), [A(x_{2}), A(x_{1})]]dx_{1}dx_{2}dx_{3}  (81)

−\frac{1}{6} ∫ ∫ ∫_{Δ_{3}}[A(x_{2}), [A(x_{3}), A(x_{1})]]dx_{1}dx_{2}dx_{3} + · · · .

The higher-order terms involve integrals of order k ≥ 4. As far as I can ascertain, this formula was first enunciated by K. Friedrichs around 1950 in his work on the foundations of Quantum Field Theory. A corollary is the Campbell-Hausdorff formula:

exp L · exp M =
exp(L + M + \frac{1}{2}[L, M] + \frac{1}{12}[L, [L, M]] + \frac{1}{12}[M, [M, L]] + ⋯). (82)

It can be derived from (80) by putting a = 0, b = 2, A(x) = M for 0 ≤ x ≤ 1 and A(x) = L for 1 ≤ x ≤ 2.

The T -exponential found lately numerous geometrical applications. If C is a curve in a space of arbitrary dimension, the line integral ∫_{C} A_{μ}(x)dx^{μ} is well-defined and the corresponding T -exponential

T exp ∫_{C} A_{μ}(x)dx^{μ}  (83)

is closely related to the parallel transport along the curve C.

3 Operational calculus

3.1 An algebraic digression: umbral calculus

We first consider the classical Bernoulli numbers. I claim that they are defined by the equation

(B + 1)^{n} = B^{n} for n ≥ 2,  (1)

together with the initial condition B^{0} = 1. The meaning is the following: expand (B + 1)^{n} by the binomial theorem, then replace the power B^{k} by B_{k}. Hence (B + 1)^{2} = B^{2} gives B^{2} + 2B^{1} + B^{0} = B^{2}, that is after lowering the indices B_{2} + 2B_{1} + B_{0} = B_{2}, that is 2B_{1} + B_{0} = 0. Treating (B + 1)^{3} = B^{3} in a similar fashion gives 3B_{2} + 3B_{1} + B_{0} = 0. We write the first equations of this kind

n = 2 2B_{1} + B_{0} = 0
n = 3 3B_{2} + 3B_{1} + B_{0} = 0
n = 4 4B_{3} + 6B_{2} + 4B_{1} + B_{0} = 0
n = 5 5B_{4} + 10B_{3} + 10B_{2} + 5B_{1} + B_{0} = 0.

Starting from B_{0} = 1 we get successively

B_{1} = −\frac{1}{2}, B_{2} = \frac{1}{6}, B_{3} = 0, B_{4} = −\frac{1}{30}, …

Using the same kind of formalism, define the Bernoulli polynomials by

B_{n}(X) = (B + X)^{n}.  (2)

According to the previous rule, we first expand (B + X)^{n} using the binomial theorem, then replace B^{k} by B_{k}. Hence we get explicitly

B_{n}(X) = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) B_{n−k}X^{k}.  (3)

Since \frac{d}{dX} (X + c)^{n} = n(X + c)^{n−1} for any c independent of X, we expect

\frac{d}{dX} B_{n}(X) = nB_{n−1}(X).  (4)

This is easy to check on the explicit definition (3). Here is a similar calculation

(B + (X + Y))^{n} = ((B + X) + Y)^{n} = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) (B + X)^{n−k}Y^{k},

from which we expect to find

B_{n}(X + Y) = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) B_{n−k}(X)Y^{k}.  (5)

Indeed from (4) we get

(\frac{d}{dX} )^{k} B_{n}(X) = \frac{n!}{(n − k)!}B_{n−k}(X)  (6)

by induction on k, hence (5) follows from Taylor’s formula B_{n}(X + Y) = ∑_{k≥0} \frac{1}{k!}(\frac{d}{dX})^{k}B_{n}(X)Y^{k}.

We deduce now a generating series for the Bernoulli numbers. Formally

(e^{S} − 1)e^{BS} = e^{S}e^{BS} − e^{BS} = e^{(B+1)S} − e^{BS}

= ∑_{n≥0} \frac{1}{n!}S^{n}((B + 1)^{n} − B^{n}) = S((B + 1)^{1} − B^{1}) = S.

Since e^{BS} = ∑_{n≥0} \frac{1}{n!}B^{n}S^{n}, we expect

∑_{n≥0} B_{n}S^{n}/n! = \frac{S}{e^{S} − 1}.  (7)

Again this can be checked rigorously.

What is the secret behind these calculations?

We consider functions F (B, X, …) depending on a variable B and other variables X, . . .. Assume that F (B, X, …) can be expanded as a polynomial or power series in B, namely

F (B, X, …) = ∑_{n≥0} B^{n}F_{n}(X, …).  (8)

Then the “mean value” with respect to B is defined by

< F (B, X, …) >= ∑_{n≥0} B_{n}F_{n}(X, …),  (9)

where the B_{n}’s are the Bernoulli numbers: this corresponds to the rule “lower the index in B^{n}”. If the function F (B, X, …) can be expanded into a series ∑_{i} F_{i}(B, X, …)G_{i}(X, …) where the G_{i}’s are independent of B, then obviously^{3}

< F (B, X, …) >= ∑_{i} < F_{i}(B, X, …) > G_{i}(X, …).  (10)

The formal calculations given above are justified by this simple rule which affords also a probabilistic interpretation (see Section 3.7).

The previous method is loosely described as “umbral calculus”. We insisted on speaking of “mean values” to keep in touch with physical applications. From a purely mathematical point of view, it is just applying a linear functional acting on polynomials in B, mapping B^{n} into B_{n} for all n’s.

3.2 Binomial sequences of polynomials

These are sequences of polynomials U_{0}(X), U_{1}(X), … in one variable X satisfying the following relations:

3 So far we considered only identities linear in the B_{n}’s. If we want to treat nonlinear terms, like products B_{m} ·B_{n}, we need to introduce two independent symbols B and B′ and use the umbral rule to replace B^{m}B′^{n} by B_{m}B_{n}. In probabilistic terms (see Section 3.7), we introduce two independent random variables and take the mean value with respect to both simultaneously.

a) U_{0}(X) is a constant;

b) for any n ≥ 1, one gets

\frac{d}{dX} U_{n}(X) = nU_{n−1}(X).  (11)

By induction on n it follows that U_{n}(X) is of degree ≤ n. The binomial sequence is normalized if furthermore U_{0}(X) = 1, in which case every U_{n}(X) is a monic polynomial of degree n, that is

U_{n}(X) = X^{n} + c_{1}X^{n−1} + … + c_{n}.

Applying Taylor’s formula as above (derivation of formula (5)), one gets

U_{n}(X + Y) = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) U_{n−k}(X)Y^{k}.  (12)

We introduce now a numerical sequence by u_{n} = U_{n}(0) for n ≥ 0. Putting X = 0 in (12) and then replacing Y by X (as a variable), we get

U_{n}(X) = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) u_{n−k}X^{k}.  (13)

Conversely, given any numerical sequence u_{0}, u_{1}, …, and defining the polynomials U_{n}(X) by (13), one derives immediately the relations

\frac{d}{dX} U_{n}(X) = nU_{n−1}(X), U_{n}(0) = u_{n}.  (14)

The exponential generating series for the constants u_{n} is given by

u(S) = ∑_{n≥0} u_{n}S^{n}/n!.  (15)

From (13), one obtains the exponential generating series

U(X, S) = ∑_{n≥0} U_{n}(X)S^{n}/n!

for the polynomials U_{n}(X), namely in the form

U(X, S) = u(S)e^{XS}.  (16)

This could be expected. Writing ∂_{X}, ∂_{S} … for the partial derivatives, the basic relation ∂_{X}U_{n} = nU_{n−1} translates as (∂_{X} −S)U(X, S) = 0 or equivalently as

∂_{X}(e^{−XS}U(X, S)) = 0.  (17)

Hence e^{−XS}U(X, S) depends only on S, and putting X = 0 we obtain the value U(0, S) = u(S).

The umbral calculus can be successfully applied to our case. Hence U_{n}(X) can be interpreted as ⟨(X +U)^{n}⟩ provided ⟨U^{n}⟩ = u_{n}. Similarly u(S) is equal to ⟨e^{US}⟩ and U(X, S) to ⟨e^{(X+U)S}⟩. The symbolic derivation of (16) is as follows

U(X, S) = ⟨e^{(X+U)S}⟩ = ⟨e^{XS} · e^{US}⟩ = e^{XS}⟨e^{US}⟩ = e^{XS}u(S).

We describe in more detail the three basic binomial sequences of polynomials:

a) The sequence I_{n}(X) = X^{n} satisfies obviously (11). In this (rather trivial) case, we get

i_{0} = 1, i_{1} = i_{2} = … = 0, I(S) = 1, I(X, S) = e^{XS}.

b) The Bernoulli polynomials obey the rule (11) (see formula (4)). I claim that they are characterized by the normalization B_{0}(X) = 1 and the further property

∫_{0}^{1} B_{n}(x)dx = 0 for n ≥ 1.  (18)

Indeed, introducing the exponential generating series

B(X, S) = ∑_{n≥0} B_{n}(X)S^{n}/n!,  (19)

the requirement (18) is equivalent to the integral formula

∫_{0}^{1} B(x, S)dx = 1.  (20)

According to the general theory of binomial sequences, B(X, S) is of the form b(S)e^{XS}, hence

∫_{0}^{1} B(x, S)dx = ∫_{0}^{1} b(S)e^{xS}dx = b(S) (\frac{e^{S} − 1}{S}) .

Solving (20) we get b(S) = S/(e^{S} − 1) and from (7) this is the exponential generating series for the Bernoulli numbers. The exponential generating series for the Bernoulli polynomials is therefore

B(X, S) = \frac{Se^{XS}}{e^{S} − 1}.  (21)

Here is a short table:

B_{0}(X) = 1

B_{1}(X) = X − \frac{1}{2}

B_{2}(X) = X^{2} − X + \frac{1}{6}

B_{3}(X) = X^{3} − \frac{3}{2}X^{2} + \frac{1}{2}X.

c) We come to the Hermite polynomials which form the normalized binomial sequence of polynomials characterized by

∫_{−∞}^{+∞} H_{n}(x)dγ(x) = 0 for n ≥ 1,  (22)

where dγ(x) denotes the normal probability law, that is

dγ(x) = (2π)^{−1/2}e^{−x^{2}/2}dx.  (23)

We follow the same procedure as for the Bernoulli polynomials. Hence for the exponential generating series

H(X, S) = ∑_{n≥0} H_{n}(X)S^{n}/n! = h(S)e^{XS}  (24)

we get

∫_{−∞}^{+∞} H(x, S)dγ(x) = 1,  (25)

that is

1/h(S) = ∫_{−∞}^{+∞} e^{xS}dγ(x).  (26)

The last integral being easily evaluated, we conclude

h(S) = e^{−S^{2}/2}.  (27)

From this relation, we can evaluate H(X, S) namely

H(X, S) = e^{XS−S^{2}/2} = e^{X^{2}/2}e^{−(X−S)^{2}/2}  (28)

and using Taylor’s expansion for e^{−(X−S)^{2}/2}, we get

H_{n}(X) = (−1)^{n}e^{X^{2}/2} (\frac{d}{dX} )^{n} e^{−X^{2}/2}.  (29)

In the spirit of operator calculus, use the identity

e^{X^{2}/2} \frac{d}{dX} e^{−X^{2}/2} = \frac{d}{dX} − X,  (30)

whence

H_{n}(X) = (X − \frac{d}{dX} )^{n} .1.  (31)

This is tantamount to a recursion formula

H_{n+1}(X) = XH_{n}(X) − \frac{d}{dX} H_{n}(X) = XH_{n}(X) − nH_{n−1}(X).  (32)

The following table is then easily derived:

H_{0}(X) = 1
H_{1}(X) = X
H_{2}(X) = X^{2} − 1
H_{3}(X) = X^{3} − 3X
H_{4}(X) = X^{4} − 6X^{2} + 3.

3.3 Transformation of polynomials

We use the standard notation C[X] to denote the vector space of polynomials in the variable X with complex coefficients. Since the monomials X^{n} form a basis of C[X], a linear operator U : C[X] → C[X] is completely determined by the sequence of polynomials U_{n}(X) defined as the image U[X^{n}] of X^{n} under U. Here are a few examples:

I identity operator I_{n}(X) = X^{n}
D derivation \frac{d}{dX} D_{n}(X) = nX^{n−1}
T_{c} translation operator T_{c,n}(X) = (X + c)^{n}.

Notice that in general T_{c} transforms a polynomial P (X) into P (X + c) and Taylor’s formula amounts to

T_{c} = e^{cD};  (33)

furthermore T_{0} = I. From the definition of the derivative, one gets

D = lim_{c→0}(T_{c} − I)/c.  (34)

We can reformulate the properties of binomial sequences:

- the definition DU_{n}(X) = nU_{n−1}(X) amounts to DU = UD;

- the exponential generating series U(X, S) is nothing else than U[e^{XS}];

- formula (12), after substituting c to Y reads as

U_{n}(X + c) = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) U_{n−k}(X)c^{k}

that is

T_{c}U[X^{n}] = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) U[X^{n−k}]c^{k} = U[(X + c)^{n}] = UT_{c}[X^{n}].

Hence this formula expresses that U commutes to T_{c}

T_{c}U = UT_{c};  (35)

- formula (13) can be rewritten as

U[X^{n}] = ∑_{k≥0} \frac{1}{k!}u_{k}D^{k}[X^{n}].  (36)

From the definition (15) of the exponential generating series, we obtain

U = u(D).  (37)

To sum up, our operators are characterized by the following equivalent properties:

a) U commutes with the derivative D;

b) U commutes with the translation operators T_{c};

c) U can be expressed as a power series u(D) in D. Furthermore, since D acts on e^{XS} by multiplication by S, the operator U = u(D) multiplies e^{XS} by u(S). Hence we recover formula (16).

3.4 Expansion formulas

As we saw before, B_{n}(X) and H_{n}(X) are monic polynomials and therefore the sequences (B_{n}(X))_{n≥0} and (H_{n}(X))_{n≥0} are two bases of the vector space C[X]. Hence an arbitrary polynomial P (X) can be expanded as a linear combination of the Bernoulli polynomials, as well as of the Hermite polynomials. Our aim is to give explicit formulas.

Consider a general binomial sequence (U_{n}(X))_{n≥0} such that u_{0} ≠ 0, with exponential generating series U(X, S) = u(S)e^{XS}. Introduce the inverse series v(S) = 1/u(S); explicitly

v(S) = ∑_{n≥0} v_{n}S^{n}/n!,

and the coefficients v_{n} are defined inductively by

v_{0} = 1/u_{0}, v_{n} = −\frac{1}{u_{0}} ∑_{k=1}^{n} ( \begin{matrix}n \\ k \end{matrix} ) u_{k}v_{n−k}.  (38)

In the spirit of umbral calculus, let us define the linear form ϕ_{0} on C[X] by ϕ_{0}[X^{n}] = v_{n}. I claim that the expansion of an arbitrary polynomial in terms of the U_{n}’s is given by

P (X) = ∑_{n≥0} \frac{1}{n!}ϕ_{0}[D^{n}P] · U_{n}(X).  (39)

Before giving a proof, let us examine the three basic examples:

a) If U_{n}(X) = X^{n}, then u(S) = 1, hence v(S) = 1. That is v_{0} = 1 and v_{n} = 0 for n ≥ 1. The linear form ϕ_{0} is given by ϕ_{0}[P] = P (0) and formula (39) reduces to MacLaurin’s expansion

P (X) = ∑_{n≥0} \frac{1}{n!}(D^{n}P)(0) · X^{n}.  (40)

b) For the Bernoulli polynomials we know that the series 1/b(S) is equal to (e^{S} −1)/S, hence v_{n} = \frac{1}{n+1}. The linear form ϕ_{0} is defined by ϕ_{0}[X^{n}] = \frac{1}{n+1}, that is

ϕ_{0}[P] = ∫_{0}^{1} P (x)dx.  (41)

Hence

P (X) = ∑_{n≥0} \frac{1}{n!} (∫_{0}^{1} D^{n}P (x)dx) · B_{n}(X).  (42)

c) In the case of Hermite polynomials, we know that 1/h(S) is equal to e^{S^{2}/2}, hence

v_{2m} = \frac{(2m)!}{m!2^{m}}, v_{2m+1} = 0.  (43)

According to (26), we get v_{n} = ∫_{−∞}^{+∞} x^{n}dγ(x), hence

ϕ_{0}[P] = ∫_{−∞}^{+∞} P (x)dγ(x) = (2π)^{−1/2} ∫_{−∞}^{+∞} P (x)e^{−x^{2}/2}dx.  (44)

In these three cases, the formula for ϕ_{0}[P] takes a similar form, namely

ϕ_{0}[P] = ∫_{a}^{b} P (x)w(x)dx  (45)

with the following prescriptions:

a = −∞, b = +∞, w(x) = δ(x) in case a),

a = 0, b = 1, w(x) = 1 in case b),

a = −∞, b = +∞, w(x) = (2π)^{−1/2}e^{−x^{2}/2} in case c). The normalization ϕ_{0}[1] = 1 amounts to ∫_{a}^{b} w(x)dx = 1, that is w(x) is the probability density of a random variable taking values in the interval [a, b] (see Section 3.7).

There is a peculiarity in case c). Namely, according to the general formula (39), an arbitrary polynomial P (X) can be expanded in a series of Hermite polynomials ∑_{n≥0} c_{n}H_{n}(X) where c_{n} is equal to \frac{1}{n!} ∫_{−∞}^{+∞} D^{n}P (x)dγ(x). Integrating by parts and taking into account the definition (29) of H_{n}(X), we obtain

c_{n} = \frac{1}{n!} ∫_{−∞}^{+∞} P (x)H_{n}(x)dγ(x).  (46)

This amounts to the orthogonality relation

∫_{−∞}^{+∞} H_{m}(x)H_{n}(x)dγ(x) = δ_{mn}n!  (47)

for the Hermite polynomials. There is no such orthogonality relation for the Bernoulli polynomials.

One final word about the proof of (39). By linearity, it suffices to consider the case P = U_{m}, that is to prove the biorthogonality relation

ϕ_{0}[D^{n}U_{m}] = n!δ_{mn}.  (48)

We first calculate ϕ_{0}[U_{m}]. From formula (13), we obtain

ϕ_{0}[U_{m}] = ∑_{k=0}^{m} ( \begin{matrix}m \\ k \end{matrix} ) u_{m−k}ϕ_{0}[X^{k}] = ∑_{k=0}^{m} ( \begin{matrix}m \\ k \end{matrix} ) u_{m−k}v_{k},

and from (38), ϕ_{0}[U_{m}] is 0 for m ≥ 1. Since D^{n}U_{m} is proportional to U_{m−n} according to the basic formula DU_{m} = mU_{m−1}, one gets ϕ_{0}[D^{n}U_{m}] = 0 for m ≠ n. Finally D^{m}U_{m} = m!u_{0}, hence ϕ_{0}[D^{m}U_{m}] = v_{0}m!u_{0} = m!.

3.5 Signal transforms

A transmission device transforms a suitable input f into an output F. Both are evolving in time and are represented by functions of time f(t) and F (t).^{4} We assume the device to be linear and in a stationary regime, that is there is a linear operator V taking f(t) into F (t) (linearity) and f(t + τ) into F (t + τ) for any fixed τ (stationarity).

Here are the main types of response:

\underline{Input}  \underline{Output}
δ(t) I(t)
e^{pt} Θ(t)e^{pt}
t^{n} V_{n}(t)

In the first case, δ(t) is a Dirac singular function, that is a pulse, and I(t) is the impulse response. By stationarity, V transforms δ(t − τ) into I(t − τ); an arbitrary input f can be represented as a superposition of pulses

f(t) = ∫_{−∞}^{+∞} f(τ)δ(t − τ)dτ  (49)

hence by linearity the output

F (t) = ∫_{−∞}^{+∞} f(τ)I(t − τ)dτ = ∫_{−∞}^{+∞} f(t − τ)I(τ)dτ.  (50)

4 For simplicity, we restrict to the case where input and output are scalars and not vectors.

In the non-anticipating case, I(t) is zero before the pulse δ(t) occurs, that is I(t) = 0 for t < 0. In this case the output is given by

F (t) = ∫_{−∞}^{t} f(τ)I(t − τ)dτ.  (51)

In the case of the exponential input f(t) = e^{pt}, the output is equal to ∫_{−∞}^{+∞} e^{pτ}I(t − τ)dτ = ∫_{−∞}^{+∞} e^{p(t−τ)}I(τ)dτ according to (50), that is, to Θ(p)e^{pt} with the spectral gain

Θ(p) = ∫_{−∞}^{+∞} e^{−pτ}I(τ)dτ.  (52)

We can give an a priori argument: f_{p}(t) = e^{pt} is a solution of the differential equation Df_{p} = pf_{p}; since the operator V is stationary, that is, commutes with the translation operators T_{c}, it commutes with D = lim_{c→0}(T_{c} − I)/c. Hence the output F_{p} corresponding to the input f_{p} is a solution of the differential equation DF_{p} = pF_{p}, hence is proportional to e^{pt}.

In a similar way, the monomials t^{n} satisfy the cascade of differential equations

D[t] = 1, D[t^{2}] = 2t, D[t^{3}] = 3t^{2}, …

Since V commutes with D and the constants are the solutions of the differential equation D(f) = 0, it follows that the images V_{n}(t) = V[t^{n}] form a binomial sequence of polynomials. Explicitly

V_{n}(t) = ∫_{−∞}^{+∞}(t − τ)^{n}I(τ)dτ = ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) v_{n−k}t^{k}

with the constants

v_{n} = (−1)^{n} ∫_{−∞}^{+∞} I(τ)τ^{n}dτ = V_{n}(0).  (53)

Comparing (52) to (53) we conclude

Θ(p) = ∑_{n≥0} v_{n}p^{n}/n!.  (54)

More generally, since e^{pt} is equal to ∑_{n≥0} \frac{1}{n!}p^{n}t^{n}, application of the linear operator V gives

V[e^{pt}] = ∑_{n≥0} \frac{1}{n!}p^{n}V [t^{n}],  (55)

that is

Θ(p)e^{pt} = ∑_{n≥0} \frac{1}{n!}p^{n}V_{n}(t).  (56)

Up to the change in notation (p for S, and t for X), the spectral gain Θ(p) is nothing else than the numerical exponential generating series associated to the binomial sequence (V_{n}(t))_{n≥0}.

Comparing with the results obtained in Section 3.3, it is tempting to write the operator V as Θ(D). According to (54), Θ(D) can be interpreted as ∑_{n≥0} v_{n}D^{n}/n!, but it is known that infinite order differential operators are not so easily dealt with. A better interpretation is obtained via Laplace or Fourier transform. Indeed since D multiplies e^{pt} by p, any function F (D) ought to multiply e^{pt} by F (p), and the rule V[e^{pt}] = Θ(p)e^{pt} is in agreement with the interpretation V = Θ(D). If the input can be represented as a Laplace transform

f(t) = ∫ e^{pt}ϕ(p)dp,  (57)

then V = Θ(D) transforms it into the output

F (t) = ∫ e^{pt}Θ(p)ϕ(p)dp.  (58)

Similarly, if the input f(t) is given by its spectral resolution (or Fourier transform)

f(t) = ∫_{−∞}^{+∞} f̂(ω)e^{iωt}dω,  (59)

then the output is given by

F (t) = ∫_{−∞}^{+∞} Θ(iω)f̂(ω)e^{iωt}dω.  (60)

This is Heaviside’s magic trilogy:

symbolic p operator D spectral iω = 2πiν (ν frequency, ω = 2πν pulsation)

Θ(p) ⟷ Θ(D) ⟷ Θ(2πiν).

Recall that in the Laplace transform, p is a complex variable, while in the Fourier transform ω and ν are real variables (see example in the next section).

3.6 The inverse problem

This is the problem of recovering the input, knowing the output. In operator terms, we have to compute the inverse U of the operator V (if it exists!). Since V is stationary, so is U, and at the level of polynomial inputs and outputs, U corresponds to a binomial sequence of polynomials U_{n}(t):

\{ \begin{matrix}input & U_{n}(t) \\ output & t^{n}. \end{matrix}

Together with the numerical sequence v_{n} = V_{n}(0), we have to consider the numerical sequence u_{n} = U_{n}(0). Introducing the exponential generating series

u(S) = ∑_{n≥0} \frac{1}{n!}u_{n}S^{n}, v(S) = ∑_{n≥0} \frac{1}{n!}v_{n}S^{n} = Θ(S),  (61)

we can write U = u(D) and V = v(D) at least when acting on polynomials. Since U and V are inverse operators, we expect the relation u(S)v(S) = 1, equivalent to the chain of relations

u_{0}v_{0} = 1, ∑_{k=0}^{n} ( \begin{matrix}n \\ k \end{matrix} ) u_{k}v_{n−k} = 0 for n ≥ 1,  (62)

to hold. Indeed, this is easily checked (see Section 3.4, formula (38)).

Since the input U_{n}(t) corresponds to the output t^{n}, a Taylor-MacLaurin expansion of the output corresponds to an expansion of the input in terms of the U_{n}(t). Fix an epoch t_{0} and use the Taylor expansion of the output

F (t) = ∑_{n≥0} \frac{1}{n!}D^{n}F (t_{0})(t − t_{0})^{n}.  (63)

Applying the operator U, we get

f(t) = ∑_{n≥0} \frac{1}{n!}D^{n}F (t_{0}) · U_{n}(t − t_{0}),  (64)

since U transforms (t − t_{0})^{n} into U_{n}(t − t_{0}) by stationarity. The reader is invited to compare this formula to formula (39).

We give one illustrating example. Let the output be a moving average of the input

F (t) = ∫_{t−1}^{t} f(s)ds,  (65)

Fig. 3. The impulse response corresponding to (65)

corresponding to the impulse response shown in Figure 3.

The spectral gain is

Θ(p) = ∫_{0}^{1} e^{−pτ}dτ = \frac{1 − e^{−p}}{p} .  (66)

The inverse series u(p) = 1/Θ(p) is given by

u(p) = \frac{p}{1 − e^{−p}} .  (67)

Since u(−p) = \frac{p}{e^{p}−1} is the exponential generating series of the Bernoulli numbers, the polynomials U_{n}(t) are easily identified,

U_{n}(t) = (−1)^{n}B_{n}(−t) = B_{n}(t) + nt^{n−1} = B_{n}(t + 1).  (68)

Notice that Θ(p) vanishes for p ≠ 0 of the form p = 2πin with an integer n; equivalently, the inverse function u(p) has poles for p ≠ 0, p = 2πin. Hence not every output is admissible since (65) entails ∑_{n} F (t + n) = ∫_{−∞}^{+∞} f(s)ds. That is, an output satisfies the necessary (and sufficient) condition

∑_{n} F (t + n) = c (c a constant)  (69)

and the input f(t) can be reconstructed from the output F (t) up to the addition of a function f_{0}(t) with

f_{0}(t) = f_{0}(t + 1), ∫_{0}^{1} f_{0}(t)dt = 0.  (70)

Exercise a) Derive from (64) the relation

f(t) = ∑_{n≥0} \frac{1}{n!}u_{n}D^{n}F (t),  (71)

for a general transmission device, where 1/Θ(p) = ∑_{n≥0} u_{n}p^{n}/n!.

b) In the particular case (65), one gets u_{n} = B_{n}(1), hence u_{n} = B_{n} if n ≥ 2, and u_{0} = 1, u_{1} = \frac{1}{2}.

c) Deduce from relation (7) that B_{n} = 0 for n ≥ 3, n odd.

d) Derive the Euler-MacLaurin summation formula

\frac{1}{2}(f(t) + f(t − 1))

= ∫_{t−1}^{t} f(s)ds + ∑_{m≥1} \frac{1}{(2m)!}B_{2m}[D^{2m−1}f(t) − D^{2m−1}f(t − 1)]. (72)

3.7 A probabilistic application

We consider a random variable ξ. In general, we denote by ⟨X⟩ the mean value of a random variable X. We want to define a probabilistic version of the so-called Wick Powers in Quantum Field Theory.

The goal is to associate to ξ a sequence of random variables : ξ^{n} : such that

a) the mean value of : ξ^{n} : is 0 for n ≥ 1;

b) there exists a normalized binomial sequence of polynomials Π_{n}(X) such that : ξ^{n} := Π_{n}(ξ).

Let w(x) be the probability density associated to ξ, hence w(x) ≥ 0 and ∫_{−∞}^{+∞} w(x)dx = 1. Moreover, for any (non random) function f(x) of a real variable x, the random variable f(ξ) has a mean value given by

⟨f(ξ)⟩ = ∫_{−∞}^{+∞} f(x)w(x)dx.  (73)

Hence the conditions a) and b) amount to

0 = ⟨Π_{n}(ξ)⟩ = ∫_{−∞}^{+∞} Π_{n}(x)w(x)dx for n ≥ 1.  (74)

Using the same method as in Section 3.2, we introduce the exponential generating series Π(X, S) = ∑_{n≥0} Π_{n}(X)S^{n}/n!, hence the relation (74) translates into

∫_{−∞}^{+∞} Π(x, S)w(x)dx = 1.  (75)

Putting π(S) = Π(0, S), hence Π(X, S) = π(S)e^{XS}, we derive

1/π(S) = ∫_{−∞}^{+∞} e^{xS}w(x)dx.  (76)

We translate these relations into probabilistic jargon: replace S by p and X by ξ to get

1/π(p) = ⟨e^{pξ}⟩  (77)
Π(ξ, p) = ∑_{n≥0} \frac{1}{n!}p^{n} : ξ^{n} :  (78)
Π(ξ, p) = π(p)e^{pξ}.  (79)

Extending the definition of : : by linearity to e^{pξ} = ∑_{n≥0} \frac{1}{n!}p^{n}ξ^{n}, we rewrite (78) as Π(ξ, p) =: e^{pξ} :. Here is the conclusion

: e^{pξ} : = \frac{e^{pξ}}{⟨e^{pξ}⟩}.  (80)

Let us specialize our results in the case of the binomial sequences considered so far:

a) If ξ = 0, then ⟨e^{pξ}⟩ = 1, hence : e^{pξ} := e^{pξ} = 1. That is : ξ^{n} := 0 for n ≥ 1.

b) Suppose that ξ is uniformly distributed in the interval [0, 1], that is w(x) = 1 if 0 ≤ x ≤ 1, and w(x) = 0 otherwise. Then

⟨e^{pξ}⟩ = ∫_{0}^{1} e^{px}dx = \frac{e^{p} − 1}{p}.  (81)

We get

∑_{n≥0} \frac{1}{n!}p^{n} : ξ^{n} : = : e^{pξ} : = \frac{pe^{pξ}}{e^{p} − 1},  (82)

40 Pierre Cartier

that is

: ξ^{n} : = B_{n}(ξ)  (83)

where B_{n}(X) is the Bernoulli polynomial of degree n. In particular

: ξ : = ξ − ⟨ξ⟩ = ξ − \frac{1}{2}

: ξ^{2} : = ξ^{2} − ξ + \frac{1}{6}

: ξ^{3} : = ξ^{3} − \frac{3}{2}ξ^{2} + \frac{1}{2}ξ, etc …

c) Assume now that ξ is normalized: ⟨ξ⟩ = 0, ⟨ξ^{2}⟩ = 1, and follows a Gaussian law. Then w(x) = (2π)^{−1/2}e^{−x^{2}/2} and

⟨e^{pξ}⟩ = (2π)^{−1/2} ∫_{−∞}^{+∞} e^{px−x^{2}/2}dx = e^{p^{2}/2}.  (84)

Reasoning as above, we obtain

: ξ^{n} : = H_{n}(ξ)  (85)

where H_{n}(X) is the Hermite polynomial of degree n. Explicitly : ξ : = ξ

: ξ^{2} : = ξ^{2} − 1
: ξ^{3} : = ξ^{3} − 3ξ.

To get a general formula, apply (80) to obtain the pair of relations

: e^{pξ} : = e^{−p^{2}/2}e^{pξ}, e^{pξ} = e^{p^{2}/2} : e^{pξ} : .  (86)

Equating equal powers of p, we derive

: ξ^{n} : = ∑_{0≤k≤n/2}(−1)^{k} \frac{n!}{2^{k}k!(n − 2k)!}ξ^{n−2k},  (87)

and conversely

ξ^{n} = ∑_{0≤k≤n/2} \frac{n!}{2^{k}k!(n − 2k)!} : ξ^{n−2k} : .  (88)

Notice that the orthogonality relation (47) for the Hermite polynomials translates in probabilistic terms into

⟨: ξ^{m} : : ξ^{n} :⟩ = m!δ_{mn},  (89)

hence the sequence 1,: ξ :,: ξ^{2} :, … is derived from the natural sequence 1, ξ, ξ^{2}, … by orthogonalization.

To conclude, we can use the reflected probability density w(−x) as an impulse response and define the input-output relation by

F (t) = ∫ f(t + τ)w(τ)dτ,  (90)

that is

F (t) = ⟨f(t + ξ)⟩  (91)

in probabilistic terms. The interpretation is that the input is spoiled by random delay in transmission. Then Π_{n}(t) is the input corresponding to the output t^{n}. Analytically this is expressed by

∫_{−∞}^{+∞} Π_{n}(t + τ)w(τ)dτ = t^{n},  (92)

and probabilistically by

⟨: (ξ + t)^{n} :⟩ = t^{n}.  (93)

3.8 The Bargmann-Segal transform

Let us consider again the input-output transformation in the Gaussian case. It is then called the Bargmann-Segal transform (or B-transform), denoted by B. According to (90), we have

Bf(z) = \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} f(x + z)e^{−x^{2}/2}dx  (94)

or

Bf(z) = \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} f(x)e^{−(z−x)^{2}/2}dx.  (95)

Comparing formulas (24) and (28), one obtains

e^{−(z−x)^{2}/2} = e^{−x^{2}/2} ∑_{n≥0} H_{n}(x)z^{n}/n!,  (96)

and by integrating term by term in the expression (95) one concludes

Bf(z) = ∑_{n≥0} Γ_{n}(f)z^{n}/n!  (97)

with

Γ_{n}(f) = ∫_{−∞}^{+∞} H_{n}(x)f(x)dγ(x).  (98)

Taking into account the orthogonality property (47) namely

∫_{−∞}^{+∞} H_{m}(x)H_{n}(x)dγ(x) = δ_{mn}n!,

one derives Γ_{n}(f) = n!c_{n} for f given by a series ∑_{n≥0} c_{n}H_{n}(x). That is, the B-transform takes ∑_{n≥0} c_{n}H_{n}(x) into ∑_{n≥0} c_{n}z^{n}.

To be more precise, we need to introduce some function spaces. The natural one is L^{2}(dγ) consisting of the (measurable) functions f(x) for which the integral ∫_{−∞}^{+∞} |f(x)|^{2}dγ(x) is finite; the scalar product is given by

⟨f_{1}|f_{2}⟩ = ∫_{−∞}^{+∞} \overline{f_{1}(x)}f_{2}(x)dγ(x).  (99)

In this space, the functions He_{n}(x) := H_{n}(x)/(n!)^{1/2} (for n = 0, 1, …) form an orthonormal basis^{5}. The B-transform takes this space onto the space of series ∑ c_{n}z^{n} with ∑_{n≥0} n!|c_{n}|^{2} < ∞.

In its original form (94), the transformation B requires z to be real, but the form (95) extends to the case of a complex number z. Indeed, from the property that ∑_{n≥0} n!|c_{n}|^{2} is finite, it follows that the series ∑_{n≥0} c_{n}z^{n} has an infinite radius of convergence, hence represents an entire function of the complex variable z. The space of such entire functions is denoted F(C) and called the Fock space (in one degree of freedom, see Section 3.9).

The elements of L^{2}(dγ) can be interpreted as the random variables of the form X = f(ξ) with ⟨|X|^{2}⟩ finite, where ξ is a normalized Gaussian random variable. We saw that B takes H_{n}(ξ) = : ξ^{n} : into z^{n}. Hence it is tempting to denote by : : the map inverse to B, so that : z^{n} : = H_{n}(x). We have a

5 The orthonormality condition ⟨He_{m}|He_{n}⟩ = δ_{mn} is nothing else than the orthogonality condition (47). But it requires a proof to show that this system is complete, that is that any function in the Hilbert space L^{2}(dγ) can be approximated by polynomials (in the norm convergence).

new kind of operational calculus

B

L^{2}(dγ) ⇀ ↽ F(C),

: :

where B transforms a function f(x) or the random variable X = f(ξ) into the entire function

BX(z) = e^{−z^{2}/2}⟨X · e^{zξ}⟩  (100)

and the inverse map takes an entire function Φ(z) = ∑_{n≥0} c_{n}z^{n} in the Fock space into the random variable : Φ(ξ) := ∑_{n≥0} c_{n} : ξ^{n} :.

According to the definition (94), B takes the function e^{px} into e^{pz+p^{2}/2}, that is it acts as e^{D^{2}/2} where D is the derivation, followed by the change of variable x into z. This result can be reformulated as follows. Using the exponential generating series

∑_{n≥0} H_{n}(x)p^{n}/n! = e^{px−p^{2}/2}  (101)

for the Hermite polynomials, and noting that e^{D^{2}/2} applied to e^{pz−p^{2}/2} gives e^{px}, we conclude that e^{D^{2}/2} takes H_{n}(x) into x^{n}, that is B coincides on the polynomials with the differential operator e^{D^{2}/2} of infinite degree.

We would like to conclude the general rule

B = e^{D^{2}/2},  (102)

hence

: : = e^{−D^{2}/2}.  (103)

One way to substantiate these claims is to consider the heat (or diffusion) equation

∂_{s}F (s, x) = \frac{1}{2}∂_{x}^{2}F (s, x)  (104)

with initial value

F (0, x) = f(x).  (105)

Since D = ∂_{x}, the solution of equation (104) can be written formally as F (s, x) = e^{sD^{2}/2}f(x), hence e^{D^{2}/2}f(x) represents the value for s = 1 of the solution of equation (104) which agrees for s = 0 with f(x). But we know an explicit solution to the heat equation

F (s, x) = \frac{1}{\sqrt{2πs}} ∫_{−∞}^{+∞} f(x + u)e^{−u^{2}/2s}du.  (106)

Comparing with (94), we obtain

Bf(x) = F (1, x),  (107)

and this relation is the true expression of B = e^{D^{2}/2}.

The operator e^{D^{2}/2} (or B) is smoothing. That is, if we simply assume that f belongs to L^{2}(dγ) (that is that the integral ∫_{−∞}^{+∞} |f(x)|^{2}e^{−x^{2}/2}dx is finite), then the function F (1, x) = e^{D^{2}/2}f(x) extends to an entire function in the complex domain. Conversely, the Wick operator : : = e^{−D^{2}/2} makes sense only for the functions g(x) (for x real) which extend in the complex domain into a function Φ(z) (for z complex) belonging to the Fock space F(C).

3.9 The quantum harmonic oscillator

Let us rewrite the definition of the B-transform as an integral operator

Bf(z) = ∫_{−∞}^{+∞} B(z, x)f(x)dx  (108)

with a kernel

B(z, x) = (2π)^{−1/2}e^{−(z−x)^{2}/2}.  (109)

It is often more convenient to replace the Hilbert space L^{2}(dγ) by the Hilbert space L^{2}(R) ^{6}. Defining the function u_{0}(x) by (2π)^{−1/4}e^{−x^{2}/4}, we get

dγ(x) = u_{0}(x)^{2}dx,  (110)

hence ∫ |f(x)|^{2}dγ(x) is finite if and only if ∫ |f(x)u_{0}(x)|^{2}dx is finite. That is, the multiplication by the function u_{0}(x) gives an isometry of L^{2}(dγ) onto

6 consisting of the (measurable) functions ϕ(x) such that ∫_{−∞}^{+∞} |ϕ(x)|^{2}dx be finite, with scalar product ⟨ϕ_{1}|ϕ_{2}⟩ = ∫_{−∞}^{+∞} \overline{ϕ_{1}(x)}ϕ_{2}(x)dx.

L^{2}(R). We can transfer the B-transform to L^{2}(R), as the isometry B′ of L^{2}(R) onto F(C) ^{7} defined by Bf = B′(fu_{0}). Explicitly

B′f(z) = ∫_{−∞}^{+∞} B′(z, x)f(x)dx  (111)

with

B′(z, x) = u_{0}(x)^{−1}B(z, x) = (2π)^{−1/4}e^{−z^{2}/2+zx−x^{2}/4}.  (112)

Many properties are easier to describe in the Fock space. For instance the function 1 is called the ground state Ω, the multiplication by z is called the creation operator, denoted by a^{∗}, and the derivation ∂_{z} = \frac{d}{dz} is the annihilation operator, denoted by a. The vectors

e_{n} = \frac{1}{\sqrt{n!}}(a^{∗})^{n}Ω,  (113)

that is the functions e_{n}(z) = \frac{1}{\sqrt{n!}}z^{n}, form an orthonormal basis of F(C) with e_{0} = Ω. An easy calculation gives

\{ \begin{matrix}ae_{n} = n^{1/2}e_{n−1} for n ≥ 1, ae_{0} = 0 \\ a^{∗}e_{n} = (n + 1)^{1/2}e_{n+1}, \end{matrix}  (114)

hence the matrices

a = ( \begin{matrix}0 & \sqrt{1} & 0 & 0 & . & . & . \\ 0 & 0 & \sqrt{2} & 0 & . & . & . \\ 0 & 0 & 0 & \sqrt{3} & . & . & . \\ 0 & 0 & 0 & 0 & . & . & . \\ . & . & . & . & . & . & . \\ . & . & . & . & . & . & . \end{matrix}) , a^{∗} = ( \begin{matrix}0 & 0 & 0 & 0 & . & . & . \\ \sqrt{1} & 0 & 0 & 0 & . & . & . \\ 0 & \sqrt{2} & 0 & 0 & . & . & . \\ 0 & 0 & \sqrt{3} & 0 & . & . & . \\ . & . & . & . & . & . & . \\ . & . & . & . & . & . & . \end{matrix})

in the basis (e_{n})_{n≥0}; it follows that a and a^{∗} are adjoint to each other. Moreover from the definitions a^{∗} = z, a = ∂_{z} follows the commutation relation

aa^{∗} − a^{∗}a = 1.  (115)

7 with the scalar product ⟨Φ_{1}|Φ_{2}⟩ = ∑_{n≥0} n!\overline{c}_{n,1}c_{n,2} for Φ_{j}(z) = ∑_{n≥0} c_{n,j}z^{n} (j = 1, 2).

Finally the number operator N = a^{∗}a is given by N = z∂_{z}, hence is diagonalized in the basis (e_{n})

Ne_{n} = ne_{n}.  (116)

In the spirit of operational calculus, we transfer these results from the Fock space model to the spaces L^{2}(dγ) and L^{2}(R). The following table summarizes these translations (where ∂_{x} is \frac{d}{dx}):

\begin{matrix}Space & L^{2}(dγ) & L^{2}(R) & F(C) \\ Ω & 1 & (2π)^{−1/4}e^{−x^{2}/4} = u_{0}(x) & 1 \\ e_{n} & (n!)^{−1/2}H_{n}(x) & (n!)^{−1/2}H_{n}(x)u_{0}(x) & (n!)^{−1/2}z^{n} \\ a^{∗} & x − ∂_{x} & x/2 − ∂_{x} & z \\ a & ∂_{x} & x/2 + ∂_{x} & ∂_{z} \\ N & x∂_{x} − ∂_{x}^{2} & −∂_{x}^{2} + x^{2}/4 − 1/2 & z∂_{z} \end{matrix}

For instance, the fact that a^{∗} corresponds to x − ∂_{x} in L^{2}(dγ) is proved as follows: from the definition of B(z, x) one gets

(x + ∂_{x})B(z, x) = zB(z, x).  (117)

Multiplying by f(x) and integrating by parts, we get

∫ B(z, x)(x − ∂_{x})f(x)dx = ∫ (x + ∂_{x})B(z, x)f(x)dx
= z ∫ B(z, x)f(x)dx,

that is B((x − ∂_{x})f) = zBf. The other cases are similar.

We apply these results to the harmonic oscillator. In classical mechanics, the harmonic oscillator is described by the Hamiltonian H = \frac{p^{2}}{2m} + \frac{Kq^{2}}{2} in canonical coordinates p,q. The equation of motion is q̈+ ω^{2}q = 0 with the pulsation ω = \sqrt{K/m}, and the momentum p = mq̇. To get the corresponding quantum Hamiltonian H, replace p by the operator p = −ih̄∂_{q} hence

H = −\frac{h̄^{2}}{2m}∂_{q}^{2} + \frac{mω^{2}q^{2}}{2} .  (118)

Introduce the dimensionless coordinate x = (2mω/h̄)^{1/2}q. Then H can be rewritten as

H = h̄ω(a^{∗}a + \frac{1}{2})  (119)

in the model L^{2}(R). From the diagonalization of N = a^{∗}a, we conclude that the energy levels of the quantum harmonic oscillator (that is, the eigenvalues of H) are given by h̄ω(n+ \frac{1}{2}) with n = 0, 1, 2, …: that is Planck’s radiation law, with the correction \frac{1}{2} giving \frac{1}{2}h̄ω for the energy of the ground state u_{0} = Ω.

4 The art of manipulating infinite series

4.1 Some divergent series

Euler claimed that S = 1−1+1−1+… is equal to \frac{1}{2}. Here is the purported proof:

S = 1 − 1 + 1 − 1 + … + S = 1 − 1 + 1 − … ———————————— 2 S = 1 + 0 + 0 + 0 + … = 1.

What is implicit is the use of two rules:

a) If S = u_{0} + u_{1} + u_{2} + …, then S = 0 + u_{0} + u_{1} + …

b) If S = u_{0} + u_{1} + u_{2} + … and S′ = u_{0} + u′_{1} + u′_{2} + …, then S + S′ = (u_{0} + u′_{0}) + (u_{1} + u′_{1}) + (u_{2} + u′_{2}) + …. These rules certainly hold for convergent series but to extend them to divergent series is somewhat hazardous.

Let us repeat the previous calculation in a slightly more general form:

S = 1 − t + t^{2} − t^{3} + … + tS = t − t^{2} + t^{3} − … ————————————— (1 + t)S = 1 + 0 + 0 + 0 + … = 1.

The result is

1 − t + t^{2} − t^{3} + … = \frac{1}{1 + t},  (1)

the classical summation of the geometric series. If t is a real number such that |t| < 1, the geometric series is convergent, and the use of rules a) and b) is justified. To get Euler’s result, take the limiting value t = 1 in (1).

What we need is the explicit description of various procedures to define rigorously the sum of certain divergent series (not all at once) and to compare these procedures. Suppose we want to define the sum

S = u_{0} + u_{1} + ….  (2)

Introduce weights p_{0,t}, p_{1,t}, … and the weighted series

S_{t} = p_{0,t}u_{0} + p_{1,t}u_{1} + ….  (3)

If the series S_{t} is convergent for each value of the parameter t, and S_{t} approaches a limit S when t approaches some limiting value t_{0}, then S is the sum for this procedure^{8}.

The previous procedure is reasonable only when lim_{t→t_{0}} p_{n,t} = 1 for n = 0, 1, . . .. Some examples:

a) p_{0,N} = p_{1,N} = … = p_{N,N} = 1, p_{n,N} = 0 for n > N and N = 0, 1, 2, . . .. Then the weighted sum amounts to the finite sum

S_{N} = u_{0} + … + u_{N}

(obviously convergent) and the convergence of S_{N} towards a limit S corresponds to the convergence of the series u_{0} + u_{1} + u_{2} + … in the standard sense, with the standard sum S.

b) Put σ_{N} = \frac{1}{N+1}(S_{0} + … + S_{N}); this corresponds to the weights

p_{n,N} = \{ \begin{matrix}1 − \frac{n}{N+1} & for 0 ≤ n ≤ N \\ 0 & for n > N. \end{matrix}  (4)

If σ_{N} converges to a limit σ, this is the Cesaro-sum of the series u_{0} + u_{1} + u_{2} + . . ..

8 See Knopp’s book [8] for this method.

c) To get the Abel summation, we introduce the weights p_{n,t} = t^{n} for n = 0, 1, 2, … and a real parameter t with 0 < t < 1. We take therefore the limit for t = 1 of the power series ∑_{n≥0} u_{n}t^{n}.

It is known that every convergent series with sum S is Cesaro-summable with the same sum σ = S. Similarly, Cesaro summation is extended by Abel summation. Euler’s example is u_{n} = (−1)^{n}, hence

S_{N} = \{ \begin{matrix}1 & if N ≥ 0 is even \\ 0 & if N ≥ 0 is odd \end{matrix}  (5)

and therefore

σ_{N} = \{ \begin{matrix}\frac{1}{2} & if N is odd \\ \frac{1}{2} + \frac{1}{2N+2} & if N is even. \end{matrix}  (6)

It follows that σ_{N} converges to σ = \frac{1}{2}. Hence the series 1 − 1 + 1 − 1 + … is Cesaro-summable to \frac{1}{2}, and a previous calculation shows that it is Abelsummable to \frac{1}{2} also.

The scope of Abel summation can be extended in various ways. For instance, if the sequence (u_{n}) is bounded, that is |u_{n}| ≤ M for n = 0, 1, 2, … with a constant M independent of n, then the power series ∑_{n≥0} u_{n}z^{n} converges for any complex number z with |z| < 1 and defines therefore a holomorphic function U(z) in the open disk |z| < 1 (see Fig. 4). If the limit lim_{r→1} U(re^{iθ})(for 0 ≤ r < 1) exists, it can be taken as an Abel sum for the series ∑_{n≥0} u_{n}e^{inθ}.

In a slightly more general way, we can assume that the sequence (u_{n}) is polynomially bounded, that is

|u_{n}| ≤ Cn^{k}

for all n = 1, 2, … and some constants C > 0 and k = 1, 2, … The radius of convergence of the series ∑_{n≥0} u_{n}z^{n} is still at least 1, and if U(1) = lim_{r→1} U(r) = lim_{r→1} ∑_{n≥0} u_{n}r^{n} exists, it is the Abel sum for u_{0}+u_{1}+u_{2}+…

We just give one example, namely u_{n} = (−1)^{n}n^{k} for k = 0, 1, 2, . . .. We calculate

U(z) = ∑_{n≥0} u_{n}z^{n} = ∑_{n≥0}(−1)^{n}n^{k}z^{n} = ∑_{n≥0} n^{k}(−z)^{n}

= ∑_{n≥0}(z∂_{z})^{k}(−z)^{n} = (z∂_{z})^{k} ∑_{n≥0}(−z)^{n} = (z∂_{z})^{k} \frac{1}{1 + z} .

50 Pierre Cartier

Fig. 4. The open unit disk

Particular cases:

k = 0, U(z) = \frac{1}{1 + z} , U(1) = \frac{1}{2}

k = 1, U(z) = \frac{−z}{(1 + z)^{2}} , U(1) = −\frac{1}{4}

k = 2, U(z) = \frac{z(z − 1)}{(1 + z)^{3}}, U(1) = 0

k = 3, U(z) = \frac{−z^{3} + 4z^{2} − z}{(1 + z)^{4}} , U(1) = \frac{1}{8},

that is

∑_{n≥0}(−1)^{n} = \frac{1}{2}

∑_{n≥0}(−1)^{n}n = −\frac{1}{4}

∑_{n≥0}(−1)^{n}n^{2} = 0

∑_{n≥0}(−1)^{4}n^{3} = \frac{1}{8}.

In general, we get U(1) = (z∂_{z})^{k} \frac{1}{1+z}|_{z=1} that is, after the change of variable z = e^{u},

∑_{n≥0}(−1)^{n}n^{k} = ∂_{u}^{k} \frac{1}{e^{u} + 1}|_{u=0}.  (7)

Using the exponential generating series for the Bernoulli numbers in the form

\frac{1}{e^{u} − 1} = \frac{1}{u} + ∑_{k≥0} \frac{B_{k+1}}{(k + 1)!}u^{k}

we obtain

\frac{1}{e^{u} + 1} = \frac{1}{e^{u} − 1} − \frac{2}{e^{2u} − 1} = ∑_{k≥0} \frac{(1 − 2^{k+1})B_{k+1}}{(k + 1)!} u^{k},  (8)

hence Euler’s result

∑_{n≥0}(−1)^{n}n^{k} = \frac{(1 − 2^{k+1})B_{k+1}}{k + 1} .  (9)

We leave it to the reader to rederive the previous cases 0 ≤ k ≤ 3 using the values for B_{1}, B_{2}, B_{3}, B_{4} given in Section 3.1. We come back to this result in Section 4.4.

Euler gave formulas for wildly divergent series like ∑_{n≥0}(−1)^{n}n!. Using the classical formula

n! = ∫_{0}^{∞} e^{−t}t^{n}dt  (10)

and assuming unrestrictedly as legitimate formal rules both term-by-term integration and summation of a geometric series, we get

∑_{n≥0}(−1)^{n}n! = ∑_{n≥0}(−1)^{n} ∫_{0}^{∞} e^{−t}t^{n}dt

= ∫_{0}^{∞} e^{−t} ∑_{n≥0}(−t)^{n} dt = ∫_{0}^{∞} \frac{e^{−t}}{1 + t}dt,

the last integral being convergent. This is just the beginning of the use of Borel transform and Borel summation for divergent series.

4.2 Polynomials of infinite degree and summation of series

It is an important principle that a polynomial can be reconstructed from its roots. More precisely, let

P (z) = c_{n}z^{n} + c_{n−1}z^{n−1} + … + c_{1}z + c_{0}  (11)

(with c_{n} ≠ 0) be a polynomial of degree n with complex coefficients. If λ_{1} is a root of P, that is P (λ_{1}) = 0, it is elementary to factorize P (z) = (z−λ_{1})P_{1}(z) where P_{1}(z) is a polynomial of degree n − 1. Continuing this process, we end up with a factorization

P (z) = (z − λ_{1}) … (z − λ_{m})Q(z),  (12)

where the polynomial Q(z) of degree n − m has no more roots. According to a highly non-trivial result, first stated by d’Alembert (1746) and proved by Gauß (1797), a polynomial without roots is a constant, hence the factorization (12) takes the form

P (z) = c_{n}(z − λ_{1}) … (z − λ_{n}).  (13)

By a well known calculation, one derives the following relations between coefficients and roots

λ_{1} + … + λ_{n} = −c_{n−1}/c_{n}

∑_{i<j} λ_{i}λ_{j} = c_{n−2}/c_{n}, etc …

For our purposes, it is better to use the inverses of the roots, assumed to be nonzero. Since the logarithmic derivative transforms product into sum and annihilates constants, we derive

DP (z)/P (z) = ∑_{i=1}^{n} \frac{1}{z − λ_{i}} .  (14)

Use of the geometric series gives

∑_{i=1}^{n} \frac{1}{z − λ_{i}} = − ∑_{i=1}^{n} ∑_{k≥0} z^{k}/λ_{i}^{k+1}.  (15)

Introducing the sums of inverse powers of roots

γ_{k} = ∑_{i=1}^{n} λ_{i}^{−k},  (16)

we conclude from this calculation

zDP (z) + P (z) ∑_{k≥1} γ_{k}z^{k} = 0.  (17)

Assuming for simplicity c_{0} = 1 and equating the coefficients of equal powers of z, we obtain the following variant of Newton’s relations

γ_{k} + c_{1}γ_{k−1} + … + c_{k−1}γ_{1} + kc_{k} = 0  (18)

for k ≥ 1. It is important to notice that the degree n of P (z) does not appear explicitly in the relation (18), which can be solved inductively

γ_{1} = −c_{1}  (19)
γ_{2} = c_{1}^{2} − 2c_{2}  (20)
γ_{3} = −c_{1}^{3} + 3c_{1}c_{2} − 3c_{3}  (21)
γ_{4} = c_{1}^{4} − 4c_{1}^{2}c_{2} + 4c_{1}c_{3} − 4c_{4} + 2c_{2}^{2}.  (22)

Around 1734, Euler undertook to calculate the sum of the series S_{2} = ∑_{n≥1} \frac{1}{n^{2}} . This series is slowly convergent, but Euler invented efficient acceleration methods for summing series and calculated the sum S_{2} = 1.64493406 …; he recognized S_{2} = π^{2}/6. He obtained also the value of S_{4} = ∑_{n≥1} 1/n^{4} to be π^{4}/90. To establish these results rigorously, he introduced the equation sin x = 0 admitting the solutions x = 0, ±π, ±2π, ±3π, … Discarding the root x = 0 and using the power series expansion of sin x, we are led to consider the equation

1 − \frac{x^{2}}{6} + \frac{x^{4}}{120} − … = 0

with roots ±π, ±2π, ±3π, . . .. With the previous notations we have

c_{1} = 0, c_{2} = −\frac{1}{6}, c_{3} = 0, c_{4} = \frac{1}{120}, …

γ_{2} = ∑_{n≥1} [\frac{1}{(πn)^{2}} + \frac{1}{(−πn)^{2}} ] = 2S_{2}/π^{2}

γ_{4} = ∑_{n≥1} [\frac{1}{(πn)^{4}} + \frac{1}{(−πn)^{4}} ] = 2S_{4}/π^{4}.

Assuming that the relations (20) and (22) still hold, we get

2S_{2}/π^{2} = γ_{2} = −2c_{2} = \frac{1}{3}

2S_{4}/π^{4} = γ_{4} = −4c_{4} + 2c_{2}^{2} = −\frac{1}{30} + \frac{1}{18} = \frac{1}{45}.

The desired relations

S_{2} = π^{2}/6, S_{4} = π^{4}/90

follow immediately.

To summarize the method used by Euler:

a) first guess the value from accurate numerical work;

b) consider the function

\frac{sin x}{x} = 1 − \frac{x^{2}}{6} + \frac{x^{4}}{120} − …

as a polynomial of infinite degree, with infinitely many roots ±π, ±2π, ±3π, …;

c) since Newton’s relations (19) to (22) do not involve explicitly the degree n of the polynomial, assume their validity in the case n = ∞ as well, and exploit them for P (x) = (sin x)/x.

4.3 The Euler-Riemann zeta function

We use Riemann’s definition and notation

ζ(s) = ∑_{n≥1} n^{−s}.  (23)

The series converges absolutely for any complex number s with real part ℜ(s) greater than 1. It has been shown by Riemann that ζ(s) can be analytically continued to the whole complex plane, the only singularity being a pole of order 1 at s = 1, that is ζ(s) − 1/(s − 1) is an entire function. Obviously ζ(1) = ∑_{n≥1} 1/n is a divergent series, but ζ(s) is defined when s ≠ 1 is an integer (positive or negative). Euler was the first to calculate ζ(s) when s is an integer.

We consider the case where s is even and strictly positive. Euler proved the formula

ζ(2k) = \frac{2^{2k−1}π^{2k}|B_{2k}|}{(2k)!},  (24)

and in particular

ζ(2) = \frac{2π^{2}|B_{2}|}{2!} = \frac{π^{2}}{6}  (25)

ζ(4) = \frac{8π^{4}|B_{4}|}{4!} = \frac{π^{4}}{90}.  (26)

Since ζ(2) = ∑_{n≥1} 1/n^{2} = S_{2} and similarly ζ(4) = S_{4}, we recover the formulas for S_{2} and S_{4}. The method we used in the previous section could be extended to cover the general case (24), but it is simpler to go back to the formula given for the logarithmic derivative in (14). For the function sin z, the logarithmic derivative is cot z = cos z/ sin z. This function is meromorphic in the whole complex plane, with simple poles of residue 1 at each integral multiple of π. Euler assumed at first that, in analogy with (14), cot z should be equal to the sum of its polar contributions, that is

cot z = ∑_{n=−∞}^{+∞} \frac{1}{z − nπ} .  (27)

Assume this relation for a moment, and derive (24). The series (27) is not absolutely convergent, but can be summed in a symmetrical way by taking ∑_{n=−∞}^{+∞} to be lim_{N→∞} ∑_{n=−N}^{+N}. Hence

cot z − \frac{1}{z} = ∑_{n=1}^{∞} \frac{2z}{z^{2} − n^{2}π^{2}} .  (28)

The right-hand side can be developed using the geometric series; for |z| < π, the series involved are absolutely convergent, hence

∑_{n=1}^{∞} \frac{2z}{z^{2} − n^{2}π^{2}} = − ∑_{n=1}^{∞} 2z ∑_{k=1}^{∞} \frac{(z^{2})^{k−1}}{(n^{2}π^{2})^{k}} = −2 ∑_{n=1}^{∞} ∑_{k=1}^{∞} \frac{z^{2k−1}}{n^{2k}π^{2k}}

= −2 ∑_{k=1}^{∞} \frac{z^{2k−1}}{π^{2k}} ∑_{n=1}^{∞} \frac{1}{n^{2k}} ,

that is

cot z = \frac{1}{z} − 2 ∑_{k≥1} \frac{ζ(2k)}{π^{2k}}z^{2k−1}.  (29)

Another use of the exponential generating series for the Bernoulli numbers yields

cot z = i\frac{e^{2iz} + 1}{e^{2iz} − 1} = \frac{2i}{e^{2iz} − 1} + i

= 2i\{\frac{1}{2iz} − \frac{1}{2} + ∑_{k≥1} \frac{B_{2k}}{(2k)!}(2iz)^{2k−1}\} + i,

hence finally

cot z = \frac{1}{z} + ∑_{k≥1} \frac{(−1)^{k}2^{2k}B_{2k}}{(2k)!} z^{2k−1}.  (30)

To establish (24), it is enough to compare the expansions (29) and (30) for cot z and to note that ζ(2k) = ∑_{n≥1} \frac{1}{n^{2k}} is the sum of a convergent series of positive numbers, hence ζ(2k) > 0.

Euler’s proof for the expansion (27) of cot z is reproduced in many textbooks. Here is a variant which seems to have been unnoticed so far. Define

Φ(z) = cot z − ∑_{n=−∞}^{+∞} \frac{1}{z − nπ} .  (31)

Examining the poles of cot z, we see that Φ(z) is an entire function of the complex variable z. A simple manipulation yields the functional equation

Φ(z) = \frac{1}{2} [Φ (\frac{z}{2}) + Φ (\frac{z + π}{2} )];  (32)

we have to prove that Φ(z) = 0 for all z.

a) The function Φ is bounded: indeed, denote by C_{n} the set of complex numbers whose modulus is at most (2^{n} + 1)π. Since C_{1} is a compact set and Φ is continuous, there exists a constant M > 0 such that |Φ(z)| ≤ M for z in C_{1}. Assuming the estimate |Φ(z)| ≤ M for z in C_{n}, we use the functional equation for z in C_{n+1},

|Φ(z)| ≤ \frac{1}{2} |Φ (\frac{z}{2})| + \frac{1}{2} |Φ (\frac{z + π}{2} )| ,  (33)

and observe that both z/2 and (z + π)/2 belong to C_{n}, hence |Φ(z/2)| ≤ M, |Φ((z + π)/2)| ≤ M; from (33) we conclude that |Φ(z)| ≤ M (for z in C_{n+1}). Every complex number belongs to some set C_{n}, hence |Φ(z)| ≤ M for all z.

b) We appeal now to Liouville’s theorem to conclude that Φ, being a bounded entire function, is a constant, hence Φ(z) = Φ(0).

c) The function Φ is odd, that is Φ(−z) = −Φ(z), hence Φ(0) = 0.

Liouville’s theorem, the main ingredient in this proof, was proved around 1850, a century after Euler worked on these questions. It is interesting to note that the d’Alembert-Gauß theorem is an easy corollary of Liouville’s theorem. [Hint: if P (z) is a polynomial without roots, the function Φ(z) = 1/P (z) is entire and bounded, hence a constant; that is, P (z) is a constant.]

4.4 Sums of powers of numbers

The other result of Euler about ζ(s) can be stated as follows:

1^{k} + 2^{k} + 3^{k} + … = −\frac{B_{k+1}}{k + 1}  (34)

for k = 1, 2, 3, … It looks at first odd, since it gives a finite value for an infinite sum of positive numbers, obviously divergent since each term is at least 1. Euler’s derivation is more or less as follows.

Formula (9) can be written as

1^{k} − 2^{k} + 3^{k} − 4^{k} + … = −(1 − 2 · 2^{k})\frac{B_{k+1}}{k + 1}.  (35)

On the other hand, multiply the left-hand side of (34) by 1 − 2 · 2^{k} and rearrange. This yields

(1 − 2 · 2^{k})(1^{k} + 2^{k} + 3^{k} + …) =
1^{k} + 2^{k} + 3^{k} + 4^{k} + 5^{k} + 6^{k} + …

−2( 2^{k} + 4^{k} + 6^{k} + …)

——————————————– = 1^{k} − 2^{k} + 3^{k} − 4^{k} + 5^{k} − 6^{k} + …

and finally (34) is obtained from (35).

This procedure is highly questionable, but can be fixed as follows. We introduce two functions

ζ(s) = ∑_{n≥1} n^{−s}, η(s) = ∑_{n≥1}(−1)^{n−1}n^{−s}.  (36)

Provided these functions can be continued analytically to the negative integers, formula (34) and (35) read respectively as

ζ(−k) = −\frac{B_{k+1}}{k + 1} (34′)

η(−k) = (2^{k+1} − 1)\frac{B_{k+1}}{k + 1} (35′)

for k = 1, 2, … Furthermore

η(s) = ∑_{n odd} n^{−s} − ∑_{n even} n^{−s} = ∑_{n≥1} n^{−s} − 2 ∑_{n even} n^{−s}

= ζ(s) − 2 ∑_{m≥1}(2m)^{−s}

and finally

η(s) = (1 − 2^{1−s})ζ(s).  (37)

Our manipulation of series is justified as long as ℜ(s) > 1, but the final formula remains valid for all s for which both ζ(s) and η(s) are regular (analytic continuation!). In particular

η(−k) = (1 − 2^{k+1})ζ(−k).  (38)

Hence formulas (34′) and (35′) are equivalent, substantiating Euler’s derivation.

Using the known values of the Bernoulli numbers, we deduce

ζ(−2) = ζ(−4) = ζ(−6) = … = 0

ζ(−1) = −\frac{1}{12}, ζ(−3) = \frac{1}{120}, ζ(−5) = −\frac{1}{252}, ….

It can also be shown that ζ(0) = −1/2. Hence we get the paradoxical results:

ζ(0) = 1 + 1 + 1 + … = −\frac{1}{2}

ζ(−1) = 1 + 2 + 3 + … = −\frac{1}{12}.

Among the many methods available to construct the analytical continuation of ζ(s), we select the following one using η(s). Indeed, from Euler’s definition of the gamma function

Γ (s) = ∫_{0}^{∞} e^{−t}t^{s−1}dt  (39)

(for ℜ(s) > 1), we deduce by a simple change of variable

Γ (s)n^{−s} = ∫_{0}^{∞} e^{−nt}t^{s−1}dt.  (40)

By summation and term by term integration

Γ (s)η(s) = ∑_{n≥1}(−1)^{n−1}Γ (s)n^{−s} = ∑_{n≥1}(−1)^{n−1} ∫_{0}^{∞} e^{−nt}t^{s−1}dt

= ∫_{0}^{∞} { ∑_{n≥1}(−1)^{n−1}e^{−nt}} · t^{s−1}dt,

that is

Γ (s)η(s) = ∫_{0}^{∞} \frac{t^{s−1}}{e^{t} + 1}dt.  (41)

All the calculations are justified as long as ℜ(s) > 1. We use now a general principle, unknown to Euler, to deal with integrals of the type

Φ(s) = ∫_{0}^{∞} F (t)t^{s−1}dt.  (42)

We split the integral as ∫_{0}^{1} + ∫_{0}^{∞}.

a) Assuming that F (t) decreases at infinity faster than any power t^{−k} (for k ≥ 0), then

Φ_{1,∞}(s) = ∫_{1}^{∞} F (t)t^{s−1}dt  (43)

extends to an entire function.

b) Assuming that F (t) is differentiable to any order on the closed interval [0, 1], then

Φ_{0,1}(s) = ∫_{0}^{1} F (t)t^{s−1}dt  (44)

extends to a meromorphic function in the complex plane C. The only singularities^{9} are at s = 0, −1, −2, … with singular part \frac{D^{k}F (0)}{k!(s+k)} around s = −k.

Applying this principle to the definition (39) of Γ (s) we recover the wellknown fact that Γ (s) extends to a meromorphic function, with poles at s = 0, −1, −2, … and singular part \frac{(−1)^{k}}{k!(s+k)} around s = −k. We use now formula

9 Hint: integrate by parts using \frac{d}{dt} t^{−s} = −st^{−s−1}.

(41) for Γ (s)η(s). Hence this function extends to a meromorphic function with poles at s = 0, −1, … and singular part \frac{c_{k}}{s+k} around s = −k, where

\frac{1}{e^{t} + 1} = ∑_{k≥0} c_{k}t^{k}.  (45)

According to (8), we get

c_{k} = \frac{(1 − 2^{k+1})B_{k+1}}{(k + 1)!} .  (46)

Dividing Γ (s)η(s) by Γ (s), the poles cancel; hence η(s) extends to an entire function, and comparing the singular parts of Γ (s)η(s) and Γ (s) around s = −k, we find

η(−k) = (−1)^{k}k!c_{k}.  (47)

We have to distinguish between several cases:

• k = 0 yields η(0) = c_{0} = −B_{1} = \frac{1}{2};

• k ≥ 2 is even yields η(−k) = 0 since B_{r} = 0 for r = k + 1 odd;

• k ≥ 1 is odd, then (−1)^{k} = −1 and

η(−k) = \frac{(2^{k+1} − 1)B_{k+1}}{k + 1} .  (48)

This is Euler’s formula (35) or formula (35′). The analytical continuation of ζ(s) can now be performed by using (37), that is we define

ζ(s) = \frac{η(s)}{1 − 2^{1−s}} .  (49)

Since η(s) is entire, the only singularity of ζ(s) is a pole at s = 1, with singular part \frac{1}{s−1}. We can now calculate ζ(−k) from η(−k) and get

ζ(−k) = −\frac{B_{k+1}}{k + 1}  (50)

for k = 1, 2, …, as expected. Furthermore, from η(0) = \frac{1}{2} we get the remaining value

ζ(0) = −η(0) = −\frac{1}{2}.  (51)

4.5 Variation I: Did Euler really fool himself?

Bourbaki wrote (in [2], page VI.29): “Mais la tendance au calcul formel est la plus forte, et l’extraordinaire intuition d’Euler lui-même ne l’empêche pas de tomber parfois dans l’absurde, lorsqu’il écrit par exemple 0 = ∑_{n=−∞}^{+∞} x^{n}”. Did Euler really fool himself?

To keep with our habits (after Cauchy!) denote by z a complex variable and try to evaluate the sum of I = ∑_{n=−∞}^{+∞} z^{n}. We break the sum into I_{+} + I_{−} − 1, where

I_{+} = ∑_{n≥0} z^{n}, I_{−} = ∑_{n≤0} z^{n}.

By the geometric series, we get I_{+} = \frac{1}{1−z} and I_{−} = \frac{1}{1−1/z} and simple algebra gives

I_{+} + I_{−} = \frac{1}{1 − z} + \frac{z}{z − 1} = \frac{1 − z}{1 − z} = 1,  (52)

hence I = 0 as claimed. What is paradoxical is that there is no complex number z ≠ 0 for which both series I_{+} and I_{−} converge simultaneously, since ∑_{n≥0} z^{n} converges for |z| < 1 and ∑_{n≤0} z^{n} converges for |z| > 1. We really need analytical continuation: I_{+} as a function of z extends from the convergence domain |z| < 1 to C − \{1\} as the rational function \frac{1}{1−z}, and one goes from I_{+}to I_{−} by inverting z (into 1/z). If both I_{+} and I_{−} are extended in this way to C − {1}, the calculation (52) is perfectly valid, hence I = 0 in this sense.

Another method to prove I = 0 is to observe that multiplication of I by z shifts z^{n} to z^{n+1}, hence rearranges the series, hence Iz = I, hence I(z − 1) = 0, and by dividing by z − 1 we get I = 0 for z ≠ 1. Nevertheless, there is some trouble. Consider the critical region |z| = 1 where both I_{+} and I_{−} diverge, and use polar coordinates z = e^{2πiu}. Then I is the series

J(u) = ∑_{n=−∞}^{+∞} e^{2πinu}.  (53)

Playing with Fourier series, introduce a test function f(u) supposed to be smooth (i.e., infinitely differentiable) and periodic, f(u + 1) = f(u). We expand it as a Fourier series

f(u) = ∑_{n=−∞}^{+∞} c_{n}e^{2πinu}  (54)

with

c_{n} = ∫_{0}^{1} e^{−2πinu}f(u)du.  (55)

From (54) we get, by putting u = 0,

f(0) = ∑_{n=−∞}^{+∞} c_{n}  (56)

hence by (55) and (53)

f(0) = ∑_{n=−∞}^{+∞} ∫_{0}^{1} e^{−2πinu}f(u)du = ∫_{0}^{1} { ∑_{n=−∞}^{+∞} e^{−2πinu}} · f(u)du,

and finally

f(0) = ∫_{0}^{1} J(u)f(u)du.  (57)

Remove now the assumption f(u + 1) = f(u) by introducing a smooth function ϕ(u) vanishing off some finite interval and by defining

f(u) = ∑_{m=−∞}^{+∞} ϕ(u + m)  (58)

(an absolutely convergent series). By an easy manipulation, one derives from (57)

∑_{m=−∞}^{+∞} ϕ(m) = ∫_{−∞}^{+∞} J(u)ϕ(u)du.  (59)

Using the standard Dirac’s function δ(u), we get by definition

ϕ(m) = ∫_{−∞}^{+∞} ϕ(u)δ(u − m)du,

hence

0 = ∫_{−∞}^{+∞} {J(u) − ∑_{m=−∞}^{+∞} δ(u − m)} ϕ(u)du.  (60)

Since the test function ϕ is arbitrary, we can omit it from (60), hence the conclusion

J(u) = ∑_{m=−∞}^{+∞} δ(u − m).  (61)

That is, by substituting e^{2πiu} to z, the series I = ∑_{n=−∞}^{+∞} z^{n} is not 0 but ∑_{m=−∞}^{+∞} δ(u−m) . So Euler was wrong, but not too much, since δ(u−m) = 0 for u ≠ m, hence ∑_{n=−∞}^{+∞} z^{n} is 0 for z ≠ 1.

Recall the other proof, using

I(z − 1) = 0;  (62)

division by z − 1 gives I = 0, provided z ≠ 1, corresponding to u ∉ Z for z = e^{2πiu}. Formula (62) is equivalent to

J(u)(e^{2πiu} − 1) = 0,  (63)

and this suggests a new proof of (61). Indeed, if f(u) is a smooth function with isolated simple zeroes u_{m}, then J(u)f(u) = 0 implies that J(u) is a linear combination of terms c_{m}δ(u − u_{m}). Here f(u) = e^{2πiu} − 1, hence u_{m} = m for m in Z, that is m = 0, ±1, ±2, … hence J(u) = ∑_{m=−∞}^{+∞} c_{m}δ(u − m) for suitable coefficients c_{m}. But J(u + 1) = J(u), hence all coefficients c_{m} are equal to some constant c and J(u) = c ∑_{m=−∞}^{+∞} δ(u − m). It remains to calculate the normalization constant c. That kind of argument could be understood by Euler, but it acquires now a rigorous meaning due to Laurent Schwartz’s theory of distributions (200 years after Euler!)^{10}.

Another version of our proof is by using a contour integral (see Fig. 5). Consider a function Φ(z) holomorphic in a domain containing the annulus r ≤ |z| ≤ R bounded by C_{+} and C_{−} (beware the orientations). The rational function R(z) = \frac{1}{z−1} is given by a convergent series ∑_{n=−∞}^{−1} z^{n} for |z| > 1, hence for z in C_{+}. It follows

∫_{C_{+}} R(z)Φ(z)dz = ∑_{n=−∞}^{−1} ∫_{C_{+}} z^{n}Φ(z)dz.  (64)

Similarly

∫_{C_{−}} R(z)Φ(z)dz = ∑_{n=0}^{∞} ∫_{C_{−}} z^{n}Φ(z)dz,  (65)

and using the residue formula,

∫_{C_{+}} ∑_{n=−∞}^{−1} z^{n} · Φ(z)dz + ∫_{C_{−}} ∑_{n=0}^{∞} z^{n} · Φ(z)dz = 2πiΦ(1).

10 Our final result can be expressed as ∑_{n=−∞}^{+∞} e^{2πinu} = ∑_{m=−∞}^{+∞} δ(u − m). It is equivalent to Poisson’s summation formula.

Fig. 5. Path for the contour integral

A shorthand would be

∑_{n=−∞}^{+∞} z^{n} = 2πiδ(z − 1)  (66)

using δ-functions in the complex domain^{11}.

Let us go back to sums of powers and Bernoulli numbers and polynomials. A classical formula reads as follows:

B_{k}(u) = −k! ∑_{n≠0} \frac{e^{2πinu}}{(2πin)^{k}} .  (67)

11 A classical formula is

δ(f(u)) = ∑_{m} \frac{1}{|f′(u_{m})|} δ(u − u_{m}),

the summation being extended over the solutions of the equation f(u_{m}) = 0 (provided f′(u_{m}) ≠ 0). In this formula f(u) is a real-valued function of a real variable u. Assuming that it remains valid for f(u) = e^{2πiu} − 1 (with complex values), we derive

2πiδ(z − 1) = ∑_{m=−∞}^{+∞} δ(u − m)

for z = e^{2πiu}. This brings together our two methods.

A complex version is the following

∑_{n≠0} \frac{z^{n}}{n^{k}} = −\frac{(2πi)^{k}}{k!} B_{k} (\frac{log z}{2πi} ) (68)_{k}

(by the change of variables z = e^{2πiu}) for k = 0, 1, 2, . . .. For k = 0, this reads as Euler’s “absurd formula”

∑_{n≠0} z^{n} = −1, (68)_{0}

since B_{0}(x) = 1. The case k = 1 is

∑_{n≠0} \frac{z^{n}}{n} = πi−log z; (68)_{1}

of course

∑_{n=1}^{∞} \frac{z^{n}}{n} = log \frac{1}{1 − z}

∑_{n=−∞}^{−1} \frac{z^{n}}{n} = − ∑_{m=1}^{∞} \frac{(z^{−1})^{m}}{m} = − log \frac{1}{1 − 1/z}

and (68)_{1} amounts to

log \frac{1}{1 − z} − log \frac{1}{1 − 1/z} = πi − log z.  (69)

Since log(−1) = −πi and \frac{1}{1−1/z} = \frac{z}{z−1} = \frac{z·(−1)}{1−z}, this relation follows from log uv = log u+log v, but some care has to be exercized with the multivalued complex logarithm (notice the ambiguity log(−1) = ±πi for instance).

For the general case, notice the following. Define

L_{k}(z) = ∑_{n≠0} \frac{z^{n}}{n^{k}} , R_{k}(z) = −\frac{(2πi)^{k}}{k!} B_{k} (\frac{log z}{2πi} ) .  (70)

We know already that L_{0}(z) = R_{0}(z) and L_{1}(z) = R_{1}(z). Furthermore, it is obvious that

z \frac{d}{dz} L_{k}(z) = L_{k−1}(z),  (71)

and from the fact that the derivative of B_{k}(x) is kB_{k−1}(x), one gets

z \frac{d}{dz} R_{k}(z) = R_{k−1}(z).  (72)

So we can easily conclude that L_{2}(z) − R_{2}(z) is a constant, which has to be shown to be 0 to prove L_{2}(z) = R_{2}(z). Then L_{3}(z) − R_{3}(z) is a constant, etc…

This line of argument can be made rigorous. Introduce the polylogarithmic functions^{12}

Li_{k}(z) = ∑_{n=1}^{∞} \frac{z^{n}}{n^{k}} .  (73)

Our formula reads now as follows:

Li_{k}(z)+(−1)^{k}Li_{k} (\frac{1}{z}) = −\frac{(2πi)^{k}}{k!} B_{k} (\frac{log z}{2πi} ) . (74)_{k}

To make sense of it, we proceed as follows:

a) We cut the complex plane along the real interval [0, +∞[, to get Ω_{0} = C − [0, +∞[.

b) In the cut plane, we choose the somewhat unusual branch of the logarithm log(re^{iθ}) = log r + iθ for 0 < θ < 2π.

c) We define the function Li_{k}(z) by the convergent series (73) for |z| < 1, and verify that

z \frac{d}{dz} Li_{k}(z) = Li_{k−1}(z)  (75)

for k = 1, 2, … and Li_{0}(z) = \frac{z}{1−z}. Since the cut plane Ω_{1} = C − [1, ∞[ is simply connected, any holomorphic function in Ω_{1} has a primitive, hence by (75), each Li_{k}(z) extends analytically to Ω_{1}.

d) For z in Ω_{0}, both z and \frac{1}{z} are in Ω_{1}, hence both Li_{k}(z) and Li_{k}(\frac{1}{z}) are defined for z in Ω_{0}, and formula (74)_{k} is asserted for z in Ω_{0}.

e) The cases k = 0 and k = 1 are settled as before.

f) From (75) and the rule for the derivative of B_{k}(x), we get that the validity of (74)_{k} for the index k implies that of (74)_{k+1} for the index k + 1

12 The dilogarithm Li_{2}(z) was known to Euler, and further developed in the 19^{th} century in connection with Lobatchevski geometry. Fifteen years ago, the subject was almost forgotten, to be resurrected by geometers and mathematical physicists alike. It is now a hot subject of research.

up to the addition of a constant. To show that it is 0 use the fact that for k ≥ 2, the series ∑_{n=1}^{∞} \frac{z^{n}}{n^{k}} converges also for |z| = 1, and study the limiting value for z → 1, using B_{k}(0) = B_{k}(1).

So after all, Euler was right!

Putting z = 1 in (74)_{k} we obtain the value of ζ(k) + (−1)^{k}ζ(k). For k odd, we get 0 = 0, but for k even, we recover the value of ζ(k) given by (24).

4.6 Variation II: Infinite products

Suppose we want to calculate ∞! = 1 · 2 · 3 · · ·. Going to logarithms we define

∞! = exp (∑_{n=1}^{∞} log n) .  (76)

Suppose we have a series ∑_{n≥1} a_{n} with na_{n} bounded. What is known as the ζ-summation procedure fits our general framework at the beginning of Section 4.1: consider the convergent series ∑_{n≥1} a_{n}n^{−ϵ} for ϵ > 0, and let ϵ tend to 0. To sum the series ∑_{n≥1} log n, we should consider ∑_{n≥1} n^{−ϵ} · log n but this converges for ϵ > 1 only and we cannot go directly to the limit ϵ → 0. What we have to do is to consider the series ∑_{n≥1} n^{−s} · log n for ℜ(s) > 1; this is obviously the derivative −ζ′(s) of the Riemann zeta function, hence it can be analytically continued to the neighborhood of 0. The regularized sum of ∑_{n≥1} log n is then −ζ′(0) and finally

∞! = e^{−ζ′(0)}.  (77)

From the formulas (37) and (41), one derives without much ado ζ′(0) = −\frac{1}{2} log 2π (a formula more or less equivalent to Stirling’s formula). Conclusion:

∞! = \sqrt{2π}.  (78)

General rule: to normalize a divergent product Π_{n≥1}a_{n}, introduce the series ∑_{n≥1} a_{n}^{−s} = Z(s), make an analytic continuation from the convergence domain ℜ(s) > σ_{0} to s = 0 and define

∏_{n≥1}^{reg} a_{n} := e^{−Z′(0)}.  (79)

Generalizing slightly (78), we can use this method to prove the identity^{13}

∏_{n≥0}^{reg}(n + v) = \frac{\sqrt{2π}}{Γ (v)}.  (80)

We can compare this to the Weierstraß product expansion for the gamma function

\frac{1}{Γ (v)} = ve^{γv} ∏_{n≥1} (1 + \frac{v}{n}) e^{−v/n}.  (81)

A careless, but nevertheless instructive, comparison of (80) and (81) is as follows:

∏_{n≥1} (1 + \frac{v}{n}) e^{−v/n} = ∏_{n≥1} \frac{n + v}{n} e^{−v/n}

= ∏_{n≥1}(n + v)( ∏_{n≥1} n)^{−1} exp ( − v ∑_{n≥1} 1/n)

and regularizing the divergent series ∑_{n≥1} \frac{1}{n} by the Euler constant γ, we are through! Notice that the two most important properties of Γ, namely

1) the functional equation Γ (v + 1) = vΓ (v);

2) the function \frac{1}{Γ (v)} of a complex variable v is entire with zeroes at 0, −1, −2, … can be read off immediately from (80).

According to our general method, the proof of (80) requires to study the function

ζ(s, v) = ∑_{n≥0}(n + v)^{−s}  (82)

known as Hurwitz zeta function (see [9] for more details). We list a few properties:

a) a particular case ζ(s, 1) = ζ(s);

b) functional equations:

ζ(s, v + 1) = ζ(s, v) − v^{−s}  (83)
∂_{v}ζ(s, v) = −sζ(s + 1, v);  (84)

13 Formally: \frac{∞!}{(∞+v)!} = Γ (v + 1).

c) analytic continuation: for fixed v, ζ(s, v) can be analytically continued to the complex plane with one singularity at s = 1, with singular part \frac{1}{s−1}; hence ζ(s, v) − ζ(s) is an entire function;

d) special values:

ζ(−k, v) = −\frac{B_{k+1}(v)}{k + 1}  (85)

for k = 0, 1, 2, … The last relation can be written, in the spirit of Euler, as

v^{k} + (v + 1)^{k} + (v + 2)^{k} + … = −\frac{B_{k+1}(v)}{k + 1}.  (86)

As a particular case we get the surprising identity

v^{0} + (v + 1)^{0} + (v + 2)^{0} + … = \frac{1}{2} − v.  (87)

5 Conclusion: From Euler to Feynman

Feynman is the modern heir to Euler. Among his many contributions to theoretical physics, the most famous one is his use of diagrams to encode in a very compact way complicated integrals with significance in experiments in high energy physics. His method of diagrams has been generalized by various authors (Cvitanovic, Penrose,…) to provide a very flexible tool for computations in tensor analysis.

His really bold discovery is the use of integrals in function spaces (see for instance [5]), the so-called Feynman path integrals. These (so far) ill-defined integrals are powerful tools to evaluate infinite series and infinite products. We give just one example. Consider the Hilbert space L^{2}(0, 2π) of functions f(x) with 0 < x < 2π and ∫_{0}^{2π} |f(x)|^{2}dx finite. The unbounded operator Δ = −d^{2}/dx^{2} can be diagonalized with eigenfunctions e_{n}(x) = e^{inx} (for n = 0, ±1, ±2, …) corresponding to the eigenvalue n^{2}. Hence the characteristic determinant det(v − Δ) is an entire function with the eigenvalues as zeroes. Using our regularized products, one now defines the regularized determinant as

det^{reg}(v − Δ) = v ( ∏_{n≥1}^{reg}(v − n^{2})^{2})  (1)

(0 is a simple eigenvalue, and 1^{2}, 2^{2}, … are eigenvalues of multiplicity 2). This can be evaluated by a formula due to Euler

sin v = v ∏_{n≥1} (1 − \frac{v^{2}}{n^{2}π^{2}} ) ,  (2)

equivalent (via logarithmic derivatives) to the formula

cot v = \frac{1}{v} + ∑_{n≥1} \frac{2v}{v^{2} − n^{2}π^{2}}  (3)

also due to Euler, and considered above.

Feynman’s bold step is as follows. From matrix calculus, we learn the following integral formula for a characteristic determinant

det(v − A) = [∫_{R^{n}} d^{n}x exp −π (v ∑_{i=1}^{n} x_{i}^{2} − ∑_{i,j} a_{i,j}x_{i}x_{j})]^{−2}  (4)

where d^{n}x is the volume element dx_{1} ⋯ dx_{n} in the Euclidean space R^{n}, and A = (a_{i,j}) is a real symmetric, positive definite, matrix of size n × n. By analogy, Feynman writes det(v − Δ) as the square of the inverse of

∫_{L^{2}(0,2π)} Dx · exp(−πS(x)),  (5)

where the so-called action S(x) is defined by

S(x) = v ∫_{0}^{2π} x(t)^{2}dt + ∫_{0}^{2π} x′(t)^{2}dt  (6)

(the variable in [0, 2π] is denoted by t, the function in L^{2}(0, 2π) by x(t), and its derivative by x′(t)). The symbol Dx is formally a volume element in the Hilbert space L^{2}(0, 2π) (infinite-dimensional generalization of the Euclidean space R^{n}), sometimes written as C ∏_{t} dx(t). Its rigorous definition is the main problem [5].

Part of these calculations have been put into a rigorous framework, but not all of them.

After all, Feynman shall be right!

Acknowledgements

Ils vont à Michel Planat et Michel Waldschmidt, pour leur impeccable organisation du colloque de Chapelle des Bois, et pour leur conviction que ce texte verrait le jour. Leurs amicales pressions et leur aide matérielle ont grandement contribué à l’heureux dénouement. Merci aussi à Christian Krattenthaler pour une relecture extrêmement soignée de ce texte.

References

1. Borel E. (1902) Leçons sur les séries à termes positifs. Gauthier-Villars.

2. Bourbaki N. (1976) Fonctions d’une variable réelle (chapitres 1 à 7). Hermann, edited by Masson since 1980.

3. Bourbaki N. (1984) Eléments d’histoire des mathématiques. Masson.

4. Cartier P. (1966) Quantum mechanical relations and theta functions, in Algebraic groups and discontinuous subgroups. Proc. Symp. Pure Maths, vol. IX, American Mathematical Society, Providence, 361-383.

5. DeWitt-Morette C., Cartier P. and Folacci P. (editors) (1997) Functional Integration, Basics and Applications, Nato Asi, Series B, vol 361. Plenum Press.

6. Euler L. (1755) Institutiones calculi differentialis cum ejus usu in analysi infinitorum ac doctrina serierum. Berlin, 1910.

7. Hardy G. H. (1924) Orders of infinity. Cambridge Univ. Press (reprinted by Hafner, New York, 1971).

8. Knopp K. (1964) Theorie und Anwendung der unendlichen Reihen. Springer (5th edition).

9. Waldschmidt M. et al (1995) From Number Theory to Physics. Springer (see Chapter 1 by P. Cartier: An introduction to zeta functions).

10. Weil A. (1984) Number Theory, An Approach Through History. Birkhäuser.