Proof: Non-negative integer matrix representations (NIM Rep) of a finite group are always permutation matrices. Let $G$ be a finite group. Definition: Non-negative integer matrix representation (NIM rep) of a finite group $G$ is a group homomorphism $\varphi : G \to GL(n,\mathbb{Z}_{\geq 0})$.
Prove that every matrix in the representation i.e. $\varphi(g)$ for every $g \in G$ is always a permutation matrix.
Phrased another way: Let $M \in  GL(n,\mathbb{Z}_{\geq 0})$. If $M^k = I$ for some $k \in \mathbb{Z}_{\geq 0}$, then prove that $M$ is a permutation matrix.
Edit: The answer by Qiaochu Yuan proves that $GL(n,\mathbb{Z}_{\geq 0}) = S_n$ which is even stronger than what is asked to be proven.
 A: The intuition (hand-waving argument) behind this is that if you have a generic invertible matrix $M$, having only non-negative integers as its elements, then multiplying it with itself more and more will make its matrix elements larger and larger. To satisfy the property $M^k=I$, for some $k\in \mathbb{Z}_+$, intuitively
means that the matrix elements can't be integers greater than $1$ otherwise taking repeated powers will make it explode. In fact, as we show below, it will indeed be a permutation matrix.
Let us define a function $f$ which sums over all the elements of a matrix $M$ i.e. $f:GL(n,\mathbb{Z}_{\geq 0})\to \mathbb{Z}_{\geq 0}$ such that $f(M) = \sum_{ij} M_{ij}$. If $M^k = I$ for some $k \in \mathbb{Z}_{\geq 0}$, $\implies det(M)\neq 0$. This means each row and each column of $M$ must have at least one non-zero element, otherwise, $det(M)$ would be $0$.  $$ \implies f(M)\geq n$$
Also, $$f(M^k) = f(I) = n$$
Claim: $f(M^m)$ is a non-decreasing function in $m \in \mathbb{Z}_{+}$.
Proof: (This is one of those proofs which are easier to visualize using matrix products rather than using index notation but which are easier to write in index notation) We will prove a general fact that if we have two invertible matrices $A$ and $B$ such that $A,B\in GL(n,\mathbb{Z}_{\geq 0}) $ then $f(AB) \geq f(A)$. $$f(AB) = \sum_{i,j,l}A_{ij} B_{jl} = \sum_{i,j}\left( A_{ij} \sum_lB_{jl} \right)$$ Since $det(B)\neq 0$ $\exists$ $l' \in \{1,2,..,n\}$ such that $B_{jl'}>0$  i.e. $B_{jl'}\in \mathbb{Z_+}$
$$\implies \sum_l B_{jl} \geq 1$$
$$\implies f(AB) = \sum_{i,j}\left(A_{ij}\sum_l  B_{jl} \right) \geq \sum_{i,j}A_{ij}$$
(Edit: There is a simpler proof at the end which uses this fact)
This means for $m>1$
$$f(M^m) = f(M^{m-1}M)\geq f(M^{m-1})$$ Therefore, $f(M^m)$ is a non-decreasing function in $m$.
But, $f(M)\geq n$ and $f(M^k) = f(I) = n$, therefore $\;f(M^m)$ must be constant in $m$.
$$f(M^m) = n \; \forall \; m \in \mathbb{Z}_+$$
In particular, this means $f(M)=n$ and since, $M \in GL(n,\mathbb{Z}_{\geq 0})$, every row and every column must contain exactly one $1$ and rest all $0$'s. This proves $M$ is a permutation matrix. $\blacksquare$
Edit: A simple proof inspired by the answer by Qiaochu Yuan:
Since we showed in proof of the claim above that for two invertible matrices $A,B \in GL(n,\mathbb{Z}_{\geq 0})$, $f(AB)\geq f(A)$, we can choose $B$ to be the inverse of $A$. This would imply
$$f(AB) = f(I) = n \geq f(A)$$. But since there is at least one positive integer in every row and every column of $A$ $(det(A) \neq 0)$ i.e. $f(A) \geq n$. The two inequalities on $f(A)$ imply $f(A) = n$. The only way to satisfy this is for $A\in GL(n,\mathbb{Z}_{\geq 0})$ to be a permutation matrix. This implies every element of $GL(n,\mathbb{Z}_{\geq 0})$ is a permutation matrix.
$$GL(n,\mathbb{Z}_{\geq 0}) = S_n$$
A: It is actually unnecessary to assume that $G$ is finite. The following simpler result suffices: if a matrix in $M_n(\mathbb{Z}_{\ge 0})$ is invertible (meaning that its inverse is another matrix in $M_n(\mathbb{Z}_{\ge 0})$), then it is a permutation matrix. In other words, $GL_n(\mathbb{Z}_{\ge 0}) = S_n$.
To see this, suppose $X, Y$ are two matrices with non-negative integer coefficients such that $XY = I$. This means that their entries $x_{ij}, y_{ij}$ satisfy
$$\sum_j x_{ij} y_{jk} = \delta_{ik}.$$
If $i \neq k$ then the RHS is $0$, and since the LHS is a sum of non-negative integers all of them must be equal to $0$. This gives that if $i \neq k$ then either $x_{ij} = 0$ or $y_{jk} = 0$.
If $i = k$ then we get $\sum_j x_{ij} y_{ji} = 1$, which gives that for fixed $i$ there is a unique $j = \sigma(i)$ such that $x_{ij} y_{ji} = x_{ij} = y_{ji} = 1$, and for all $j \neq \sigma(i)$ we have that either $x_{ij} = 0$ or $y_{ji} = 0$. Since $x_{i \sigma(i)} \neq 0$ it follows by the above that $y_{\sigma(i) k} = 0$ for all $k \neq i$, and similarly since $y_{\sigma(k) k} \neq 0$ it follows by the above that $x_{i \sigma(k)} = 0$ for all $i \neq k$.
So we've learned that $X$ has exactly one nonzero entry in each row (and $Y$ has exactly one nonzero entry in each column). Since $X$ is an invertible matrix, these entries must be distinct, so $\sigma$ is a permutation and $X$ and $Y$ are permutation matrices as desired.
