I have found this identity by computing the same quantity in two different ways but I have not found a way to prove it explicitly. Taking $s \in \mathbb{R}^+$, one has $$ F(s) \equiv \int_0^\infty du_1 \int_0^\infty du_2 \ e^{-2 s (u_1+u_2+1)} \left((u_1+u_2+1) \operatorname{arctanh}\left(\sqrt{\frac{u_1}{u_1+1}} \sqrt{\frac{u_2}{u_2+1}}\right)\right)= \frac{K_0(s) K_2(s)-K_1(s){}^2}{16 s} $$ Do you have any suggestions about how this can be simplified?
Update 1
One can make the change of variable $r = \sqrt{u / (u + 1)}$, leading to the expression $$ F(s) = \int_0^1 dr_1 \int_0^1 dr_2 P(r_1) P(r_2) \left(\frac{1}{r_1 r_2}-r_1 r_2\right) \operatorname{arctanh}(r_1 r_2) $$ where we introduced the unnormalised probability distribution $$ P(r) = \frac{2 r^2 \exp (s) \exp \left(-\frac{2 s}{1-r^2}\right)}{\left(1-r^2\right)^3} $$ One can express the (unnormalised) even moments of this probability distribution as $$ R_k = \int_0^1 dr \; P(r) r^{2k} = e^{-s} \Gamma \left(k+\frac{3}{2}\right) U(k+\frac{3}{2},3,2 s)$$ where $U$ is the Tricomi's (confluent hypergeometric) function. From the power series of the $\operatorname{arctanh}$, we thus obtain the representation $$ F(s) = -\sum_{k=0}^\infty \frac{R_{k+1}^2 - R_{k}^2}{2k + 1} $$ Applying summation by parts, we can recast this into $$ F(s) = -\frac{K_1(s)^2}{16 s^2} + \sum_{k=0}^\infty R_k^2 \left(\frac{1}{2k+1} - \frac{1}{2k-1} \right) $$ Note that the expression to be proved can be equivalently rewritten as $$ F(s) = -\frac{1}{4} \pi e^{-2 s} s \left(U\Bigl(\frac{3}{2},3,2 s\Bigr)^2-U\Bigl(\frac{1}{2},1,2 s\Bigr) U\Bigl(\frac{5}{2},5,2 s\Bigr)\right) $$
Update 2
Another approach is to introduce the variables $p = u_1 u_2$ and $q = u_1 + u_2$. One obtains the alternative representation $$ F(s) = 2\int_0^\infty dq \int_0^{q^2/4} dp \frac{(q+1) e^{-2 (q+1) s} \operatorname{arctanh}\left(\sqrt{\frac{p}{p+q+1}}\right)}{\sqrt{q^2-4 p}} $$ Integrating by parts over the variable $p$, it can be recast into $$ F(s) = \frac{1}{2} \int_0^\infty dq \int_0^{q^2/4} dp \; (q+1) \sqrt{\frac{q^2-4 p}{p (p+q+1)}} e^{-2 (q+1) s} $$ The integral in $p$ can now be performed leading to $$ F(s) = \int_0^\infty dq \frac{1}{2} \sqrt{q+1} (q+2)^2 K\left(-\frac{q^2}{4 (q+1)}\right) e^{-2 (q+1) s}-2(q+1)^{3/2} e^{-2 (q+1) s} E\left(-\frac{q^2}{4 (q+1)}\right) $$ where $K(k)$ and $E(k)$ are the complete elliptic integrals.