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?

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]


