Diagonal invariants of the symmetric group on $k[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n]$ This sounds like something that must have been answered long ago, but for some reason I can find nothing on it in the internet. (There has been lots of recent activity in diagonal covariants, related to the $n!$ conjecture, but invariants seem to have become a stepchild in this process.)
Let $k$ be a field of characteristic $0$. Let $n\in\mathbb N$. The group $S_n\times S_n$ acts on the polynomial ring $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]$ by
$\left(\sigma,\tau\right)\left(P\right) = P\left(X_{\sigma\left(1\right)},X_{\sigma\left(2\right)},...,X_{\sigma\left(n\right)},Y_{\tau\left(1\right)},Y_{\tau\left(2\right)},...,Y_{\tau\left(n\right)}\right)$.
Thus, the symmetric group $S_n$ also acts on $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]$ due to the diagonal embedding $S_n\to S_n\times S_n$.
Since the action of $S_n$ is not generated by pseudoreflections, it follows from the converse of the Chevalley-Shephard-Todd theorem that $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]$ is not a free $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]^{S_n}$-module. But it is easy to see that $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]$ is a free $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]^{S_n\times S_n}$-module of rank $n!^2$.
Question: Is $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right]^{S_n}$ a free $k\left[X_1,X_2,...,X_n,Y_1,Y_2,...,Y_n\right] ^ {S_n\times S_n}$-module?
 A: 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.
A: I believe the answer to your question is yes.
The Reynolds operator 
$$
  R: k[X_1,\ldots,X_n,Y_1,\ldots,Y_n] \longrightarrow k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n}
$$
is $k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n}$-equivariant, and therefore in particular 
$k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n\times S_n}$-equivariant. Also, $R$ is a projection, and therefore the image of $R$ is a direct summand of of its domain of definition.
It follows therefore that the $k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n\times S_n}$-module $ k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n}$ is a direct summand of the free $k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]^{S_n\times S_n}$-module $k[X_1,\ldots,X_n,Y_1,\ldots,Y_n]$.
Now since the action of $S_n\times S_n$ is generated by reflections, its invariant ring is a polynomial ring, and therefore any direct summand of a free module over this ring is again free (this is the Quillen–Suslin theorem).

Edit (concerning the rank, which is $n!$): Let $G$ be a finite group acting on $R=k[x_1,\ldots,x_n]$ and let $H\leq G$ be a subgroup. Then 
$$
   {\rm frac}(R^G)\otimes_{R^G} R^H \cong (R^G-\{0\})^{-1}R^H \cong {\rm frac}(R^H)
$$
where ${\rm frac}$ denotes the fraction field (this is easily seen from the fact that by a construction similar to the Reynols operator, the denominator of a fraction over $R^H$ can always be made $R^G$-invariant). Let $Q={\rm frac}(R)$, then one can show in a similar fashion $Q^G = {\rm frac}(R^G)$. We conclude that
$$
   {\rm rank}_{R^G} R^H = {\rm dim}_{{\rm frac}(R^G)} {\rm frac} (R^H) = [Q^H:Q^G] = [G:H]
$$ 
the last equality being due to Galois theory.
