Warning: I have not carefully checked the argument below.
First note that $\text{SL}_n(\mathbb{Z}_p)$ contains row operations of determinant $1$ (adding a multiple of a row to another row, transposing two rows and negating one of them, and multiplication by a diagonal matrix of determinant $1$). Given $M \in \mathcal{M}_n(\mathbb{Z}_p)$ of determinant $p$, we will row reduce it: if the entries in the first column are not all divisible by $p$, one of them is invertible, so by transposing and multiplying by a diagonal matrix we can place it in the first row and set it equal to $1$ to pivot on it. We can continue row reducing until arriving at a column where the only possible pivots are divisible by $p$ (this must eventually occur by the determinant condition), and then we just pivot on the one with minimal $p$-adic valuation (it must be precisely divisible by $p$ by the determinant condition). By multiplying by a diagonal matrix we can set this pivot to $p$, and proceed with row reduction.
The final result should be a matrix all of whose diagonal entries are $1$ except for a single diagonal entry of $p$ somewhere. In every column except the one containing $p$ all of the other entries are zero, and in the column containing $p$ the entries above it range from $0$ to $p-1$ and the entries below it are also zero. I think a tedious calculation will now confirm that if $A, B$ are distinct matrices in this form then $AB^{-1} \not \in \text{SL}_n(\mathbb{Z}_p)$. For $n = 2$ the representatives this algorithm gives are
$$\left[ \begin{array}{cc} p & 0 \\ 0 & 1 \end{array} \right], \left[ \begin{array}{cc} 1 & b \\ 0 & p \end{array} \right], b = 0, 1, ... p-1.$$
(If I'm not horribly mistaken, what you've written down are representatives for the orbits of the action by right multiplication.) This should generalize to a notion of Hermite normal form for matrices over a PID.
