Integrate $\sqrt{1+x^2+y^2}$ Calculate $$I=\int_{-1}^1\int_{-1}^1\sqrt{1+x^2+y^2}\,\mathrm{d}y\,\mathrm{d}x.$$
It's a problem from a book about calculus. My attempt:
$$\begin{align}
I &= \int_{-1}^1\int_{-1}^1\sqrt{1+x^2+y^2}\,\mathrm{d}y\,\mathrm{d}x \\
  &= \int_{-1}^1\left.\frac{x}{2}\sqrt{1+x^2+y^2}+\frac{1+y^2}{2}\log\left(x+\sqrt{1+x^2+y^2}\right)\right|_{-1}^{1}\,\mathrm{d}y \\
  &= \int_{-1}^1\sqrt{2+y^2}+\frac{1+y^2}{2}\left(\log \left(\sqrt{2+y^2}+1\right)-\log\left(\sqrt{2+y^2}-1\right)\right)\,\mathrm{d}y \\
  &= \sqrt{3}+2\operatorname{arsinh}\frac{1}{\sqrt{2}}+\int_{-1}^1(y^2+1) \operatorname{arsinh} \frac{1}{\sqrt{y^2+1}}\,\mathrm{d}y \\
  &= \color{red}\ldots \\
  &= -\frac{2}{9} (\pi + 12 \log 2 - 6 \sqrt{3} - 24 \log (1+ \sqrt{3}))
\end{align}$$
(answer taken from the solutions, no idea how to reach it).
[edit] Here is an attempt with polar coordinates. Due to symmetry, it's enough to integrate over $0 \le x \le 1$ and $0 \le y \le x$, 1/8th of the initial square.
$$\begin{align}
I &= 8\int_0^{\pi/4} \int_0^{1/\cos \theta}r \sqrt{1+r^2}\,\mathrm{d}\theta\\
&= 8\int_0^{1/\cos \theta} \frac{(1+1/\cos^2\theta)^{3/2}-1}{3}\,\mathrm{d}\theta\\
&= {?}
\end{align}$$
 A: Not an answer, but a start. Might not work.
In polar coordinates, it is $$4\int_{0}^{\sqrt2} f(r)\sqrt{1+r^2}\,dr$$ where $f(r)=\pi r/2$ for $0\leq r\leq 1.$
When $1<r\leq \sqrt2,$ we need $f(r)$ to be the length of the arc of radius $r$ around $(0,0)$ contained in the square $[0,1]^2.$
This is $$f(r)=r\left(\frac\pi2-2\arctan\left(\sqrt{r^2-1}\right)\right).$$
So the integral is:
$$4\frac{\pi}2\int_0^1 r\sqrt{1+r^2}\,dr+4\int_1^{\sqrt2}r\left(\frac\pi2-2\arctan\left(\sqrt{r^2-1}\right)\right)\sqrt{1+r^2}\,dr$$
The first integral is easy. Don't know what to do about the second, however.
When $x\leq 1,$ $2\arctan x =\arctan\frac{2x}{1-x^2}$ and $\frac{\pi}2-\arctan(y)=\arctan(1/y).$
So $$\frac{\pi}2-2\arctan\left(\sqrt{r^2-1}\right)=\arctan\left(\frac{2-r^2}{2\sqrt{ r^2-1}}\right)$$
If we set $r=\sec\theta,$ then $$\arctan\left(\frac{2-r^2}{2\sqrt{ r^2-1}}\right)=\frac\pi2-2\theta.$$
So the integral becomes:
$$\int_{0}^{\pi/4}\left(\frac\pi2-2\theta\right)\sec^2\theta\tan(\theta)\sqrt{1+\sec^2\theta}\,d\theta$$
A: Here is another approach using the Leibniz Integral Rule. This process gets ugly, so I will skip some steps.
For this solution, I will assume
$$I = \sqrt{3}+2\operatorname{arcsinh}\left(\frac{1}{\sqrt{2}}\right)+\int_{-1}^{1}\left(x^{2}+1\right)\operatorname{arccsch}\left(\sqrt{1+x^{2}}\right)dx$$
is correct and evaluate
$$J:= \int_{-1}^{1}\left(x^{2}+1\right)\operatorname{arccsch}\left(\sqrt{1+x^{2}}\right)dx.$$
(Solution) Letting $x \to \tan{(x)}$, we get
$$J = \int_{-\frac{\pi}{4}}^{\frac{\pi}{4}}\sec\left(x\right)^{4}\operatorname{arccsch}\left(\sec\left(x\right)\right)dx = 2\int_{0}^{\frac{\pi}{4}}\frac{\operatorname{arcsinh}\left(\cos\left(x\right)\right)}{\cos\left(x\right)^{4}}dx.$$
Let
$$J(a) = 2\int_{0}^{\frac{\pi}{4}}\frac{\operatorname{arcsinh}\left(a\cos\left(x\right)\right)}{\cos\left(x\right)^{4}}dx.$$
Differentiating with respect to $a$ yields
$$J'(a) = 2\int_{0}^{\frac{\pi}{4}}\frac{\sec^{3}\left(x\right)}{\sqrt{a^{2}\cos^{2}\left(x\right)+1}}dx.$$
I used an integral calculator to evaluate the antiderivative of the integrand. The integral becomes
$$
\eqalign{
& 2\left[\dfrac{\left(1-a^2\right)\ln\left(\left|\sqrt{\tan^2\left(x\right)+a^2+1}+\tan\left(x\right)\right|\right)+\tan\left(x\right)\sqrt{\tan^2\left(x\right)+a^2+1}}{2}\right]_{0}^{\frac{\pi}{4}} \cr
=& \sqrt{2+a^{2}}+\left(1-a^{2}\right)\left(\ln\left(1+\sqrt{2+a^{2}}\right)-\ln\left(\sqrt{1+a^{2}}\right)\right).
}
$$
We integrate both sides with respect to $a$ and over the interval $[0,1]$. Using the same integral calculator again, we get
$$
\eqalign{
\int_{0}^{1}J'\left(a\right)da &= \int_{0}^{1}\left(\sqrt{2+a^{2}}+\left(1-a^{2}\right)\left(\ln\left(1+\sqrt{2+a^{2}}\right)-\ln\left(\sqrt{1+a^{2}}\right)\right)\right)da \cr
&= -\dfrac{8\arcsin\left(\frac{\sqrt{2}\sqrt{3}+\sqrt{2}}{4}\right)-10\ln\left(2\sqrt{3}+2\right)-4\ln\left(\sqrt{3}+1\right)-6\operatorname{arsinh}\left(\frac{1}{\sqrt{2}}\right)+17\ln\left(2\right)-2{\pi}-2\sqrt{3}}{6}.
}
$$
We also know from the Fundamental Theorem of Calculus that
$$\int_{0}^{1}J'\left(a\right)da = J(1) - J(0) = J.$$
By transitivity, we get
$$J = -\frac{8\arcsin\left(\frac{\sqrt{2}\sqrt{3}+\sqrt{2}}{4}\right)-10\ln\left(2\sqrt{3}+2\right)-4\ln\left(\sqrt{3}+1\right)-6\operatorname{arsinh}\left(\frac{1}{\sqrt{2}}\right)+17\ln\left(2\right)-2\pi-2\sqrt{3}}{6}.$$
Putting everything together and doing a lot of simplifying, we conclude the integral $I$ to be
$$\int_{-1}^{1}\int_{-1}^{1}\sqrt{1+x^{2}+y^{2}}dydx = -\frac{2}{9}\left(\pi+12\ln\left(2\right)-6\sqrt{3}-24\ln\left(1+\sqrt{3}\right)\right).$$
If there is a doubt about the computations, I have added my Desmos link here. Let me know if there are any questions.
