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 A008292 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.)