How to integrate the multinomial over a ball in $\mathbb{R}^{n}$? I got an interesting question. Consider this integral:
$$ \int_{B(0,1)}\bigg(\sum_{j=1}^{n}a_{j}x_{j}^2\bigg)^m \mbox{d}x, \quad m,n\in \mathbb{N}, \  a_{i}>0, \ i=1,2,\ldots,n.$$
It is clear that if we choose all $a_{j}=1$, then it will be simple! Futhermore, I want to know if the answer is connected with some special functions like hyperbolic functions, please do not use a finite summation in the result.
If we just expand this multinomial directly, that is,
$$ \bigg(\sum_{j=1}^{n}a_{j}x_{j}^2\bigg)^m =\sum_{\lambda_{1}+\lambda_{2}+\ldots+\lambda_{n}=m}\frac{m!}{\lambda_{1}!\lambda_{2}\ldots\lambda_{n}!}a_{1}^{\lambda_{1}}a_{2}^{\lambda_{2}}\ldots a_{n}^{\lambda_{n}}x_{1}^{2\lambda_{1}}x_{2}^{2\lambda_{2}}\ldots x_{n}^{2\lambda_{n}},
$$
then, after a simple calculation, in Table of Integrals I found a related formula,
$$\frac{\bigg(\sum_{i=1}^{r}a_{i}^{2}\bigg)^{\frac{n}{2}}}{n!}H_{n}\left(\frac{\sum_{i=1}^{r}a_{i}x_{i}}{\sqrt{\sum_{i=1}^{r}a_{i}^{2}}}\right)=\sum_{\substack{m_{1}+\ldots+m_{r}=n,\\{m_{i}\geq 0}}}\prod_{k=1}^{r}\bigg\{\frac{a_{k}^{m_{k}}}{m_{k}!}H_{m_{k}}(x_{k})\bigg\},
$$
where $H_{n}(x)$ denotes the Hermite polynomials; I don't know if it works. If this result is a little hard, one can just find the asymptotic representation as $m\leq n$ and $m\rightarrow \infty$.
 A: $\newcommand\R{\mathbb R}\newcommand{\Ga}{\Gamma}$The integral in question is
\begin{equation*}
\begin{aligned}
    \int_0^1 dr\,\int_{rS_{n-1}}dx\,f(x)
    &=\int_0^1 dr\,r^{2m+n-1}\int_{S_{n-1}}du\,f(u) \\ 
    &=\frac1{2m+n}\,|S_{n-1}|Ef(U_n), 
\end{aligned}
\end{equation*}
where
\begin{equation*}
    f(x):=\Big(\sum_{j=1}^n a_j x_j^2\Big)^m 
\end{equation*}
(so that $f$ is $2m$-homogeneous), $|S_{n-1}|$ is the surface area of the unit sphere $S_{n-1}$ in $\R^n$, and $U_n$ is a random vector uniformly distributed on $S_{n-1}$.
So, the problem reduces to finding the asymptotics of $Ef(U_n)$. Without loss of generality,
\begin{equation*}
    U_n=\frac G{\|G\|},
\end{equation*}
where $G=(G_1,\dots,G_n)$ is a standard Gaussian vector in $\R^n$ and $\|\cdot\|$ is the Euclidean norm.
Let us normalize the $a_j$'s by assuming that
\begin{equation*}
    \sum_{j=1}^n a_j=1 \tag{0}
\end{equation*}
and consider the following two extreme cases:
Case 1: $a_1=1$, $a_2=\cdots=a_n=0$ and
Case 2: $a_1=\cdots=a_n=1/n$.
In Case 1,
\begin{equation*}
    f(U_n)=\Big(\frac{G_1^2}{\|G\|^2}\Big)^m, 
\end{equation*}
and $\dfrac{G_1^2}{\|G\|^2}$ has the beta distribution with parameters $1/2,(n-1)/2$. So, here
\begin{equation*}
    Ef(U_n)=\frac{\Ga(n/2)\Ga(m+1/2)}{\Ga(1/2)\Ga(m+n/2)}. 
\end{equation*}
So, if e.g. $m=n\to\infty$, then the integral in question is
\begin{equation*}
    \sim\frac1{3n}\,\sqrt6\,\Big(\frac2{3\sqrt3}\Big)^n\,|S_{n-1}|. \tag{1}
\end{equation*}
So, here the coefficient of $|S_{n-1}|$ decreases exponentially.
In the trivial Case 2, $f(U_n)=1/n^m$, and hence for $m=n\to\infty$ the integral in question is
\begin{equation*}
    =\frac1{3n^{n+1}}\,|S_{n-1}|, 
\end{equation*}
which is very different from the asymptotics for Case 1 given in (1).
Between these two extreme cases, any intermediate asymptotics (or absense of any asymptotics) should be possible. As was noted in comments, the asymptotics in general (if any) will very much depend on how much the $a_j$'s differ from one another.
Addendum 1:
To be more specific, consider the following setting, intermediate between Cases 1 and 2:

$a_1=\cdots=a_k=1/k$, $a_{k+1}=\cdots=a_n=0$ for natural $k\in[2,n-1]$, where $k$ is allowed to vary with $n$ and $m$.

Then
\begin{equation*}
    f(U_n)=\frac1{k^m}\,\Big(\frac{G_1^2+\cdots+G_k^2}{\|G\|^2}\Big)^m, 
\end{equation*}
and $\dfrac{G_1^2+\cdots+G_k^2}{\|G\|^2}$ has the beta distribution with parameters $k/2,(n-k)/2$. So, here
\begin{equation*}
    Ef(U_n)=\frac1{k^m}\,\prod_{r=0}^{m-1}\frac{k/2+r}{n/2+r}. \tag{2}
\end{equation*}
So, if $k/m^2\to\infty$, then
\begin{equation*}
    Ef(U_n)\sim\frac1{n^m}, \tag{3}
\end{equation*}
which is the same asymptotics as for $m$ fixed.
However, if $n/m^2\to\infty$ and $k\sim cm^2$ for some $c\in(0,\infty)$, then
\begin{equation*}
    Ef(U_n)\sim\frac{e^{1/c}}{n^m}, 
\end{equation*}
now with the additional factor $e^{1/c}$ as compared to (3).
If we now let $n/m^2\to\infty$ and $k\sim cm$ for some $c\in(0,\infty)$, then an additional exponentially growing (with $m$) factor (as compared to (3)) will appear.
Addendum 2: Letting $g(a_1,\dots,a_n):=Ef(U_n)$, we see that the function $g$ is symmetric and convex, and hence Schur convex -- see e.g. Theorem A on p. 258. So, in view of condition (0),
\begin{equation*}
    Ef(U_n)\ge g(1/n,\dots,1/n)=\frac1{n^m}. \tag{$\clubsuit$}
\end{equation*}
That is, the smallest value of $Ef(U_n)$ is attained in Case 2 ($a_1=\cdots=a_n=1/n$), considered above.
Addendum 3: For completeness, consider also the case when $m$ is fixed (even though $m\to\infty$ in the OP).
Then, by (2),
\begin{equation*}
    Ef(U_n)\sim\frac1{n^m}\,\prod_{r=0}^{m-1}\Big(1+\frac{2r}k\Big) 
\end{equation*}
if $a_1=\cdots=a_k=1/k$ and $a_{k+1}=\cdots=a_n=0$ for a fixed natural $k\in[2,n-1]$ and $n\to\infty$. We see that, even when $m$ is fixed, the asymptotics depends on $k$ and, more generally, on how much the $a_j$'s differ from one another.
Addendum 4: Here we complement the lower bound on $Ef(U_n)$ given by ($\clubsuit$) in Addendum 2 by providing a matching upper bound on $Ef(U_n)$ that implies the following:

If the $a_i$'s are uniformly small in the sense that
\begin{equation*}
    a:=\max_{i=1}^n a_i\to0 
\end{equation*}
and, moreover, $m$ is at most moderately large in the sense that
\begin{equation*}
    m^3 a\to0,  \tag{4}
\end{equation*}
then
\begin{equation*}
    Ef(U_n)\sim\frac1{n^m}  \tag{$\heartsuit$}
\end{equation*}
(as $n\to\infty$).

Indeed, the random point $(x_1^2,\dots,x_n^2)$ has the Dirichlet distribution with parameters $1/n,\dots,1/n$, and the Dirichlet distribution has the negative association (NA) property. So, by (say) Theorem 2,
\begin{equation*}
    Ef(U_n)\le E\Big(\sum_{j=1}^n a_j Y_j\Big)^m,  \tag{5}
\end{equation*}
where the $Y_j$'s are iid random variables each with the beta distribution with parameters $1/2,(n-1)/2$. Denoting now the $L^m$ norm by $\|\cdot\|_m$ and using Minkowski's inequality, we get
\begin{equation*}
    (Ef(U_n))^{1/m}\le\Big\|\sum_{j=1}^n a_j Y_j\Big\|_m
    \le \Big\|\sum_{j=1}^n a_j EY_j\Big\|_m +\Big\|\sum_{j=1}^n a_j Z_j\Big\|_m,
\end{equation*}
where $Z_j:=Y_j-EY_j$. Since $EY_j=1/n$, we have
\begin{equation*}
    \Big\|\sum_{j=1}^n a_j EY_j\Big\|_m
=   \sum_{j=1}^n a_j EY_j=\frac1n. \tag{6}
\end{equation*}
Note also that $Var\,Y_j\sim2/n^2$ and $\|Z_j\|_m\le2\|Y_j\|_m\ll m/n$; here and in what follows, $A\ll B$ means $A\le CB$ for some universal real constant $C>0$.
Using now an appropriate version of Rosenthal's inequality (see e.g. Theorem 6.1), we get
\begin{equation*}
\begin{aligned}
\Big\|\sum_{j=1}^n a_j Z_j\Big\|_m
&\ll\frac1n\,(m^2 a^{1-1/m}+m^{1/2} a^{1/2}) \\ 
&=\frac1n\,\frac1m\,((m^3 a)^{1-1/m}m^{3/m}+(m^3 a)^{1/2})
=o\Big(\frac1n\,\frac1m\Big),
\end{aligned}
\end{equation*}
by (4).
So, by (5) and (6),
\begin{equation*}
    (Ef(U_n))\le\frac1{n^m}\Big(1+\frac{o(1)}m\Big)^m
    =\frac{1+o(1)}{n^m}. 
\end{equation*}
Now ($\heartsuit$) follows, in view of ($\clubsuit$).
A: Let me separate the radial integration from the angular integration,
$$\int_{|\mathbf{x}|\leq 1}f(\mathbf{x})d\mathbf{x}=\frac{2\pi^{n/2}}{\Gamma(n/2)}\int_0^1 r^{n-1}\bar{f}(r)\,dr,$$
where $\bar{f}(r)$ is the average of $f$ over the surface of the $n$-dimensional hypersphere of radius $r$. In our case
$$f(\mathbf{r})=\bigg(\sum_{j=1}^{n}a_{j}x_{j}^2\bigg)^m.$$
For $n\gg 1$ at fixed $m$ the concentration of measure allows us to replace $x_j^2$ by $r^2/n$, so
$$\bar{f}(r)\approx (r^2/n)^m\bigg(\sum_{j=1}^{n}a_{j}\bigg)^m,$$
and thus we estimate
$$\int_{B(0,1)}\bigg(\sum_{j=1}^{n}a_{j}x_{j}^2\bigg)^m \,d\mathbf{x}\approx
\frac{2\pi^{n/2}}{n^m(2m+n)\Gamma(n/2)}\bigg(\sum_{j=1}^{n}a_{j}\bigg)^m,\;\;n\gg 1.$$

Obviously, this estimate is exact if all $a_j$'s are equal to each other. To check for the other extreme, let me try (following the suggestion by Iosif Pinelis) the simple case $m=1$, $a_1=1$ and $a_j=0$ for $j=2,3,\ldots n$. The exact integration gives
$$\int_{B(0,1)}x_1^2\, d\mathbf{x}=\frac{2 \pi ^{\frac{n-1}{2}}}{\Gamma \left(\frac{n-1}{2}\right)}\int_0^1 r^{n+1}dr\int_0^\pi\cos^2\phi\sin^{n-2}\phi\, d\phi=  \frac{2\pi^{n/2}}{n(n+2)\Gamma(n/2)},$$
which again equals the large-$n$ estimate.
