I have a function $I:\mathbb{N}\setminus\{0\}\to\mathbb{R}$ such that: $$ I(n) = \int_{[0,2\pi]^n}f_n(\theta_1,\dots,\theta_n) d\theta_1\dots d\theta_n $$ where $f_n:\mathbb{R}^{n}\to\mathbb{R}$ is a bounded real function of class $C^{\infty}(\mathbb{R}^n)$. In particular in my specific case: $$ f_n(\theta_1,\dots,\theta_n) = \Bigl(\sum_{i<j}^{n}\cos(\theta_i-\theta_j)\Bigr)^k $$ I would like to compute an asymptotic formula for $I(n)$ as $n\to\infty$. In my particular case, I should obtain something like:
$$ I(n) \sim !k\cdot (2\pi)^n\cdot \Bigl(\frac{n}{2}\Bigr)^k\mbox{ as $n\to\infty$} $$
where $!k$ is the subfactorial of $k$ and $k\in\mathbb{N}$,$k\geq 2$ is the parameter of the function $f_n$.
I was trying to expand the function $f_n$ using the "multidimensional" taylor series, but I was wondering about how many terms I should consider before truncating the series (or I have to consider all the infinite terms of the series to apply this idea?)
Are there any other "famous" methods to approach this computation? Searching online and in some analysis books I didn't find anything useful, unfortunately. Every suggestion/idea is appreciated! Thank you in advance.
Example
For $k = 2$ I found out that: $$ \int_{[0,2\pi]^n}\Bigl(\sum_{i<j}^{n}\cos(\theta_i-\theta_j)\Bigr)^2 d\theta_1\dots d\theta_n = \frac{n!}{(n-2)!} 2^{n-2} \pi^n = \frac{n(n-1)}{4}(2\pi)^n $$
and in fact, for $k = 2$: $$ \lim_{n\to\infty} \frac{I(n)}{!2\cdot(2\pi)^n\cdot\Bigl(\frac{n}{2}\Bigr)^2} = \lim_{n\to\infty} \frac{\frac{n(n-1)}{4}(2\pi)^n}{1\cdot(2\pi)^n\cdot\Bigl(\frac{n}{2}\Bigr)^2} = 1 $$
where we have used the fact that $!2 = 1$. I prove a similar result for $k=3$ but for $k>3$ it seems difficult to find a closed formula for $I(n)$.
UPDATE
I was trying to prove this asymptotic behavior using induction on $k$. In particular I prove that:
$$ I_{k+1}(n) = \int_{[0,2\pi]^n}\Bigl(\sum_{i<j}^{n}\cos(\theta_i-\theta_j)\Bigr)^{k+1} d\theta_1\dots d\theta_n = $$ $$ = \frac{n(n-1)}{2}\int_{[0,2\pi]^n}\cos(\theta_1-\theta_2)\Bigl(\sum_{i<j}^{n}\cos(\theta_i-\theta_j)\Bigr)^{k} d\theta_1\dots d\theta_n $$
but from this point I don't know how to proceed. In particular I don't know how to write the last integral in terms of $I_k(n)$ although it's very similar to it. Do you have any idea?
(Too long for a comment)
As a quick thought, I would advise you against expanding the multidimensional Taylor series. The reazon behind this is that you are trying to integrate your function $f(x)$ and thgough it holds that $f(x_1,\dots, x_n) = \lim_{j\to \infty} \sum_{i=0}^{j} p_i(x1,\dots,x_n)$ where each $p_i$ is a polynomial, it may not hold that $\int_U f(x_1,\dots, x_n)dx = \lim_{j\to \infty} \sum_{i=0}^{j} \int_U p_i(x1,\dots,x_n)dx$. This happens because interchanging a limit operation and an integral (as in $\int_U \lim_{j\to \infty} p_j(x1,\dots,x_n)dx = \lim_{j\to \infty} \int_U p_j(x1,\dots,x_n)dx$ ) is a nontrivial step which needs to be cautiously analyzed, and there are conditions which assure this operation can be done. So if you can't find any other way, you can try this, but make sure the limit operation can be extracted from the integral in your case.