I finally got around to writing the promised details. I tried to make this a bit instructive, I hope you still find it useful.
First I will expand a bit on my comment above. A good reference is Stanley's article "Invariants of finite groups and their applications to combinatorics". There is a folklore theorem (which first appeared in print in M. Hochster and J. A. Eagon, Cohen-Macaulay rings, invariant theory, and the generic perfection of determinantal loci, Amer. J. Math. 93 (1971), 1020-1058) which says that for a finite subgroup $G$ of $Gl_n(\mathbb C)$ the algebra of invariants $\mathbb C[x_1,x_2,\dots,x_n]^G$ is Cohen-Macaulay. Therefore if $G$ is a subgroup of $G'$ and $G'$ is generated by pseudoreflections we get as a corollary of the Chevalley-Shephard-Todd theorem that $\mathbb C[x_1,\dots,x_n]^G$ is free over $\mathbb C[x_1,\dots,x_n]^{G'}$. In particular, this holds for $G'=S_n\times S_n$ and $G$ its diagonal subgroup.
Now, from a combinatorics perspective, we aren't simply satisfied by calculating the dimension of a polynomial algebra over another, but we would also like to exhibit a nice basis. I think it's worth spending sometime understanding the case of $R=\mathbb C[x_1,\dots,x_n]$ over $\mathbb C[x_1,\dots,x_n]^{S_n}$, first. Because the multivariate cases are similar in nature.
It's not hard to prove that the dimension of $R$ over $R^{S_n}$ is $n!$ and moreover there are two standard bases one learns about:
- The Artin Basis, consisting of monomials $x_1^{a_1}x_2^{a_2}\cdots x_n^{a_n}$, with $0\le a_i\le n-i$ for all $i$.
- Schubert polynomials.
Schubert polynomials are a nice basis, because they are indexed over combinatorial objects, and satisfy many combinatorial and geometric properties. The Artin basis, on the other hand, makes it easy to see that the Hilbert series is $(1+q)(1+q+q^2)\cdots (1+q+\cdots+q^{n-1})=[n]_q!$, yet it somewhat hides the presence of the symmetric group.
We know that $[n]_q!$ is the generating function of any Mahonian statistics on the symmetric group, so it would be nice to have a basis to reflect that. We can do this by the so called "descent" basis, and we will see that this is a construction that generalizes to the bivariate case (this is essentially the content of Bergeron and Lamontagne's paper).
The most famous Mahonian statistic is the Major index. Our basis is modeled after this statistic and is indexed over permutations. Since we have
$$\operatorname{maj} (\sigma) = \sum _{\sigma _{i+1} < \sigma _i} i ,$$
the most natural thing to try is the collection of monomials
$$b _{\sigma} = \prod _{i = 1} ^{n-1} (x _{\sigma_1} \cdots x _{\sigma _i}) ^{\chi (\sigma _i > \sigma _{i+1})}.$$
To prove that the $b_{\sigma}$'s form a basis it is enough to show that the polynomials $m_{\lambda}b_\sigma$ where $m_\lambda$ ranges over all symmetric monomials parametrized by partitions $\lambda$ are linearly independent. And the proof goes like this: We will construct a bijection $(\lambda,\sigma)\leftrightarrow \lambda$ where $\mu,\lambda$ are partitions with at most $n$ parts and $\sigma\in S_n$, and then use a Grobner type argument, i.e. show that the leading monomial in $m_{\lambda}b_{\sigma}$ is precisely $X^{\mu}=x_1^{\mu_1}\cdots x_n^{\mu_n}$. You will find this argument spelled out in detail in section 2 of Allen's "The Descent Monomials and a Basis for the Diagonally Symmetric Polynomials".
To guess a basis for the bivariate case Bergeron and Lamontagne, play a similar game, where they first calculate the Hilbert series of $R^{S_n}$ over $R^{S_n\times S_n}$. They actually calculate the Frobenius series and obtain the expression $$(q;q)_n(t;t)_n h_n\left[\frac{1}{(1-q)(1-t)}\right]$$
in plethystic notation. But this is a well known generating function over $S_n$. Namely it is
$$\sum _{\sigma\in S _n} q^{\operatorname{maj}(\sigma)}t^{\operatorname{maj}(\sigma^{-1})}.$$
So they construct a basis $$B _{\sigma}=\rho b _{\sigma}(X)b _{\sigma^{-1}}(Y)$$
similar to the construction above, where $\rho$ is the Reynolds operator. To be able to use a Grobner type argument, they construct a bijection $(\lambda _1,\lambda _2,\sigma)\leftrightarrow (\mu _1,\mu _2)$ (section 12), and they are able to show that the polynomials $m _{\lambda _1}(X)m _{\lambda _2}(Y)B _{\sigma}(X,Y)$ are linearly independent (section 13).
Finally, a word on the case of diagonal coinvariants. It is a big open problem to exhibit a basis for the space of diagonal coinvariants, $\mathbb C[X,Y]$ over $\mathbb C[X,Y]^{S_n}$. Even though, there is a conjectured form for the Hilbert series as a generating function of two statistics over parking functions (the dimension here is no longer $n!$, rather $(n+1)^{n-1}$, proved by Haiman in 2001), these statistics are not natural enough to let one guess what the corresponding basis will be.