plotting - How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters)

I've been working on this research problem for months with some success by no global/proper solution to the problem.

So I have square, Hermitian matrices of various sizes ($8,20,38$ dimensions) with functional dependence on two variables, call them x and y. I am interested in plotting all the eigenvalues as functions of the $2$ variables x and y.

Let's call M the $8 \times 8$ hermitian matrix.

Plot3D[ Eigenvalues[M], {x,-1,1}, {y,-1,1}]

gives $8$ sheets corresponding to the eigenvalues of M as functions of x and y.

However, there will be many "jumps" between the sheets which totally obscure the plots. Basically, I think that during the plotting of the eigenvalues, for certain regions of the parameter space, Mathematica organizes the $8$ values differently, so the different sheets swap places, so we get tons of random garbage vertical bits between sheets.

I have spent months trying to resolve this in different ways, mostly involving methods to force Mathematica into a consistent sorting method. Here are some examples of failed methods:

Plot3D[ Sort[Eigenvalues[M]], {x,-1,1}, {y,-1,1}]


Plot3D[ Sort[Eigenvalues[M], Abs[#1] <= Abs[#2]&], {x,-1,1}, {y,-1,1}]

Block[{set=Ordering[ Eigenvalues[ M/.x->0 /. y->0], All, Abs[#1] <= Abs[#2]&]},
  Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]

Block[{set=Ordering[ Ordering[Eigenvalues[ M/.x->0 /.y->0], All, 
        Abs[#1] <= Abs[#2]&]]},
  Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]


Block[{set = Ordering[Eigenvalues[ M /.x->0 /.y->0], All, Abs[#1] <= Abs[#2]&]},
  Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]

Block[{set = Ordering[ Ordering[ Eigenvalues[ M /.x->0 /.y->0], All,
       Abs[#1] <= Abs[#2]&]]},
  Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]

as well as other permutations of similar ideas and more. Some of these methods have worked somewhat, reducing the jumping or working in certain regimes of the parameter domain. But nothing works consistently.

Please help me understand the proper solution to this problem!

Answer

Thanks to everyone for the help. Here's my working solution:

======
First, the spirit of the underlying problem and the nature of the ultimate solution. Simply sorting by the function values (executed with the default one-argument Sort command) is enough to separate the functions and prevent the jumping. However, it does not give the user control over the functions and it is difficult to do anything further (e.g. remove several of the functions from the plot, apply colors in a predictable way, etc.). If you sort instead by the two argument Sort with a pure absolute value function as

Sort[(*function*),Abs[#1]<=Abs[#2]&] 

there is a pairwise degeneracy problem when two functions cross planes equidistant above and below the zero plane, and you get tons of pairwise swapping at those places and a garbage plot. The spirit of this attempt, however, would be to gain control over the smallest absolute value function bands, for example.

The solution that follows uses a different technique. We instead use a reference point (the origin) to compute a single list of eigenvalues, perform simple operations there to compute which bands in the final plot will correspond to meeting a desired criteria (such as the smallest n bands in terms of absolute value) and then plot the functions sorted by a simple Sort[] with these indices
======

"lownumber" is the desired number of plotted eigenvalues, the smallest lownumber bands in terms of absolute value. lowb and lowt are the lower and upper bound indices for the plotted functions

lowb=Min[Drop[Ordering[

     Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All, 
     Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];  

lowt=Max[Drop[Ordering[
     Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All, 
     Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];

After this, there are two possible solutions. I have tested them both to be of the same speed to within about 3%, for matrices up to size 38. For much larger matrices, I am not sure, one might be preferable. Here they are:

Plot3D[
 Table[Sort[Eigenvalues[hamiltonian]][[s]], {s, lowb, lowt}] 

or

Plot3D[
 Table[RankedMin[Eigenvalues[hamiltonian],s], {s, lowb, lowt}]