Number of permutations in $S_{a+b}$ with $\operatorname{maj}(\pi)=a$ and $\operatorname{maj}(\pi^{-1})=b$ $\DeclareMathOperator\maj{maj}\DeclareMathOperator\inv{inv}$Major index, $\maj$, of a permutation on $1,2,\dotsc,n$ is defined as
$$
\maj(\pi) \mathrel{:=} \sum_{i=1}^{n-1} i \cdot \chi(\pi(i)\gt \pi(i+1))
$$
where $\chi$ is 1 if the statement inside is true, 0 otherwise.
Let $t_{a,b}$ be the numbers
$$
 t_{a,b} \mathrel{:=} \lvert\{ \pi \in S_{a+b} : \maj(\pi)=a \text{ and }  \maj(\pi^{-1})=b \}\rvert.
$$
Here, $S_{a+b}$ denotes the set of permutations of $1,2,\dotsc,a+b$.
By a result of Foata, one can also look at the pair of statistics $(\maj, \inv)$, and a few other combinations — these pairs of statistics will produce the same numbers.
Now, according to the OEIS entry A090806, it is proved by Garsia–Gessel that
\begin{equation}
 \sum_{a,b} t_{a,b} q^a t^b = \prod_{i,j \geq 1} \frac{1}{1-q^i t^j}
\qquad (\ast)
\end{equation}
I cannot see exactly where in their paper one can deduce this.
I have tried to prove this myself (mainly by resorting to RSK, the Cauchy identity, and some symmetric function identities).
This leads to the following (which appears in Stanley's EC2):
$$
 \sum_{n \geq 0} \frac{z^n}{(1-q)^n[n]_q!(1-t)^n [n]_t!} \sum_{\pi \in S_n} t^{\maj(\pi)} q^{\maj(\pi^{-1})} 
 =
 \prod_{i,j \geq 0} \frac{1}{1-z q^i t^j},
$$
where $[n]_q! \mathrel{:=} [1]_q [2]_q \dotsm [n]_q$, and $[n]_q = 1+q+q^2+\dotsb + q^{n-1}$.
However, I do not see some short way to deduce the above generating function from this.
Question: Is there some alternative (more recent?) reference where $(\ast)$ is
stated and easily referenced? Alternatively, someone who can see exactly where in the paper $(\ast)$ is proven?
Garsia, A. M.; Gessel, I., Permutation statistics and partitions, Adv. Math. 31, 288-305 (1979). ZBL0431.05007.
 A: Here is a derivation of $(\ast)$ from the displayed equation
$$\sum \frac{z^n}{(1-q)^n [n]_q! (1-t)^n [n]_t!} \sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \prod_{i,j \geq 0} \frac{1}{1-zt^i q^j}.$$
Taking the coefficient of $z^n$ on both sides, we have
$$\frac{\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})}}{\prod_{i=1}^n (1-q^i) \prod_{j=1}^n (1-t^j)} = h_n(\{t^i q^j : i,j \geq 0\}).$$
Here the RHS is the complete homogenous symmetric function evaluated on the set of monomials $\{t^i q^j : i,j \geq 0\}$.
Now, $h_n(1, u_1, u_2, u_3, \ldots) = \sum_{k=0}^n h_k(u_1, u_2,\ldots).$ So we can rewrite this RHS to get
$$\frac{\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})}}{\prod_{i=1}^n (1-q^i) \prod_{j=1}^n (1-t^j)} = \sum_{k=0}^n h_k(\{ q^i t^j : i,j \geq 0,\ (i,j) \neq (0,0) \}). \qquad (\clubsuit)$$
According to the OEIS entry, the quantity $\# \{ \pi \in S_n : \mathrm{maj}(\pi) = a,\ \mathrm{maj}(\pi^{-1}) = b \}$ stabilizes at $t_{ab}$ as $n \to \infty$; the OP's mention of $n=a+b$ is just a particular value in the stable range. So we can take the limit of both sides of $(\clubsuit)$ as $n \to \infty$ to get
$$\frac{\sum_{a,b \geq 0} t_{ab} q^a t^b}{\prod_{i=1}^{\infty} (1-q^i) \prod_{j=1}^{\infty} (1-t^j)} = 
\sum_{k=0}^{\infty} h_k(\{ q^i t^j : i,j \geq 0,\ (i,j) \neq (0,0) \}) = \prod_{\substack{i,j \geq 0 \\ (i,j) \neq (0,0)}} \frac{1}{1-q^i t^j} . \qquad (\diamondsuit)$$
Now cancel common factors from both sides of $(\diamondsuit)$ to get the claim.

Here is another approach, which suggests that something more interesting may be going on. Recall that the RS correspondence is a bijection between the symmetric group $S_n$ and pairs of SYT $(T,U)$ of the same shape $\lambda$, where $|\lambda| = n$. We define $\mathrm{maj}(T)$ for an SYT $T$ to be the sum of those $i$ such that $i$ occurs in a strictly higher row of $T$ than $i+1$ does. Then, if RS maps $\pi$ to $(T,U)$, we have
$\mathrm{maj}(\pi)= \mathrm{maj}(T)$ and $\mathrm{maj}(\pi^{-1})= \mathrm{maj}(U)$.
Define $f^{\lambda}(q)$ to be the sum, over $SYT$ of shape $\lambda$, of $q^{\mathrm{maj}(T)}$. So
$$\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \sum_{|\lambda| = n} f^{\lambda}(q) f^{\lambda}(t).$$
Given a permutation $\mu$, and $n > |\mu| + \mu_1$, let $\mu[n]$ be the partition $(n-|\mu|, \mu_1, \mu_2, \ldots, \mu_k)$. So every partition of $n$ is of the form $\mu[n]$ for a unique $\mu$ and that, for any $\mu$, the partition $\mu[n]$ is well defined for $n$ large enough.
It is easy to see that $\lim_{n \to \infty} f^{\mu[n]}(q)$ exists. Set $f^{\mu[\infty]}(q)$ to be $\lim_{n \to \infty} f^{\mu[n]}(q)$.
So we have
$$\lim_{n \to \infty} \sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \sum_{\mu} f^{\mu[\infty]}(q) f^{\mu[\infty]}(t).$$
On the other hand, the Cauchy identity gives
$$\prod_{i,j \geq 1} \frac{1}{1-q^i t^j} = \sum_{\mu} s_{\mu}(q,q^2, \cdots) s_{\mu}(t,t^2, \cdots).$$
So here is the weird thing: It turns out that $f^{\mu[\infty]}(q)$ and $s_{\mu}(q)$ are the same thing! This seems like it should have a combinatorial proof. We can clearly think of $f^{\mu[\infty]}(q)$ as a generating function for tableau of shape $\mu$ with distinct entries, counted by a variant of major index. And $s_{\mu}(q)$ is the generating function for semistandard tableau of shape $\mu$, counted by weight. It feels like there should be an easy bijection here.
Well, I couldn't find one. But it isn't hard to prove the equality using hook length formulas. Let $m = |\mu|$ and let $h_1$, $h_2$, ..., $h_m$
be the hook lengths of $\mu$. Set $M = \sum (i-1) \mu_i$. We have
$$s_{\mu}(q) = \frac{q^M}{\prod_i (1-q^{h_i})} .$$
On the other hand, the hook lengths of $\mu[n]$ are $h_1$, $h_2$, ..., $h_m$ together with an $n-m$ additional hook lengths $S$. The exact set $S$ doesn't matter; what is important is  $\{1,2,\ldots, n-\mu_1-m \} \subseteq S \subseteq \{1,2,\ldots, n \}$. So, by a formula of Stanley (EC2, Cor. 7.21.5), we have
$$f^{\mu[n]}(q) = \frac{q^M (1-q)(1-q^2) \cdots (1-q^N)}{\prod_{s \in S} (1-q^s) \prod_i (1-q^{h_i})}.$$
In the limit as $n \to \infty$, both $(1-q)(1-q^2) \cdots (1-q^N)$ and $\prod_{s \in S} (1-q^s)$ approach $\prod_{k=1}^{\infty} (1-q^k)$, so they cancel and we are left with
$$ f^{\mu[\infty]}(q) = \frac{q^M}{\prod_i (1-q^{h_i})} = s_{\mu}(q,q^2, \cdots).$$
