On the integral $\int_{(0,1)^n}\frac{\prod\sin\theta_k}{\sum\sin\theta_k}d\mu$

This question is a followup to MSE2732980, where it is shown that $$ \mathcal{J}_2=\iint_{(0,1)^2}\frac{dx\,dy}{\sqrt{1-x^2}+\sqrt{1-y^2}}=\frac{\pi(4-\pi)}{4}.\tag{$n=2$}$$

It comes natural to wonder if there is a simple closed form, in terms of elementary functions and polylogarithms, for the generalized integral $$\begin{eqnarray*}\mathcal{J}_n&=&\int_{(0,1)^n}\frac{d\mu}{\sum_{k=1}^{n}\sqrt{1-x_k^2}}\\&=&\int_{0}^{+\infty}\left(\int_{0}^{+\infty}\sin(\theta)e^{-x\sin\theta}\,d\theta\right)^n\,dx\\&=&\left(\frac{\pi}{2}\right)^n\int_{0}^{+\infty}\left(L_{-1}(x)-I_1(x)\right)^n\,dx. \end{eqnarray*}$$

Remark 1: from the second integral representation for $\mathcal{J}_n$ we have that $\{\mathcal{J}_n\}_{n\geq 1}$ is a log-convex sequence, and it might be related to some probabilistic problem, like the probability that a random walk in $\mathbb{Z}^n$ returns to the origin.
Remark 2: $$ \int_{0}^{1}\frac{du}{A+\sqrt{1-u^2}} = \frac{\pi}{2}-\frac{\text{arctanh}\sqrt{1-A^2}}{\sqrt{1-A^2}} A=\frac{\pi}{2}-\sum_{k\geq 0}\frac{A(1-A^2)^k}{2k+1}$$ hence the computation of $\mathcal{J}_2$ is equivalent to the computation of $\phantom{}_2 F_1\left(\tfrac{1}{2},\tfrac{1}{2};\tfrac{3}{2};1\right)$.
Remark 3: by the tangent half-angle substitution $$\mathcal{J}_3=8\iiint_{(0,1)^3}\left(\sum_{k=1}^{3}\frac{1-u_k^2}{1+u_k^2}\right)^{-1}\prod_{k=1}^{3}\frac{1-u_k^2}{(1+u_k^2)^2}\,d\mu$$ and $\mathcal{J}_3$ might be related to Watson's triple integrals.

This is too long for a comment but I believe some insight has been achieved so I decided to post it . Maybe it will inspire somebody to further ideas.

Let us define:

\begin{eqnarray} Q_+(u_1,u_2)&=& (1+u_1)^2+(1+u_2)^2-1+u_1 u_2(2 u_1+2 u_2+u_1 u_2)\\ Q_-(u_1,u_2)&=& -(1+u_1)^2-(1+u_2)^2+1+u_1 u_2(2 u_1+2 u_2-u_1 u_2) \end{eqnarray} Then we have: \begin{eqnarray} {\mathcal J}_3 &=& 2^5 \int\limits_{[0,1]^3} \prod\limits_{\xi=1}^3\left(\frac{u_\xi}{1+u_\xi^2}\right) \cdot \frac{1}{\sum\limits_{\xi=1}^3 u_\xi \prod\limits_{\eta\neq \xi} (1+u_\eta^2)} \prod\limits_{\xi=1}^3 d u_\xi\\ &=&\frac{\pi}{2} + 2^6 \int\limits_{[0,1]^2} \prod\limits_{\xi=1}^2 \left(\frac{u_\xi}{(1+u_\xi^2)^2}\right) (u_1+u_2)(1+u_1 u_2) \cdot \frac{\arctan(\frac{(1-u_2)^2(1+u_1^2)-2 u_1(1+u_2^2)}{\sqrt{Q_+(u_1,u_2)Q_-(u_1,u_2)}})}{\sqrt{Q_+(u_1,u_2)Q_-(u_1,u_2)}} du_1 du_2 \end{eqnarray} In the first line we used the tangent half-angle substitution, ie.e $u_\xi = \tan(\theta_\xi/2)$ for $\xi=1,2,3$ and in the second line we carried out the integral over $u_3$ by decomposing the integrand into partial fractions and integrating term by term and then simplifying the result using inverse trigonometric identities.

Unfortunately now we are not sure how to proceed. It would make sense to substitute for the $Q_\pm$ quantities however we don't know how to change variables .