How to determine the limiting distribution of $Y=(\bar{X})^2$, where $X=(X_1,...,X_n)$ and $X_i \overset{\mathrm{iid}}{\sim} \textrm{Unif}(0,1)$? For my introductory statistics class, I've been trying to determine the limiting distribution of $Y=(\bar{X})^2$, where $X_i \overset{\mathrm{iid}}{\sim}\textrm{ Unif}(0,1)~\forall ~i \in \{1,\ldots,n\}$.
At first, I've tried to obtain some data by doing experiments with R, just to get an idea of what my result should look like. After plotting a good amount of samples on a density histogram, I concluded that the shape pretty much looked like a normal distribution, so I fitted the data on that model.
This is the R code I wrote:
means = c()

for(a in 1:5000) {
  x = runif(5000, 0, 1)
  means = c(means, mean(x)**2)
}

library(fitdistrplus)
normal_ = fitdist(means, "norm")

plot(normal_)
print(normal_$estimate)

This code outputs something along these lines:
       mean          sd 
0.249928052 0.004100903 

The model seems to fit well as the QQ-plot is almost a perfect diagonal.
This being said, I couldn't find a way to express these empirical findings into something more formal. How can I calculate the asymptotic distribution in this case, from a theoretical point of view?
 A: Hint: Have a look at the mean and variance of the uniform distribution and use this to define the statistic:
$$Z_n \equiv \sqrt{12n} (\bar{X}_n - \tfrac{1}{2}),$$
giving you:
$$\bar{X}_n^2 = \frac{(Z_n + \sqrt{3n})^2}{{12n}}.$$
Now, have a think about how you can approximate the distribution of $Z_n$ and what the corresponding approximating distribution for $\bar{X}_n^2$ would be.  (Further hint: Using a normal distribution for the latter is asymptotically valid, but not a great option for good performance for lower values of $n$.  There is a much better approximating distribution that can be used here.)
A: In the comments you are close to a correct answer. You already have that $\sqrt{n}\left(\bar{X}_n-\mu\right) \overset{d}{\rightarrow} \mathop{\mathcal{N}}\left(0, \sigma^2\right)$ by the Lindeberg–Lévy version of the CLT. You have also correctly calculated $\mu = \frac{1}{2}$ and $\sigma^2 = \frac{1}{12}$.
Now (with the suitably scaled sequence), using the delta method as described in this Wikipedia article with $\mathop{g}\left(x\right) \mathrel{:=} x^2$ will work.1

1 You can ignore the assumption $g'(x) \neq 0$ if you allow for singular normal limits.
A: This is just an extended comment in that a display is produced showing how close the normal approximation is to the true probability density function.  The density for $Y$ for any sample size can be found as a function of the Irwin-Hall distribution (the sum of independent uniform random variables).
I've programmed this in Mathematica but I'm sure it can be easily done in R, Matlab, Maple, etc.
p = ConstantArray[0, 15];
Do[
 dist = TransformedDistribution[(x/n)^2, x \[Distributed] UniformSumDistribution[n]];
 pdf = PDF[dist, y];
 pdfN = PDF[NormalDistribution[Mean[dist], StandardDeviation[dist]], y];
 pdfA = PDF[NormalDistribution[1/4, Sqrt[1/(12 n)]], y];
 p[[n]] = Plot[{pdf, pdfN, pdfA}, {y, 0, 1},
   PlotLegends -> {"pdf of \!\(\*SuperscriptBox[OverscriptBox[\(X\), \(_\)], \(2\)]\)", 
     "Normal approximation\nwith same mean and variance", "Asymptotic normal"},
   PlotStyle -> {Blue, {Black, Dotted}, Red},
   PlotRange -> {{0, 1}, {0, 6}}, Frame -> True, 
   FrameLabel -> (Style[#, Bold, 14] &) /@ {"Y", "Probability density function"},
   PlotLabel -> Style["n = " <> ToString[n], Bold, 18]],
 {n, 1, 15}]

Export["display.gif", Flatten[{p, Reverse[p]}], 
 "DisplayDurations" -> 1]


