# 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 $\frac{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 $\frac{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. Apr 11, 2011 at 2:56
• 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/… Apr 11, 2011 at 4:42
• @Alex, I take it that by "simply connected" you mean "connected in a simple way". Apr 11, 2011 at 8:32
• I suppose I don't agree with the premise. Bernoulli numbers, it seems clear to me, do not appear everywhere. They're relatively rare objects that come up in fairly particular circumstances. If they appeared everywhere, I imagine I would be seeing them more often. Dec 10, 2014 at 4:37
• Basically, where ever you see the derivative you could replace it with the Bernoulli polynomials and make some sense of it, but you are right, Ryan. I've never seen even one under my bed. Never seen an exponential there either. What was the OP thinking of?! Dec 10, 2014 at 20:07

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. Apr 11, 2011 at 3:42
• You can do this in higher dimensions too, integrating over polytopes. Apr 11, 2011 at 12:44
– aorq
Apr 12, 2011 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. Sep 18, 2012 at 2:56
• From the previous comment, the Bernoulli polynomials and numbers as defined by the operator relation $$B_n(x+1)-B_n(x)=(B.(x+1))^n-(B.(x))^n = dx^n/dx = n x^{n-1}$$ are directly related to differentiation of analytic functions through $$f(B.(x+1))-f(B.(x))= f'(x)$$ and to expressing the power sums that crop up in the Riemann sum approximations (e.g., as sum of the areas of n rectangles with bases of equal length) to integrals of analytic functions expressible as power series through $$\sum_{k=1}^n k^p = [B_{p+1}(n+1)-B_{p+1}(1)]/(p+1).$$ Dec 26, 2018 at 19:27

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 ... Apr 11, 2011 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?? :) Apr 11, 2011 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 Apr 12, 2011 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)$. Sep 18, 2012 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/… . Nov 25, 2014 at 18:42

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?" Apr 11, 2011 at 14:44
• See also Vector Bundles and K-Theory by Allen Hatcher. Dec 11, 2014 at 0:11
• @Neil Strickland, do you have some particularly good Web-accessible refs sketching these ideas? Jan 8, 2017 at 17:10

I have been told (by J.H. Conway) that Bernoulli numbers were first discovered by Faulhauber. See the 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/… Nov 25, 2014 at 19:10

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. Nov 20, 2014 at 9:36
• How about for penance saying two Hail Marys and finding two new articles? My answer is way too long to add more, and I'd be interested in some more perspectives. (I'll delete my comments.) Nov 20, 2014 at 11:13
• OK. Saying two Hail Marys produced the following two new (I hope) articles: maths.ed.ac.uk/~aar/papers/hirzrem.pdf (which explains why Bernoulli numbers do appear in topology) and sciencedirect.com/science/article/pii/S0022247X07008694 (which you probably already know). And a special prayer for the Dodgers resulted in emis.de/journals/INTEGERS/papers/o17/o17.pdf (wich probably will give some more perspectives). Nov 21, 2014 at 10:58
• Yep, thanks. I was aware of those, but through the sciencedirect suggested references I spotted a very important connection to solitons (and the Ricatti equation) I hadn't noted before. Thanks very much for leading me near that trail. Nov 22, 2014 at 1:10

(Revamped. Majority of entry at my mini-arxiv now.)

The natural, unique extension of the Bernoulli numbers is the Bernoulli Appell polynomial sequence, which can be operationally defined by action on polynomials and analytic functions, such as the exponential and logarithm, when convergent, or order by order for a formal power series, through the umbral relation $f(B.(x+1))-f(B.(x))={f}'(x)=D_x \; f(x)$, where D is the derivative. Action on $exp(xt)$ gives the e.g.f., of the Bernoulli polynomials without any resort to numerical values of the Bernoulli numbers. Action on the Mellin transform of $exp(-xt)$ defines the Bernoulli numbers in terms of the Riemann zeta values. The e.g.f. of their umbral compositional inverses is the reciprocal of the one for the Bernoullis, which gives the "reciprocal polynomials", based on reciprocal natural numbers, very naturally associated with not only the exponential divided by its argument, but also the logarithm.

From grafting the Bernoulli and reciprocal polynomials together stem a pair of Lie operator derivatives for powers of the state number, or Euler op, and associated normal ordered ops. The Lie ops are related to the compositional inverse pair of functions that are the log and exp functions for the multiplicative formal group law associated to the Todd class. The matrix reps are conjugates of the infinitesimal generator of the Pascal triangular matrix by the mutually orthogonal Stirling number matrices and encode the combinatorics of simplices. The multiplicative, compositional, and umbral compositional inversions are inextricably bound together and reveal myriad associations to combinatorics, Lie theory, and topology.

The interplay of the Bernoulli and reciprocal polynomials reveal this. It also provides easy proofs of many, if not most, identities involving the Bernoullis and a way of looking at the Riemann zeta function that can not be readily achieved from the perspective of the e.g.f. operators of the Euler-Maclaurin expansion. For example, regarding the Mellin transform as a means of interpolation, the natural extension of the Bernoulli polynomials is the Hurwtz 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)$.

For relations to Hirzebruch genera / Todd class, see this MOQ.

• I don't see how any of this explains why Bernoulli numbers show up in topology. Nov 18, 2014 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. Nov 19, 2014 at 2:34
• The series formula for zeta above involves the base number sequence for the generalized Bernoulli polynomials of the second order, whose asymptotics show the the series above must be an asymptotic expansion. ----- The MOQ link above for the Hirzebruch genera no longer contains my answers. I transferred them to my website, where I can refine and expand on them. Jan 15, 2015 at 7:33

Occurence of these numbers in the formula for the Todd class is related with Campbell-Hausdorff, look for instance http://arxiv.org/abs/math/0610553

I also always wondered about the appearance of the quite arbitrary looking function $$\dfrac{x}{\exp(x)-1}$$. Here is one way of not letting it come out of the blue but getting it by a line of argument. This is an answer to the p.s. of the questioner, so only a partial answer to the initial question, since it will not contribute to the question why the Bernoulli numbers appear everywhere (which, as far as I can see, has not been satisfactorily answered in this thread, anyway), but will refer to the aspect that they are natural things to consider. As, by the way, are the Bernoulli polynomials, which on our way of developing things also make a natural appearance.

The point of departure is the determination of the power sum polynomial $$S_p(x)$$, $$p = 0,1,2,\dots$$. This is the unique polynomial with $$\forall{k \in \mathbb{N}:\quad S_p(k) = 0^p+1^p+2^p+\cdots+(k-1)^p},$$ where $$0^p := [p=0]$$. It cann alternatively be characterized as the unique polynomial $$S_p(x)$$ satisfying $$\tag{1} \Delta S_p(x) = x^p \quad,\quad S_p(0) = 0.$$ where $$\Delta$$ is the forward difference operator acting on polynomials via $$\Delta f(x) := f(x+1) - f(x).$$ Differentiating (1) gives, since $$\Delta$$ commutes with differentiation, $$\Delta S'_p(x) = px^{p-1} = p\Delta S_{p-1}(x)$$ or $$\Delta(S'_p(x) - pS_{p-1}(x)) = 0$$ or $$\tag{1a} S'_p(x) - pS_{p-1}(x) = B_p$$ for some constant $$B_p$$, (dependent on $$p$$). To steer this ship into familiar waters, define polynomials $$\tag{2} B_p(x) := S'_p(x) \quad,\quad p=0,1,2,\dots$$ of degree $$p$$ as to have, on differentiating (2): $$\tag{3} B'_p(x) - pB_{p-1}(x) = 0,$$ or $$\tag{4} \fbox{B'_p(x) = pB_{p-1}(x)\quad,\quad p=0,1,2,\dots}$$ so that we have done away succesfully with the unwanted (for the time being) constants $$B_p$$. The condition $$S_p(0)=0$$ can be transformed, since $$S_p(1) = [p=0]$$ by definition of $$S_p(x)$$, into the equivalent condition $$\tag{5} S_p(1) - S_p(0) = \int_0^1 B_p(t) dt = [p=0].$$ Thus we have reformulated the problem of finding $$S_p(x)$$ of degree $$p+1$$ as finding $$B_p(x)$$ of degree $$p$$ such that (4) and (5) hold; then we recover $$S_p(x)$$ as $$\tag{5a} S_p(x) = \int_0^x B_p(t) dt \overset{(3)}{=} \dfrac{1}{p+1} \int_0^x B'_{p+1}(t) dt = \dfrac{1}{p+1}(B_{p+1}(x)-B_{p+1}(0)),$$ and so we have to take care of the $$B_p(x)$$, in particular to show that they are well-defined by the conditions (4) and (5).

As said before, the intention behind the introduction of the $$B_p(x)$$ was to reach familiar waters, and, indeed, we have done so, because (4) shows that the sequence of the $$B_p(x)$$ is an Appell sequence, so named after Paul Émile Appell (1855-1930):

Definition A sequence of polynomial $$A_p(x), p=0,1,2,\dots$$ is called an Appell sequence if

(i) $$A_0(x)$$ is a nonzero constant;

(ii) $$\forall{p=1,2,\dots}\quad\quad A'_p(x) = pA_{p-1}(x)$$.

Since in an Appell sequence $$A_p(x)$$ is, by integration, determined by $$A_{p-1}(x)$$ up to a constant, $$a_p$$, say, via $$\tag{6} A_p(x) = p\int_0^x A_{p-1}(t) dt + a_p,$$ where $$a_0 \not= 0$$ because of (i), the Appell sequences over $$\mathbb{R}$$ -- or any field of characteristic 0 -- correspond in a one-to-one fashion to invertible power series $$a(x) = \sum_{p=0}^{\infty} a_p\dfrac{t^p}{p!} \in \mathbb{R}[[t]].$$ If the Appell sequence $$(A_p(x))_{p \in \mathbb{N}}$$ corresponds to $$a(x)$$ via (6), $$a(x)$$ is called its characteristic series.

To identify the Appell sequence $$(B_p(x))_{p \in \mathbb{N}}$$, our task is to identify its characteristic series $$b(x)$$, and for this we need more theory. So consider, for a given Appell sequence $$(A_p(x))_{p \in \mathbb{N}}$$ its Maclaurin series $$A_p(x) = \sum_{k=0}^p A^{(k)}_p(0) \dfrac{x^k}{k!}.$$ Because of (ii), one has $$A^{(k)}_p(0) = p(p-1)\cdots(p-k+1)A_{p-k}(0) = p(p-1)\cdots(p-k+1)a_{p-k},$$ so that $$\tag{7} A_p(x) = \sum_{k=0}^p \binom{p}{k} a_{p-k}x^k.$$ This is a binomial convolution, so it suggests to introduce the exponential generating function (EGF) of any sequence $$(f_p(x))_{p \in \mathbb{N}}$$ of polynomials as $$f(x,t) := \sum_{p=0}^{\infty} f_p(x) \dfrac{t^p}{p!} \in \mathcal{P}[[t]]$$ with $$\mathcal{P} := \mathbb{R}[t]$$ the polynomial ring. For example, if $$f(x) = x^p$$, then $$f(x,t) = \sum_{p=0}^{\infty} x^p \dfrac{t^p}{p!} = \exp(xt).$$ Then, for an Appell sequence $$(A_p(x))_{p \in \mathbb{N}}$$, with $$\tag{8} A(x,t) := \sum_{p=0}^{\infty} A_p(x) \dfrac{t^p}{p!},$$ one has $$\tag{9} A(0,t) = a(t) \quad,\quad A(x,t) = A(0,t)\exp(xt) = a(t)\exp(xt).$$ Namely, if I put $$x:=0$$ in (8): $$A(0,t) = \sum_{p=0}^{\infty} A_p(0) \dfrac{t^p}{p!} =:a(t),$$ which is the first claim of (9), while (7) implies the second claim of (9). These relations show explicitely how the Appell sequence $$(A_p(x))_{p \in \mathbb{N}}$$ and its characteristic series $$a(t)$$ determine each other: given $$A(x,t)$$, one has $$a(t) = A(0,x)$$, and given $$a(x)$$, one has $$A(x,t) = a(t)\exp(xt)$$, so that the $$A_p(x)$$ are just given by the binomial convolution of the sequences $$(a_0,a_1,a_2,\dots)$$ and $$(1,x,x^2,\dots)$$.

By (6), we see that at any stage the polynomial $$A_p(x)$$ of an Appell sequence is determined by its predecessor up to a constant which is nonzero for $$p=0$$, so this triggers the idea of fixing an Appell sequence uniquely by putting a nondegenerate linear constraint on each member; in fact, this is exactly what happens in (5), since $$f(x) \mapsto \int_0^1 f(t) dt$$ is a linear form on $$\mathcal{P}$$. To explicate this, let $$L:\mathcal{P} \rightarrow \mathbb{R}$$ a linear form on $$\mathcal{P}$$. We denote the value of $$L$$ on the polynomial $$f(x) \in \mathcal{P}$$ by $$L[f(x)]$$. Let $$F:=(f_p(x))_{p \in \mathbb{N}}$$ be a basic sequence, i.e. having $$\deg f_p(x) = p$$ for all $$p$$. Since a linear form is uniquely defined by the sequence of values it takes on a basic sequence, the linear form $$L$$ is uniquely defined by, and uniquely defines, the power series $$c_{L,F}(t) := \sum_{p=0}^{\infty} L[f_p(x)] \dfrac{t^p}{p!}$$ In particular, if $$M:=(x^p)_{p \in \mathbb{N}}$$ is the standard monomial basis, this is the case for $$c_L(t) := c_{L,M}(t) = \sum_{p=0}^{\infty} L[x^p] \dfrac{t^p}{p!}$$ which we call the canonical series of $$L$$, and we call $$L$$ nondegenerate if $$L \not= 0$$ (this terminology is not standard). Then

Theorem Let $$L \in \mathcal{P}^*$$ be a nondegenerate linear form on $$\mathcal{P}$$. Then there is a unique Appell sequence $$A:=(A_p(x))_{p \in \mathbb{N}}$$ with $$\tag{10} \forall p \in \mathbb{N}:\quad L[A_p(x)] = [p=0],$$ namely the unique Appell sequence with characteristic series $$\tag{11} a(t) = (c_{L}(t))^{-1},$$ i.e. with EGF $$\tag{12} A(x,t) = (c_{L}(t))^{-1} \exp(xt).$$

Proof $$\,\,$$ By (7) we have, for any Appell sequence $$A:=(A_p(x))_{p \in \mathbb{N}}$$: $$\ L[A_p(x)] = \sum_{k=0}^p \binom{p}{k} a_{p-k}L[x^k]$$ which just says $$\ c_{L,A} = a(t) c_L(t).$$ But the conditions (10) are just equivalent to $$c_{L,A} = 1$$, and so we are done. QED

So, to determine the Appell sequence $$B = (B_p(x))_{p \in \mathbb{N}}$$, all we have to do is to determine its characteristic series $$b(t)$$ and so, by (11) and (12), the canonical series $$c_L(t)$$ for $$L$$ given by (5): $$c_L(t) = \sum_{p=0}^{\infty} \left(\int_0^1 x^p dx \right)\dfrac{t^p}{p!} = \sum_{p=0}^{\infty} \dfrac{1}{p+1} \dfrac{t^p}{p!} = \sum_{p=0}^{\infty} \dfrac{t^p}{(p+1)!} = \dfrac{\exp(t)-1}{t},$$ which gives us the EGF of its constant coefficients $$B_p:=B_p(0)$$ as $$\tag{13} b(t) = \left(\sum_{p=0}^{\infty} \dfrac{t^p}{(p+1)!}\right)^{-1} = \dfrac{t}{\exp(t)-1},$$ and the EFG of the $$B_p(x)$$ as $$\tag{14} B(x,t) = b(t)\exp(xt) = \dfrac{t\exp(xt)}{\exp(x)-1}.$$ hence so we have recovered the formerly unwanted constants $$B_p$$ of (1a) and call them Bernoulli numbers; they are given, but not defined, by their EGF(13). And so have we the polynomials $$B_p(x)$$ and call them Bernoulli polynomials; they are given, but not defined, by their EGF(14). And, in passing, we have proved, via (5a), Jakob Bernoulli's celebrated formula $$\ S_p(x) = \dfrac{1}{p+1} \sum_{k=1}^{p+1} \binom{p+1}{k}B_{n-k+1} x^k,$$ so to speak the cradle of the Bernoulli numbers.

Final Remarks It is an annoying fact that it has become common practice to define the Bernoulli numbers by their EGF (13) (this EGF goes back to Euler, by the way), since this ad hoc definition comes out of the blue and has motivating value approaching zero exponentially. The same holds for the common definition of the Bernoulli polynomials via (14), which comes staggering along with the same motivational value. The computation above makes them appear, at least, as a result of a comprehensible line of argument. Unfortunately, however, this line of argument appears to add nothing to the clearance of the question why the Bernoulli numbers are so ubiquitous.

• Thank you, my mistake. Aug 23, 2020 at 7:35

You can look for the insight here.

In short, it we introduce a non-Archimedean numerical system that interprets generalized summations of infinite series (in Ramanujan's sense or Zeta function regularization) as "standard part" of the sum of the series, then in this system the quantity of all natural numbers $\omega_-$ has the standard part $-1/2$ and the quantity of non-negative integers $\omega_+$ is greater by 1 (for zero), so has the standard part $1/2$.

Consequently, from Faulhaber's formula for Ramanujan's summation,

$$\operatorname{st}\omega_-^n=B_n$$

$$\operatorname{st}\omega_+^n=B^*_n$$

where Where $B_n$ are the first Bernoulli numbers and $B^∗_n$ are the second Bernoulli numbers.

So the Bernoulli numbers are the standard part of the powers of the quantity of naturals.

Similarly,

$$\operatorname{st} e^{z\omega_-}=\frac{z}{e^{z}-1}$$