Double Bubble

Manipulate[
 With[{d = Sqrt[1 + r2^2 - 2 r2 Cos[\[Pi]/3]], R = 1/(1/r2 - 1)^2},
  With[{y2 = Sqrt[(4 d^2 - (d^2 - r2^2 + 1)^2)/(4 d^2)],
    x2 = (d^2 - r2^2 + 1)/(2 d)}, Show[{
     ContourPlot3D[x^2 + y^2 + z^2 == 1,
      {x, -1.2, 3.}, {y, -2.1, 2.1}, {z, -2.1, 2.1},
      RegionFunction ->
       Function[{x, y, z}, (x - d)^2 + y^2 + z^2 >= r2^2],
      ColorFunction ->
       Function[{x, y, z, f},
        Directive[Opacity[0.4], ColorData["LakeColors"][0.5 + z/2]]],
      Mesh -> 2, Boxed -> False, Axes -> None
      ],
     ContourPlot3D[(x - d)^2 + y^2 + z^2 == r2^2,
      {x, d - 1, d + 1}, {y, -1, 1}, {z, -1, 1},
      RegionFunction -> Function[{x, y, z},
        x^2 + y^2 + z^2 >= 1
        ],
      ColorFunction ->
       Function[{x, y, z, f},
        Directive[Opacity[0.4], ColorData["Aquamarine"][0.5 + y/2]]],
      Mesh -> 2],
     ContourPlot3D[(x - x2 - Sqrt[R - y2^2])^2 + y^2 + z^2 == R,
      {x, x2 + Sqrt[R - y2^2] - R^2,
       x2 + Sqrt[R - y2^2] + R^2}, {y, -1, 1}, {z, -1, 1},
      RegionFunction -> Function[{x, y, z},
        x^2 + y^2 + z^2 <= 1 && (x - d)^2 + y^2  + z^2 <= r2^2
        ],
      ColorFunction ->
       Function[{x, y, z, f}, Directive[Opacity[1.0], Gray]], Mesh -> 2
      ]
     }, ImageSize -> 500]]], {{r2, .53,
   "\!\(\*SubscriptBox[\(r\), \(B\)]\)"}, .25, .95}]