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.
Is there anything wrong with the following argument?
First of all, by scaling all $x_i$ and all $y_i$ by a positive constant, we may safely assume that $\prod_i x_i = \prod_i y_i =1$.
The result would now follow from the following more general conjecture.
$\bf Conjecture:$ For ${\bf a}=(a_1,\ldots,a_{n-1})\in (\mathbb R_{>0})^{n-1}$ consider
$h_{\bf a}(z) = z^n+a_{n-1}z^{n-1}+\cdots+a_1z+1$. Define the function $f:(\mathbb R_{>0})^{n-1}\to {\mathbb R}$ by
$$
f(a_1,\ldots,a_{n-1}) = \sum_{z\vert h_{\bf a}(z)=0} (\log(-z))^2
$$
where the branch of $\log$ is picked as usual on the complement of ${\mathbb R}_{\leq 0}$ in $\mathbb C$. Then all partial derivatives $\frac{\partial f}{\partial a_k}$ are positive.
First, a few comments. The function $f$ is well-defined because none of the roots
of $h_a(z)$ are positive reals. In addition, $f$ is real-valued, because roots of $h_a(z)$ come in complex conjugate pairs for which values of $\log(-z)$ are complex conjugates.
The consequence of this conjecture is that if ${\bf a}\leq {\bf b}$ coordinate-wise,
then $f({\bf a})\leq f({\bf b})$. Applied to the case when roots of $h_{\bf a}(z)$
and $h_{\bf b}(z)$ are real numbers $-x_i$ and $-y_i$, we get the original conjecture.
Now, let me describe what I think is the proof of this new conjecture.
First of all, as standard, I can write the function $f(\bf a)$ by a contour integral. In a neighborhood of fixed ${\bf a}$ for any large enough $R$ and small enough $\epsilon>0$ there holds
$$
f({\bf a}) = \frac 1{2\pi i}\int_{C} (\log(-z))^2 \frac {h_{\bf a}'(z)}{h_{\bf a}(z)}\,dz
$$
over the contour $C$ on the union of $\mathbb C$ with two copies of $\mathbb R_{>0}$
which I will call "north" and "south" shores (so that $\log(-t)=\log t -\pi i$ on the north shore and $\log(-t)=\log t +\pi i$ on the south shore).
The contour $C$ is the union of the following four pieces $C_\epsilon$, $C_R$, $C_+$, $C_-$.
$C_\epsilon$ is a circle of radius $\epsilon$ around $z=0$ traveled from $\epsilon$-south
to $\epsilon$-north clockwise.
$C_R$ is a circle of radius $R$ around $z=0$ traveled from $R$-north to $R$-south counterclockwise.
$C_+$ is the line segment $[\epsilon,R]$-north.
$C_-$ is the line segment $[R,\epsilon]$-south.
The derivative $\frac{\partial f({\bf a})}{\partial a_k}$ is the integral of the
derivative, so we get:
$$
\frac{\partial f({\bf a})}{\partial a_k}= \frac 1{2\pi i}\int_{C} (\log(-z))^2 \frac \partial{\partial a_k}\frac {h_{\bf a}'(z)}{h_{\bf a}(z)}\,dz
$$
$$
=
\frac 1{2\pi i}\int_{C} (\log(-z))^2 \Big(\frac {z^k}
{h_{\bf a}(z)}
\Big)'\,dz
=
-\frac 1{2\pi i}\int_{C} \Big((\log(-z))^2 \Big)'\frac {z^k}
{h_{\bf a}(z)}
\,dz
$$
$$
=-\frac 1{\pi i}\int_{C} \log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz = -\frac 1\pi {\rm Im}\Big(\int_{C} \log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz\Big)
$$
We can take a limit as $R\to +\infty$ and $\epsilon\to 0$. Since $k\leq {n-1}$ the integral over $C_R$ goes to zero (the length is $2\pi R$ and the size of the function is
$O(R^{-2}\log R)$). The integral over $C_\epsilon$ also goes to zero, because the $k>=1$,
so the function is $O(\log \epsilon)$ and the length is $2\pi\epsilon$.
So we get
$$
\frac{\partial f({\bf a})}{\partial a_k}
= -\lim_{\epsilon\to 0^+}\frac 1\pi {\rm Im}\Big( \int_{[\epsilon,+\infty]-{\rm north}\cup [+\infty,\epsilon]-{\rm south}}
\log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz\Big)
$$
$$
=-\lim_{\epsilon\to 0^+}\frac 1\pi
\int_{\epsilon}^{+\infty}{\rm Im}\Big(
(\log(t)-\pi i)\frac {t^{k-1}}
{h_{\bf a}(t)} -(\log(t)+\pi i)\frac {t^{k-1}}
{h_{\bf a}(t)}\Big)\,dt
$$
$$
=2\lim_{\epsilon \to 0^+}
\int_{\epsilon}^{+\infty}
\frac {t^{k-1}}{h_{\bf a}(t)}\,dt
>0.$$
This finishes the proof of the new conjecture, and consequently of the old conjecture.
Remark:
I am guessing that there is a simpler argument for
$$
\frac{\partial f({\bf a})}{\partial a_k} = 2\int_{0}^{+\infty}
\frac {t^{k-1}}{h_{\bf a}(t)}\,dt
$$
but I am just writing the first thing that came to my mind.
Best Answer
I heard of three relatively recent works in that direction, taking different routes and arriving to interesting information about diagonal invariants.
First, there is the paper of Vaccarino that Darij mentioned in his comment (http://arxiv.org/abs/math/0205233).
Second, there are results of Domokos (http://arxiv.org/abs/0706.2154).
Third, there is a by-product of works of Buchshtaber and Rees on Frobenius n-characters, which also leads to new insight in that direction (http://www.ams.org/mathscinet-getitem?mr=2069166).
Hope that helps!