$$\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

and confirms the obtained formulas for the antiderivative.