# Why do Bernoulli numbers arise everywhere?

I have seen Bernoulli numbers many times, and sometimes very surprisingly. They appear in my textbook on complex analysis, in algebraic topology, and of course, number theory. Things like the criteria for regular primes, or their appearance in the Todd class, zeta value at even numbers looks really mysterious for me. (I remember in Milnor's notes about characteristic class there is something on homotopy group that has to do with Bernoulli numbers, too, but I don't recall precisely what that is. I think they also arise in higher K-theory.)

The list can go on forever. And the wikipedia page of Bernoulli number is already quite long.

My question is, why do they arise everywhere? Are they a natural thing to consider?

==========================================

p.s.----(maybe this should be asked in a separate question)

Also, I've been wondering why it is defined as the taylor coefficient of the particular function x/(e^x-1), was this function important? e.g. I could have taken the coefficient of the series that defines the L-genus, namely $\dfrac{\sqrt{z}}{\text{tanh}\sqrt{z}}$, which only amounts to change the Bernoulli numbers by some powers of 2 and some factorial. I guess many similar functions will give you the Bernoulli numbers up to some factor. Why it happen to be the function x/(e^x-1)?

-
According to history of mathematics, at the beginning of Bernoulli's original idea, he did not begin defined as the Taylor coefficients of the function x/(e^x-1) at 0. –  yaoxiao Apr 11 '11 at 2:56
That's a good point. So do you know how did he define it and what's his motivation? –  36min Apr 11 '11 at 4:19
Just a minor remark: the space of the areas you have mentioned might be simply connected. For example, the criteria for regular primes, the zeta values at even integers, and higher K-theory - these are all closely related areas. If you are willing to take on board everything we know and everything we conjecture to be true, then it follows that Bernoulli numbers appear in one of these if and only if they appear in all of them. Somewhat related is my answer here: mathoverflow.net/questions/45376/… –  Alex B. Apr 11 '11 at 4:42
@Alex, I take it that by "simply connected" you mean "connected in a simple way". –  gowers Apr 11 '11 at 8:32
This has become a list of reasons without one being obviously definitive, so I am making this Community Wiki. –  Todd Trimble yesterday

I don't know of a universal theory of all places where Bernoulli numbers arise, but Euler-Maclaurin summation explains many of their more down-to-earth occurrences.

The heuristic explanation (due to Lagrange) is as follows. The first difference operator defined by $\Delta f(n) = f(n+1)-f(n)$ and summation are inverses, in the same sense in which differentiation and integration are inverses. This just amounts to a telescoping series: $\sum_{a \le i < b} \Delta f(i) = f(b) - f(a)$.

Now by Taylor's theorem, $f(n+1) = \sum_{k \ge 0} f^{(k)}(n)/k!$ (under suitable hypotheses, of course). If we let $D$ denote the differentation operator defined by $Df = f'$, and $S$ denote the shift operator defined by $Sf(n) = f(n+1)$, then Taylor's theorem tells us that $S = e^D$. Thus, because $\Delta = S-1$, we have $\Delta = e^D - 1$.

Now summing amounts to inverting $\Delta$, or equivalently applying $(e^D-1)^{-1}$. If we expand this in terms of powers of $D$, the coefficients are Bernoulli numbers (divided by factorials). Because of the singularity at "$D=0$", the initial term involves antidifferentiation $D^{-1}$, i.e., integration. Thus, we have expanded a sum as an integral plus correction terms involving higher derivatives, with Bernoulli number coefficients.

Specifically, $$\sum_{a \le i < b} f(i) = \int_a^b f(x) \, dx + \sum_{k \ge 1} \frac{B_k}{k!} (f^{(k-1)}(b) - f^{(k-1)}(a)).$$ (Subtracting the values at $b$ and $a$ just amounts to the analogue of turning an indefinite integral into a definite integral.)

This equation isn't literally true in general: the infinite sum usually won't converge and there's a missing error term. However, it is true when $f$ is a polynomial, and one can bootstrap from this case to the general one using the Peano kernel trick.

So from this perspective, the reason why $t/(e^t-1)$ is a natural generating function to consider is that we sometimes want to invert $e^t-1$ (the factor of $t$ is just to make it holomorphic), and the most important reason I know of to invert it is that we want to invert $\Delta = e^D-1$.

-
If I recall correctly, the infinite sum converges to the correct answer also for polynomials times exponentials. –  Allen Knutson Apr 11 '11 at 3:42
You can do this in higher dimensions too, integrating over polytopes. –  Steve Huntsman Apr 11 '11 at 12:44
See also this question: mathoverflow.net/questions/10667 –  aorq Apr 12 '11 at 0:34
Also umbrally $p(B(x)+1)-p(B(x))={p}'(x)$ where $B(x)$ are the Bernoulli polynomials and $p(x)$ any polynomial, so you might expect them to pop up in approximation and linearization problems. –  Tom Copeland Sep 18 '12 at 2:56
Define the Bernoulli polynomials as the Appell sequence, i.e., $(B.(0)+x)^n=B_n(x)$, such that $f(B.(x+1))-f(B.(x))={f}'(x)$ when convergent. Then $e^{B.(x+1)t}-e^{B(x)t}=te^{xt}$ implies $e^{B.(x)t}(e^t-1)=te^{xt}$ and the e.g.f. and makes the appearance of $t$ natural. Relates the Bernoulli polynoms to the tangent space. –  Tom Copeland 2 days ago

The main reason I know for the appearance of Bernoulli numbers is the one Henry Cohn already explained: we'd like to invert the difference operator $e^D - 1$, so we'd like to expand $1/(e^D - 1)$ as a Taylor series. But $1/(e^x - 1)$ doesn't have a Taylor series, because it has a pole at the origin. It has a perfectly nice Laurent series, but just to make things more obscure people prefer to discuss the Taylor series of $x/(e^x - 1)$. And the coefficients of this are called Bernoulli numbers.

I understand how Bernoulli numbers are used to compute $\sum_{i=1}^n i^k$ and how they show up in formulas for the Riemann zeta function.

However, Alain Connes loses me here:

• Alain Connes, Andre Lichnerowicz and Marcel Paul Schutzenberger, A Triangle of Thoughts, AMS, Providence, 2000.

He points out that if $H$ is the Hamiltonian for some sort of particle in a box and $\beta$ is the inverse temperature,

$$1/(1 - e^{-\beta H}) = 1 + e^{-\beta H} + e^{-2 \beta H} + \cdots$$

is the operator you take the trace of to get the partition function for a collection of an arbitrary number of particles of this sort. And he claims that pondering this explains all the appearances of $x/(1 - e^x)$ and the Bernoulli numbers in topology!

Does anyone understand that? I imagine he's hinting at some relation between characteristic classes, the heat equation, the Laplacian on differential forms, and things like that. But I've never understood how these pieces are supposed to fit together.

And here's something that remains more mysterious to me. The paper by Kervaire and Milnor has a cool formula for the order of the group of smooth structures on the $(4n-1)$-sphere for $n > 1$. It's:

$$2^{2n-4} (2^{2n-1} - 1) P(4n-1) B(n) a(n) / n$$

where:

$P(k)$ is the order of the $k$th stable homotopy group of spheres

$B(k)$ is the $k$th Bernoulli number, in the sequence 1/6, 1/30, 1/42, 1/30, 5/66, 691/2730, 7/6, ...

$a(k)$ is 1 or 2 according to whether k is even or odd

How do the Bernoulli numbers weasel their way into this game?

-
I think I have a rough idea of how Bernoulli numbers get into exotic spheres. There's this "J-homomoprhism" Z x BO --> Pic(S), which sends a (stable) real vector bundle to the (stable) spherical fibration given by one-point compactifying the fibers. The Spivak normal fibration over S^n gives a canonical map S^n --> Pic(S), and finding manifold structures on S^n is tied to lifting this map along J (so saying that the normal fibration actually came from a normal vector bundle). This lifting problem comes down to understanding what J does on homotopy groups. This in turn can be understood ... –  Dustin Clausen Apr 11 '11 at 16:19
in terms of Pontryagn-Thom theory: a class in the kernel of J on pi_n is represented by an "exotic" null-bordism of S^n. These correspond to almost parallelizable manifolds (manifolds parallelizable away from a fixed point), and once you've got one of those most of your characteristic classes vanish, and (for instance) the integrality of the A-hat genus will tell you something about just one particular Bernoulli number. Well OK, but what does this have to do with Euler-Maclaurin summation?? :) –  Dustin Clausen Apr 11 '11 at 16:23
John discussed the relevant passage of 'A Triangle of Thoughts' at greater length here golem.ph.utexas.edu/category/2008/02/metric_spaces.html#c014936 –  David Corfield Apr 12 '11 at 8:27
Also Harer and Zagier have shown that the orbifold Euler characteristic of the moduli space $\mathcal M_{g,n}$ is $(-1)^{n-1}\frac{(2g+n-3)!}{(2g-2)!} \zeta(1-2g)$. –  Tom Copeland Sep 18 '12 at 2:46
For the Harer-Zagier formula the Bernoulli numbers weasel their way in through an asymptotic expansion of the digamma fct. See last page of ocw.mit.edu/courses/mathematics/… . –  Tom Copeland 2 days ago

In algebraic topology one key point is as follows. The complex $K$-theory spectrum has homotopy groups $KU_{\ast}=\mathbb{Z}[u,u^{-1}]$, with $u$ in degree two. This ring maps in an obvious way to $\pi_{\ast}(H\wedge KU)$, and it is not hard to calculate that the resulting map $\mathbb{Z}[u,u^{-1}]\to\pi_{\ast}(H\wedge KU)$ induces an isomorphism $\mathbb{Q}[u,u^{-1}]\to\pi_{\ast}(H\wedge KU)$. The ring $(H\wedge KU)^0(\mathbb{C}P^\infty)$ can be described as $\mathbb{Q}[[ux]]$ or as $\mathbb{Q}[[y]]$, where $x$ comes from $H^2(\mathbb{C}P^\infty)$ and $y$ comes from $KU^0(\mathbb{C}P^\infty)$. Specifically, if we let $L$ denote the tautological line bundle, then $y$ can be taken to be the $K$-theory class of the virtual bundle $L-1$. It then works out that $y=e^x-1$, so $x/y$ is the Bernoulli series. The Bernoulli numbers occur as coefficients of $x^k/k!$ rather than $x^k$ itself, which suggests that one should work with $\Omega S^3$ rather then $\mathbb{C}P^\infty$: there is a canonical map $\Omega S^3\to\mathbb{C}P^\infty$ using which we can identify $H^*(\Omega S^3)$ with the ring of all series of the form $\sum_ka_kx^k/k!$ with $a_k\in\mathbb{Z}$. All this is of course linked with Adams's treatment of the $J$-homomorphism and the $e$-invariant. However, I think that much of this is still mysterious, at least to me. I think there are some missing ingredients involving the relationship between $R$ and $gl_1(R)$ for various $E_\infty$ ring spectra $R$, particularly those related to surgery and the $J$-homomorphism. There is a lot of literature about this kind of thing from the 1970s but I have not managed to extract the answers that I wanted.

-
I remember this question coming up one tea time when I was sitting next to Frank Adams, and he said the reasons for the occurrence of Bernoulli numbers in topology were "trivial". Makes me wish I'd asked "compared to what?" –  Charles Matthews Apr 11 '11 at 14:44

I have been told (by J.H. Conway) that Bernoulli numbers were first discovered by Faulhauber. See the http://en.wikipedia.org/wiki/Faulhaber's_formula wikipedia article for details. This reference hardly answers your question, but one possible characterization is that they are useful in summing nth powers. That fact alone indicates that their ubiquity is quite natural to expect.

-

Another way where they show up is in Lie theory. If you want to compute the derivative of the exponential map (of a Lie group) you encounter the function $x/(e^x - 1)$ quite inevitably. The already posted questions can partly be viewed as incarnations of this. THis results also in the appearence of the Bernoulli numbers in the BCH series, which is of course of fundamental importance far beyond the usage in Lie algebra theory...

-
Hausdorff (1906) via Iserles, "Expansions that grow on trees": Matrix diff. eqn. ${Y}'= A(t)\; Y$ with change of variables $Y(t)=e^{\Omega(t)}$ becomes ${\Omega}'=exp[B.\;ad_{\Omega}]A=\frac{ad_{\Omega}}{e^{ad_{\Omega}-1}}A\;.$ citeseerx.ist.psu.edu/viewdoc/… –  Tom Copeland 2 days ago

The Bernoulli numbers naturally enter so many domains because their escorts are the simple, but elegant, reciprocal integers.

The formalism that links the two and reveals the associations is that of the grand mathmages Blissard and Appell. First we make the simple complicated and then the complicated becomes simple. Appell polynomial sequences are an extension of the Kronecker delta base sequence $\delta.=\delta.(0)=(1,0,0,0, \cdots.)$, i.e., $(\delta.)^n=\delta_n=\delta_n(0)$. Translating the sequence to a higher plane through the enchanted binomials, we encounter the power basis

$$(\delta.(0)+x)^n= \sum_{k=0}^{n} \binom{n}{k}\delta_j(0)x^{n-j}=x^n=\delta_n(x)=(\delta_.(x))^n$$

with lowering (creation) operator $L=D=\frac{d}{d(\delta.(x))}=\frac{d}{dx}$ and raising (destruction / annihilation) operator $R=\delta_.(x) =x$, i.e.,

$$D\; \delta_n(x)=n \cdot \delta_{n-1}(x), \;\;\; R \; \delta_n(x) = \delta_{n+1}(x)\;,$$

with commutator $[L,\;R] = [D,x] = 1$ (Graves, 1857, contemporary of Blissard and Sylvester, the compulsive neologian who coined the phrase umbral calculus).

(So already we see shadows of the Heisenberg-Weyl Lie algebra in the iconic Appell sequence, and it's no surprise that the probabilist's Hermite (Appell) polynomials appear for the harmonic oscillator in quantum mechanics and the Bernoulli's for other Q.M. domains, and the BCH theorem.)

The exponential generating fct. for the Appell sequence is then $$E_{\delta}(x,t) = e^{\delta.(x)t}=e^{(\delta.(0)+x)t}=e^{\delta.(0)t}e^{xt}=e^{xt}$$ and an ordinary generating fct. $O_{\delta}(x,t)$is the formal Borel-Laplace transform of the e.g.f.

$$\int_{0}^{\infty } e^{\delta.(x)u}e^{-\frac{u}{t}}du=\frac{t}{1- \delta.(x)t}=\sum_{k \ge 0}\delta_{n}(x)t^{n+1}=t [1-(\delta.(0)+x)t]^{-1}.$$

We rise to an even higher plane with the greatest of all the intuitive mathmages Ramanujan and use his master formula (the Mellin transform) to make the indices continuous

$$\int_{0}^{\infty } e^{-\delta.(x)u}\frac{u^{s-1}}{(s-1)!}du=\frac{1}{ (\delta.(x))^{s}}= \delta_{-s}(x)$$ $$=[1-[1-(\delta.(0) + x)]]^{-s}$$

so we have the Mellin transform or a Newton-Gauss interpolator for extending (and analytically continuing to the complex domain) the base sequence. For the Kronecker base sequence this is $x^{-s}$.

Now simply apply the formalism to the Bernoulli numbers $B.(0)$ and out pops the Bernoulli polynomials and the Hurwitz zeta function, which specializes to the Riemann zeta for $x = 1$, for which $B.(x=1)=-B.(x=0)$.

But looking at the umbral compositional inverse of an Appell sequence really allows us to understand how the reciprocal integers and the Bernoullis are mated to give their offspring in zeta functions, Lie algebras / matrix operator reps, convolution algebras, Euler-Maclaurin operator expressions, Faulhaber and other identities, K-theory, and interesting combinatorics, such as that of permutahedra and associahedra through reciprocation and functional compositional inversion.

Define the umbral compositional inverse $\bar{\delta}.(x)$ by

$$\delta_n(\bar{\delta}.(x))= x^n = \bar{\delta}_n(\delta.(x)).$$

Then use the translation property twice to give

$$\delta_n(\bar{\delta}.(x))= (\delta.(0)+\bar{\delta}.(0)+x)^n =x^n,$$ and setting $x=0$ defines the base sequence of the umbral inverse as

$$(\delta.(0)+\bar{\delta}.(0))^n =\delta_n.$$

Exponentiating helps us to readily interpret this as

$$e^{(\delta_.(0) + \bar{\delta}.(0))t} = e^{\delta.(0)}e^{\bar{\delta}.(0)t}=e^{\delta.t}=1.$$

The e.g.f.s of the base sequences are reciprocals of each other. This means the base sequences (and these could be almost any abelian numbers, operators, matrices, etc.) are connected by the combinatorics of surjections and permutahdera A133314 (A049019), among other important far-reaching implications. The algebra can be mapped to finite and infinite matrix algebras (Lie op algebras) with infinitesimal generators (nilpotent for finite rank). But that's for another night.

Back to the Bernoullis extended to polynomials defined by

$$e^{B.(x)t}=e^{(B.(0)+x)t}=e^{B.(0)}e^{xt}=\frac{t}{e^t-1}e^{xt}$$

with the umbral inverse polynomials, their escorts, the elegant reciprocal integers,

$$e^{\bar{B}.(x)t}=e^{(\bar{B}.(0)+x)t}=e^{\bar{B.}(0)}e^{xt}=\frac{e^t-1}{t}e^{xt}$$

with $$\bar{B}_n=\frac{(x+1)^{n+1}-x^{n+1}}{n+1}$$ and

$$\bar{B}_n(0)=\frac{1}{n+1}.$$

(The Pascal matrix nudges its way in here with all the combinatorial import (see OEIS-A074909, A135278), introducing the simplices.)

The e.g.f.s morphed into operators give you the Euler-Maclaurin expansion (and more since the two e.g.f.s for the base sequence are inverse by construction, independent of their interpretation as shift operators):

For an analytic function $f$,

$$\frac{D_y}{e^{D_y}-1}e^{xD_y}f(y)=e^{B.(x)D_y}f(y)= f(B.(x+y)),$$

and

$$\frac{e^{D_y}-1}{D_y}e^{xD_y}f(y)=e^{\bar{B}.(x)D_y}f(y)=f(\bar{B}.(x)+y)= f(\bar{B}.(x+y))$$

$$= D_y^{-1} [f(x+y+1) - f(x+y)],$$

where $D_y^{-1} y^n/n!= y^{n+1}/(n+1)!$.

The operators are clearly an inverse pair from the umbral inverse properties and commute, so

$$\frac{D_y}{e^{D_y}-1}e^{xD_y}f(\bar{B}.(y))=f(\bar{B}.(B.(x+y)))=f(x+y)$$

$$=\frac{D_y}{e^{D_y}-1}e^{xD_y}D^{-1} [f(y+1) - f(y)]=D^{-1}[f(B.(x+y+1)) - f(B(x+y))],$$

and

$$f(B.(x+y+1)) - f(B(x+y))= f'(x+y).$$

Using these properties and expanding (usually with asymptotic results, see Hardy, Divergent Series),

$$\frac{D_y}{e^{D_y}-1}e^{xD_y}=-\sum_{k\ge 0}e^{(n+x)D_y}D_y = \sum_{k\ge 1}e^{-(n-x)D_y}D_y,$$

the Euler-MacLaurin series can be generated and Faulhaber's formula as well.

Now to the raising operators, for the Bernoullis

$$R_B \;B_n(x) = e^{B.(0)D_x}x\;e^{\bar{B}.(0)D_x}B_n(x)= e^{B.(0)D_x}x\;B_n(\bar{B}.(x))$$ $$= e^{B.(0)D_x}x^{n+1}=(B.(0)+x)^{n+1}=B_{n+1}(x).$$

Likewise for the umbral inverse, $$R_{\bar{B}} = e^{\bar{B}.(0)D_x}x\;e^{B.(0)D_x},$$ and we are conjugating the basic raising op for the Kronecker base sequence. There's more hidden here, and we can reveal it by invoking a commutator and the Pincherle derivative:

$$R_B = x-x + e^{B.(0)D_x}x \; e^{\bar{B}.(0)D_x}= x - e^{B.(0)D_x}[e^{\bar{B}.(0)D_x},x].$$

For a general pair of lowering and raising ops, the Pincherle derivative is

$$[f(L),R]=\frac{d}{dL}f(L)=f'(L)=f(B.(L+1))-f(B.(L)),$$

so we expect the Bernoullis to pop up in all of these algebras one way or another, and we have several further interesting relations (recall the e.g.f.s are reciprocals):

$$R_B = x - e^{B.(0)D}\frac{d}{dD} e^{\bar{B}.(0)D} = x - \frac{d}{dD}ln[e^{\bar{B}.(0)D}] = x + \frac{d}{dD}ln[e^{B.(0)D}]$$ and so, with a simple change of sign,

$$R_{\bar{B}} = x + e^{B.(0)D}\frac{d}{dD} e^{\bar{B}.(0)D} = x + \frac{d}{dD}ln[e^{\bar{B}.(0)D}] = x - \frac{d}{dD}ln[e^{B.(0)D}],$$

which hold for general Appell sequences. For any Sheffer sequence, of which the Appells are a subgroup, the raising op separates into $x+g(D_x)$, so the commutator remains invariant upon substitution of any Sheffer raising op for $x$ in the commutator. The Pincherle derivative is conjugated by the basic underlying e.g.f. of each Appell sequence.

More specifically for the Bernoulli couple, working out the Pincherle derivative and bouncing between e.g.f.s leads to a pairing and back to our familiar Riemann zeta

$$R_B = x + B.(0) e^{-B.(0)\bar{B}.(0)D}=x + \sum_{k \ge 0}(-1)^n\frac{B_{n+1}(0)}{n+1}\frac{D^n}{n!}=x+exp[\zeta(-n)D].$$

($\zeta(-n)$ in the exponential here is meant to be shorthand for $(\zeta(-.))^n=\zeta(-n)$.)

Reprising,

$$R_B = x + exp[\zeta(-n)D], \;\;\;\;\;\;\; R_{\bar{B}} = x - exp[\zeta(-n)D].$$ Or,

$$R_B = x - \frac{1}{2}+ \sum_{n \ge 1} \zeta(1-2n) \frac{D^{2n-1}}{(2n-1)!}= x- \frac{1}{2}+ \sum_{n \ge 1} (-1)^n \frac{2 \zeta(2n)}{(2 \pi)^{2n}} D^{2n-1}$$

So we can see how deeply entwined the reciprocals of the integers, the Riemann zeta, and the Bernoullis are with each other and important families of operator algebras.

Using the Mellin-Riemann-Ramanujan interpolation, the natural extension of the the Bernoullis is the Hurwitz zeta function

$$B_{-s}(x)=s \sum_{n \ge 0}\frac{1}{(n+x)^{s+1}},$$ which with $x=1$ becomes $s\cdot zeta(s+1)$, and for the reciprocal integers,

$$\bar{B}_{-s}(x)=\frac{(x+1)^{1-s}-x^{1-s}}{1-s}.$$ The two are related through umbral composition and inversion, so that the pole singularities are reflected in each other and, in fact, the Gauss-Newton series and umbral composition leads to

$$\zeta(s)=\sum_{n \ge 0}(-1)^{n+1}\;\frac{(-s)!}{n!(2-s-n)!} \frac{2^{2-s}-2^n}{2^n}C_n,$$

for $Re(s)<1$ (but gives good results with just eight terms over the range of reals $-6 \le s \le 2$--it's capturing the dependence of zeta on the singularity, the falling factorials of $s$, and zeta's first three simple zeroes up to that approx.--ten terms captures the dependence on the next zero) where $C.=(1,1,5/6,1/2,1/10,-1/6,-5/42,1/6,...)$ are determined by $C_n=(1-G.)^n$ and $G.=(1,0,-1/6,0,1/10,0,-5/42,0,...)$ come from the umbral composition of the Bernoulli polynomials with the Bernoulli numbers $B_n(1)=(-1)^nB_n(0)$.

I think this makes a good case for the explanatory and constructive power of viewing the Bernoullis' richness as a result of their intimate association with their companion sequence the reciprocal integers through the umbral Appell formalism and the operator algebra it entails. If I had to propose one property of the Bernoulli polynomials that seems to make them unique, it would be their ability to generate the tangent of a function through umbral composition, but they really go hand-in-hand with the reciprocal integers.

As for the o.g.f.s, well, that's a rich story too (take a look at A074909), not to mention cumulants, partition functions, etc., but best saved for another night.

@Qiaochu, see God for why. Google for how and associations, e.g., Coates and Givental, "Quantum cobordisms and formal group laws" (pg. 15), introduce the Bernoullis through Euler-Maclaurin. You already know, I'm sure, how FGLs are related to functional inversion and simple translation, i.e., the derivative operator (note the relation of the Bernoullis to derivatives if you read my answer). In another paper, (Dunne and Schubert, "Bernoulli numbers identities from quantum field theory and topological string theory"), the Bernoullis slip in through differentiation (pg 3) and the digamma fct, which is another raising operator for another Appell sequence (the gamma genus) involving the zeta and cycle index polynomials which have all the dressings of Chern characteristic classes (MOQ-111165) ($\zeta(n+)$, but it seems some MO users don't know there is a reflection formula for the zeta fct, pg. 3 D & S , MOQ-112062). (Follow the OEIS entries also related to o.g.f.s.) Given the relation of the Bernoullis to differentiation and differentiation to the Pincherle commutator, the question to me is really how the Bernoullis could not appear.

Just found this reference rife with the Bernoullis (thanks to MOQ-16169 and Zoran Skoda) "A universal formula for representing Lie algebra generators as formal power series with coefficients in the Weyl algebra" by Durov, Meljanec, Samsarov, and Skoda in which the authors give an alternate briefer proof of their results using a dual approach (coderivations) that runs parallel to the Appell formalism above, no doubt unaware of the connection. Page 43 Eqn.37 is a g.f. for $a.\bar{B}.(T)$, umbrally evaluated, with $a.$ and $T$ in the paper. They also encounter in their original approach the Fibonnaci polynomials A011973 (pg. 16), with no remark on their identity, which are also the Pascal rows read along anti-diagonals (no doubt connected to the duality) and therefore connected to the e.g.f. for $\bar{B}.$--they are also the coefficients for the characteristic polynomials of the Coxeter adjacency matrix for $A_n$, related to the Chebyshev polynomials of the second kind, and to Cartan matrices, and the shifted version, well, that has exciting connections to crossing partitions, positroids, and a general Appell sequence related to compositional inversion and ... . The Appell formalism should bring out sharper connections between all these structures.

And, after following up on David Speyer's answer to the same MO question, I believe the association of the normalized Bernoulli numbers (and Euler and Gennochi numbers) to alternating permutations may give a combinatorial perspective on their presence in the analysis. The normalized versions pop up when considering polylogarithms, trigonometric, and hyberbolic functions that are related to the exponentials. The o.g.f.

$$O_{\bar{B}.}(x,t) = ln\left [ \frac{1-xt}{1-(1+x)t)} \right ]=ln\left [ 1 + L[t,-(x+1)]\right ]$$

where $L(t,x)=t/(1+tx)$ with inverse $L(t,-x)$, a special Mobius transformation (which is the iconic o.g.f.), has the compositional inverse

$$O_{\bar{B}.}^{(-1)}(x,t)=L\left [ e^t-1, (1+x) \right ]=\frac{e^t-1}{1+(1+x)(e^t-1)}.$$ The forms $L(g(t),x)$ are related to colored compositions. You can read how the inverse o.g.f. is related to Eulerians, permutahedra, probablity theory, a Weierstrass elliptic function, and a formal group law, related to a generalized cohomology, through comments and references in A008282 and A074909. It can also be rewritten in terms of the e.g.f.s of the Bernoullis and their umbral inverses. (The combinatorics that underlie reciprocation and compositional inversion are those of the permutahedra, associahedra, crossing partitions, and the myriad combinatoric structures related to them, so no surprise that they make an appearance in all this.)

So, we have this interplay among the Mobius transformation, reciprocation and mutiplicative inversion, umbral and regular composition, and umbral and regular compositional inversion of the logarithm and the exponential that accounts in my mind for the prevalence of the Bernoullis and pairing with the reciprocal integers, and of course the royal binomials (Pascal matrix). The Bernoullis are intimately related to differentiation and, therefore, clearly to the exponential; $e^{B.(x+1)}-e^{B.(x)}=De^x=e^x$ defines them, and the connection of the exponential operator to Lie theory, movement around a manifold, is the fundamental action translation, which is also at the heart of the Appell formalism, in both the indices and the independent variable.

(Edit, Nov 21 2014)

$O_{\bar{B}}^{(-1)}(x,t)$ is an e.g.f. for signed reverse face polynomials of the permutahedra and has the infinitesimal generator $$g(x,u)\frac{d}{du} = [(1-xu)(1-(1+x)u)]\;\frac{d}{du},$$ i.e., $$exp\left [ t\;g(x,u)\;\frac{d}{du} \right ]\;u\; |_{u=0} = O_{\bar{B}}^{(-1)}(x,t).$$ G. Rzadowski in "Bernoulli numbers and solitons revisited" explicitly shows the links between derivatives of $g(x,u)$ to solutions of the Ricatti differential equation, soliton solns. of the KdV equation, and the Eulerian and Bernoulli numbers. In addition, the comp. inversion formula A145271 connects products of derivatives of $g(x,u)$ and the refined Eulerian numbers to $O_{B.}(x,t)$, which gives the normalized, reverse face polynomials of the simplices (A135278, divided by n+1). Or, apply the inversion method of A134264 (intimately related to Appell polynomials in general and associated interpolated families of polynomials spanning the Coxeter group $A_n$) to $$h(x,t) = \frac{t}{O_{\bar{B}.}^{(-1)}(x,t)} = (1+x)t+\frac{t}{e^t-1} = 1 + (1+x)t + 2! B_2 + 3!B_3 + \;...$$ and you get a relation between noncrossing partitions or Dyck lattice paths weighted by the normalized Bernoullis and the face polynomials of the simplices. (Now morph all of this into totally non-negative grassmannians, positroids, binary trees, operads, and computations of characteristic classes of genera (Hirzebruch) and you have can have an exciting math weekend.)

Relation to Hirzebruch genera

R. Lu in his thesis "Regularized equivariant Euler classes and gamma functions" states (pg 43 & 44), "The crucial observation in Hirzebruch's theory [on multiplicative sequences and multiplicative genera] is that every multiplicative sequence gives rise to a multiplicative genus, which is the evaluation of a polynomial of characteritic classes against the fundamental class of the manifold." And, just prior to that, "The key point is that, to every power series with constant term equal to unity, we can associate a multiplicative sequence of polynomials." Lu wasn't aware that these sequences are Appell polynomials in the linear coefficient and that they may be generated by the raising operator of an Appell sequence with substitution of $B_n s_n$ for $B_n$ before any final numerical values are assigned to any of the indeterminates, where $s_n$ is the $n$-th power symmetric polynomial. This is tantamount to substituting $B_{n+1}(0) /(n+1)$ for $\zeta_{n+1}$ in the Libgover gamma multiplicative  polynomials on pg. 51 of Lu, or $s_{n+1}B_{n+1}(0)/(n+1)$ in the polynomials of my MOQ on an Appell sequence based on Riemann zeta values or for the indeterminates of the signed cycle index polynomials for the symmetric group.The cycle index polynomials A036039 themselves are an Appell sequence in the indeterminate $x_1$, and its umbral compositional inverse is itself, mod signs. So now we see that the Appell formalism provides an overarching context for the presence of the Bernoullis in so many apparently disparate domains along with the operational definition $f[B.(x+1)]-f[B.(x)]={f}'(x)$ (and Appell translation property $(B.(0)+x)^n=B(x)$), giving it a special relation to the tangent space. (In addition,the defining relation for the umbral compositional inverse in terms of the Kronecker delta above is identical in form to Grothendieck's axiomatic formula (renormalized by the factorials)  for the Chern class in the Wikipedia article on the topic.)

-
I don't see how any of this explains why Bernoulli numbers show up in topology. –  Qiaochu Yuan Nov 18 at 20:06
@Qiaochu, I really don't understand the point of your comment. Is it a judgement or question? The original question is about everywhere, not just topology. Cohn stated that he doesn't know of a comprehensive answer. This answer presents the Bernoullis in a broader context, which I'm sure many users are not familiar with. The Euler-Maclarin view is very limited, conservative and well-established by the likes of Hardy, so naturally more attractive to the conservative cliques here and at mathstack exchange. I find it lacking. –  Tom Copeland Nov 19 at 2:34
The Bernoullis result from a dance of the reciprocals across the permutahedra, and vice versa. See oeis.org/A049019 and oeis.org/A133314 for relation to surjections, matrix reps for reciprocation, and a weighted graphs interpretation. See Buchstaber and Panov's "Toric Topology" for connections of polytopes to topology. –  Tom Copeland Nov 20 at 20:51
Erratum: For the soliton, from A145271, $dz/dt=df(t)/dt=g(f(t))=g(z)$ here with $f(x,t)= O^{(-1)}_{\bar{B}}(x,t)=z$, so connects the derivatives of $g$ to the f-vectors of the simplicial duals of permutahedra (Eulerians are the h-vectors), not those of the simplices. A134264 gives the critical link between the Bernoullis and normalized, reverse f-vectors of the self-dual simplices $\bar{B}_n$. Also, another link related to Todd class: mathoverflow.net/questions/10630/… . –  Tom Copeland Nov 22 at 20:51
For the Todd-Chern-Hirzebruch connection see mathoverflow.net/questions/60478/… –  Tom Copeland Nov 23 at 16:51

One of the nicest facts about Bernoulli numbers is the formula for zeta evaluated at an even positive integer:
$\zeta(2n) = \dfrac{(-1)^{n+1}(2\pi)^{2n}B_{2n}}{2(2n)!}$
A lot of contexts give rise to this particular kind of series. A simple example, which isn't restricted to even arguments for zeta, is how the set of positive integers that are $k$th-power-free (they are not multiples of $d^{k}$ for any $d > 1$) has density $\frac{1}{\zeta(k)}$ for any $k > 1$ (actually true for $k=1$ too, in a stupid way).
Another example is where Eisenstein series' normalized versions come from (the presence of Bernoulli numbers in the defining formula can look quite mysterious), though I do not have the time or the clarity of thought to work this out here.

-

Apart from pure mathematics, the Bernoulli numbers appear quite often in quantum ﬁeld theory computations due to their relation with the Riemann zeta function. The fundamental reason for this is explained in http://arxiv.org/abs/math/0406610 (Bernoulli Number Identities from Quantum Field Theory and Topological String Theory, by Gerald V. Dunne and Christian Schubert -- the reference indicated in the Tom Copeland's answer) as follows:

"This comes about at a very basic level: perturbative loop calculations in quantum ﬁeld theory generally involve traces of inverse powers of derivatives of functions deﬁned on a circle. Since the spectrum of the ordinary derivative operator $\partial_P$ with periodic boundary conditions consists of the integer numbers, one has $$\rm{tr}(\partial_P^{-2n})\sim\sum\limits_{k=1}^n\frac{1}{k^{2n}}=\zeta(2n).$$ But $\zeta(2n)$ is related to the Bernoulli numbers through Euler’s identity".

Anyway, it seems rather misterious that the wisdom of quantum ﬁeld theory can be used (as described in the Dunne and Schubert's paper cited) to simplify the proof of very non-trivial Miki's identity $$\sum\limits_{k=2}^{n-2}\beta_k\beta_{n-k}=\sum\limits_{k=2}^{n-2}\binom{n}{k}\beta_k\beta_{n-k}+2H_n\beta_n,$$ where $n>2$, $\beta_n=(-1)^nB_n/n$ and $H_n=1+\frac{1}{2}+\ldots+\frac{1}{n}$ is nth harmonic number (see http://people.mpim-bonn.mpg.de/zagier/files/doi/10.1007/978-4-431-54919-2/curious-bernoulli.pdf -- Curious and Exotic Identities for Bernoulli Numbers, by Don Zagier).

-
No, I have not use the reference from you answer. I'm sorry that I have not noticed that the reference appeared in your answer. I first found Zagier's paper and then the search of literature for alternative proofs of Miki's identity resulted in the Dunne and Schubert's paper. –  Zurab Silagadze Nov 20 at 9:36