A random sphere containing the center of the unit cube

As the probability is the same in all cubes, we can calculate it in the cube $[-1,1]^3$.

We can limit ourselves to the sixth of the cube where $z$ is positive and has the greatest absolute value of the three coordinates. Then the radius of the sphere is $1-z$, and the centre of the cube is in the sphere if $x^2+y^2+z^2\le(1-z)^2$. Thus the admissible area of $(x,y)$ is the intersection of the square $[-z,z]^2$ with the circle $x^2+y^2=1-2z$. A corner of the square lies on the circle if $3z^2=(1-z)^2$, that is, $z=\frac{\sqrt3-1}2$, and a midpoint of the square lies on the circle if $2z^2=(1-z)^2$, that is, $z=\sqrt2-1$.

Thus, for $0\le z\le\frac{\sqrt3-1}2$, the entire square lies within the circle, so the area is $4z^2$.

For $\frac{\sqrt3-1}2\le z\le\sqrt2-1$ the circle and square intersect. The four segments of the circle that extend beyond the square each have area $(1-2z)\arccos\frac z{\sqrt{1-2z}}-z\sqrt{1-2z-z^2}$, so the area is $\pi(1-2z)-4\left((1-2z)\arccos\frac z{\sqrt{1-2z}}-z\sqrt{1-2z-z^2}\right)$.

For $\sqrt2-1\le z\le\frac12$ the entire circle lies within the square, so the area is $\pi(1-2z)$; and for $z\gt\frac12$ the area is $0$.

Thus the desired probability is

$$ \frac68\left(\int_0^\frac{\sqrt3-1}24z^2\mathrm dz+\int_\frac{\sqrt3-1}2^{\sqrt2-1}\left((1-2z)\left(\pi-4\arccos\frac z{\sqrt{1-2z}}\right)+4z\sqrt{1-2z-z^2}\right)\mathrm dz+\int_{\sqrt2-1}^\frac12\pi(1-2z)\mathrm dz\right)\;. $$

The first and last integral evaluate to $\frac43\left(\frac{\sqrt3-1}2\right)^3=\sqrt3-\frac53$ and $\frac\pi4\left(1-2\left(\sqrt2-1\right)\right)^2=\pi\left(\frac{17}4-3\sqrt2\right)$, respectively. Wolfram|Alpha evaluates the indefinite form of the second integral to

$$ -\pi z^2+\pi z+4\sqrt{1-2z-z^2}\left(\frac{z^2}3+\frac z6-\frac56\right)+(6-z)\sqrt{1-2z-z^2}+\frac{15}2\arctan{\frac{1+z}{\sqrt{1-2z-z^2}}}+\frac12\arctan\frac{1-3z}{\sqrt{1-2z-z^2}}-4\arcsin\frac{1+z}{\sqrt2}+4(z-1)z\arccos\frac z{\sqrt{1-2z}} $$

but refuses to evaluate it with limits. Substituting the limits by hand yields

$$ -\pi\left(3-2\sqrt2\right)+\pi\left(\sqrt2-1\right)+\frac{15}2\cdot\frac\pi2-\frac12\cdot\frac\pi2-4\cdot\frac\pi2=\left(3\sqrt2-\frac52\right)\pi $$

at the upper limit and

$$ -\pi\left(1-\frac{\sqrt3}2\right)+\pi\cdot\frac{\sqrt3-1}2+\frac23-\sqrt3+\frac72\sqrt3-4+\frac{15}2\cdot\frac{5\pi}{12}+\frac12\left(-\frac\pi{12}\right)-4\cdot\frac{5\pi}{12}+4\cdot\frac{\sqrt3-3}2\cdot\frac{\sqrt3-1}2\cdot\frac\pi4=-\frac{10}3+\frac52\sqrt3+\frac{17}{12}\pi $$

at the lower limit, so the second integral evaluates to

$$ \frac{10}3-\frac52\sqrt3+\left(3\sqrt2-\frac{47}{12}\right)\pi\;. $$

Thus, the desired probability is

$$ \frac34\left(\sqrt3-\frac53+\frac{10}3-\frac52\sqrt3+\left(3\sqrt2-\frac{47}{12}\right)\pi+\pi\left(\frac{17}4-3\sqrt2\right)\right)\\=\boxed{\frac\pi4+\frac54-\frac98\sqrt3\approx0.086841}\;, $$

in agreement with Aaron’s calculation and simulation.

I'll try an approach using order statistics. Suppose the coordinates of the chosen random point are $(X_1, X_2, X_3)$. We can assume that the coordinates of the chosen point are all positive. (If they are not, we can reflect the point into the first octant.) Hence, $X_1, X_2, X_3$ are independent uniform random variables on $(0, 1)$.

Now, we define the order statistics $Y_1, Y_2, Y_3$ so that $Y_1$ is the smallest of the $X_i$ values, $Y_2$ is the middle value, and $Y_3$ is the largest value. Note that the $Y_i$ variables are neither uniform on $(0, 1)$ nor independent of one another.

There are two variables of interest: $R$, the radius of the sphere, and $D$, the distance from the chosen point to the origin. Note that $R = \min\{1 - X_1, 1 - X_2, 1 - X_3\} = 1 - Y_3 $ and that $D = \sqrt{X_1^2 + X_2^2 + X_3^2} = \sqrt{Y_1^2 + Y_2^2 + Y_3^2}$. The operative question is: what is $\mathbb P(D < R)$? That is, what is $$\mathbb P \left(\sqrt{Y_1^2 + Y_2^2 + Y_3^2} < 1 - Y_3 \right)?$$

First, a bit of algebra to clean this up: \begin{align*} \mathbb P \left(\sqrt{Y_1^2 + Y_2^2 + Y_3^2} < 1 - Y_3 \right) &= \mathbb P \left( Y_1^2 + Y_2^2 + Y_3^2 < (1 - Y_3)^2 \right) \\ &= \mathbb P \left(Y_1^2 + Y_2^2<1-2 Y_3 \right) \end{align*} We know the joint pdf of these order statistics to be $$f(y_1, y_2, y_3) = \begin{cases} 3!, & 0 < y_1 < y_2 < y_3 < 1 \\ 0, & \text{otherwise} \end{cases}$$ so we just need to integrate that density over the set $\{y_1^2 + y_2^2 < 1 - 2 y_3\}$ in the cube $[0, 1]^3$. Note that this requires in particular that $y_3 \leq 1/2$. I claim that this triple integral can be expressed as \begin{align*} \int_0^{1/2} \int_0^{\min\{y_3, \sqrt{1 - 2 y_3}\}} \int_0^{\min\{y_2, \sqrt{1 - 2y_3 - y_2^2}\}} 6 \, \textrm d y_1 \, \textrm d y_2 \, \textrm d y_3. \end{align*}

I don't know the analytic value of that integral (I haven't tried very hard yet), but Wolfram Alpha estimates it to be $\fbox{0.0868}$.

When I do a really long probability computation like this, I always assume I screwed up somewhere and verify my work with a Monte Carlo simulation. Here's that work in R:

spherecube <- function(){
  center <- runif(3, min=-1, max=1)
  radius <- min(abs(1 - center), abs(1 + center))
  sum(center * center) < radius^2
}

mean(replicate(100000, spherecube()))
# 0.08674

Aside from being an interesting problem, it's a great advertisement for the power of Monte Carlo simulations!