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}$?
 A: 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.
A: 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?
A: (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.
A: 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).
A: 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
A: 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.  
A: 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.
A: 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$.
A: 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...
A: 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[1] \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.
A: 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}$$
