If $X_n \sim \text{Beta}(n, n)$ Show that $[X_n - \text{E}(X_n)]/\sqrt{\text{Var}(X_n)} \stackrel{D}{\longrightarrow} N(0,1)$ Let $X_n \sim \mathbf{B}(n,n)$ (Beta distribution), with pdf
$$
f_n(x) = \frac{1}{\text{B}(n,n)}x^{n-1}(1 - x)^{n-1},~~ x \in (0,1).
$$
Knowing that $\text{E}(X_n) = 1/2$ and that $\text{Var}(X_n) = 1/[4(2n+1)]$, prove that
$$
2\sqrt{2n + 1}(X_n - \small{\frac{1}{2}}) \stackrel{D}{\longrightarrow} N(0,1).
$$

Attempt
Definition. A sequence of random variables $X_1, X_2, ...$, converges in distribution to a random variable X if
$$
\text{lim}_{n \to \infty} F_{X_n}(x) = F_X(x)
$$
So we have to prove that
$$
\text{lim}_{n \to \infty} F_{Y_n}(x) = \int_{-\infty}^{x} \frac{1}{ \sqrt{2\pi}} e^{-y^2/2}dy
$$
Where $Y_n = 2\sqrt{2n + 1}(X_n - \small{\frac{1}{2}}) $.
Now,
$$
\begin{align}
P(Y_n \leq x) & = P(2\sqrt{2n + 1}(X_n - \small{\frac{1}{2}}) \leq x) \\
              & = P(X_n - 1/2 \leq \frac{x}{2\sqrt{2n+1}} \\
              & = P(X_n \leq \frac{x}{2\sqrt{2n+1}} + 1/2) \\
              & = F_{X_n} \Bigl( \frac{x}{2\sqrt{2n+1}} + \frac{1}{2} \Bigr) \\
              & = \frac{1}{B(n,n)}\int_{0}^{ \frac{x}{2\sqrt{2n+1}} + 1/2 } t^{n-1}(1 - t)^{n-1}dt
\end{align}
$$
We use Stirling's approximation to $\text{B}(n,n)$:
$$
B(a, b) \approx \sqrt{2\pi} \frac{a^{a - 1/2}b^{b - 1/2}}{(a + b)^{a + b - 1/2}}
$$
So $\text{B}(n, n) \approx \frac{\sqrt{\pi}}{2^{2n - 1}} \frac{1}{\sqrt{n}} $, after simplification.
Substituting the Stirling approximation (we do this because it converges asymptotically and we're taking the limit), we obtain
$$
\frac{1}{\frac{\sqrt{\pi}}{2^{2n - 1}} \frac{1}{\sqrt{n}}}\int_{0}^{ \frac{x}{2\sqrt{2n+1}} + 1/2 } t^{n-1}(1 - t)^{n-1}dt.
$$
So what's left to do is prove that
$$
\text{lim}_{n \to \infty}  \frac{1}{\frac{\sqrt{\pi}}{2^{2n - 1}} \frac{1}{\sqrt{n}}}\int_{0}^{ \frac{x}{2\sqrt{2n+1}} + 1/2 } t^{n-1}(1 - t)^{n-1}dt = \int_{-\infty}^{x} \frac{1}{ \sqrt{2\pi}} e^{-y^2/2}dy.
$$
I don't know how to do this last step, finishing the proof. I asked my professor for guidance on how to finish the last step. All he said was "apply the limit theorem to solve directly".
 A: I was pondering how to formulate the simplest possible elementary solution to this problem and it occurred to me we can avoid any consideration of Beta functions (no Stirling's approximation needed; indeed, even information about the moments of Beta distributions is unnecessary).  The result is extremely general and, I hope, interesting.
Here, for the record, is what I will show:

Let $f$ be a positive multiple of any probability density function that is
bounded, unimodal, and twice differentiable in a neighborhood of
its mode.  Let the second derivative at the mode equal $-a$.  Then any sequence
of random variables $X_n$ with distribution functions proportional to
$$t\to f^n\left(\frac{t}{\sqrt{an}}\right)$$ converges in distribution to the
Standard Normal distribution.


Notation, assumptions, and preliminary simplifications
Permit me to use $n+1$ rather than $n$ as the index, so that $$f_n(t)\ \propto\ t^n(1-t)^n = (t(1-t))^n = f(t)^n$$ (for $0\le t\le 1$), thereby avoiding writing "$n-1$" too often.  In the question $f(t) = t(1-t)$ for $0\le t \le 1$ (and otherwise equals zero).  However, this formula is a distracting, irrelevant detail.
Here's all we need to assume about $f:$

*

*There is a constant $c$ for which $cf$ is a probability density function.  This means it is defined almost everywhere on all real numbers, integrable, with unit integral.  Obviously $c^{-1}=\int f(t)\,\mathrm{d}t.$

*$f$ is bounded and unimodal.  That is, $f$ has a unique finite maximum value.

*$f$ has a second derivative in a neighborhood of its mode.
These are clearly true of the $f$ in the question.
Letting $\mu$ be the mode, we may with no loss of generality analyze the function $t\to f(t-\mu),$ which has all the properties assumed of $f$ and whose mode is $0.$
Writing
$$f(t) = 1 - \frac{a}{2}\left(1 + g(t)\right)t^2,$$
the third assumption implies
$$\lim_{t\to 0} g(t) = 0$$
and there is some positive number $\epsilon$ for which whenever $|t|\le \epsilon,$ $g(t) \ge 0.$  Moreover, since $0$ is the unique mode, $a$ must be positive.
Without any loss of generality, replace $f$ by the function $t\to f(t)/f(0),$ making the largest value of $f$ exactly $1,$ attained at its mode $0.$
We are going to consider a sequence of probability density functions determined by powers of $f.$  First we need to normalize those powers, so let
$$c_n^{-1} = \int f^n (t)\,\mathrm{d}t.$$
This is always possible because
$$\int f^n(t)\,\mathrm{d}t \le \sup(f)\int f^{n-1}(t)\,\mathrm{d}t\  = \int f^{n-1}(t)\,\mathrm{d}t$$
shows recursively that the integrals of $f^n$ cannot increase and therefore are bounded.
A final preliminary manipulation is to standardize $f^n:$ we are going to analyze the sequence
$$f_n(t) = f\left(\frac{t}{\sqrt{an}}\right)^n.$$
The next few steps will show why this is effective at producing just the right cancellation of factors in the calculation.  First, though, let's look at an example.

As $n$ grows, $f$ spreads out from its mode, pushing all "satellites" out and dampening them, leaving a graph that rapidly approaches a multiple of a Normal pdf.  (The plot of $f$ in the upper left corner has not yet been rescaled to a height of $1$ at its mode.  The next plot of $f_1$ has been so scaled and is plotted on an $x$ axis expanded by a factor of $\sqrt{a}$ to show detail.)
Analysis
Let $t$ be any real number.  Once $n$ exceeds $N(t)=t^2 / (a\epsilon^2),$ $|t|/\sqrt{an}\le \epsilon$ puts this value into the neighborhood where $f$ behaves nicely.  From now on take $n\gt N(t).$
We are going to estimate the value of $f^n(t)$ by using logarithms.  This is the crux of the matter and it is where all the algebra is done.  Fortunately, it's easy:
$$\begin{aligned}
\log\left(f^n(t)\right) &= n \log(f(t)) \\
&= n \log f\left(\frac{t}{\sqrt{an}}\right) \\
&= n \log \left(1 - \frac{a}{2}\left(\frac{t}{\sqrt{an}}\right)^2\left(1 + g\left(\frac{t}{\sqrt{an}}\right) \right) \right) \\
&= n\log\left(1 - \frac{t^2}{2n}\left(1 + g\left(\frac{t}{\sqrt{an}}\right)\right)\right)
\end{aligned}$$
Because $g$ shrinks to $0$ for small arguments, a sufficiently large value of $n$ assures that the argument of the logarithm in that last expression is of the form $1-u$ for an arbitrarily small value of $u.$  This permits us to approximate the logarithm using Taylor's Theorem (with remainder), giving
$$\begin{aligned}
n\log\left(f^n(t)\right) &= -\frac{t^2}{2}\left(1 + g\left(\frac{t}{\sqrt{an}}\right)\right) + \frac{R}{n}\, \tilde{t}^4 \left(1 + g\left(\frac{\tilde t}{\sqrt{an}}\right)\right)^2
\end{aligned}$$
where $0\le |\tilde{t}| \le |t|$ and $R$ is some number (related to the remainder term in the Taylor expansion).  Taking the limit as $n\to\infty$ makes the remainder and all the $g()$ terms disappear, leaving
$$\lim_{n\to\infty} \log\left(f(t)^n\right) = -\frac{t^2}{2},$$
whence
$$\lim_{n\to\infty} f(t)^n = \exp\left(-\frac{t^2}{2}\right).$$
It follows (requiring only an intuitive, elementary proof) that the sequence of normalizing constants $c_n$ must approach the normalizing constant for the right hand side--which exists and, as is well known, equals $\sqrt{2\pi}.$  Consequently
$$\lim_{n\to\infty} f_n(t) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{t^2}{2}\right),$$
which is the standard Normal density $\phi.$
Conclusions
When $X_n$ is a sequence of random variables having densities $f_n,$ for every number $t$ the limit of their densities is $\phi(t).$  It follows easily that the limit of their distribution functions is $\Phi,$ the standard Normal distribution.
In the case of the Beta$(n,n)$ distributions, $f(t)=t(1-t)$ has a unique mode at $\mu=1/2,$ where it can be expressed (up to a constant multiple) as
$$4f(t) = 1 - \frac{8}{2}(t-1/2)^2.$$
From this we can read off the value $a=8.$  Following our preliminary simplifications, this says the distribution of $\sqrt{an}(X_n - \mu) = \sqrt{8n}(X_n-\mu)$ converges to the standard Normal distribution.  Because asymptotically the ratio of $\sqrt{8n}$ and $2\sqrt{2n+1}$ becomes unity, the statement in the original question is proven.
