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?

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.)
