On the closed-form of $\int_0^1\int_0^1\int_0^1\frac{dxdydz}{1-\frac{z}{3}(x+\sqrt{xy}+y)}$

Question

I would like to know if it is possible to calculate in closed-form, or well what work can be done about it, the definite integral $$\int_0^1\int_0^1\int_0^1\frac{3dxdydz}{3-z(x+\sqrt{xy}+y)},\tag{1}$$ where I was inspired in a well-known integral representation for the Apéry constant involving the volume $x\cdot y\cdot z$ in the denominator, and in the formula for the volume of a square frustum of basis $\sqrt{x}$ and $\sqrt{y}$ and height $z$, as reference for all I add the Wikipedia Heronian mean.

It is easy to check the integration of the logartihm $$\int_0^1 \frac{1}{3-z(x+\sqrt{xy}+y)}dz=\frac{\log 3-\log(3-x-\sqrt{xy}-y)}{x+\sqrt{xy}+y}$$ but the computations using a CAS (and standard time of computation with my computer), that I tried after this step, seems to me very tedious to evaluate. I suspect that now an important key to evaluate it should be to exploit symmetry or a suitable change of variable.

Question. I would like to know if it is possible to evaluate in closed-form (in terms of well-known constants and/or particular values of special functions) previous definite integral $(1)$. If isn't feasible to get the closed-form that I evoke explain why or add what work can be done. Many thanks.

If this integral is in the literature feel free to refer the literature in your answer or comment and I can to read the result from the literature.

Edit: (see comments, please). I don't know if from this step one can to get the integral in closed-form. Feel free to do more feedback, many thanks.

Using previous hints in comments I can to write that it is possible to reduce the integral over $xy$ to one-dimensional integral by using polar coordinates in $xy$-plane

$$\int_0^1\int_0^1\frac{dxdy}{1-\frac{z}{3}(x+\sqrt{xy}+y)}=2\int_0^{\frac{\pi}{4}}\int_0^{\sec \theta}\frac{drd\theta}{\frac{3}{r}-z(\cos\theta+\sqrt{\cos\theta\sin\theta}+\sin\theta)},$$ where the inner integral is equals to $$\int_0^{\sec\theta}\frac{dr}{\frac{3}{r}-z(\cos\theta+\sqrt{\cos\theta\sin\theta}+\sin\theta)}=\frac{-A\sec\theta-3\log(3-A\sec\theta)+3\log 3}{A^2},$$ being $A=z(\cos\theta+\sqrt{\cos\theta\sin\theta}+\sin\theta)$.

---

Since the integral seems very difficult I am going to accept an answer showing what work can be done (see the Question) as soon as expires the bounty.

Answer 1

Let me expand the integrand in powers of $z$ and integrate over $z$, $$I=\int_0^1 dx\int_0^1 dy\int_0^1 dz\;\frac{3}{3-z(x+\sqrt{xy}+y)}$$ $$\qquad\qquad=\sum_{n=0}^\infty\,\frac{3^{-n}}{n+1} \int_0^1 dx\int_0^1 dy\;(x+\sqrt{xy}+y)^n.$$ The integral over $x$ and $y$ is an element $c_{n}\in\mathbb{Q}$, $$I=\sum_{n=0}^\infty\,\frac{3^{-n}}{n+1}c_n,\;\;c_n=\int_0^1 dx\int_0^1 dy\;(x+\sqrt{xy}+y)^n,$$ $$\{c_0,c_1,c_2,\ldots\}=\{1,\tfrac{13}{9},\tfrac{149}{60},\tfrac{1667}{350},\tfrac{18623}{1890},\tfrac{69667}{3234},\ldots\}.$$ I have not yet succeeded in finding a closed-form expression for $c_n$ valid for all integer $n$, but that seems doable. At least we would then have the desired integral as a sum over rational coefficients.

---

_{With FusRoDah's input for $c_n$, and carrying out the sum over $n$, I find $$I=\sum_{a,b=0}^\infty\frac{12 \Gamma (a+b+1) \left(3^{-a-\frac{b}{2}} B_{\frac{1}{3}}\left(\frac{b}{2}+1,-a-b\right)-B_{\frac{1}{3}}(a+b+1,-a-b)\right)}{(2 a+b) (2 a+b+2) \Gamma (a+1) \Gamma (b+1)},$$ (with $B$ the incomplete beta function).}