Formula for number of $(3,k)$ magic squares

Question

Let $n,k\in\mathbb{N}.$ By a $(n,k)$ magic square, we mean a $n×n$ matrix containing non-negative integer entries such that the sum of entries of any given row or column is $k.$ Note that we don't need the diagonals to add up to $k.$ Prove that the number of $(3,k)$ magic squares is ${{k+4}\choose{4}}+{{k+3}\choose{4}}+{{k+2}\choose{4}}.$

I tried induction on $k,$ but that doesn't seem to work. In an attempt to work backwards, I realised that if we fill the $4$ top left entries of the matrix, then the remaining elements are uniquely determined. Maybe this is where the $4$ in the ${{k+4}\choose{4}}+{{k+3}\choose{4}}+{{k+2}\choose{4}}$ comes from. I didn't get far with this either.

Other than this, I've noticed that the number of $(2,k)$ magic squares is $k+1.$ Also, if a $(3,k)$ magic square contains $k$ somewhere, that row and column must be filled with zeroes. Then, the remaining $4$ elements of the matrix form a $(2,k)$ magic square.

The formula makes it seems like we may use stars and bars somewhere. But it is not obvious to me.

I haven't gotten much else. How to solve this?

Answer 1

Let me begin by invoking some heavy machinery. All this theory can be intimidating, but it will pay off in the end. I will put the part it's okay to skim inside a quote (though I'm not quoting anyone in particular).

If we take an $(3,k)$ magic square, and divide all its entries by $k$, we get a doubly stochastic matrix: a nonnegative matrix whose rows and columns all add up to $1$.

Geometrically, we can think of the doubly stochastic matrices as points in $\mathbb R^9$ (where the $9$ entries of the matrix are coordinates), but the row and column sum conditions define an affine subspace of $\mathbb R^9$ which has dimension $4$. (This is essentially the observation already stated in the question: the $4$ top left entries of a magic square determine the other $5$.) Additionally, we don't get the whole affine subspace: we carve off a bounded region of the subspace by the $9$ nonnegativity constraints.

This $4$-dimensional polytope in $\mathbb R^9$ with is called the Birkhoff polytope $B_3$. By a case of the Birkhoff–von Neumann theorem, it is the convex hull of its $6$ corner points, which are the $3\times 3$ permutation matrices.

In terms of $B_3$, we can describe the $(3,k)$ magic squares as follows: they are the integer points of $k B_3$, the polytope we get by scaling every point in $B_3$ by a factor of $k$. The problem of counting integer points in a scaled-up version of a polytope is a generally studied one. In particular, it is known that in the case of an integral polytope (a polytope that, like $B_3$, has corner points with integer coordinates), the number of integer points is given by a polynomial function. That polynomial function is called the Ehrhart polynomial of the polytope.

Here's the reason for bringing in all this high-powered discrete geometry. One of the facts that we know about the Ehrhart polynomial is that it has degree $d$, where $d$ is the dimension of the polytope. In other words: we know that the number of $(3,k)$ magic squares is some degree-$4$ polynomial in $k$.

To prove that it's actually $\binom{k+4}{4} + \binom{k+3}{4} + \binom{k+2}{4}$, it's enough to check that this formula gives the correct number for $5$ different values of $k$: one more than the degree. Two degree-$4$ polynomials that agree on $5$ different values must be equal, or else their difference would be a polynomial of degree at most $4$ with $5$ roots, which violates the fundamental theorem of algebra.

Three values of $k$ can be checked without too much suffering:

• When $k=0$, the formula gives us $\binom 44 + \binom 34 + \binom 24 = 1$, and in fact there is only one $(3,0)$ magic square: the zero matrix.
• When $k=1$, the formula gives us $\binom 54 + \binom 44 + \binom 34 = 6$, and in fact there are six $(3,1)$ magic squares: the six permutation matrices.
• When $k=2$, the formula gives us $\binom 64 + \binom 54 + \binom 44 = 21$. The $(3,2)$ magic squares come in two types: $6$ of them are obtained by doubling a permutation matrix, and $\binom 62$ of them are obtained by adding two permutation matrices together. (It can be checked that these sums are all distinct.) $6+\binom62 = 21$, so the formula still works.

I really don't want to check $k=3$ and $k=4$, where the formula gives $55$ and $120$, though that would be sufficient to solve the problem. So I will invoke some more heavy machinery: Ehrhart–Macdonald reciprocity.

This is the statement that plugging in negative values into the Ehrhart polynomial is also meaningful! Up to a factor of $(-1)^d$ (which in our case is $(-1)^4 =1$), evaluating the polynomial at $-k$ gives the number of interior points of the polytope scaled up by $k$. A point is an interior point if it does not lie on any boundary of the polytope. In our case, the boundaries are the nonnegativity constraints, so a $(3,k)$ magic square is an interior point if all its entries are positive.

For $k=1$ and $k=2$, there are no such magic squares, so the Ehrhart polynomial should be $0$ at $k=-1$ and $k=-2$. In fact, that's also true of our conjectured formula: $$\binom34 + \binom 24 + \binom 14 = \binom 24 + \binom14 + \binom04 = 0.$$ So our conjectured formula agrees with the Ehrhart polynomial at five values: $k=-2,-1,0,1,2$. Since both of them are degree-$4$ polynomials, they actually agree at all values of $k$, proving that the formula is always correct!

As a bonus, we can consider the formula at $k=-3$, where $\binom{k+4}{4} + \binom{k+3}{4} + \binom{k+2}{4}$ simplifies to $0 + 0 + \binom{-1}{4} = 1$. This agrees with the fact that there's only one $(3,3)$ magic square whose entries are all positive: the magic square in which all entries are $1$.