How do you add up the first $n$ numbers? The first $n$ squares? The first $n$ cubes?

These sums of powers — $S_k(n) = 1^k + 2^k + 3^k + \cdots + n^k$ — have been studied for over 2,500 years. Each one has a closed-form formula, and the pattern that connects them leads directly to some of the deepest objects in number theory.

The classical formulas

$$S_0(n) = \underbrace{1 + 1 + \cdots + 1}_{n} = n$$

$$S_1(n) = 1 + 2 + 3 + \cdots + n = \frac{n(n+1)}{2}$$

The story goes that Gauss discovered this as a schoolboy, but the result was known to the Pythagoreans around 550 BC. The trick: pair $1$ with $n$, $2$ with $n-1$, etc. — each pair sums to $n+1$, and there are $n/2$ pairs.

$$S_2(n) = 1^2 + 2^2 + \cdots + n^2 = \frac{n(n+1)(2n+1)}{6}$$

Known to Archimedes (c. 250 BC), who needed it for computing the area under a parabola.

$$S_3(n) = 1^3 + 2^3 + \cdots + n^3 = \frac{n^2(n+1)^2}{4} = \left(\frac{n(n+1)}{2}\right)^2$$

Discovered by Aryabhata (476 AD). Notice: the sum of cubes equals the square of the sum of the first $n$ numbers. This beautiful identity has no analogue for higher powers.

The general formula and Bernoulli numbers

Is there a pattern? Fermat, Pascal, and Jakob Bernoulli all worked on the general case. Bernoulli found it:

$$S_k(n) = \frac{1}{k+1} \sum_{j=0}^{k} \binom{k+1}{j} B_j \, n^{k+1-j}$$

where $B_j$ are the Bernoulli numbers — a sequence of rational numbers that appears throughout number theory:

$$B_0 = 1, \quad B_1 = -\frac{1}{2}, \quad B_2 = \frac{1}{6}, \quad B_3 = 0, \quad B_4 = -\frac{1}{30}, \quad B_5 = 0, \quad B_6 = \frac{1}{42}, \ldots$$

The odd-indexed Bernoulli numbers (except $B_1$) are all zero. The even-indexed ones grow rapidly and alternate in sign.

From finite to infinite

What happens when $n \to \infty$? For $k \geq 1$, $S_k(n)$ grows without bound — you’re adding bigger and bigger terms. But for negative powers, the terms shrink:

$$S_{-2}(n) = \frac{1}{1^2} + \frac{1}{2^2} + \cdots + \frac{1}{n^2}$$

As $n \to \infty$, this converges to $\pi^2/6$ — the Basel problem.

More generally, $S_{-k}(\infty) = \zeta(k)$ — the Riemann zeta function. And Euler showed:

$$\zeta(2m) = \frac{|B_{2m}|}{2 \cdot (2m)!} (2\pi)^{2m}$$

The Bernoulli numbers connect finite sums of powers to infinite sums of inverse powers. The same sequence that gives you a formula for $1^k + 2^k + \cdots + n^k$ also tells you the exact value of $1/1^{2m} + 1/2^{2m} + 1/3^{2m} + \cdots$.

The $S_{-1}$ gap

For $k = -1$, we get the harmonic series:

$$S_{-1}(n) = 1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n} \approx \ln n + \gamma$$

This diverges — it’s the singular case, the pole of the zeta function. Every other negative power converges. The harmonic series sits exactly at the boundary.