Expectation of square root of sum of independent squared uniform random variables Let $X_1,\dots,X_n \sim U(0,1)$ be independent and identicallly distributed standard uniform random variables.
$$\text{Let }\quad Y_n=\sum_i^nX_i^2  \quad \quad \text{I seek: } \quad \mathbb{E}\big[\sqrt{Y_n } \big]$$

The expectation of $Y_n$ is easy:
$$\begin{align}
\mathbb{E}\left[X^2\right] &=\int_0^1\frac{y}{2\sqrt{y}}=\frac{1}{3}\\
\mathbb{E}\left[Y_n\right] &=\mathbb{E}\left[\sum_i^nX_i^2\right] = \sum_i^n\mathbb{E}\left[X_i^2\right]=\frac{n}{3}
\end{align}$$
Now for the boring part. To apply LOTUS, I would need the pdf of $Y_n$. Of course the pdf of the sum of two independent random variables is the convolution of their pdfs. However, here we have $n$ random variables and I guess the convolution would lead to a...convoluted expression (horrible pun intended). Is there a smarter way?
I would prefer to see the correct solution, but if it's impossible or too complicated, an asymptotic approximation for large $n$ could be acceptable. By Jensen's inequality, I know that 
$$\sqrt{\mathbb{E}[Y_n]}=\sqrt{\frac{n}{3}}\geq\mathbb{E}\left[\sqrt{Y_n}\right]$$
But this doesn't help me much, unless I can find also a non-trivial lower bound. Note that the CLT doesn't directly apply here, because we have the square root of the sum of independent RVs, not just the sum of independent RVs. Maybe there could be other limit theorems (which I ignore) that may be of help here.
 A: As an extended comment: it seems clear here that $E\left[\sqrt{Y_n}\right]= E\left[\sqrt{\sum_i X_i^2}\right]$ starts with $E\left[\sqrt{Y_n}\right] =\frac12=\sqrt{\frac{n}{3}-\frac{1}{12}}$ when $n=1$ and then approaches $\sqrt{\frac{n}{3}-\frac{1}{15}}$ as $n$ increases, related to the variance of $\sqrt{Y_n}$ falling from $\frac{1}{12}$ towards $\frac{1}{15}$. My linked question which S.Catterall answered provides a justification for the $\sqrt{\frac{n}{3}-\frac{1}{15}}$ asymptotic result based on each $X_i^2$ having mean $\frac13$ and variance $\frac{4}{45}$, and for the distribution being approximately and asymptotically normal. 
This question is effectively about the distributions of distances from the origin of random points in an $n$-dimensional unit hypercube $[0,1]^n$.  It is similar to a question on the distribution of distances between points in such a hypercube, so I can easily adapt what I did there to show the densities for various $n$ from $1$ to $16$ using numerical convolution.  For $n=16$, the suggested normal approximation shown in red is a good fit, and from $n=4$ you can see a bell curve appearing.  

For $n=2$ and $n=3$ you get a sharp peak at the mode of $1$ with what looks like the same density in both cases. Compare this with the distribution of $\sum_i X_i$, where the bell curve appears with $n=3$ and where the variance is proportional to $n$      
A: One approach is to first calculate the moment generating function (mgf) of $Y_n$ defined by $Y_n = U_1^2 + \dotsm + U_n^2$ where the $U_i, i=1,\dotsc, n$ is independent and identically distributed standard uniform random variables.
When we have that, we can see that
$$ \DeclareMathOperator{\E}{\mathbb{E}}
    \E \sqrt{Y_n} 
$$
is the fractional moment of $Y_n$ of order $\alpha=1/2$.  Then we can use results from the paper Noel Cressie and Marinus Borkent: "The Moment Generating Function has its Moments", Journal of Statistical Planning and Inference 13 (1986) 337-344, which gives fractional moments via fractional differentiation of the moment generating function.
First the moment generating function of $U_1^2$, which we write $M_1(t)$.
$$
   M_1(t) = \E e^{t U_1^2} = \int_0^1 \frac{e^{tx}}{2\sqrt{x}}\; dx 
$$
and I evaluated that (with help of Maple and Wolphram Alpha) to give
$$ \DeclareMathOperator{\erf}{erf}
 M_1(t)= \frac{\erf(\sqrt{-t})\sqrt{\pi}}{2\sqrt{-t}}        
$$ where $i=\sqrt{-1}$ is the imaginary unit.
(Wolphram Alpha gives a similar answer, but in terms of the Dawson integral.)  It turns out we will mostly need the case for $t<0$. Now it is easy to find the mgf of $Y_n$:
$$
    M_n(t) = M_1(t)^n
$$
Then for the results from the cited paper.  For $\mu>0$ they define the $\mu$th order integral of the function $f$ as
$$
   I^\mu f(t) \equiv \Gamma(\mu)^{-1} \int_{-\infty}^t (t-z)^{\mu-1} f(z)\; dz
$$
Then, for $\alpha>0$ and nonintegral, $n$ a positive integer, and $0<\lambda<1$ such that $\alpha=n-\lambda$.  Then the derivative of $f$ of order $\alpha$ is defined as
$$
   D^\alpha f(t) \equiv \Gamma(\lambda)^{-1}\int_{-\infty}^t (t-z)^{\lambda-1} \frac{d^n f(z)}{d z^n}\; dz.
$$
Then they state (and prove) the following result, for a positive random variable $X$: Suppose $M_X$ (mgf) is defined.  Then, for $\alpha>0$,
$$
   D^\alpha M_X(0) = \E X^\alpha < \infty
$$
Now we can try to apply these results to $Y_n$.  With $\alpha=1/2$ we find
$$
  \E Y_n^{1/2} = D^{1/2} M_n (0) = \Gamma(1/2)^{-1}\int_{-\infty}^0 |z|^{-1/2} M_n'(z) \; dz
$$
where the prime denotes the derivative.  Maple gives the following solution:
$$
  \int_{-\infty}^0 \frac{n\cdot\left(\erf(\sqrt{-z})\sqrt{\pi}-2e^z\sqrt{-z} \right)e^{\frac{n(-2\ln 2 +2 \ln(\erf(\sqrt{-z}))-\ln(-z)  +\ln(\pi))}{2}}}{2\pi(-z)^{3/2}\erf(\sqrt{-z})} \; dz
$$
I will show a plot of this expectation,  made in maple using numerical integration, together with the approximate solution $A(n)=\sqrt{n/3-1/15}$ from some comment (and discussed in the answer by @Henry).  They are remarkably close:

As a complement, a plot of the percentage error:

Above about $n=20$ the approximation is close to exact.  Below the maple code used:
int( exp(t*x)/(2*sqrt(x)), x=0..1 ) assuming t>0;
int( exp(t*x)/(2*sqrt(x)), x=0..1 ) assuming t<0;
M := t -> erf(sqrt(-t))*sqrt(Pi)/(2*sqrt(-t))
Mn := (t, n) -> exp(n*log(M(t)))
A  :=  n -> sqrt(n/3 - 1/15)
Ex :=  n ->   int( diff(Mn(z, n), z)/(sqrt(abs(z))*GAMMA(1/2) ), 
                   z=-infinity..0 , numeric=true)

plot([Ex(n), A(n)], n=1..100, color=[blue, red], legend= 
      [exact, approx], labels=[n, expectation], 
      title="expectation of sum of squared uniforms")
plot([((A(n)-Ex(n))/Ex(n))*100], n=1..100, color= 
      [blue], labels=[n, "% error"], 
      title="Percentage error of approximation")

