# Remove the blank of branch cuts

H = {{I*(x + I*y), 1}, {1, -I*(x + I*y)}};
{x1, x2} = Eigenvalues[H];

Plot3D[{Im[x1], Im[x2]}, {x, 0, 2}, {y, -1, 1}, PlotStyle -> {LightRed, LightRed}, Mesh -> 0] There is a white blank on the surface because of the branch cut. How to remove this "white line" to make the two surfaces as a whole part?

The desired imaginary parts are represented by v in algebraic formulation below. Put u before v in the variables for GroebnerBasis to get the first polynomial to define v in terms of x and y. Reverse (v before u, both before x, y) to get the real part u in terms of x and y.

branches = u + I v - Eigenvalues[H];
surf = Times @@ branches // Simplify;
parts = ReIm[surf] // ComplexExpand;
gb = GroebnerBasis[parts, {u, v, x, y}];
ContourPlot3D[
First[gb] == 0 // Evaluate, {x, 0, 2}, {y, -1, 1}, {v, -2, 2},
Mesh -> None, ContourStyle -> LightRed, MaxRecursion -> 6,
PlotPoints -> 20] Another way that happens to remove the branch-cut jump (the use of polar coordinates is similar to Henrik's, but it's the combination with ComplexExpand that moves the branch cut to the boundary of the plot):

imEVs = FullSimplify[
ComplexExpand[Im@Eigenvalues[H], TargetFunctions -> {Re, Im}],
x > 0 && -1 < y < 1];
ParametricPlot3D[
{x, 0, 2}, {y, -1, 1}, PlotStyle -> LightRed, Mesh -> None,
Exclusions -> None, BoxRatios -> {1, 1, 1}
] Using polar coordinates helps a little

plot = With[{
opts =
Sequence[PlotStyle -> {LightRed, LightRed}, Mesh -> 0,
PlotPoints -> {73, 20},
BoundaryStyle -> None,
PlotRange -> {{0, 2}, {-1, 1}, {-2, 2}}
],
ϵ = 1. 10^-6,
aa = {x, y, Im[x1]} /. {x -> 1 + r Cos[t], y -> r Sin[t]},
bb = {x, y, Im[x2]} /. {x -> 1 + r Cos[t], y -> r Sin[t]}
},
Show[
ParametricPlot3D[aa, {t, ϵ, 2 Pi - ϵ}, {r, 0, 2},
opts],
ParametricPlot3D[bb, {t, ϵ, 2 Pi - ϵ}, {r, 0, 2},
opts]
]
] However, the slit can still be seen from different view points. But one can try and glue the two GraphicsComplexes within plot together in order to fix this.

Somewhat inelegant, but works?

Plot3D[{Max[Im[x1], Im[x2]], Min[Im[x1], Im[x2]]}, {x, 0, 2}, {y, -1, 1},
PlotStyle -> {LightRed, LightRed}, Mesh -> 0, ExclusionsStyle -> LightRed] Just for fun note the eigenvalue expression is: $$\lambda^2=1-(x^2+2ixy-y^2)=1-(x+iy)^2$$. Letting $$z=x+iy$$ we have the expression $$\lambda=\sqrt{1-z^2}$$

So basically, $$\lambda$$ is just the square root function offset to x=1. And one easy way to plot this is:

myf[z_] := Sqrt[1 - z^2];
p1 = ParametricPlot3D[{Re[z], Im[z], Im@myf[z]} /.
z -> 1 + r Exp[I t], {r, 0, 1}, {t, 0, 2 \[Pi]}]
p2 = ParametricPlot3D[{Re[z], Im[z], -Im@myf[z]} /.
z -> 1 + r Exp[I t], {r, 0, 1}, {t, 0, 2 \[Pi]}]
Show[{p1, p2}] 