Consider $C_0, v \in \mathbb{R}^n$ with $\|v\| = 1$ and $\epsilon > 0$, then $C_1 = C_0 + \epsilon \cdot v$ and $C_2 = C_0 - \epsilon\cdot v$.
What is the smallest volume ellipsoid containing $\mathcal{B}(C_1, R) \cup \mathcal{B}(C_2,R)$ for some $R> 0$.
I guess its centre should be in $C_0$ so one looks for $(x-C_0)^T\cdot P \cdot (x-C_0) \leq 1$. But how to proceed?
Should one consider the two centers as focal points ?
Turn the problem into a 2D problem, by defining the x-axis along the two spheres center, and the y-axis being arbitrary and the origin in the midpoint between the spheres.
Then the solution would be defined by the following sketch
Given the sphere radius $r$ at a distance $c$ from the origin, find the ellipse semi-major and semi-minor axis as
$$ \begin{aligned} a & = c + r & & \text{semi-major} \\ b & = \sqrt{r (c+r)} & & \text{semi-minor} \\ f & = \sqrt{c\,a} & & \text{focal point} \\ \epsilon & = \sqrt{c/a} & & \text{eccentricity} \\ \end{aligned} \tag{1}$$
Then the ellipsoid has the same semi-minor axis in z direction as it has in the y direction, $b$ since the problem is symmetric about the x-axis.
The development of the above equations is as follows.
Given any ellipse defined by $a$ and $b$, and a circle of radius $r$, you can fit ther circle inside the ellipse if the center of the circle is located at
$$ c = \sqrt{ \left(a^2-b^2\right) \left(1 - \left( \tfrac{r}{b} \right)^2 \right)} \tag{2}$$
which obviously only produces a result if $b<a$ and $r<b$. This describes the family of curves that fit the spheres.
Now since $c$ is fixed, and $a$ and $b$ vary, I solved the above for $a=b \sqrt{ \frac{b^2+c^2-r^2}{b^2-r^2} }$ and then set $\tfrac{{\rm d}a}{{\rm d}b} = 0$ to get the solution
$$b = \sqrt{r ( c+r)} \tag{3}$$
This describes the optimal solution which minimizes the semi-major axis $a$.
here is another example with rounded numbers and $a=20+6$ and $b=\sqrt{6×26}$
Just some hints: I think that without loss of generality, you can consider $C_0=0$ and $v=e_1$ (the first vector of any orthonormal basis). Via a roto-translation (in particular, an isometry, conserving volumes) will recover the general case.
Moreover note that by radial symmetry of the ellipsoid, the right sphere is included if and only if the left sphere is included...so only 1 sphere is needed, I think.
Then I would write it as a semidefinite optimization problem ($P\succ 0$ the variable minimizing the volume of $\{x^\top Px \leq 1\}$, constraint that one sphere $\{(x-c_1)^\top I(x-c_1)\leq R\}$ is included).
;tldr
Without loss of generality we can assume the circles have Radius R and are centered on the $x$ axis at $(\pm p, 0)$, $a$ is semi-major axis, $b$ is semi-minor axis and $\sqrt{a^2-b^2}=\pm c$ is where the foci are on the x-axis. $A$ is area.
$$b^2 = \frac{4R^2+p^2\sqrt{1+8R^2/p^2}-p^2}{4}\approx 2R^2$$
$$a^2= \frac{5p^2+4R^2+3p^2\sqrt{1+8R^2/p^2}}{4} \approx 2(p^2+2R^2 ) $$
$$c^2=\frac{6p^2+2p^2\sqrt{1+8R^2/p^2}}{4} \approx 2(p^2+R^2)$$
$$A=\pi ab \approx 2\pi R\sqrt{p^2+2R^2}$$
Let circle be $(x-p)^2+y^2=R^2$ and the ellipse be $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$.
Solve for $y$ in terms of $x$. Requiring the discriminant of the resulting quadratic in $x$ to be $0$ establishes a condition for the circle being tangent to the ellipse:
$$\frac{p^2}{c^2}+\frac{R^2}{b^2}=1$$
Where $a^2=b^2+c^2$. $a$ is the semi-major axis, $b$ is the semi-minor, and $c$ is the distance of the focus from the center.
Minimizing the squared area of the ellipse will minimize the Area of the ellipse so :
$$A^2 = \pi^2 a^2b^2$$
The tangency condition suggests we parmeterise with:
$$c^2=p^2\sec^2\theta$$ $$b^2=R^2\csc^2 \theta$$
So :
$$A^2/(\pi^2R^2)=(p^2\sec^2\theta + R^2\csc^2\theta)csc^2\theta$$
So:
$0= (2p^2\sec^2\theta \tan\theta - 2R^2\csc^2\theta \cot\theta)\csc^2\theta-2(p^2\sec^2\theta + R^2 \csc^2\theta)\csc^2\theta\cot\theta$
$(p^2\sec^2\theta + R^2 \csc^2\theta) \cot \theta= (p^2\sec^2\theta \tan \theta - R^2\csc^2\theta\cot \theta)$
$2R^2\csc^2\theta \cot\theta = p^2\sec^2\theta(\tan \theta - \cot \theta)$
$2R^2 \cot \theta = p^2 \tan^2\theta (\tan \theta - \cot \theta)$
$2R^2=p^2\tan^2 \theta (\tan^2\theta -1)$
$p^2\tan^4 \theta - p^2\tan^2\theta -2R^2=0$
$$\tan^2\theta = \frac{p + \sqrt{p^2+8R^2}}{2p}$$
Where we only take the plus sign since $\tan ^2 \theta \ge 0$ and it's assumed the radius is non-zero.
$\sec^2\theta = \frac{3p + \sqrt{p^2+8R^2}}{2p}$
$\cos^2\theta = \frac{2p}{3p + \sqrt{p^2+8R^2}}$
$\sin ^2 \theta = \frac{p +\sqrt{p^2+8R^2}}{3p + \sqrt{p^2+8R^2}}$
$\csc^2\theta = \frac{3p + \sqrt{p^2+8R^2}}{p +\sqrt{p^2+8R^2}}=1+ \frac{2p}{p +\sqrt{p^2+8R^2}}$
$$c^2 = \frac{6p^2 + 2p^2\sqrt{1+8R^2/p^2}}{4}$$
$$b^2=\frac{3R^2 + R^2\sqrt{1+8R^2/p^2}}{1 +\sqrt{1+8R^2/p^2}}=R^2+ \frac{2R^2}{1 +\sqrt{1+8R^2/p^2}}=\frac{R^2(\sqrt{1+8R^2/p^2}+3)}{\sqrt{1+8R^2/p^2}+1}$$
$$b^2 = \frac{4R^2+p^2\sqrt{1+8R^2/p^2}-p^2}{4}\approx 2R^2$$
$$a^2= \frac{5p^2+4R^2+3p^2\sqrt{1+8R^2/p^2}}{4} \approx 2(p^2+2R^2 ) $$
$A\approx 2\pi R\sqrt{p^2+2R^2} $
To minimize $a$:
$$a^2=(p^2\sec^2\theta + R^2\csc^2\theta) $$
$p^2 \tan^4 \theta = R^2 $
$\tan^2\theta = R/p$
$\sec^2\theta = (p+R)/p$
$\csc^2\theta = (p+R)/R$
$$a=p+R$$
$$b^2= R(p+R)=aR$$
$$c^2=p(p+R)=ap$$
After converting the problem to $2D$, the smallest ellipse is one that has the same curvature at its common vertex with the circle. The curvature of the ellipse at the vertex is given by
$ \kappa = \dfrac{b^2}{a} $
Now from the description of the problem,
$ a = \epsilon + R $
Therefore, we want
$ \dfrac{1}{R} = \dfrac{b^2}{\epsilon + R} $
Hence
$ b =\sqrt{\dfrac{\epsilon + R}{R}} $
The ellipsoid will be an ellipsoid of revolution about the vector $v$
Its equation is
$ (r - C)^T Q (r - C) = 1 $
where $C$ is the given center, and $r = [x, y, z]^T $
and
$Q = V D V^T$
with $D = \begin{bmatrix} \dfrac{1}{b^2} && 0 && 0 \\ 0 && \dfrac{1}{b^2} && 0 \\ 0 && 0 && \dfrac{1}{a^2} \end{bmatrix} $
And
$V = [ u_1 , u_2, v ] $
where $u_1$ is any vector perpendicular to the given vector $v$, and $u_2 = v \times u_1 $.