z1 = x + I y;
z2 = z + I w;
n = 2;
eq = ComplexExpand[ReIm[z1^n + z2^n] - {1, 0}];
wsol = Solve[eq[[2]] == 0, w, Reals];
eq3D = Expand[z^2 eq[[1]] /. wsol[[1]]];
R = 5;
ContourPlot3D[eq3D == 0, {x, -R, R}, {y, -R, R}, {z, -R, R},
 RegionFunction -> Function[{x, y, z}, x^2 + y^2 + z^2 <= 25],
 PlotPoints -> 60,
 Mesh -> 25
 ]