How to find $\lim_{n \to \infty} \int_0^1 \cdots \int_0^1 \sqrt{x_1+\sqrt{x_2+\sqrt{\dots+\sqrt{x_n}}}}dx_1 dx_2\dots dx_n$ Here I mean the limit of the following sequence:
$$p_1=\int_0^1 \sqrt{x} ~dx=\frac{2}{3}$$
$$p_2=\int_0^1 \int_0^1 \sqrt{x+\sqrt{y}} ~dxdy=\frac{8}{35}(4 \sqrt{2}-1) =  1.06442\dots$$
$$p_3=\int_0^1 \int_0^1 \int_0^1 \sqrt{x+\sqrt{y+\sqrt{z}}} ~dxdydz =  1.242896586866\dots$$
$$p_4 \approx 1.314437693607766$$
$$p_5 \approx  1.34186271753784$$
Here the approximate values are computed by Mathematica. In principle every one of these integrals can be evaluated in closed form, but it becomes very complicated (see $p_3$ at the bottom of the post).
How can we find the limit at $n \to \infty$? It should be finite because of the range of variables chosen.

$$\lim_{n \to \infty}p_n=\lim_{n \to \infty} \int_0^1 \cdots \int_0^1 \sqrt{x_1+\sqrt{x_2+\sqrt{\dots+\sqrt{x_n}}}}dx_1 dx_2\dots dx_n=?$$

I find it very likely that $\lim_{n \to \infty}p_n=\phi$ (the Golden Ratio), but I'm not sure (this is not correct, see the comments).

Edit: With the help of Wolfram Alpha I tackled $p_3$ (see the updated numerical value above):

$$p_3=\frac{64}{135135} (2 \sqrt{3244081+2294881 \sqrt{2}}-664\sqrt{2}-1092\cdot 2^{3/4}+305)$$

This confirms my suspicions that there is no hope for apparent pattern in the first few $p_k$. Now an interesting challenge is to see how many $p_k$ can be realistically computed in closed form.
 A: Let us define:
\begin{equation}
I_3(a) := \int\limits_{[0,1]^3} \sqrt{x+\sqrt{y+\sqrt{z+a}}} dx dy dz
\end{equation}
Then by using elementary integration we have the following result:
\begin{eqnarray}
I_3(a) &=& \frac{32}{31} \left(
\sum\limits_{k=0}^3 |[\begin{array}{r} 3 \\ k \end{array}]|(-1)^{3-k} \frac{(u^+)^{\frac{9}{2}+k} - (u^-)^{\frac{9}{2}+k}}{\frac{9}{2}+k}
\right)+\\
&-&
\frac{32}{15}\left(
\sum\limits_{k=0}^3 \binom{3}{k}(-1)^{3-k} \frac{(u_1^+)^{\frac{7}{2}+k} - (u_1^-)^{\frac{7}{2}+k}}{\frac{7}{2}+k}
\right)+\\
&-&\frac{16}{21}\left(
\frac{(u^+_2)^{\frac{15}{4}} - (u^-_2)^{\frac{15}{4}}}{\frac{15}{4}}-\frac{(u^+_2)^{\frac{11}{4}} - (u^-_2)^{\frac{11}{4}}}{\frac{11}{4}}-\frac{(u^+_3)^{\frac{15}{4}} - (u^-_3)^{\frac{15}{4}}}{\frac{15}{4}}
\right)
\end{eqnarray}
where
\begin{eqnarray}
(u^+, u^-) &:=& (\sqrt{1+\sqrt{1+a}}+1,\sqrt{1+\sqrt{0+a}}+1) \\
(u_1^+,u_1^-) &:=& (\sqrt{0+\sqrt{1+a}}+1,\sqrt{0+\sqrt{0+a}}+1)\\
(u_2^+,u_2^-)&:=&(\sqrt{1+a}+1,\sqrt{0+a}+1)\\
(u_3^+,u_3^-)&:=&(\sqrt{1+a}+0,\sqrt{0+a}+0)
\end{eqnarray}
Now, clearly the next integral we need to compute is $I_4(a) := \int\limits_0^1 I_3(\sqrt{\xi+a}) d\xi$. All the integrals are can be expressed through elementary functions by substituting for the respective $u_j^{\pm}$ for $j=0,\cdots,3$ and integrating power functions. Therefore with some effort it is easy to compute higher elements of this sequence. Yet it is unclear for me at this stage if I might be able to find a neat expression for arbitrary elements.
