Expand $\binom{xy}{n}$ in terms of $\binom{x}{k}$'s and $\binom{y}{k}$'s

The identity $\binom{x+y}{n} = \sum_{k=0}^{n} \binom{x}{k} \binom{y}{n-k}$ is well-known, it is called the Vandermonde identity. The answer for $\binom{xy}{n}$ can be explained using the notion of a $\lambda$-ring, where here we consider the binomial ring $\mathbb{Z}$ with $\lambda^n(x)=\binom{x}{n}$. The main theorem on symmetric polynomials enables us to write $$\sum_{k=0}^{n} P_k(\sigma_1,\dotsc,\sigma_k,\tau_1,\dotsc,\tau_k) \,t^k = \prod_{i,j=1}^{n} (1+ t x_i y_j)$$ for some polynomials $P_k \in \mathbb{Z}[x_1,\dotsc,x_k,y_1,\dotsc,y_k]$, where $\sigma_1,\dotsc,\sigma_n$ are the elementary symmetric polynomials in $x_1,x_2,\dotsc,x_n$ and $\tau_1,\dotsc,\tau_n$ are the elementary symmetric polynomials in $y_1,\dotsc,y_n$. Then, one has (by definition) $$\lambda^n(xy)=P_n\bigl(\lambda^1(x),\dotsc,\lambda^n(x),\lambda^1(y),\dotsc,\lambda^n(y)\bigr).$$ For example: $$\lambda^2(xy) = x^2 \lambda^2(y) + \lambda^2(x) y^2 - 2 \lambda^2(x) \lambda^2(y)$$ $$\lambda^3(xy)=x^3 \lambda^3(y) + \lambda^3(x) y^3 + x \lambda^2(x) y \lambda^2(y) - 3 x \lambda^2(x) \lambda^3(y) - 3 \lambda^3(x) y \lambda^2(y) + 3 \lambda^3(x) \lambda^3(y)$$ Of course, one has to prove that every binomial ring becomes a $\lambda$-ring. See for instance Darij Grinberg's notes, Theorem 7.2 (which is a corollary of Theorem 7.1).

A beautiful combinatorial answer was found by Gjergji Zajmi. For an expanded version, see the solution to Exercise 3.9 in my Notes on the combinatorial fundamentals of algebra. (At least, it is Exercise 3.9 in the version of 10 January 2019. In future versions, the numbering can shift.)

Unlike Martin's answer, this one directly gives a formula for $\dbinom{xy}{n}$ as a nonnegative linear combination of products $\dbinom{x}{k}\dbinom{y}{\ell}$, as opposed to a polynomial in the $\dbinom{x}{k}$ and the $\dbinom{y}{\ell}$. This, of course, does not work for arbitrary $\lambda$-rings.