$\def\des{\operatorname{des}}$Let $\des(\pi)$ be the number of descents of the permutation $\pi$. Then for any permutation $\pi$ in $S_k$, we have $$\binom{xy+k-\des(\pi)-1}{k} =\sum_{\sigma\tau=\pi}\binom{x+k-\des(\tau)-1}{k} \binom{y+k-\des(\sigma)-1}{k}.\tag{$*$}\label{star}$$
If $\pi$ is the identity then $\des(\pi)=0$, so the left side of \eqref{star} is $\binom{xy+k-1}{k}$ and on the right, $\tau=\sigma^{-1}$. The OP's identities correspond to the fact that if $k\le 3$ then the number of descents of $\pi$ is the same as the number of descents of $\pi^{-1}$, but this is not true for $k>3$. This also explains the occurrence of the Eulerian numbers as coefficients for $k\le3$.
As far as I know, the earliest proof of \eqref{star} is in Bogdan Mielnik and Jerzy Plebański, Combinatorial approach to Baker–Campbell–Hausdorff exponents, Annales de l’I. H. P., section A, tome 12, no 3 (1970), pp. 215–254 (see equation (11.10)); further references can be found in Jason Fulman and T. Kyle Petersen, Card shuffling and P-partitions, arXiv:2004.01659 [math.CO].