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?

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.
