Démonstration. Notons que clairement $B_1=1$. Soit $n\in {\mathbb{N}}^{\ast }$, exprimons $B_{n+1}$ en fonction des termes précédents. Pour tout $k\le n$, on note $E_k$ l’ensemble des partitions $P$ de $[[1,n+1]]$ tel que la partie de $P$ qui contient l’entier $n+1$ est de taille $k+1$. Choisir une partition dans $E_k$, c’est choisir $k$ entiers de $[[1,n]]$ (ceux de l’ensemble qui contient $n+1$ dans la partition), puis construire une partition des $n-k$ éléments restants. Donc $|E_k|=\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_{n-k}$.

Comme $E_0,\dots ,E_n$ forment une partition de l’ensemble des partitions de $[[1,n+1]]$, on obtient : $$B_{n+1}=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_{n-k}=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{n-k}\right)B_k=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k\text{\hspace{1em}(∗)}$$ À toute partition $P$ de $[[1,n]]$, on peut associer une permutation ${\sigma }_P\in S_n$, qui est le produit des cycles de supports les éléments de $P$. On construit ainsi une application $P\mapsto {\sigma }_P$ injective. D’où : $$B_n\le |S_n|=n!$$ On en déduit en particulier que $\frac{B_n}{n!}\le 1$. En vertu du lemme d’Abel, le rayon de convergence $R$ de la série entière $\sum \frac{B_n}{n!}{x}^n$ est supérieur ou égal à $1$. On peut donc définir $$B:\begin{array}{ccc}]-R,R[& \to & \mathbb{R}\\ x& \mapsto & {\sum }_{n=0}^{+\infty }\frac{B_n}{n!}{x}^n\end{array}$$ et en dérivant, $\forall x\in ]-R,R[$ : $$\begin{array}{rl}{B}^{\prime }(x)& =\sum _{n=0}^{+\infty }\frac{B_{n+1}}{n!}{x}^n\\ & \stackrel{(\ast )}{=}\sum _{n=0}^{+\infty }\frac{1}{n!}\left(\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k\right)x^n\\ & =\sum _{n=0}^{+\infty }\left(\sum _{k=0}^n\frac{B_k}{k!}\frac{1}{(n-k)!}\right)x^n\end{array}$$ On reconnaît là le produit de Cauchy suivant : $$B^{\prime }(x)=\sum _{n=0}^{+\infty }\left(\sum _{k=0}^n\frac{B_k}{k!}\frac{1}{(n-k)!}\right)x^n=\left(\sum _{n=0}^{+\infty }\frac{B_n}{n!}{x}^n\right)\left(\sum _{n=0}^{+\infty }\frac{x^n}{n!}\right)=B(x)e^x$$ Reste à résoudre cette équation différentielle linéaire homogène d’ordre $1$ : $$B(x)=\lambda e^{e^x}$$ Or, $B(0)=B_0=1=\lambda e^1$. D’où $B(x)=\frac{1}{e}{e}^{e^x}$.

La série entière définissant l’exponentielle a un rayon de convergence infini. On peut écrire, pour tout $z\in \mathbb{C}$ : $$e^{e^z}=\sum _{n=0}^{+\infty }\frac{e^{nz}}{n!}=\sum _{n=0}^{+\infty }\sum _{k=0}^{+\infty }\underset{u_{n,k}(z)}{\underbrace{\frac{(nz{)}^k}{n!k!}}}$$ On va appliquer le théorème de Fubini-Lebesgue à $u_{n,k}(z)$ (où $z\in \mathbb{C}$ est fixé) : $$\sum _{n=0}^{+\infty }\sum _{k=0}^{+\infty }|u_{n,k}(z)|=\sum _{n=0}^{+\infty }\frac{e^{n|z|}}{n!}=e^{e^{|z|}}<+\infty$$ Donc on peut intervertir les signes de sommations. Pour tout $x\in ]-R,R[$, $$\begin{array}{rl}B(x)& =\frac{1}{e}{e}^{e^x}\\ & =\frac{1}{e}\sum _{n=0}^{+\infty }\sum _{k=0}^{+\infty }{u}_{n,k}(x)\\ & =\frac{1}{e}\sum _{k=0}^{+\infty }\sum _{n=0}^{+\infty }{u}_{n,k}(x)\\ & =\frac{1}{e}\sum _{k=0}^{+\infty }\left(\sum _{n=0}^{+\infty }\frac{n^k}{n!}\right)\frac{x^k}{k!}\end{array}$$ Par unicité du développement en série entière d’une fonction, on en déduit, par identification des coefficients : $$\forall k\in {\mathbb{N}}^{\ast },B_k=\frac{1}{e}\sum _{n=0}^{+\infty }\frac{n^k}{n!}$$ ◻