Nested root integral $\int_0^1 \frac{dx}{\sqrt{x+\sqrt{x+\sqrt{x}}}}$

$$\color{brown}{\textbf{Simple example ($\mathrm{k=2}$).}}$$

Let us consider the integral $$N_2(x) = \int \dfrac{\mathrm dx}{\sqrt{x+\sqrt x}}.$$ Substitution $$y = \sqrt{1+\frac1{\sqrt x}},\quad x=\dfrac1{(y^2-1)^2},\quad \mathrm dx = -\dfrac {4y}{(y^2-1)^3}\,\mathrm dy,\quad x>0,\quad y>1, \tag{i1}$$ gives $$\sqrt{x+\sqrt x} = \sqrt{\dfrac1{(y^2-1)^2}+\dfrac1{y^2-1}} = \dfrac y{y^2-1},\tag{i2}$$ so $$N_2(x)=Q_2\left(\sqrt{1+\frac1{\sqrt x}}\right),\tag{i3}$$ where $$Q_2(y) = -4\int\dfrac{\mathrm dy}{(y^2-1)^2} = \int\dfrac{-4+4y^2}{(y^2-1)^2}\,\mathrm dy - \int \dfrac{4y^2}{(y^2-1)^2}\mathrm dy$$ $$ = 4\int\dfrac{\mathrm dy}{y^2-1} + 2\int y\,\mathrm d\left(\dfrac1{y^2-1}\right) = \dfrac{2y}{y^2-1}+ 2\int\dfrac{\mathrm dy}{y^2-1},$$ $$Q_2(y) = \dfrac {2y}{y^2-1} + \log\left(\frac{y-1}{y+1}\right) + \mathrm{const}.\tag{i4}$$ Taking in account $(\mathrm{i2})$ and $(\mathrm{i4}),$ solution $(\mathrm{i3})$ can be written as $$N_2(x) = 2\sqrt{x + \sqrt x} + \ln\dfrac{\sqrt{1+\frac1{\sqrt x}}-1}{\sqrt{1+\frac1{\sqrt x}}+1}+\mathrm{const} = 2\sqrt{x + \sqrt x} + \ln\dfrac{\sqrt{1+\sqrt x}-\sqrt[4]x}{\sqrt{1+\sqrt x}+\sqrt[4]x}+\mathrm{const},$$ or, eliminating the numerator, $$\boxed{N_2(x) = 2\sqrt{x + \sqrt x} - \ln\left(1 + 2\sqrt x + 2\sqrt{x+\sqrt x}\right) + \mathrm{const}}.$$

$$\color{brown}{\textbf{OP antiderivative ($\mathrm{k=3}$).}}$$

$\color{brown}{\underline{\textrm{Transformation to single radical.}}}$

Substitution.

Taking in account $(\mathrm i1)-(\mathrm i4),$ can be written $$N_3(x) = \int \dfrac{\mathrm dx}{\sqrt{x+\sqrt {x+\sqrt x}}} = Q_3\left(\sqrt{1+\frac1{\sqrt x}}\right),\tag{1}$$ where $$\begin{align} &\sqrt{x+\sqrt{x+\sqrt x}} = \sqrt{\dfrac1{(y^2-1)^2}+\dfrac y{y^2-1}},\\[4pt] &Q_3(y) = \int\dfrac{-\dfrac {4y}{(y^2-1)^3}\,\mathrm dy}{\sqrt{\dfrac1{(y^2-1)^2}+\dfrac y{y^2-1}}} = \int\dfrac{-\dfrac {4y}{(y^2-1)^3}-\dfrac{y^2+1}{(y^2-1)^2}}{\sqrt{\dfrac1{(y^2-1)^2}+\dfrac y{y^2-1}}}\,\mathrm dy + \int\dfrac{\dfrac{y^2+1}{(y^2-1)^2}\,\mathrm dy }{\sqrt{\dfrac1{(y^2-1)^2}+\dfrac y{y^2-1}}}\\[4pt] & = \sqrt{\dfrac1{(y^2-1)^2}+\dfrac y{y^2-1}} + \int \dfrac{y^2+1}{y^2-1} \dfrac{\mathrm dy}{\sqrt{y^3-y+1}}. \end{align}$$

Then $$N_3(x) = 2\sqrt{x+\sqrt {x+\sqrt x}} +Q_{30}\left(\sqrt{1+\frac1{\sqrt x}}\right),\tag2$$ where $$Q_{30}(y) = \int \dfrac{y^2+1}{y^2-1} \dfrac{\mathrm dy}{\sqrt{y^3-y+1}}.\tag3$$

In this way, substitution $(\mathrm i1)$ has transformed given nested root integral to elliptic type.

Algebraic transformations.

Denote $$\begin{cases} P_3(y) = y^3-y+1\\[4pt] r = \sqrt[3]{\dfrac{9+\sqrt{69}}{18}} + \sqrt[3]{\dfrac{9-\sqrt{69}}{18}}\approx 1.32471\,79572\,44746\\[4pt] \lambda = \sqrt{3r^2-1} = \dfrac12\sqrt{4+\sqrt[3]{800+96\sqrt69}+\sqrt[3]{800-96\sqrt69}} \approx 2.06509\,87866\,78274\\[4pt] w = \dfrac12+\dfrac34\dfrac r\lambda\approx 0.98110\,94143\,9836556,\\[4pt] p = r+\lambda \approx 3.38981\,67439\,23020. \end{cases}\tag4$$

Then \begin{align} &P_3(-r)=-\left(\dfrac{9+\sqrt{69}}{18}+\dfrac{9-\sqrt{69}}{18}\right) -3\left(\sqrt[3]{\dfrac{9+\sqrt{69}}{18}} \cdot \sqrt[3]{\dfrac{9-\sqrt{69}}{18}}\right)\\ &\times\left(\sqrt[3]{\dfrac{9+\sqrt{69}}{18}} + \sqrt[3]{\dfrac{9-\sqrt{69}}{18}}\right) +\sqrt[3]{\dfrac{9+\sqrt{69}}{18}} + \sqrt[3]{\dfrac{9-\sqrt{69}}{18}}+1\\ &=-3\sqrt[3]{\dfrac{12}{324}}\left(\sqrt[3]{\dfrac{9+\sqrt{69}}{18}} + \sqrt[3]{\dfrac{9-\sqrt{69}}{18}}\right)+\sqrt[3]{\dfrac{9+\sqrt{69}}{18}} +\sqrt[3]{\dfrac{9-\sqrt{69}}{18}}=0, \end{align} i.e. $(-r)$ is the root of $P_3(y),$ $$r^3=r+1\tag5.$$

According to the Bezou theorem, $P_3(y)$ allows the decomposition in the form of $$P_3(y) = y^3-y+1 = (y+r)(y^2-ry+r^2-1) = (y+r)\big((y+r)^2-3r(y+r)+3r^2-1\big),$$ or $$P_3(y) = (y+r)\big((y+r)^2-3r(y+r)+\lambda^2).\tag6$$

On the other hand, for any constant $p$ can be obtained the fractional decomposition of $$R_3(y) = \dfrac{y^2+1}{y^2-1} = A + B\dfrac{y+p}{y-1} + C\dfrac{y+p}{y+1},$$ where $$B = \lim\limits_{y\to1}\dfrac{y-1}{y+p}R_3(y) = \dfrac1{p+1},$$ $$C = \lim\limits_{y\to-1}\dfrac{y+1}{y+p}R_3(y) = -\dfrac1{p-1},$$ $$A = \lim\limits_{y\to\infty}\left(R_3(y)-B\dfrac{y+p}{y-1} - C\dfrac{y+p}{y+1}\right) = \dfrac{p^2+1}{p^2-1}.$$

Therefore, $$\begin{align} &Q_{30}(y) = \dfrac{p^2+1}{p^2-1}I_0(y) + \dfrac1{p+1}I_1(y,p+1) - \dfrac1{p-1} I_1(y,p-1),\tag7\\[4pt] &\text{where}\\[4pt] &\begin{cases} I_0(y) = \int\dfrac{\mathrm dy}{\sqrt{y^3-y+1}}\\ I_1(y,q) = \int\dfrac{y+p}{y+p-q}\dfrac{\mathrm dy}{\sqrt{y^3-y+1}}\\ \end{cases}\tag8 \end{align}$$

$\color{brown}{\underline{\textrm{The first elliptic integral.}}}$

Substitution $(\mathrm i1)$ transforms the interval $x\in(0,1)$ to the interval $y\in(\infty,\sqrt2).$ Should be found the real-valued antiderivative for this interval.

Substitution $$y=g(t) = -r+\lambda\tan^2t,\quad \mathrm dy = 2\lambda\dfrac{\tan t}{\cos^2t}\,\mathrm dt, \quad t=\arctan\sqrt{\dfrac{y+r}\lambda}\tag{f1}$$

for the presentation $(6)$ gives \begin{align} &I_0(g(t)) = \int\dfrac{2\lambda\dfrac{\tan t}{\cos^2t}}{\sqrt{\lambda\tan^2t\left(\lambda^2\tan^4t -3r\lambda\tan^2t+\lambda^2\right)}}\mathrm dt\\[4pt] &=\int\dfrac{2\mathrm dt}{\sqrt{\lambda(\cos^2t+\sin^2t)^2-(2\lambda+3r)\sin^2t\cos^2t}}\\[4pt] &=\dfrac1{\sqrt \lambda}\int \dfrac{\mathrm d(2t)}{\sqrt{1-\dfrac{2\lambda+3r}{4\lambda}\sin^2 2t}} = \dfrac1{\sqrt \lambda}\mathrm F(2t\ |\ w)+\mathrm{const}, \end{align}

where

$$\quad F(\varphi\ |\ w) = \int\limits_0^\varphi \dfrac{\mathrm d\varphi}{\sqrt{1-w\sin^2\varphi}}\tag{f3}$$ is the elliptic integral of the first kind.

Assuming the domain of $x\in(0,\infty)$ and taking in account $(\mathrm f1),$ one can get $$\left(y\in(1,\infty)\right)\wedge(y+r\in(\lambda,\infty)) \Rightarrow t\in\left(\dfrac\pi4,\dfrac\pi2\right),$$ $$\sin^2 2t = \dfrac{4\tan^2 t}{(1+\tan^2 t)^2},$$ then $$\sin^2 2t = \sin^2(\pi-2t) = \dfrac{4\lambda(y+r)}{(y+r+\lambda)^2},\quad \cos 2t = -\cos(\pi-2t) = \dfrac{\lambda - y - r}{\lambda+y+r}\tag{f4}$$ $$\pi - 2t = 2\arctan\sqrt{\dfrac\lambda{y+r}} = \arccos\dfrac{y+r-\lambda}{y+r+\lambda}.\tag{f5}$$

At the same time,

$$\mathrm F(\pi-\varphi\ |\ w) = \mathrm F\left(\dfrac\pi2\ \bigg|\ w\right) - \mathrm F(\varphi\ |\ w).\tag{f6}$$

Therefore, $$\color{green}{\mathbf{\dfrac{p^2+1}{p^2-1}I_0(y) = -\dfrac1{\sqrt \lambda}\dfrac{(\lambda+r)^2+1}{(\lambda+r)^2-1} F\left(\arccos\dfrac{y+r-\lambda}{y+r+\lambda}\ \Bigg|\ w\right) + const. \tag{f7}}}$$

Note that the choice $\mathrm{const}=0$ provides the condition $I_0(\infty)= 0.$

Substitution $(\mathrm f1)$ effectively transforms the expression under the radical and provides simple form of solution for $I_0(y)$.
Effective integration of $I_1(y,q)$ looks more serious problem. Decomposition in the form $(7),(8)$ for choosen value $(4)$ of parameter $p$ allows to simplify integration.

$\color{brown}{\underline{\textrm{The second elliptic integral.}}}$

From $(4),(8)$ should $$I_1(y,q) = \int\dfrac{y+r+\lambda}{y+r+\lambda-q}\dfrac{\mathrm dy}{\sqrt{y^3-y+1}}.\tag{s1}$$

Applying substitution $(\mathrm f1)$ and taking in account $(4),$ one can get \begin{align} &I_1(g(t),q)=\int \dfrac{\lambda(1+\tan^2 t)}{\lambda(1+\tan^2 t) -q}\cdot \dfrac{2\mathrm dt}{\sqrt \lambda\sqrt{1-w\sin^2 2t}} =\int \dfrac1{\lambda-q\cos^2 t}\cdot \dfrac{\mathrm 2\sqrt \lambda\,\mathrm dt}{\sqrt{1-w\sin^2 2t}}\\[4pt] &= \int \dfrac1{2\lambda-q-q\cos 2t}\cdot \dfrac{4\sqrt \lambda\,\mathrm dt}{\sqrt{1-w\sin^2 2t}} = \int \dfrac{2\lambda-q+q\cos 2t}{(2\lambda-q)^2-q^2\cos^2 2t}\cdot \dfrac{4\sqrt \lambda\,\mathrm dt}{\sqrt{1-w\sin^2 2t}}\\[4pt] &= -\int \dfrac{2\lambda-q+q\cos 2t}{4\lambda(q-\lambda)-q^2\sin^2 2t}\cdot \dfrac{4\sqrt \lambda\,\mathrm dt}{\sqrt{1-w\sin^2 2t}}\\[4pt] &= -\int \dfrac{2\lambda-q-q\cos(\pi-2t)}{4\lambda(q-\lambda)-q^2\sin^2(\pi-2t)}\cdot \dfrac{4\sqrt \lambda\,\mathrm dt}{\sqrt{1-w\sin^2(\pi-2t)}}\\[4pt] &= \int \dfrac{\dfrac{2\lambda-q}{2\sqrt \lambda(q-\lambda)} -\dfrac{q\cos(\pi-2t)}{2\sqrt \lambda(q-\lambda)}} {1-\dfrac{q^2}{4\lambda(q-\lambda)}\sin^2(\pi-2t)}\cdot \dfrac{\mathrm d(\pi-2t)}{\sqrt{1-w\sin^2(\pi-2t)}}, \end{align} $$\begin{align} &\mathbf{\dfrac1q I_1(g(t),q)= A(q)\Pi\left(u(q),\pi-2t\ |\ w \right) -B(q)P\left(u(q),\pi-2t\ |\ w \right) + const}, \end{align}\tag{s2}$$

where $$\begin{cases} A(\lambda+r\pm1)=\dfrac1{2\sqrt\lambda}\dfrac{1\pm(r-\lambda)}{(r\pm1)(\lambda+r\pm1)} \approx \dbinom{0.00885\,15573\,11770}{-0.78032\,04227\,56824}\\[4pt] B(\lambda+r\pm1) = \dfrac1{2\sqrt\lambda(r\pm1)}\approx\dbinom{0.14966\,81250\,93454}{1.07150\,27311\,21398}\\[4pt] u(\lambda+r\pm1)=\dfrac{(\lambda+r\pm1)^2}{4\lambda(r\pm1)}\approx\dbinom{1.00350\,99620\,67355}{2.12922\,75171\,98979}, \end{cases}\tag{s3}$$

$$\Pi(u,\varphi\ |\ w) = \int\limits_0^\varphi \dfrac{\mathrm d\varphi}{(1-u\sin^2\varphi)\sqrt{1-w\sin^2\varphi}}\tag{s4}$$ is the elliptic integral of the third kind, and $$\mathrm P(u,\varphi\ |\ w) = \int\limits_0^\varphi \dfrac{\cos\varphi\,\mathrm d\varphi}{(1-u\sin^2\varphi)\sqrt{1-w\sin^2\varphi}} =\dfrac1{2\sqrt{u-w}}\ln\left|\dfrac{\sqrt{u-w}+\sqrt{\csc^2\varphi-w}} {\sqrt{u-w}-\sqrt{\csc^2\varphi-w}}\right|\tag{s5}$$ (see also Wolfram Alpha test) can be considered as "pseudo-elliptic" integral, linked with known rational integral in the form of $$\int \dfrac{\mathrm ds}{(as^2+b)\sqrt{fs^2+g}}= {\small \begin{cases} \dfrac 1{\sqrt{b(ag-bf)}}\arctan \dfrac{s\sqrt{ag-bf}}{\sqrt{b(fs^2+g)}},\quad\text{if}\quad ag-bf > 0\\[4pt] \dfrac1{2\sqrt{b(bf-ag)}}\ln \left|\dfrac{\sqrt{b(fs^2+g)}+s\sqrt{bf-ag}}{\sqrt{b(fs^2+g)}-s\sqrt{bf-ag}}\right| \quad\text{if}\quad ag-bf < 0. \end{cases}}\tag{s6} $$

$\color{brown}{\underline{\textrm{Transformations of the second integral solution.}}}$

Solving of the convergency problem.

From $(4),(\mathrm s3)$ should $w < 1 < u(q),$ and the first inequality provides real-valued final expression of $(\mathrm s3)$.

On the other hand, the second inequality presents hyperbolic case of $\Pi-$function. This leads to the convergency problem, because the denominator equals to zero if $u(q)\sin^2(\pi-2t) = 1.$

Taking in account $(\mathrm f4)$ and $(\mathrm s3),$ one can get \begin{align} &u(p+1)\sin^2(\pi-2t) = \dfrac{y+r}{(\lambda+y+r)^2}\left(\dfrac{1+r}{(\lambda+1+r)^2}\right)^{-1}\in(0,1),\\[4pt] &u(p-1)\sin^2(\pi-2t) =\dfrac{y+r}{(\lambda+y+r)^2}\left(\dfrac{r-1}{(\lambda+r-1)^2}\right)^{-1}\in(0,3.897), \end{align} $$u(p+1)\sin^2(\pi-2t)\bigg|_{y=1}=1,\tag{t1}$$ $$u(p-1)\sin^2(\pi-2t)\bigg|_{y=\dfrac{2r^2+r-1}{r+1}\approx 11.8} = 1.\tag{t2}$$

To avoid this problem, the integral $(\mathrm s2)$ should be transformed.

Let $$u(q)v(q)=w,\quad b=\sqrt{(u(q)-1)(1-v(q))},\tag{t3}$$

then $$\dfrac1{1-u\sin^2 \varphi}+\dfrac1{1-v\sin^2\varphi} = 1+\dfrac{1-w\sin^4\varphi}{1-(u+v)\sin^2\varphi+w\sin^4\varphi},$$ \begin{align} &\dfrac {\mathrm d}{\mathrm d\varphi}\left(\ln\left(\sqrt{1-w\sin^2 \varphi}+b\tan\varphi\right) - \ln\left(\sqrt{1-w\sin^2\varphi}-b\tan\varphi\right)\right)\\[4pt] &=\dfrac{2b(w\sin^2\varphi\tan^2\varphi-\sec^2\varphi)}{\sqrt{1-w\sin^2\varphi}(b^2\tan^2\varphi+w\sin^2\varphi-1)}\\[4pt] &=\dfrac{2b(1-w\sin^4\varphi)}{\sqrt{1-w\sin^2\varphi} (-b^2\sin^2\varphi+(1-w\sin^2\varphi)(1-\sin^2\varphi))}\\[4pt] &=\dfrac{2b(1-w\sin^4\varphi)} {\sqrt{1-w\sin^2\varphi}(1-(b^2+w+1)\sin^2\varphi+w\sin^4\varphi)}\\[4pt] &=\dfrac{2b(1-w\sin^4\varphi)} {\sqrt{1-w\sin^2\varphi}(1-(u+v)\sin^2\varphi+w\sin^4\varphi)}. \end{align}

So $$\mathbf{\mathrm \Pi(u(q),\varphi,w) = \dfrac1{2b}\ln \dfrac{\sqrt{\csc^2\varphi-w}\cos\varphi+b}{\sqrt{\csc^2\varphi-w}\cos\varphi-b} +\mathrm F(\varphi,w) - \mathrm \Pi(v(q),\varphi,w),\tag{t4}}$$ wherein $$v(\lambda+r\pm1) = \dfrac{(3r+2\lambda)(r\pm1)}{(\lambda+r\pm1)^2}\approx\dbinom{0.97767\,78023\,97853}{0.46078\,18593\,70775}.\tag{t5}$$

Simplifications.

Firstly, the factor near $\mathrm F(\varphi,w)$ is \begin{align} &-\dfrac1{\sqrt\lambda}\dfrac{p^2+1}{p^2-1}+A(p+1) - A(p-1)\\[4pt] &= -\dfrac1{\sqrt \lambda((\lambda+r)^2-1)}\bigg((\lambda+r)^2+1+\dfrac{(\lambda+r-1)(r-\lambda+1)}{2(r+1)}+\dfrac{(\lambda+r+1)(\lambda-r+1)}{2(r-1)}\bigg)\\[4pt] &= -\dfrac1{\sqrt \lambda((\lambda+r)^2-1)}\bigg((\lambda+r)^2+1+\dfrac{(r^2-(\lambda-1)^2)(r-1)+((\lambda+1)^2-r^2)(r+1)}{2(r^2-1)}\bigg)\\[4pt] &= -\dfrac1{\sqrt \lambda((\lambda+r)^2-1)}\bigg((\lambda+r)^2+1 +\dfrac{-r^2+2\lambda r+\lambda^2+1}{r^2-1}\bigg)\\[4pt] &= -\dfrac1{\sqrt\lambda((\lambda+r)^2-1)}\bigg((\lambda+r)^2-1 +\dfrac{r^2+2\lambda r+\lambda^2-1}{r^2-1}\bigg)\\[4pt] &= -\dfrac1{\sqrt\lambda((\lambda+r)^2-1)}((\lambda+r)^2-1+r((r+\lambda)^2-1)) = -\dfrac{r+1}{\sqrt\lambda}. \end{align}

At the second, applying formulas $(\mathrm f4),(\mathrm s3),(\mathrm t1),(\mathrm t5),$ one can get: $$\begin{cases} u(p\pm1)-1 = \dfrac{(r\pm1+\lambda)^2}{4\lambda(r\pm1)}-1 = \dfrac{(r\pm1-\lambda)^2}{4\lambda(r\pm1)},\\[4pt] 1-v(p\pm1) = \dfrac{u(p\pm1)-w}{u(p\pm1)} =\dfrac1{(r\pm1)(\lambda+r\pm1)^2},\\[4pt] b = \sqrt{(u-1)(1-v)} = \dfrac1{2\sqrt\lambda(r\pm1)}\dfrac{1\pm(r-\lambda)} {\lambda+r\pm1} = A(q), \end{cases}\tag{t6}$$ so the factor near the logarithm in formulas $(\mathrm t4)$ is $\frac12.$

Besides, \begin{align} &\csc^2(\pi-2t) - w = \dfrac{(y+r+\lambda)^2}{4\lambda(y+r)}-\dfrac12 - \dfrac{3r}{4\lambda} = \dfrac{(y+r)^2 + \lambda^2-3r(y+r)}{4\lambda(y+r)} = \dfrac{y^3-y+1}{4\lambda(y+r)^2}, \end{align}

then $$2\sqrt\lambda\sqrt{\csc^2\varphi-w}\cos\varphi = \dfrac{\sqrt{y^3-y+1}}{y+r}\dfrac{y+r-\lambda}{y+r+\lambda} = \sqrt{y^3-y+1}\left(\dfrac2{y+r+\lambda}-\dfrac1{y+r}\right).\tag{t7}$$

At last, \begin{align} &u(p\pm1)-w = \dfrac{(\lambda+r\pm1)^2}{4\lambda(r\pm1)}-\dfrac12-\dfrac{3r}{4\lambda}\\[4pt] &= \dfrac{3r^2-1+2\lambda(r\pm1)+(r\pm1)^2-2\lambda(r\pm1)-3r(r\pm1)}{4\lambda(r\pm1)} =\dfrac{r^2\mp r}{4\lambda(r\pm1)} =\dfrac{r^3-r}{4\lambda(r\pm1)^2},\\[4pt] \end{align} $$\begin{cases} u(p\pm1)-w = \dfrac1{4\lambda(r\pm1)^2}\\[4pt] \sqrt{\dfrac{\csc^2(\pi-2t) - w}{u(p\pm1)-w}} = \dfrac{y\pm1}{y+r}\sqrt{y^2-y+1}\\[4pt] \dfrac{B(p\pm1)}{2\sqrt{u(p\pm1)-w}} =\dfrac12\dfrac{2\sqrt \lambda(r\pm1)}{2\sqrt \lambda(r\pm1)} = \dfrac12,\\[4pt] \end{cases}\tag{t8}$$

so $$\mathbf{B(p\pm1)\mathrm P\left(u(p\pm1),\arccos\dfrac{y+r-\lambda}{y+r+\lambda}\ \Big|\ w \right) = \dfrac12\ln\left|\ \dfrac {y+r + (r\pm1)\sqrt{y^3-y+1\phantom{\Big|}}} {y+r - (r\pm1)\sqrt{y^3-y+1\phantom{\Big|}}}\ \right|. \tag{t9}}$$

Closed form of the antiderivative.

Finally, the closed form of antiderivative can be presented by the formula $(2),$ where $$\color{green}{\boxed{\mathbf{\begin{align} &Q_{30}(\infty)-Q_{30}(y) = \dfrac{r+1}{\sqrt\lambda} \mathrm F\left(\arccos\dfrac{y+r-\lambda}{y+r+\lambda}\ \Bigg|\ w\right)\\[4pt] &+\dfrac12\ln\left|\dfrac{\dfrac{\sqrt{y^3-y+1}}{y+r}\dfrac{y+r-\lambda}{y+r+\lambda}+C_1} {\dfrac{\sqrt{y^3-y+1}}{y+r}\dfrac{y+r-\lambda}{y+r+\lambda}-C_1}\right| +\dfrac12\ln\left|\dfrac{\dfrac{\sqrt{y^3-y+1}}{y+r}\dfrac{y+r-\lambda}{y+r+\lambda}+C_2} {\dfrac{\sqrt{y^3-y+1}}{y+r}\dfrac{y+r-\lambda}{y+r+\lambda}-C_2}\right|\\[4pt] &-\dfrac{C_1}{2\sqrt\lambda}\Pi\left(\dfrac{(3r+2\lambda)(r+1)}{(\lambda+r+1)^2},\arccos\dfrac{y+r-\lambda}{y+r+\lambda}\ \Bigg|\ w \right)\\[4pt] &-\dfrac{C_2}{2\sqrt\lambda}\Pi\left(\dfrac{(3r+2\lambda)(r-1)}{(\lambda+r-1)^2}, \arccos\dfrac{y+r-\lambda}{y+r+\lambda}\ \Bigg|\ w \right)\\[4pt] & +\dfrac12\ln\left| \dfrac {y + r + (r+1)\sqrt{y^3-y+1\phantom{\Big|}}} {y + r - (r+1)\sqrt{y^3-y+1\phantom{\Big|}}}\right| -\dfrac12\ln\left|\dfrac {y + r + (r-1)\sqrt{y^3-y+1\phantom{\Big|}}} {y + r - (r-1)\sqrt{y^3-y+1\phantom{\Big|}}}\right|,\\[4pt] &C_1 = \dfrac{r-\lambda+1}{(r+1)(r+\lambda+1)},\quad C_2 = \dfrac{\lambda-r+1}{(r-1)(r+\lambda-1)}, \end{align}}}\tag{*}}$$ wherein the parameters $r,\,\lambda,\,w$ are given by $(4),$ and $Q_{30}(\infty) = 0.$

Since \begin{align} &\int_{\sqrt2}^\infty \dfrac {y^2+1}{y^2-1} \dfrac{\mathrm dy}{\sqrt{y^3-y+1}} =Q_{30}(\infty) - Q_{30}(\sqrt2) \end{align} $\approx4.01502\,83666\,96210\,78140$ (term1) $+0.33133\,29488\,62457\,95551$ (term2)

$+0.03549\,45749\,48612\,01694$ (term3) $-0.13827\,66191\,08129\,76996$ (term4)

$-2.80594\,91779\,29332\,44409$ (term5) $+0.80579\,17865\,34765\,72804$ (term6+term7)

$\approx 2.24342\,18800\,04584\,26784,$

then $$\int_0^1 \dfrac{\mathrm dx}{\sqrt{x+\sqrt{x+\sqrt x}}} = 2\sqrt{1+\sqrt2}-(Q_{30}(\infty)-Q_{30}(\sqrt2)) \approx 0.86412\,60680\,55497.$$

This result corresponds with the numeric calculations

Numeric calculations of the issue integral

and confirms the obtained formulas for the antiderivative.