Let $M(n,k)$ be the set of $n\times n$ matrices of nonnegative integers such that every row and every column sums to $k$. Let $P(n,k)$ be the fraction of such matrices which have no zero entries, equivalently the probability that a random matrix from the uniform distribution on $M(n,k)$ has no zero entries.
One thing to note is that $$P(n,k)=\frac{|M(n,k-n)|}{|M(n,k)|}$$
(think about subtracting 1 from every entry). Also note that $P(n,k)$ is the fraction of integer points in the $k$-dilated Birkhoff polytope that lie in the interior.
It seems "obvious" that $P(n,k)$ is a non-decreasing function of $k$. For large enough $k$ it is strictly increasing by Ehrhart theory, but I'd like to see a proof for all $k$. So the problem is:
Prove that $P(n,k)$ is a non-decreasing function of $k$ for fixed $n$.
COMMENT: By a result of Stanley (see David Speyer's answer for refs) there are non-negative integers $\{h_i\}$ such that
$$|M(n,k)| = \sum_{i=0}^d h_i \binom{k+i}{d}$$
where $d=(n-1)^2$. I'm wondering if this is enough. Does every polynomial of that form have the desired properties? [Gjergji showed not.]
COMMENT2: The reciprocity theorem for Ehrhart series provides a formula for the number of points in the interior in terms of the number in the whole (closed) polytope. Making use of the above, we find that if the $H_n(k)$ is the polynomial equal to $|M(n,k)|$ for positive integers $k$, then the number of interior points (already identified as $H_n(k-n)$) equals $(-1)^{n+1}H_n(-k)$. So what we have to prove is that
$$(-1)^{n+1}\frac{H_n(-k)}{H_n(k)}$$
is non-decreasing for integer $k\ge n$. Experimentally, it is not increasing for real $k$ until $k$ is larger.
Best Answer
This is frustrating because there is a lot of study of this sort of question in the combinatorial commutative algebra literature, and the exact example you discuss is a favorite example in this field, but I can't find an answer to your question. So here are some pointers to the relevant background.
Fix $n$. Let $R$ be the semigroup ring corresponding to the semigroup of $n \times n$ nonnegative integer matrices all of whose row and column sums are equal. So the Hilbert series of $R$ is $$h(x) := \sum_k |M(n,k)| x^k$$ By the Ehrhart theorem, $|M(n,k)|$ is a polynomial in $k$, so we can write $$h(x) = \frac{\delta(x)}{(1-x)^{(n-1)^2+1}}$$ for some polynomial $\delta(x)$ of degree $\leq (n-1)^2$. Richard Stanley pioneered the study of the relation between commutative algebra properties of $R$ and combinatorial properties of $\delta$. In particular, the ring $R$ is Gorenstein (i.e. the canonical module is generated by the all ones matrix) and this implies that $\delta$ is palindromic with positive coefficients. The standard reference for this material is Stanley's book "Combinatorics and Commutative Algebra".
It might be worth pausing for an example: According to this paper, $$|M(3,k)| = 1 + \frac{9 k}{4} + \frac{15 k^2}{8} + \frac{3 k^3}{4} + \frac{k^4}{8}=\frac{(k + 1) (k + 2) (k^2 + 3 k + 4)}{8}$$ and we can compute $$h(x) = \frac{1+x+x^2}{(1-x)^5}.$$ The polynomial $\delta$ is $1+x+x^2$.
Now, we have the following implications:
(1) $\delta(x)$ has all real roots $\implies$
(2) Writing $\delta(x) = \sum \delta_k x^k$, we have $\delta_k^2 \geq \delta_{k-1} \delta_{k+1}$ $\implies$
(3) We have $M(n,k)^2 \geq M(n,k-1) M(n,k+1)$ $\implies$
(4) $M(n,k) M(n,k+n-1) \geq M(n,k-1) M(n,k+n)$, which is the relation you want.
The implication $(3) \implies (4)$ is elementary; the others are discussed in Stanley's superb survey "Log-concave and Unimodal sequences in Algebra, Combinatorics, and Geometry".
In that survey, Stanley made Conjecture 4, claiming that the Hilbert series of any Cohen-Macaulay domain should obey (2). $R$ is a Cohen-Macualay domain (any normal semi-group ring is, by a result of Hochster), but Conjecture 4 turned out to be false: according to "Log-concave and Unimodal sequences in Algebra, Combinatorics, and Geometry: an update" by Brenti (scroll down, about the tenth reference from the bottom of the page), a counterexample can be found in Niesi and Robbiano; I haven't checked this reference. Stanley and Brenti both suggest that the conjecture is more plausible for Gorenstein rings, and a quick skim through the Mathscinet papers which cite them suggest that no Gorenstein counterexample is known.
I went over to Dennis Pixton's article which tabulates the Ehrhart polynomials for this exact problem. The largest value he lists is $n=9$. I computed the corresponding $\delta$ (coefficients available on request) and found that it violated (1) but obeyed (2).
Specifically, for $n=9$, the polynomial $\delta$ has degree $56$. Of its roots, $52$ are real and clearly isolated. The four remaining roots are $(-170629.9 \pm 70111.4 i)^{\pm 1}$. (This all assumes you trust Mathematica's numerical algorithms.)
Finally, you should probably know a few keywords: The polytope of $n \times n$ matrices with nonnegative entries and row and column sums equal to $1$ is the "Birkhoff polytope". The more general case where you just fix $(a_1, a_2, \ldots, a_n)$ and $(b_1, b_2, \ldots, b_n)$ with $\sum a_i = \sum b_i$ and ask for row sums $a_i$ and column sums $b_i$ is the "transportation polytope".