Central Limit Theorem for symmetric Beta Distribution 
Prove that for $X_n\sim Beta(n,n)$ $$2\sqrt{2n}\left( X_n-\frac{1}{2}\right)\Rightarrow\mathscr N(0,1)$$ in distribution as $n\rightarrow\infty$.

Let $Y_n$ and $Z_n$ be $Gamma(n,1)$ independent random variables. Then using the facts, that $$X_n=\frac{Y_n}{Y_n+Z_n},$$ and that both $Y_n$ and $Z_n$ are sum of independent $Exp(1)$ random variables, I can apply CLT to those, so both $Y_n$ and $Y_n+Z_n$ are normally distributed. Then I can use the delta method to show, that $X_n$ is normally distributed also.
That is the part where I get confused. Shall I apply the delta method to $Y_n$ with $g(x)=\frac{x}{x+Z_n}$?
 A: Maybe not the method you need/prefer, but if you define $T_n=2\sqrt{2n}(X_n-1/2)$, the pdf of $T_n$ is
$$
f_{T_n}(t)=\frac{1}{B(n,n)}\int_0^1 dx\ x^{n-1}(1-x)^{n-1}\delta(t-2\sqrt{2n}x+\sqrt{2n})\ ,
$$
where $B(\alpha,\beta)$ is a Beta function. Therefore
$$
f_{T_n}(t)=\frac{\mathbf{1}_{-\sqrt{2n}<t<\sqrt{2n}}}{2\sqrt{2n}B(n,n)}\left(\frac{t+\sqrt{2n}}{2\sqrt{2n}}\right)^{n-1}\left(1-\frac{t+\sqrt{2n}}{2\sqrt{2n}}\right)^{n-1}\to \frac{e^{-\frac{t^2}{2}}}{\sqrt{2 \pi }}\ ,
$$
as $n\to\infty$.
A: One can use the multivariate version of the delta method: 
Assume that $a_n \to \infty$ and $a_n(T_n - \mu) \stackrel{d}{\to} U$ where $T_n,\mu,U \in \mathbb R^d$. Then, under mild regularity assumptions on $f:\mathbb R^d \to \mathbb R$, we have 
\begin{align*}
 \sqrt{n} \big(f(T_n) - f(\mu)\big) \stackrel{d}{\to} J_\mu U
\end{align*}
where $J_\mu$ is the Jacobian  matrix of $f$ evaluated at $\mu$, that is $ [J_\mu]_{ij} = \partial f_i / \partial x_j \; \big|_{x = \mu}$.
In this example, let $\bar Y_n = Y_n / n$ and $\bar Z_n = Z_n/n$. Then, by the multivariate CLT (or arguing by independence from univariate CLT), we have 
\begin{align*}
 \sqrt{n} \Big[ 
  \begin{pmatrix}
  \bar Y_n \\ \bar Z_n
 \end{pmatrix} - 
 \begin{pmatrix}
  1 \\ 1
 \end{pmatrix} 
  \Big]
 \stackrel{d}{\to}  U := 
\begin{pmatrix}
  U_1 \\ U_2
 \end{pmatrix} \sim N(0,I_2).
\end{align*}
Let $f(y,z) = y/(y+z)$. The Jacobian can be identified with a row vector $$\Big[\frac{z}{(y+z)^2},\; \frac{-y}{(y+z)^2}\Big].$$ Since $X_n = f(\bar Y_n, \bar Z_n)$, we have by the multivariate delta method
\begin{align*}
  \sqrt{n} \big( X_n - f((1,1))\big) \stackrel{d}{\to} 
\begin{bmatrix}
1/4 & -1/4
 \end{bmatrix}
\begin{bmatrix}
  U_1 \\ U_2
 \end{bmatrix} = \frac{U_1-U_2}{4} \sim N(0, \frac18).
\end{align*}
Rescaling, we have $\sqrt{8n} (X_n - 1/2) \stackrel{d}{\to} N(0,1)$, the desired result.
