The two spheres have equations
x² + y² + z² = r²
(x - 1)² + y² + z² = 1 - r²
By combining the equations, eliminating y and z and solving for x, we get
(x - 1)² - x² = (1 - r²) - r²
-2x + 1 = 1 - 2r²
x = r²
which is to say the two spheres meet in the plane x = r².
If we solve the two equations above for x, we can get two equations for the surfaces x = f(y, z) bounding the region of interest:
x² + y² + z² = r² ⇒ x = ± √(r² - y² - z²)
(x - 1)² + y² + z² = 1 - r² ⇒ x = 1 ± √(1 - r² - y² - z²)
There are actually four surfaces here, but two of these will not touch each other. For the sphere centered at (0, 0, 0), we want the "upper" hemisphere (in the region x > 0), and for the sphere at (1, 0, 0), we want the "lower" hemisphere (in x < 1). So our choice boils down to
x = √(r² - y² - z²)
x = 1 - √(1 - r² - y² - z²)
When x = r², the two equations reduce to the equation of circle,
(r²)² + y² + z² = r² ⇒ y² + z² = r² - r⁴
Solving for y gives an upper and lower bound for y in terms of z :
y² + z² = r² - r⁴ ⇒ y = ± √(r² - r⁴ - z²)
If y = 0, then we can find constant upper and lower bounds for z :
z² = r² - r⁴ ⇒ z = ± √(r² - r⁴) = ± r √(1 - r²)
If R is the region of interest, then we can describe it in Cartesian coordinates by the set
[tex]R = \left\{ (x, y, z) : 1 - \sqrt{1-r^2-y^2-z^2} \le x \le \sqrt{r^2 - y^2 - z^2} \right. \\\\ \text{ and } -\sqrt{r^2 - r^4 - z^2} \le y \le \sqrt{r^2 - r^4 - z^2} \\\\ \left. \text{ and } -r \sqrt{1-r^2} \le z \le r \sqrt{1-r^2} \right\}[/tex]
and so the volume we want is
[tex]V(r) = \displaystyle \iiint_R dV[/tex]
[tex]V(r) = \displaystyle \int_{-r\sqrt{1-r^2}}^{r\sqrt{1-r^2}} \int_{-\sqrt{r^2-r^4-z^2}}^{\sqrt{r^2-r^4-z^2}} \int_{1-\sqrt{1-r^2-y^2-z^2}}^{\sqrt{r^2-y^2-z^2}} dx \, dy \, dz[/tex]
To make computation easier, convert to cylindrical or spherical coordinates. Cylindrical is the better choice, so I'll do that. Set
x(ξ, ρ, θ) = ξ … … … (ξ = "xi")
y(ξ, ρ, θ) = ρ cos(θ) … … … (ρ = "rho")
z(ξ, ρ, θ) = ρ sin(θ)
The Jacobian for this transformation is
[tex]J = \begin{bmatrix}x_\xi & x_\rho & x_\theta \\ y_\xi & y_\rho & y_\theta \\ z_\xi & z_\rho & z_\theta\end{bmatrix} = \begin{bmatrix}1&0&0\\0&\sin(\theta)&\rho\cos(\theta) \\ 0&\cos(\theta)&-\rho\sin(\theta)\end{bmatrix}[/tex]
with |det(J)| = |ρ|. Then the volume element is
dV = dx dy dz = |ρ| dξ dρ dθ
Under this transformation, the bounding surfaces' equations become
ξ = √(r² - ρ² cos²(θ) - ρ² sin²(θ)) = √(r² - ρ²)
ξ = 1 - √(1 - r² - ρ² cos²(θ) - ρ² sin²(θ)) = 1 - √(1 - r² - ρ²)
ρ represents the distance in the plane x = r² from the center of the circular intersection of the two spheres, (r², 0, 0), to the edge of this circle. The maximum distance is the radius of this circle, √(r² - r⁴). (ρ is thus a positive number, so |ρ| = ρ.)
θ is the angle made by a vector pointing from this center to the circle's edge, running from 0 to 2π, relative to some arbitrary line in the plane x = r².
Putting everything together, we can describe R in cylindrical coordinates by the set
[tex]R = \left\{ (\xi, \rho, \theta) : 1 - \sqrt{1-r^2-\rho^2} \le \xi \le \sqrt{r^2-\rho^2} \right. \\\\ \text{ and } 0 \le \rho \le \sqrt{r^2 - r^4} \\\\ \left. \text{ and } 0 \le \theta \le 2\pi \right\}[/tex]
The volume is then
[tex]\displaystyle V(r) = \iiint_R dV[/tex]
[tex]\displaystyle V(r) = \int_0^{2\pi} \int_0^{\sqrt{r^2-r^4}} \int_{1-\sqrt{1-r^2-\rho^2}}^{\sqrt{r^2-\rho^2}} \rho \, d\xi \, d\rho \, d\theta[/tex]
The integral with respect to θ is just a constant:
[tex]\displaystyle V(r) = 2\pi \int_0^{\sqrt{r^2-r^4}} \int_{1-\sqrt{1-r^2-\rho^2}}^{\sqrt{r^2-\rho^2}} \rho \, d\xi \, d\rho[/tex]
Integrate with respect to ξ :
[tex]\displaystyle V(r) = 2\pi \int_0^{\sqrt{r^2-r^4}} \rho \left(\sqrt{r^2-\rho^2} - 1 + \sqrt{1-r^2-\rho^2}\right) \, d\rho[/tex]
Integrate with respect to ρ; easily done with substitutions:
[tex]\displaystyle \int \rho \sqrt{r^2 - \rho^2} \, d\rho = -\frac12 \int (-2) \rho \sqrt{r^2 - \rho^2} \, d\rho \\ = -\frac12 \int \sqrt{r^2-\rho^2} \, d(r^2-\rho^2) \\ = -\frac12 \cdot \frac23 (r^2-\rho^2)^{3/2} \\ = -\frac13 (r^2-\rho^2)^{3/2}[/tex]
[tex]\displaystyle \int \rho \sqrt{1 - r^2 - \rho^2} \, d\rho \\ = -\frac12 \int (-2) \rho \sqrt{1-r^2 - \rho^2} \, d\rho \\ = -\frac12 \int \sqrt{1-r^2-\rho^2} \, d(1-r^2-\rho^2) \\ = -\frac12 \cdot \frac23 (1-r^2-\rho^2)^{3/2} \\ = -\frac13 (1-r^2-\rho^2)^{3/2}[/tex]
which lead us to
[tex]\displaystyle V(r) = -2\pi \left(\left(\frac{r^6}3 + \frac{r^2-r^4}2 + \frac{(1-r^2)^3}3\right) - \left(\frac{r^3}3 + \frac{(1-r^2)^{3/2}}3\right)\right)[/tex]
[tex]\displaystyle V(r) = \boxed{-2\pi \left(\frac{r^4}2 - \frac{r^3}3 - \frac{r^2}2 + \frac13 - \frac13(1-r^2)^{3/2}\right)}[/tex]
I've included below a plot of V(r) showing confirming its value is positive for 0 < r < 1. The volume is maximized at r = (√2/3 - 5/12) π ≈ 0.172 with a maximum value of 1/√2 ≈ 0.707.
