Are the “generalized Catalan numbers” of Dumitrescu–Mulase the "moments" of some "multivariate Wigner semicircle distribution"? The classical Catalan numbers
$$ C_n = \frac{1}{n+1} \binom{2n}{n}, $$
well-known for their numerous combinatorial interpretations (the second volume of Stanley's Enumerative Combinatorics famously lists a total of 66), essentially arise as the even moments $E[X^{2n}] = \int x^{2n} d\mu$ of the Wigner semicircle distribution $d\mu$: this distribution depends on a parameter $R$, the radius of the semicircle, and in general one has $E[X^{2n}] = (R/2)^{2n} C_n$, so choosing $R=2$ gives precisely the Catalan numbers. (The odd moments are all zero.)
Now, associated to this distribution are a family of very classical orthogonal polynomials: the Chebyshev polynomials of the second kind, $U_n(x)$. (Up to a simple change of variables, $U_n(x)$ coincides with the matching polynomial of the path graph on $n$ vertices; for other interpretations and discussion see this blog post and this MO question.)
In their paper Lectures on the topological recursion for Higgs bundles and quantum curves, O. Dumitrescu and M. Mulase define a sequence of numbers $C_{g,n}(\mu_1,\ldots,\mu_n)$ which they call generalized Catalan numbers. The definition (found on p. 26) is as follows:
Definition: Let $\vec{\Gamma}_{g,n}(\mu_1,\ldots,\mu_n)$ denote the set of arrowed cell graphs drawn on a closed, connected, oriented surface of genus $g$ with $n$ labeled vertices of degrees $\mu_1,\ldots,\mu_n$. Then
$$ C_{g,n}(\mu_1,\ldots,\mu_n) := |\vec{\Gamma}_{g,n}(\mu_1,\ldots,\mu_n)|. $$
Below this, the authors justify the terminology by pointing out that $C_{0,1}(2m) = C_m = \frac{1}{m+1} \binom{2m}{m}$.
I don't know much about the story of multivariate orthogonal polynomials, and there are quite a few parameters here, but I can't help but wonder: might there perhaps be some multivariate analogue of the Wigner semicircle distribution, having these quantities $C_{g,n}(\mu_1,\ldots,\mu_n)$ as its (appropriately defined) moments – and with it, a corresponding family of orthogonal polynomials?
I also find it curious that the above paper discusses Laplace transforms, since the (bilateral) Laplace transform is exactly how one passes from a probability density function to its moment-generating function (as one learns in a first course on probability theory.) Admittedly, this is probably just a coincidence.
 A: I'm new here and not sure if this constitutes an answer, but I saw this was open for a while and figure I can shed some light (please remove, edit, etc. if needed).
You say that the $k$-th Catalan number is equal to $\int x^{2k}d\mu$ where $\mu$ is the Wigner SC distribution. Let me rephrase this in terms of the density $\rho(x)$ so that $C_k=\int x^{2k}\rho(x)dx$ (i.e., $\rho(x)$ is just the Radon-Nikodym derivative of $\mu$). The underlying story here is that we have an $n\times n$ matrix, belonging to the Gaussian unitary ensemble, whose eigenvalues have density $\rho_1(x)$. Moreover, the full joint p.d.f. for these eigenvalues is actually given by the $n$-point correlation function $\rho_n(x_1,\ldots,x_n)$; note that the $m$-point correlation function is defined by
$$
\rho_m(x_1,\ldots,x_m):=\frac{(n-1)!}{(n-m)!}\int\rho_n(x_1,\ldots,x_n)dx_{m+1}\cdots dx_n.
$$
If you haven't encountered this before, an example to keep in mind is that $\int_a^b\int_c^d\rho_2(x,y)dxdy$ is the joint probability that there is an eigenvalue in $(a,b)$ while there is another eigenvalue in $(c,d)$.
Notice that I refer to the Wigner SC law as $\rho(x)$ and the spectral density of the $n\times n$ GUE matrix as $\rho_1(x)$. The distinction here is that
$$
\rho(x):=\lim_{n\to\infty}\rho_1(x),
$$
which is an important part of the story (this is where the genus parameter $g$ comes into play). You mention that the Wigner SC $\rho(x)$ is related to the Chebyshev orthogonal polynomials of the second kind; one should instead consider the Hermite orthogonal polynomials, which relate to $\rho_1(x)$ in the same way. Then, it is possible to express the $m$-point correlation functions as determinantal processes with kernel given in terms of the Hermite orthogonal polynomials. You'll want to look up chapter 5 of Mehta's book, 'Random Matrices' and maybe the Christoffel-Darboux formula. I hope this partially answers your question on orthogonal polynomials (there are of course studies on multivariate orthogonal polynomials, but I'm not sure how relevant they are to this question).
Moving on, note that it is important to distinguish the eigenvalue density $\rho_1(x)$ from its large $n$ limit $\rho(x)$. Now introduce the smoothed density,
$$
\tilde{\rho}_1(x):=\rho(x)+\sum_{g=1}^{\infty}\frac{\rho_{1,g}(x)}{n^{2g}},
$$
which is defined such that its moments coincide with those of $\rho_1(x)$. However, the formal series defining $\tilde{\rho}_1(x)$ and the true density $\rho_1(x)$ differ by oscillatory contributions which aren't reflected in the moments (thus, the two densities agree only in the large $n$ limit). The point here is that the generalised Catalan numbers of Dumitrescu and Mulase are such that
$$
C_{g,1}(\mu)=\int x^\mu\rho_{1,g}(x)dx,
$$
with $\rho_{1,0}(x):=\rho(x)$ given by the Wigner SC law. This can be expressed in a multivariate way by defining
$$
\tilde{\rho}_n(x_1,\ldots,x_n):=\sum_{g=0}^{\infty}\frac{\rho_{n,g}(x_1,\ldots,x_n)}{n^{n-1+2g}}
$$
so that the moments of $\tilde{\rho}_n$ and the cumulants of $\rho_n$ agree (the moment-cumulant relation is why $\tilde{\rho}_n$ is O$(n^{1-n}$). Then we have that
$$
C_{g,1}(\mu)=\int\left(\frac{1}{n}\sum_{i=1}^nx_i^\mu\right)\rho_{n,g}(x_1,\ldots,x_n)dx_1\cdots dx_n
$$
(use symmetry arguments to compare with the earlier equation if you doubt this), and more generally
$$
C_{g,m}(\mu_1,\ldots,\mu_m)=\int\prod_{i=1}^m\left(\frac{1}{n}\sum_{j=1}^nx_j^{\mu_i}\right)\rho_{n,g}(x_1,\ldots,x_n)dx_1\cdots dx_n.
$$
I hope this is helpful, though I warn that I may have some proportionality constants missing, and I've left out many details as this is quite a vast topic. I'll end with pointing out that the presence of the Laplace transform is not coincidental. Transitioning between p.d.f.s and moment generating functions is key to making things rigorous.
