How prove this integral limit $=f(\frac{1}{2})$

We can take advantage of a probabilistic interpretation. Let $X_1,\ldots,X_n$ be uniformly distributed random variables on $[0,1]$ and $\bar X=\frac{1}{n}(X_1+\cdots+X_n)$. Then $$\lim\limits_{n\to \infty}\int_0^1\cdots\int_0^1f\left(\frac{x_1+\cdots+x_n}{n}\right)dx_1\cdots dx_n =\lim\limits_{n\to\infty}E[f(\bar X)]=E\left[f\left(\frac12\right)\right]=f\left(\frac12\right)$$ since $\bar X$ converges in distribution to $\frac12$ as $n\to \infty$.


Starting with OP's hint,

$\displaystyle \lim\limits_{n\to\infty}\int_{0}^{1}\cdots\int_{0}^{1} \left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right)dx_{1}dx_{2}\cdots dx_{n}=\dfrac{1}{2}$,

Now, we try to show, $\displaystyle \lim\limits_{n\to\infty}\int_{0}^{1}\cdots\int_{0}^{1} \left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right)^k dx_{1}dx_{2}\cdots dx_{n}=\dfrac{1}{2^k}$

If we count the number of terms in the multinomial expansion of $\bigg(\sum\limits_{i=1}^n x_i\bigg)^k$, that contains all variables with power not exceeding $1$, is $n(n-1)(n-2)\cdots(n-k+1) = n^k + O(n^{k-1})$,

and the number of terms that has atleast one $x_i$ term with powers exceeding $1$ is not more than $n.n^{k-2} = n^{k-1}$.

So, combining all the terms we get, $\displaystyle \int_{0}^{1}\cdots\int_{0}^{1} \left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right)^k dx_{1}dx_{2}\cdots dx_{n}=\dfrac{1}{2^k} + O\bigg(\frac{1}{n}\bigg)$

Thus the above limit is true for polynomials.

Using Weierstrass Approximation theorem, we can choose a polynomial $P$, such that $|f(x) - P(x)| < \epsilon/3$ for all $x \in [0,1]$.

Thus, we can find an $N \in \mathbb{N}$, such that $\forall n >N$,

$\displaystyle \left| \int_{0}^{1}\cdots\int_{0}^{1} P\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) dx_{1}dx_{2}\cdots dx_{n} - P\left(\frac12\right)\right| < \epsilon/3$

Therefore, for $n > N$,

$\displaystyle \left| \int_{0}^{1}\cdots\int_{0}^{1} f\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) dx_{1}dx_{2}\cdots dx_{n} - f\left(\frac12\right)\right| < \left| \int_{0}^{1}\cdots\int_{0}^{1} f\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) - P\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) dx_{1}dx_{2}\cdots dx_{n}\right| + \left| \int_{0}^{1}\cdots\int_{0}^{1} P\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) dx_{1}dx_{2}\cdots dx_{n} - P\left(\frac12\right)\right| + \left|f\left(\frac12\right) - P\left(\frac12\right)\right| $

$\le \displaystyle \int_{0}^{1}\cdots\int_{0}^{1} \left| f\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) - P\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right)\right| dx_{1}dx_{2}\cdots dx_{n} + 2\epsilon/3 < \epsilon$

Thus, $\displaystyle \lim\limits_{n\to\infty} \int_{0}^{1}\cdots\cdots\int_{0}^{1}f\left(\dfrac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right)dx_{1}dx_{2}\cdots dx_{n}=f\left(\dfrac{1}{2}\right)$.


$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ $\ds{\bbox[5px,#ffd]{\lim_{n \to \infty} \int_{0}^{1}\!\!\!\!\cdots\!\!\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n} = \fermi\pars{\half}}:\ {\Large ?}}$


\begin{align} &\bbox[#ffd,5px]{% \lim_{n \to \infty}\int_{0}^{1}\cdots\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n}} \\[5mm] = &\ \lim_{n \to \infty}\ \int_{0}^{1}\cdots\int_{0}^{1} \int_{-\infty}^{\infty} \\[2mm] &\ \phantom{\lim_{n \to \infty}\,\,\,} \tilde{\fermi}\pars{k} \exp\pars{\ic k\,{x_{1} + \cdots + x_{n} \over n}} \,{\dd k \over 2\pi}\dd x_{1}\ldots\dd x_{n} \end{align} where $\ds{% \tilde{\fermi}\pars{k} \equiv \int_{-\infty}^{\infty}\fermi\pars{x} \expo{-\ic k x}\,\dd x}$ is the $\ds{\fermi\pars{x}}$ Fourier Transform.

\begin{align} &\bbox[#ffd,5px]{% \lim_{n \to \infty} \int_{0}^{1}\cdots\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n}} \\[5mm] = & \lim_{n \to \infty}\int_{-\infty}^{\infty}\tilde{\fermi}\pars{k} \pars{\int_{0}^{1}\expo{\ic kx/n}\,\dd x}^{n}\,{\dd k \over 2\pi} \\[5mm] = &\ \lim_{n \to \infty}\int_{-\infty}^{\infty}\tilde{\fermi}\pars{k} \pars{\expo{\ic k/n} - 1 \over \ic k/n}^{n}\,{\dd k \over 2\pi} \\[5mm] = &\ \lim_{n \to \infty}\int_{-\infty}^{\infty} \tilde{\fermi}\pars{k} \exp\pars{\ic k \over 2} \braces{\sin\pars{k/\bracks{2n}} \over k/\bracks{2n}}^{n}\,{\dd k \over 2\pi} \end{align}


\begin{align} &\bbox[#ffd,5px]{\lim_{n \to \infty} \int_{0}^{1}\cdots\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n}} \\[2mm] = &\ \int_{-\infty}^{\infty} \fermi\pars{x}\lim_{n \to \infty} \operatorname{K}_{n}\pars{x - \half}\,\dd x \end{align} where $$ \operatorname{K}_{n}\pars{x} \equiv \int_{-\infty}^{\infty} \exp\pars{-\ic k x} \braces{\sin\pars{k/\bracks{2n}} \over k/\bracks{2n}}^{n} \,{\dd k \over 2\pi} $$ Since $\ds{\lim_{n \to \infty}{\rm K}_{n}\pars{x} = \delta\pars{x}}$, we'll have: \begin{align} &\bbox[#ffd,5px]{\lim_{n \to \infty} \int_{0}^{1}\cdots\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n}} \\[5mm] = &\ \int_{-\infty}^{\infty} \fermi\pars{x}\delta\pars{x - \half}\,\dd x \end{align}

$\ds{\delta}$ is the Dirac Delta Function.


Finally, \begin{align} &\bbox[#ffd,5px]{\lim_{n \to \infty} \int_{0}^{1}\cdots\int_{0}^{1} \fermi\pars{x_{1} + \cdots + x_{n} \over n} \,\dd x_{1}\ldots\dd x_{n}} \\[5mm] = &\ \bbox[5px,border:1px groove navy]{\fermi\pars{\half}} \\ & \end{align}