A beautiful result due to Evans and Gillis is that the function
$$A(n_1,n_2,\cdots,n_r)=\int_0^\infty L_{n_1}(x)L_{n_2}(x)\cdots L_{n_r}(x)e^{-x}dx$$
counts the number of generalized derangements (up to a sign, I believe). Here, $L_{n}(x)$ is a Laguerre polynomial There's a nice application of it here. For example the usual derangements are counted by:
$$(-1)^n\int_0^\infty L_1^n(x)e^{-x}dx.$$
Messing around with Mathematica, it seems like one can try using other polynomials with their respective weights to produce integer valued $A(n_1,n_2,\cdots,n_r)$. For example, Hermite Polynomials:
$$\int_{-\infty}^\infty H_{n_1}(x)\cdots H_{n_r}(x)e^{-x^2}dx$$
which tends to give either an integer, or an integer multiple of $\sqrt{\pi}$. For example:
$$\int_{-\infty}^\infty H_{1}(x)H_{2}(x)H_{3}(x)H_{4}(x)e^{-x^2}dx=4224\sqrt{\pi},$$ $$\int_{-\infty}^\infty H_{1}(x)H_{2}(x)H_{3}(x)H_{5}(x)e^{-x^2}dx=19776.$$
$H_1(x)=2x$ is the simplest case and
$$\frac{(2n)!}{n!}=\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty H_1(x)^{2n}e^{-x^2}dx,$$
which corresponds to this OEIS sequence. The only permutation related thing on that sequence is a "downgrade permutation" which is a permutation that sends another permutation to a descending sequence. Then this integral counts the number of self-downgrade permutations. This seems semi-promising as it's some kind of order-derangement condition.
I didn't have much luck with Jacobi polynomials, which gave rational numbers.
So here's a general question. Has this kind of phenomena, as in the case of Laguerre polynomials, been studied for other families? In other words, for a family of orthogonal polynomials $P_n(x)$, with orthogonality weight $w(x)$ on a space $A$, is there a combinatorial interpretation for:
$$\int_A P_{n_1}(x)\cdots P_{n_r}(x)w(x)dx?$$
Are there any references on this? As far as I'm aware, Askey was able to reprove Evan's and Gillis result using analytic techniques (Evan's proof was actually combinatorial).
Naively, all orthogonal polynomials fall into some category of the Askey Scheme, meaning that really these identities are somehow connected to the right set of hypergeometric functions and their weights. Unfortunately, that's not exactly illuminating.
