How to ListContourPlot an eigenvalue spectrum without jumping?

Question

MWE:

Block[{dim = 4, grain = 250, matrix, domain, spectrum, locus},

matrix = DiagonalMatrix[Table[Sqrt[(x + 1/2 - (1 + Exp[-(i - 1)])^(-1))^2 +
   (y (1 + Exp[-2 (i - 1)])^(-1))^2], {i, 1, dim}]];
domain = Flatten[Table[{x, y}, {x, -1, 1, 2./(grain - 1)}, {y, -1, 1,2./(grain - 1)}], 1];
spectrum = Table[Append[σ, Sort@Re@Eigenvalues[ReplaceAll[matrix, {x -> σ[[1]], 
    y -> σ[[2]]}]]], {σ, domain}];
locus[energy_] := Table[Join[spectrum[[σ, 1 ;;2]], Take[SortBy[spectrum[[σ, 3]], 
    Abs[# - energy] &],1]], {σ, 1, Length@domain}]; 

ListContourPlot[locus[0.1], Contours -> {0.1}, InterpolationOrder -> 4] /. 
    _Polygon -> Sequence[]

]   

This code block functions roughly as desired, but the quality of the plot is very bad and cannot be satisfactorily improved by increasing grain and InterpolationOrder. Clearly, in this example matrix, the contours should be smooth ellipses and we need to recover that to a higher fidelity.

The matrix here is a toy, the only real constraint on the matrix I need to preserve is that it must be Hermitian and depend on two independent, real variables.

My solution attempt:

1. discretizes the domain into a mesh of (x,y) points
2. collects into array spectrum a dim-dimensional eigenvalue-list (of real numbers) at every point in the domain
3. Defines a function locus which reduces the spectrum array to include only a single eigenvalue (rather than dim-many) at every point in the domain.
4. Perform a ListContourPlot of the locus

However, I believe that the Take of the "correct" eigenvalue in step 3 eigenvalue at every point in the domain is likely the source of the problem.

In my MWE, I Take the eigenvalue which is absolutely nearest to the function argument energy. However, I can envision several pathological function landscapes which break this selection and, indeed, the plots illustrate this problem.

The goal of the code is to plot a level set (contour) of all points where there is any eigenvalue of the matrix

Answer 1

To see what is happening here, first plot a blow-up of the 2D region to show clearly the contours.

ListContourPlot[locus[0.1], Contours -> {0.1}, InterpolationOrder -> 1, 
    ContourShading -> None, PlotRange -> {{-0.2, 0.6}, {-0.4, 0.4}}]

The plot appears to consist of four ellipses plus two ragged curves. Note that InterpolationOrder -> 1 is used instead of 4, as in the question, to avoid any possible oscillations in the interpolation, which could occur near discontinuities. Note also, that ContourShading -> None is used to display contours only. The question used /. _Polygon -> Sequence[] instead, but the latter may be less transparent to some readers.

To clarify the nature of the ragged curves, next plot a slice through the 2D region, for instance, y == 0.

ListPlot[Delete[#, 2] & /@ Select[locus[0.1], Abs[#[[2]]] < .01 &], 
    PlotRange -> {{-0.2, 0.6}, {0, .2}}]

With this plot aligned with the first, it is evident that the two ragged curves are associated with discontinuities in locus[0.1]. J.M. suggested in a comment a related 1-D problem in which the solution was to sort eigenvalues to produce continuous arrays. The following seems more straightforward in 2D, however: Plot contours for each of the four sets of eigenvalues and superimpose them. This works well, because Eigenvalues sorts the eigenvalues it produces by size, and those eigenvalues vary smoothly with {x, y} except where the eigenvalues intersect.

eig[n_] := Extract[#, {{1}, {2}, {3, n}}] & /@ spectrum    
Show @@ (ListContourPlot[eig[#], Contours -> {0.1}, InterpolationOrder -> 1, 
    ContourShading -> None, PlotRange -> {{-0.2, 0.6}, {-0.4, 0.4}}] & /@ Range[1, 4])

as desired

Answer 2

I want to give a big thanks to user bbgodfrey for the elegant solution to the problem posed.

However, for variety and extension, I'm posting this answer in addition. Basically, I modified his/her solution by

1. exchanging Append to Join in spectrum definition
2. exchanging Map,Extract to Part in what would be the definition of eig
3. forgoing definition of eig in favor of building that functionality into the ListContourPlot
4. exchanging Show, Apply to Show, Table in the ListContourPlot
5. including plotting qualities in the function

The majority of these changes are aesthetic, simply because I don't have a lot of experience coding with Map and Apply, so they aren't intuitive. The last point is a simple extension that is beyond the scope of the question I posed here but I include only for other interested readers.

Block[{dim = 4, grain = 250, matrix, domain, range, cuts, spectrum, colors},
matrix = DiagonalMatrix[Table[Sqrt[(x + 1/2 - (1 + Exp[-(i - 1)])^(-1))^2 + (y (1 + Exp[-2 (i - 1)])^(-1))^2], {i, 1, dim}]];
domain = Flatten[Table[{x, y}, {x, -1, 1, 2./(grain - 1)}, {y, -1, 1,2./(grain - 1)}], 1];  
range = Range[0, 0.15, 0.15/(cuts - 1)];  
spectrum = Table[Join[\[Sigma],Sort@Re@Eigenvalues[ReplaceAll[matrix, {x -> \[Sigma][[1]],y -> \[Sigma][[2]]}]]], {\[Sigma], domain}];

discreteConPlot[energy_, quality_] := Show@Table[
ListContourPlot[spectrum[[All, {1, 2, \[Sigma] + 2}]], 
 ContourStyle -> {quality[[1]], Thickness[quality[[2]]]}, 
 InterpolationOrder -> quality[[3]], ImageSize -> quality[[4]], 
 Contours -> {energy}, ContourShading -> None, 
 PlotRange -> {{-0.2, 0.6}, {-0.4, 0.4}}], {\[Sigma], 1, dim}];

colors = Table[Hue[\[Sigma], 1, 0.7], {\[Sigma], Range[0, 0.8, 0.8/(cuts -1)]}];

ListAnimate[
Table[discreteConPlot[range[[s]], {colors[[s]], 0.005, 1, 400}], {s,1, Length@range}]
, AnimationRunning -> False]

]  

This block produces cuts-many contours across the range of values chosen, with a hue-spread across Hue[0,1,0.7] to Hue[0.8,1,0.7] and other qualities that can be adjusted readily.

Several of the Table throughout this algorithm can likely be improved by Parallelization but the optimal choice may depend on the user's machine and the size of dim and grain, amongst others.

Example output, by exchanging ListAnimate here for Show.