Find $\lim\limits_{n\to\infty}\left(\frac{(2n+1)!}{(n!)^2}\right)^2\int_0^1 \int_0^1 (xy(1-x)(1-y))^n f(x,y)dxdy$ 
Note: even though this is technically a duplicate of this post, I'd like further justification on what exactly the claim about "convergence in distribution" in @JackD'Aurizio's answer means, and, if possible, a proof as to why that claim holds.


Let $f\colon [0,1]^2\to\mathbb{R}$ be a continuous function. Find $$\lim\limits_{n\to\infty} \left(\dfrac{(2n+1)!}{(n!)^2}\right)^2\int_0^1 \int_0^1 (xy(1-x)(1-y))^n f(x,y)\ \mathrm dx\mathrm dy.$$

I think it might be useful to first consider the case where $f$ is a polynomial. By the linearity of integrals, it suffices to consider the case where $f(x,y)=x^ky^l$ for all x,y. Also, Stirling's formula might be useful. There's probably a general formula for evaluating the double integral $I(n,k,l) := \int_0^1 \int_0^1 x^{n+k} y^{n+l} (1-x)^n (1-y)^n dxdy.$ It seems possible but tedious to evaluate the latter integral using the Binomial theorem. Alternatively it might be possible to use induction if one can guess the general formula for the integral. For instance in the specific case that $n=0,$ we have the integral $\int_0^1 \int_0^1 x^k y^l dxdy = \dfrac{1}{(k+1)(l+1)}$. By the linearity of the integral we also have $I(n,k,l) = I(n-1,k+1,l+1) -I(n-1,k+2,l+1)-I(n-1,k+1,l+2)+I(n-1,k+2,l+2).$  Also a special case of the Stone-Weierstrass theorem says that every continuous function $f:[0,1]^2\to\mathbb{R}$ can be uniformly approximated by polynomials. But I'm not sure if one can reduce to the case of polynomials.
 A: Another way is to substitute $x=(1+u/\sqrt{n})/2,\ y=(1+v/\sqrt{n})/2$ and apply DCT.
We get $\lim\limits_{n\to\infty}A_n\iint_{[-\sqrt{n},\sqrt{n}]^2}f_n(u,v)\,du\,dv$, where $$A_n=\frac{(2n+1)!^2}{2^{4n+2}n!^4 n}\underset{n\to\infty}{\longrightarrow}\frac1\pi$$ by Stirling's asymptotics, and $$f_n(u,v)=\left(1-\frac{u^2}{n}\right)^n\left(1-\frac{v^2}{n}\right)^n f\left(\frac12+\frac{u}{2\sqrt{n}},\frac12+\frac{v}{2\sqrt{n}}\right)$$ satisfies the premises of DCT: $$\lim_{n\to\infty}f_n(u,v)=f\left(\frac12,\frac12\right)e^{-u^2-v^2},\quad\big|f_n(u,v)\big|\leqslant e^{-u^2-v^2}\sup_{0\leqslant x,y\leqslant 1}\big|f(x,y)\big|.$$ Hence, by the theorem, the limit equals $\frac1\pi f\left(\frac12,\frac12\right)\iint_{[-\infty,\infty]^2}e^{-u^2-v^2}\,du\,dv=f\left(\frac12,\frac12\right)$.
A: As requested, here is a justification of Jack D'Aurizio's claim about convergence in distribution in the duplicate discovered by TheBestMagician (thanks to him!).
First note that by Stirling's formula,
$$C_n:=\frac1{\mathrm B(n+1,n+1)}=\frac{(2n+1)!}{(n!)^2}\sim\frac{2^{2n+1}\sqrt n}{\sqrt\pi}.$$
This is why the  probability measure on $[0,1]$ with density
$$\varphi_n(t):=C_n\,t^n(1-t)^n$$
converges in distribution to the Dirac measure $\delta_{\frac12},$ namely:
$$\forall x\in\left[0,\frac12\right)\quad\int_0^x\varphi_n(t)\,\mathrm dt=\int_{1-x}^1\varphi_n(s)\,\mathrm ds\to0$$
because
$$\forall t\in[0,x]\quad\varphi_n(t)\le C_nx^n(1-x)^n\sim\frac{2\sqrt n}{\sqrt\pi}q^n,\quad\text{where}\quad q:=4x(1-x)<1.$$
As a consequence, the probability measure on $[0,1]^2$ with density $\varphi_n(u)\varphi_n(v)$ converges in distribution to $\delta_{\left(\frac12,\frac12\right)},$ namely:
$$\iint_A\varphi_n(u)\varphi_n(v)\,\mathrm du\mathrm dv\to\delta_{\left(\frac12,\frac12\right)}(A)$$
for every Borel set $A\subset[0,1]^2$ not containing $\left(\frac12,\frac12\right).$
Since convergence in distribution implies weak convergence , we conclude that for any continuous (hence bounded) function $f:[0,1]^2\to\mathbb R,$
$$\iint f(u,v)\varphi_n(u)\varphi_n(v)\,\mathrm du\mathrm dv\to f\left(\frac12,\frac12\right).$$
A: The expression
\begin{align*}
\left(\frac{(2n+1)!}{(n!)^2}\right)^2(xy(1-x)(1-y))^n= \frac{x^n(1-x)^n\cdot y^n (1-y)^n}{B(n+1, n+1)^2} = f_{n+1, n+1}(x) f_{n+1, n+1}(y)
\end{align*}
where $f_{\alpha, \beta}(x)$ represents the density of a Beta distribution with shape parameters $\alpha, \beta$, and $B(\cdot, \cdot)$ is the Beta function. Changing to a probabilistic interpretation, we can represent the limit as
\begin{align*}
\lim_{n\rightarrow \infty} \mathbb{E}\left[f\begin{pmatrix} X_n \\ Y_n \end{pmatrix}\right]
\end{align*}
where $X_n, Y_n$ both independently follow $\text{Beta}(n+1, n+1)$. By representation, we have $X_n \overset{\mathcal{D}}{=} \frac{\frac{1}{n+1}\sum_{k=1}^{n+1} A_k}{\frac{1}{n+1}\sum_{k=1}^{n+1} (A_k + B_k)}$ and $Y_n \overset{\mathcal{D}}{=} \frac{\frac{1}{n+1}\sum_{k=1}^{n+1} C_k}{\frac{1}{n+1}\sum_{k=1}^{n+1} (C_k + D_k)}$ for $A_k, B_k, C_k, D_k$ all independently following $\text{Exp}(1)$. By Law of Large Numbers and Sluksky's, $X_n, Y_n \overset{\mathcal{P}}{\rightarrow} \frac{1}{2}$.
Since $f$ is continuous on a closed, bounded interval, $f$ is also bounded. We have
\begin{align*}
\lim_{n\rightarrow \infty} \mathbb{E}\left[f\begin{pmatrix} X_n \\ Y_n \end{pmatrix}\right] &=  \mathbb{E}\left[\lim_{n\rightarrow \infty}f\begin{pmatrix} X_n \\ Y_n \end{pmatrix}\right] && \text{(Dominated Convergence)} \\
&=  \mathbb{E}\left[f\begin{pmatrix} \lim_{n\rightarrow \infty} X_n \\ \lim_{n\rightarrow \infty} Y_n \end{pmatrix}\right] && \text{(Continuity of $f$)} \\
&=  \mathbb{E}\left[f\begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \end{pmatrix}\right] && \text{(LLN and Slutsky's)}\\
&=  f\begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \end{pmatrix}
\end{align*}
as desired.
A: Here is a probabilistic short proof of the statement. Let $(X_n:n\in\mathbb{N})$ an i.i.d. sequence of exponential distributions with parameter $\theta$ (the value of $\theta$ is not relevant as we will see, and can be taken as $1$ for example). The following are basic facts in probability:

*

*For any $m$ and $k$, $X_1+\ldots + X_k$ and $X_{k+1}+\ldots + X_{k+m}$ are independent random variables with $\operatorname{Gamma}(k,\theta)$ and $\operatorname{Gamma}(m,\theta)$ distributions respectively (easily seen by taking Fourier transform of the convolution of exponentials).

*If $U$ and $V$ are idependent random variables with distributions $\operatorname{Gamma}(\alpha,\theta)$ and $\operatorname{Gamma}(\beta,\theta)$, then
$\frac{U}{U+V}$ is a $\operatorname{Beta}(\alpha,\beta)$ distributed random variable (an application of change of variables for an intergral over $(0,\infty)\times(0,\infty)$).

It follows that
$$Y_n=\frac{X_1+\ldots + X_n}{(X_1+\ldots + X_n)+(X_{n+1}+\ldots + X_{2n})}$$
is a $\operatorname{Beta}(n,n)$ random variable. By the law of large numbers,
$$\lim_{n\rightarrow\infty}Y_n=\lim_{n\rightarrow\infty}\frac{\tfrac1n(X_1+\ldots+ X_n)}{\tfrac1n(X_1+\ldots+X_n)+\tfrac1n(X_{n+1}+\ldots+X_{2n})}=\frac12$$
almost surely. Consequently, $Y_n\Longrightarrow\delta_{1/2}$ in law (simple application of dominated convergence).
To conclude, consider two independent sequences $(X_n)$ and $(X'_n)$ of i.i.d exponential random variable with parameter $1$, and define $Y_n$ and $Y'_n$ as above. Then $Y_n\stackrel{d}{=} Y'_n$
The joint distribution of $(Y_n,Y'_n)$ is $f_n(u,v)=\Big(\tfrac{\Gamma(2n)}{\Gamma(n)\Gamma(n)}\Big)^2u^{n-1}(1-u)^{n-1}v^n(1-v)^{n-1}\mathbb{1}_{[0,1]^2}(u,v)$. Thus, $$(Y_n,Y'_n)\stackrel{n\rightarrow\infty}{\Longrightarrow}\delta_{1/2,1/2}$$
