Start with $n = 2$. An element of $M(2,k)$ is determined by a single $a \in \{0,\ldots,k\}~$ for $k \geq 2$. Since we wish to avoid $a=0,k~$ it is clear that $P(2,k) = \frac{k-1}{k+1}~$ which is certainly non-decreasing in $k$. Keep in mind the valid choices of $a$: these come from $[0,k] \cap \mathbb{Z}~$. Observe that $[0,k]$ is the standard simplex $S$ of side $k$ in dimension $1$ and that the $0$ entries correspond to the boundary of this simplex.

The case $n=3$ is much more interesting: the matrices in $M(3,k)$ correspond to four integers $a_{11},a_{12},a_{21},a_{22}~$ from $0,\ldots,k~$ such that

1. For $i = 1,2~$ the sum $\sum_j a_{ij} \leq k$
2. For $j = 1,2~$ the sum $\sum_i a_{ij} \leq k$
3. And finally, $\sum_{ij} a_{ij} \geq k$.

We get a zero entry if either one of the $a_{ij}$ is zero or if one of the inequalities above is an equality. Each of these situations corresponds to the point with coordinates $a_{ij}$ lying on a boundary face of the four dimensional "prism" which consists of the cube $[0,k]^4$ intersected with the half space from the third inequality.

So bottom line: I think that the ratio $$\frac{\text{integral points on the boundary of this prism}}{\text{ the total number of integral points in this prism}}$$ decreases as one increases $k$, which leads to your desired result. The process generalizes to higher dimensions.

Okay, so this is nowhere near a complete solution, but this is as far as I got and hopefully someone else sees it from here:

It is clear that $P(2,k) = \frac{k-1}{k+1}~$ which is certainly non-decreasing in $k$.

When $n=3$, the matrices in $M(3,k)$ correspond to four integers $a_{11},a_{12},a_{21},a_{22}~$ from $0,\ldots,k~$ such that

1. For $i = 1$ or $2~$, $\sum_j a_{ij} \leq k$
2. For $j = 1$ or $2~$, $\sum_i a_{ij} \leq k$
3. And finally, $\sum_{ij} a_{ij} \geq k$.

These hyperplane inequalities carve out a convex region $C \subset \mathbb{R}^4$ from the cube $[0,k]^4$ and the "zero entry" cases of $M(3,k)$ are precisely the bounding faces of this region.

So, if a theorem establishes that the ratio $$\frac{\text{integral points on the boundary of } C}{\text{ total number of integral points in }C}$$ decreases as one increases $k$, then we obtain your desired result. I don't know enough about convex polytopes to cite something here but it sounds reasonable just from dimension considerations...

Ideally, this process would generalize to higher dimensions. An element of $M(n,k)$ has zero entries if and only if the vector of entries in the first $(n-1) \times (n-1)$ block lies in the boundary of convex polytope carved from the cube $[0,1]^{(n-1)^2}$ by $2n-1$ hyperplanes.

Start with $n = 2$. An element of $M(2,k)$ is determined by a single $a \in \{0,\ldots,k\}~$ for $k \geq 2$. Since we wish to avoid $a=0,k~$ it is clear that $P(2,k) = \frac{k-1}{k+1}~$ which is certainly non-decreasing in $k$. Keep in mind the valid choices of $a$: these come from $[0,k] \cap \mathbb{Z}~$. Observe that $[0,k]$ is the standard simplex $S$ of side $k$ in dimension $1$ and that the $0$ entries correspond to the boundary of this simplex.

The case $n=3$ is much more interesting: the matrices in $M(3,k)$ correspond to four integers $a_{11},a_{12},a_{21},a_{22}~$ from $0,\ldots,k~$ such that

1. For $i = 1,2~$ the sum $\sum_j a_{ij} \leq k$
2. For $j = 1,2~$ the sum $\sum_i a_{ij} \leq k$
3. And finally, $\sum_{ij} a_{ij} \geq k$.

We get a zero entry if either one of the $a_{ij}$ is zero or if one of the inequalities above is an equality. Each of these situations corresponds to the $a_{ij}$ lying on a face of the four dimensional simplex which consists of the cube $[0,k]^4$ minus the open simplex which comes from reversing the third inequality above.

So bottom line: I think that if the ratio $$\frac{\text{integral points on the boundary of this simplex}}{\text{ the total number of integral points in this simplex}}$$ decreases as one increases $k$, which leads to your desired result. The process generalizes to higher dimensions.

Start with $n = 2$. An element of $M(2,k)$ is determined by a single $a \in \{0,\ldots,k\}~$ for $k \geq 2$. Since we wish to avoid $a=0,k~$ it is clear that $P(2,k) = \frac{k-1}{k+1}~$ which is certainly non-decreasing in $k$. Keep in mind the valid choices of $a$: these come from $[0,k] \cap \mathbb{Z}~$. Observe that $[0,k]$ is the standard simplex $S$ of side $k$ in dimension $1$ and that the $0$ entries correspond to the boundary of this simplex.

The case $n=3$ is much more interesting: the matrices in $M(3,k)$ correspond to four integers $a_{11},a_{12},a_{21},a_{22}~$ from $0,\ldots,k~$ such that

1. For $i = 1,2~$ the sum $\sum_j a_{ij} \leq k$
2. For $j = 1,2~$ the sum $\sum_i a_{ij} \leq k$
3. And finally, $\sum_{ij} a_{ij} \geq k$.

We get a zero entry if either one of the $a_{ij}$ is zero or if one of the inequalities above is an equality. Each of these situations corresponds to the point with coordinates $a_{ij}$ lying on a boundary face of the four dimensional "prism" which consists of the cube $[0,k]^4$ intersected with the half space from the third inequality.

So bottom line: I think that the ratio $$\frac{\text{integral points on the boundary of this prism}}{\text{ the total number of integral points in this prism}}$$ decreases as one increases $k$, which leads to your desired result. The process generalizes to higher dimensions.