definite integralsintegration
How to solve the integral
$$\int_0^a\sqrt{\frac{1}{x}-\frac{1}{a}}\cosh\left(\sqrt{x}\right)\textrm{d}x=\pi I_1(\sqrt{a})$$
where $a>0$ and $I_1$ is the modified Bessel function of $1^\textrm{st}$ kind and order $1$? Mathematica gives the solution but Rubi not, so it seems there is no simple integration procedure.
Edit: Related question concerning Mathematica
$$\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.
So it looks like Gary's right, and the same integral was also identified as an incomplete modified Bessel function in a different answer here by Paul. Some relevant papers are: Jones, Cicchetti and Faraone, and Shu and Shastri.
Best Answer
From the Digital Library of Mathematical Functions we have the integral formula $$ I_{\nu}(z) =\frac{\left(\frac{1}{2}z \right)^{\nu}}{\pi^{\frac{1}{2}}\Gamma\left(\nu + \frac{1}{2}\right)}\int_{-1}^{1}\left(1-t^2\right)^{\nu - \frac{1}{2}}e^{\pm zt} \ dt, \qquad \Re(\nu) > -\frac{1}{2} $$ Which for the case of $\nu=1$ reduces to$$ I_{1}(z) =\frac{z}{\pi}\int_{-1}^{1}\sqrt{1-t^2}e^{\pm zt} \ dt $$ Notice that because we can choose either the "plus" or "minus" symbol on the exponential and still get the same $I_{1}(z)$, we can say that \begin{align*} I_{1}(z) &= \frac{I_{1}(z) + I_{1}(z)}{2} \\ & =\frac{z}{\pi}\int_{-1}^{1}\sqrt{1-t^2}\frac{e^{\color{blue}{+}zt}}{2} \ dt + \frac{z}{\pi}\int_{-1}^{1}\sqrt{1-t^2}\frac{e^{\color{blue}{-}zt}}{2} \ dt\\ & =\frac{z}{\pi}\int_{-1}^{1}\sqrt{1-t^2}\frac{e^{\color{blue}{+}zt}+e^{\color{blue}{-}zt}}{2} \ dt \\ & =\frac{z}{\pi}\int_{-1}^{1}\sqrt{1-t^2}\cosh(zt) \ dt \\ & =\color{purple}{2}\frac{z}{\pi}\int_{\color{purple}{0}}^{1}\sqrt{1-t^2}\cosh(zt) \ dt \end{align*} where on the last step you notice that $\sqrt{1-t^2}\cosh(zt)$ is an even function in $t$. Lastly, we just need to make the change of variable $zt = \sqrt{x}\implies zdt =\frac{1}{2\sqrt{x}} dx$. We thus get \begin{align*} \pi I_{1}(z)& =2 z\int_{0}^{z^2}\sqrt{1-\frac{x}{z^2}}\cosh(\sqrt{x})\frac{1}{2z \sqrt{x}} \ dx =\int_{0}^{z^2}\sqrt{\frac{1}{x}-\frac{1}{z^2}}\cosh(\sqrt{x}) \ dx \end{align*} which gives your desired integral with $z^2 = a$.