Solution (source code)

= Solution

A <differential form> of degree $k$ is integrated over an oriented $k$-dimensional <smooth manifold> by integrating its coefficient in orientation-preserving coordinates. For an oriented parametrization $F:U\to M$, this means integrating the <pullback> $F^*\alpha$ over $U$; a <partition of unity> combines charts. The <change of variables formula> makes the result independent of the charts. Reversing the <orientation of a smooth manifold> reverses the integral. The <Generalized Stokes theorem> says, for a compact oriented $k$-dimensional <smooth manifold> with boundary and a smooth $(k-1)$-form $\beta$,
$$
\int_M d\beta=\int_{\partial M}\beta,
$$
where the right-hand side includes the boundary <pullback>, and the boundary has the <outward-normal-first boundary orientation>. Compact support suffices on a noncompact manifold. This unifies the fundamental theorem of calculus, circulation and flux identities.

For example, take the unit disk $D$ oriented by $dx\wedge dy$ and the <differential one-form> $\beta=\tfrac12(x\,dy-y\,dx)$. Its <exterior derivative> is $d\beta=dx\wedge dy$, so the area integral is $\pi$. Parametrizing the positively oriented circle by $(x,y)=(\cos t,\sin t)$ gives the <pullback> $\beta=\tfrac12dt$ and the boundary integral $\int_0^{2\pi}\tfrac12dt=\pi$, directly illustrating <Generalized Stokes theorem>.

For the distance calculation, write $u=|\mathbf r_1|$, $v=|\mathbf r_2|$, $x=u+v$ and $y=u-v$; these are scalar radii, not position vectors. The <triangle inequality> and its reverse give $r\leq u+v=x$ and $|y|=|u-v|\leq r$. Also,
$$
x+|y|=2\max(u,v)\leq2R.
$$
Thus \b[$r\leq x\leq2R-|y|$ and $|y|\leq r$]. These conditions, with $u,v\geq0$, are also sufficient for a triangle with side lengths $u,v,r$. Endpoints correspond to collinear configurations.

Assume the two positions are <independent random variables>, each having constant volume <probability density> inside the radius-$R$ <ball>. Independence is needed: uniform marginal distributions alone do not determine the distance distribution. Put $C=3/(4\pi R^3)$ and $\mathbf s=\mathbf r_2-\mathbf r_1$. The <change of variables> $(\mathbf r_1,\mathbf r_2)\mapsto(\mathbf r_1,\mathbf s)$ has unit <Jacobian determinant>, so the joint probability volume <differential form> is the product of the two normalized volume forms. In spherical coordinates for $\mathbf r_1$ and for $\mathbf s$ relative to its axis, it is
$$
C^2u^2\sin\theta_1\,du\wedge d\theta_1\wedge d\phi_1\wedge r^2\sin\theta\,dr\wedge d\theta\wedge d\chi.
$$
The polar coordinates fail on axes and at zero radii, which are sets of zero volume and do not affect the <probability>. Here $\theta$ is the angle from $\mathbf r_1$ to $\mathbf s$. The <law of cosines> becomes
$$
v^2=u^2+r^2+2ur\cos\theta.
$$
Differentiating and taking the <wedge product> with $du\wedge dr$ eliminates the terms involving $du$ and $dr$:
$$
du\wedge dr\wedge(v\,dv)=-ur\sin\theta\,du\wedge dr\wedge d\theta.
$$
Consequently the positive integration density transforms by $|\sin\theta\,d\theta|=v\,dv/(ur)$ at fixed $u,r$. The minus sign is accounted for by reversal of limits: $v$ decreases as $\theta$ increases. The radial variable retained is $r$, while the polar angle is replaced by $v$; equivalently one can start with the angle between the two position vectors and replace that angle by $r$. The wording of the printed hint conflates these two coordinate choices, but its volume form gives the stated transformation directly.

Integrating $\sin\theta_1\,d\theta_1\,d\phi_1$ over the first orientation gives $4\pi$, and integrating $d\chi$ gives $2\pi$. The <probability density function> of $r$ is therefore
$$
f_R(r)=8\pi^2C^2r\iint_{D_r}uv\,du\,dv=\frac{9r}{2R^6}\iint_{D_r}uv\,du\,dv,
$$
where $D_r$ has $0\leq u,v\leq R$ and $|u-v|\leq r\leq u+v$. There are no such configurations for $r>2R$.

To evaluate this integral explicitly, the <Jacobian determinant> of $(u,v)=((x+y)/2,(x-y)/2)$ has absolute value $1/2$. Hence $uv\,du\,dv=(x^2-y^2)\,dx\,dy/8$ as a positive density. For $0\leq r\leq2R$, set $a=\min(r,2R-r)$. The domain is $-a\leq y\leq a$, $r\leq x\leq2R-|y|$, and evenness in $y$ gives
$$
I(r)=\iint_{D_r}uv\,du\,dv=\frac14\int_0^a\left[\frac{(2R-y)^3-r^3}{3}-y^2(2R-y-r)\right]dy.
$$
An antiderivative, zero at $y=0$, yields
$$
I(r)=\frac1{12}\left[(8R^3-r^3)a-6R^2a^2+ra^3+\frac{a^4}{2}\right].
$$
Substituting $a=r$ for $r\leq R$ and $a=2R-r$ for $r\geq R$ gives the same polynomial in both intervals:
$$
I(r)=\frac{2R^3r}{3}-\frac{R^2r^2}{2}+\frac{r^4}{24}.
$$
Thus the <distance between two uniform points in a three-dimensional ball> has the final density
$$
\boxed{dP=\left(\frac{3r^2}{R^3}-\frac{9r^3}{4R^4}+\frac{3r^5}{16R^6}\right)dr,\qquad 0\leq r\leq2R.}
$$
It is zero outside this interval. Its nonnegativity also follows from $f_R(r)=3r^2(4R+r)(2R-r)^2/(16R^6)$, and direct integration gives $\int_0^{2R}f_R(r)dr=1$. The derivation uses the full six-dimensional probability volume <differential form>, rather than treating the three scalar distances as independent.