Important combinatorial and algebraic interpretations of the coefficients in the polynomial $[n]!_q = (1+q)(1+q+q^2) \ldots (1+q+\cdots + q^{n-1})$ 
What are some important combinatorial and algebraic interpretations of the coefficients in the polynomial 
$$[n]!_q = (1+q)(1+q+q^2) \ldots (1+q+\cdots + q^{n-1})?$$

As motivation, I will give three interpretations, ask for a fourth, and raise a related question about the unimodality. I would be particularly interested in answers using the RSK-correspondence or subspaces of $\mathbb{F}_q^n$.


*

*Given a permutation $\sigma \in \mathrm{Sym}_n$, let $\mathrm{inv}(\sigma)$ denote the number of inversions of $\sigma$; that is, pairs $(x,y)$ with $x < y$ and $\sigma(x) > \sigma(y)$. Then $[n]!_q = \sum_{\sigma \in \mathrm{Sym}_n} q^{\mathrm{inv}(\sigma)}$.

*An element $x \in \{1,\ldots, n-1\}$ is a descent of $\sigma \in \mathrm{Sym}_n$ if $\sigma(x) > \sigma(x+1)$. The major index $\mathrm{maj}(\sigma)$ is the sum of the descents of $\sigma$. Then $[n]!_q = \sum_{\sigma \in \mathrm{Sym}_n} q^{\mathrm{maj}(\sigma)}$. I think this is due to MacMahon.

*In the 'inside-out' version of the Fisher–Yates shuffle on an $n$-card deck, at step $j-1$, card $j-1$ from the top is swapped with one of cards in positions $0, 1, \ldots, j-1$ from the top, chosen uniformly at random. These choices are enumerated by $1 + q + \cdots + q^{j-1}$. After $n$ steps (starting with $j=1$), each permutation has equal probability.  (This is essentially coset enumeration in the symmetric group by the chain $\mathrm{Sym}_1 \le \mathrm{Sym}_2 \le \ldots \le \mathrm{Sym}_n$.) Hence $[n]!_q$ enumerates permutations according to the sum of the positions chosen at each stage. 

Does the normal Fisher–Yates shuffle have a similar combinatorial interpretation? Is there a more natural interpretation of the $q$-power, still using the inside-out Fisher–Yates shuffle?

Finally, (1) makes it easy to see that $[n]!_q$ is symmetric, i.e. the coefficients of $q^m$ and $q^{\binom{n}{2}-m}$ are the same: use the Coxeter involution, thinking of $[n]_q!$ as the Poincaré series of the Coxeter group $\mathrm{Sym}_n$. This can also be seen in a similar way from (2). But it does not seem to be obvious from (3).

Which interpretation is the best way to show that the coefficients in $[n]!_q$ are unimodal, i.e. first increasing then decreasing?

 A: In FindStat, the database of combinatorial statistics, the major index of a permutation is http://www.findstat.org/StatisticsDatabase/St000004/. By clicking "search for distribution" there, you will find other statistics that (FindStat thinks) are equidistributed with major index. I get 41 results. As mentioned in a comment above, these statistics are called "Mahonian" in honor of MacMahon.
A: This answer concerns an interpretation of $[n]!_q$ coming from invariant theory.
There is a natural action of $S_n$ on $V=\mathbb{C}^n$, and hence on polynomial functions on $V$. This action preserves degree. The coinvariant algebra is the quotient of the polynomial functions on $V$ by (the ideal generated by) the $S_n$-invariant functions of positive degree. This (finite-dimensional!) algebra inherits an $S_n$ representation. As an ungraded $S_n$ module it is isomorphic to the regular representation (this was first proved by Chevalley). But it also inherits a grading, giving a natural grading to the regular representation of the symmetric group. Various fundamental results in invariant theory (the Chevalley-Shephard-Todd theorem, Molien's theorem, etc.) imply that the Hilbert series of the coinvariant algebra is $[n]!_q$. Hence the coefficient of $q^i$ in $[n]!_q$ is the dimension of the coinvariant algebra in degree $i$.
There is a connection here to RSK, and since the question asked about RSK, let me explain it.
Let me use $R_n$ to denote the coinvariant algebra. Since as an ungraded algebra, $R_n$ is the regular representation of $S_n$, we have $R_n \simeq \bigoplus_{\lambda \vdash n} V_{\lambda}^{f^{\lambda}}$, where $V_{\lambda}$ is the $S_n$-irrep indexed by the partition $\lambda$, and $f^{\lambda} := \mathrm{dim}(V_{\lambda})$. Recall that also $f^{\lambda} = \#$ Standard Young Tableaux (SYTs) of shape $\lambda$, and that we have the famous hook-length formula $f^{\lambda} = \frac{n!}{\prod_{u\in \lambda}h(u)}$, where $h(u)$ is the hook-length of the box $u\in \lambda$.
But we again get more by considering the grading on $R_n$. Namely, for a partition $\lambda\vdash n$, define the fake degree polynomial $f^{\lambda}(q)$ by
$$f^{\lambda}(q) := \sum_{i \geq 0} \# \textrm{ (copies of $V_{\lambda}$ in the degree $i$ component of $R_n$)} \cdot q^i.$$
Note that $f^{\lambda}(q)$ is a $q$-analog of $f^{\lambda}$. Stanley (see Section 4 of the paper cited below and Enumerative Combinatorics, Vol. 2, 7.21.5) proved a $q$-analog of the hook-length formula: $f^{\lambda}(q) = q^{b(\lambda)} \frac{[n]!_q}{\prod_{u\in \lambda}[h(u)]}$, where $b(\lambda) := 0\lambda_1 + 1\lambda_2 + 2\lambda_3 + \cdots$. Furthermore, Lusztig gave a statistical formula for $f^{\lambda}(q)$. Namely, for an SYT $T$, a  descent of $T$ is an entry $i$ for which $i+1$ is in a strictly lower row, and the major index of $T$, denoted $\mathrm{maj}(T)$, is a sum over all descents $i$ of $T$ of the value $i$. Then in unpublished work (again see Stanley's paper cited below), Lusztig showed that $f^{\lambda}(q) = \sum_{\substack{T \textrm{ a SYT}, \\ \mathrm{sh}(T) =\lambda}}q^{\mathrm{maj}(T)}$.
Now note that the decomposition $R_n \simeq \bigoplus_{\lambda \vdash n} V_{\lambda}^{f^{\lambda}}$ implies $[n]!_q = \sum_{\lambda\vdash n}f^{\lambda} \cdot f^{\lambda}(q)$, and hence due to Lusztig's result that the coefficient of $q^i$ in $[n]!_q$ is the sum over all SYT $T$ with $n$ boxes and with $\mathrm{maj}(T)=i$ of $f^{\mathrm{sh}(T)}$.
Here's where the connection to RSK comes in. Recall that the usual Robinson-Schensted (no Knuth!) algorithm is a bijection between permutations $w \in S_n$ and pairs $(P,Q)$ of n-boxed SYTs of the same shape. Usually $P$ is called the insertion tableau and $Q$ the recording tableau. It is not hard to see from the usual "insertion and bumping" description of R-S that the descents of $w$ are exactly the same as the descents of its recording tableau $Q$. In particular, their major indices are the same. Thus R-S proves bijectively that
$$ \sum_{w \in S_n} q^{\mathrm{maj}(w)} = \sum_{\substack{T \textrm{ a SYT} \\ \textrm{with $n$ boxes}}} f^{\mathrm{sh}(T)} \cdot q^{\mathrm{maj}(T)},$$
where the RHS is the interpretation of $[n]!_q$ we just mentioned, and the LHS is the well-known one due to MacMahon mentioned in the original question.
Stanley, Richard P., Invariants of finite groups and their applications to combinatorics, Bull. Am. Math. Soc., New Ser. 1, 475-511 (1979). ZBL0497.20002.
A: For a simple but noncombinatorial proof of the unimodality of the coefficients of $[n]!_q$, see 
George E. Andrews, 
A theorem on reciprocal polynomials with applications to permutations and compositions.  Amer. Math. Monthly 82 (1975), no. 8, 830–833.
A: There is a nice generalization to Catalan objects as follows.
Acyclic orientations of unit-interval graphs
Note that for the complete graph on $n$ vertices,
the number of acyclic orientations weighted by ascents (edges oriented from smaller to larger label) is given by $[n]_q!$.
Now, there are a Catalan number of unit interval graphs, which can be indexed by area sequences. The complete graph has area sequence $(0,1,2,\dotsc,n-1)$.
The weighted sum over acyclic orientations of an unit interval graph with area sequence $(a_1,\dotsc,a_n)$ is given by the product $[a_1+1]_q\dotsb [a_n+1]_q$.
Summing all these polynomials over all $(n+1)^{-1}\binom{2n}{n}$ area sequences of length $n$ gives the $n$th Touchard-Riordan polynomial, which is defined as the sum over all perfect matchings of $2n$ vertices on a circle, weighted by $q$ to the number of crossings.
Perfect matchings
In particular, $[n]_q!$ is the sum over all perfect matchings on $2n$ vertices, where vertices $1,2,\dotsc,n$ are matched with $n+1,n+2,\dotsc,2n$.
Rook placements on Ferrers board
Let $a$ be an area sequence. 
A Ferrers board with row lengths given by
$a_i+n-i$ for $i=1,\dotsc,n$ allows for $[a_1+1]_q\dotsb [a_n+1]_q$ different ways to place $n$ non-attacking rooks,
and the weight is now given by rook inversions.
I have some more info and references here.
A: After posting the question, I remembered an interpretation using formal characters of representations of $\mathrm{SL}_2(\mathbb{C})$ that gives the shortest proof I know that $[n]!_q$ is unimodal.
Write $1+q+\cdots + q^{j-1}$ as $q^{(j-1)/2}(q^{-(j-1)/2} + q^{-(j-3)/2} + \cdots + q^{(j-1)/2})$ and substitute $q= Q^2$ to get $Q^{j-1}(Q^{-(j-1)} + Q^{-(j-3)} + \cdots + Q^{j-1})$. Hence
$$[n]!_q = Q^{\binom{n}{2}}( Q^{-1} + Q) ( Q^{-2} + 1 + Q^2) \ldots (Q^{-(n-1)} + Q^{-(n-3)} + \cdots + Q^{n-1}).$$
Up to the power of $Q$, the right-hand side is the formal character of the tensor product of the (unique) irreducible representations of $\mathrm{SL}_2(\mathbb{C})$ of dimensions $2, 3, \ldots, n$. Since any formal character is a sum of the characters of irreducible representations, and each irreducible representation has a unimodal character, this proves that $[n!]_q$ is unimodal in $Q$, and hence in $q$.
A: This answer concerns a geometric/Lie-theoretic interpretation of $[n]!_q$.
$[n]!_q$ gives the number of points in the full flag variety of full flags of subspaces in an $n$-dimensional vector space $\mathbb{F}_q^n$ over the finite field $\mathbb{F}_q$.
Recall that the full flag variety (over any field) has a natural stratification, the Bruhat stratification. Due to the above point-counting remark, it follows that the coefficient of $q^i$ in $[n]!_q$ is the number of $i$-dimensional cells in the Bruhat stratification.
There is also a way to deduce the unimodality of the coefficients from this geometric perspective. Namely, the partial order on the Bruhat cells whereby $C \leq C'$ if $C$ is contained in the closure of $C'$ is called the Bruhat order or strong order. The strong order can be viewed as an order on the symmetric group because the Bruhat cells are naturally labeled by permutations. Strong order is graded, and the rank sizes are precisely the coefficients of $[n]!_q$ (i.e., the number of permutations with given inversion number). Richard Stanley showed in the "Weyl groups..." paper cited below that in this situation (when you have a complex projective variety with a cellular decomposition satisfying certain conditions), the poset in question is necessarily graded, rank-symmetric, rank-unimodal, and strongly Sperner, which in particular implies the unimodality of the coefficients of $[n]!_q$. His proof employed the hard Lefschetz theorem and so can hardly be called elementary, but it is conceptual.
Stanley, Richard P., Weyl groups, the hard Lefschetz theorem, and the Sperner property, SIAM J. Algebraic Discrete Methods 1, 168-184 (1980). ZBL0502.05004.
A: The polynomials have nonnegative, symmetric, unimodal coefficients. So the the product of polynomials also has this property.
