n-point correlations of Gaussian random variables can be simplified with Wick expansion.
$$ \langle x_{i_1} x_{i_2} \dots x_{i_{2n-1}} x_{i_{2n}} \rangle = \int_{\mathbb{R}^n} x_{i_1} \dots x_{i_{2n}} e^{- \frac{1}{2}\sum x_i^2} = \sum_{\sigma \in \text{matchings}} \; \prod_{\{i,j\}\in \sigma} \langle x_i x_j \rangle$$
A *matching* is a partition of $\{ 1,2,\dots ,2n\}$ into two-element sets, e.g. $\{ \{1,2\},\dots, \{2n-1,2n\}\}$.

The Chebyshev polynomials are orthogonal with respect to $\mu(dx) = \sqrt{4-x^2} dx$ but there probably isn't a way to simplify n-point correlations, \[ \int_{S^n} x_{i_1}\dots x_{i_{2n}} \sqrt{1 - (x_1^2+\dots + x_{n}^2)} \; dx \]

Is there an analogue of Wick's formula for other orthogonal polynomials? It seems rather unlikely since the measure can be arbitrary.

Maybe look at Wick expansion as the sum over matchings. Is there a version of the Wick expansion for non-crossing matchings or for Schroder paths?