Find all degenerate eigenvalues of a cubic equation

I have an equation that is cubic in w. The three solutions correspond to bands in a bandstructure, and are a function of wavevector {x, y}. Additionally it is parameterised by 0 <= θ <= π and 0 <= φ <= 2 π. I need a procedure that takes θ and φ as inputs, and can reliably calculate if there are any points at all in wavevector space where either two or all of the bands touch (i.e. are there any values of {x, y} which give rise to a degenerate solution w). Ideally the procedure would then be able to calculate the coordinates in momentum space of these touching points. The code to generate the cubic equation is given below.

M = {{s, ab, ac}, {ab, s, bc}, {ac, bc, s}};
f[e_] := 2 Norm[e]^-3 (1 - 3 Sin[θ]^2 Cos[φ - ArcTan[e[[1]], e[[2]]]]^2)Cos[{x, y}.e];
ab = f[{1, 0}];
ac = f[{0, 1}];
bc = f[{1, 1}] + f[{-1, 1}];
s = f[{2, 0}] + f[{0, 2}];


A slight variation on this question has already been solved here by @bbgodfrey and others and can serve as a more detailed reference. However that question called for a solution when all three bands touch, as opposed to two or more, which is the current question.

One of my ultimate goals is to be able to generate a plot in {θ, φ} space showing what values give rise to a bandstructure with touching points. As an example, in previous work a procedure for a different physical system generated the plot below, where purple is a gapless bandstructure and cream is a gapped bandstructure:

• @bbgodfrey thanks for suggesting I create a new question. Here it is. – Tom Jul 14 '15 at 10:00
• Any reason not to use Eigenvalues[] here? – J. M. is away Jul 14 '15 at 10:08
• It always spits out something in the form of Root[] which I didn't really know how to handle, so I did it this way. – Tom Jul 14 '15 at 10:44
• That's its default behavior; if you want radicals in your output, use ToRadicals[]. – J. M. is away Jul 14 '15 at 10:46
• Thanks. Have modified question to incorporate this. – Tom Jul 14 '15 at 10:55

For a numerical approximation you may try something like:

M = {{s, ab, ac}, {ab, s, bc}, {ac, bc, s}};
disc = Discriminant[CharacteristicPolynomial[M, x], x] // FullSimplify
f[e_] := 2 Norm[e]^-3 (1 - 3 Sin[θ]^2 Cos[φ - ArcTan[e[[1]], e[[2]]]]^2) Cos[{x, y}.e];
ab = f[{1, 0}];
ac = f[{0, 1}];
bc = f[{1, 1}] + f[{-1, 1}];
s = f[{2, 0}] + f[{0, 2}];

Quiet@RegionPlot[With[{θ = t, φ = p},
FindMinValue[disc, {x, y}, AccuracyGoal -> Infinity]] < 10^-6,
{t, 0, Pi}, {p, 0, 2 Pi}, PlotPoints -> 60]


You should take care of the numeric subtleties. Please note:

Quiet@RegionPlot[
FindMinValue[disc /. {θ -> t, φ -> p}, {{x, 0, Pi}, {y, 0, Pi}}] < .001,
{t, 0, Pi}, {p, 0, 2 Pi}, PlotPoints -> 30]


• Oooo interesting thank you very much. I'll have a play with this. – Tom Jul 14 '15 at 20:39