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$.
