# Finding the conditions for degenerate solutions to a cubic equation

I have a cubic equation (corresponding to the band structure of a physical system) given by w in the code below.

M = {{s - w, ab, ac}, {ab, s - w, bc}, {ac, bc, s - w}};
w = w /. Solve[Det[M] == 0, w];
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}];


With 0<=θ<=Pi and 0<=ϕ<=2Pi tuneable parameters. I would like to develop a code that finds the coordinates {x,y} at which the three different solutions (the three bands) touch for a given {θ,ϕ}. For example, for {θ,ϕ}={0,0} the band structure is given below.

I already know the touching point here is {x,y}={Pi/2,Pi/2} from previous studies. However I would really like to be able to find the touching points (if there are any) for any arbitrary {θ, ϕ}.

There are three basic ways I can think to approach this are (though none of them have worked for me):

1) Use the known properties of cubic equations. When the discriminant of the cubic equation equals zero then the equation has a multiple root and all of its roots are real. Thus some sort of code such as

c = ab^2 + ac^2 + bc^2;
d = 2 ab ac bc;
discriminant = Chop[4 c^3 + 27 d^2];
{θ, ϕ} = {0, 0};
Solve[discriminant == 0 && 0 <= θ <= Pi &&  0 <= ϕ <= 2 Pi, {x, y}]


However this process spits out some horrendously complicated expression that I can't seem to simplify.

2) Using the explicit expression for the three solutions w to find where they touch.

{θ, ϕ} = {0, 0};
Solve[w[[1]] == w[[2]] == w[[3]], {x, y}]


3) Use some kind of numerical approach with FindRoot. Either with the discriminant==0 condition or explicit touching condition. Start with a known solution, for example {x,y}={Pi/2,Pi/2} when {θ,ϕ}={0,0}. Then slowly varying theta and phi using the previous found touching point as the starting point to look for a solution and then in this way map out the touching point for the whole parameter space of theta and phi.

Any suggestions greatly appreciated

• I appreciate this is quite a hefty question and not overly general, however I'm feeling very stuck and so often a small nudge from SE has sent me in the right direction, so thanks for any help! – Tom Jul 3 '15 at 11:48
• Discriminant[] is built-in. ;) – J. M. will be back soon Jul 3 '15 at 14:28

Further edited to simplify results

It is not difficult to show that the three eigenvalues, w, of M are equal if and only if M is diagonal; i.e., ab = ac = bc = 0.

f[e_] := 2 Norm[e]^-3 (1 - 3 Sin[θ]^2 Cos[ϕ - ArcTan[e[[1]], e[[2]]]]^2) Cos[{x, y}.e]
abeval = f[{1, 0}]
(* 2 Cos[x] (1 - 3 Cos[ϕ]^2 Sin[θ]^2) *)
aceval = f[{0, 1}]
(* 2 Cos[y] (1 - 3 Sin[θ]^2 Sin[ϕ]^2) *)
bceval = Simplify[TrigExpand[f[{1, 1}] + f[{-1, 1}]]]
(* (Cos[x] Cos[y] (1 + 3 Cos[2 θ]) + 6 Sin[x] Sin[y] Sin[θ]^2 Sin[2 ϕ])/(2 Sqrt[2]) *)


Solutions can be obtained in an orderly manner as follows. For {Cos[x] -> 0, Cos[y] -> 0}, abeval and aceval are zero, and bceval reduces to

Sin[θ] Sin[2 ϕ] == 0


For {Sin[x] -> 0, Cos[y] -> 0}, aceval and bceval are zero, and abeval reduces to

1 - 3 Cos[ϕ]^2 Sin[θ]^2 == 0


Finally, for {Cos[x] -> 0, Sin[y] -> 0}, abeval and bceval are zero, and aceval reduces to

(1 - 3 Sin[θ]^2 Sin[ϕ]^2) == 0


These curves are plotted below and indicate where solutions are located in θ and ϕ.

ContourPlot[{Sin[θ] Sin[2 ϕ] == 0,
(1 - 3 Cos[ϕ]^2 Sin[θ]^2) == 0,
(1 - 3 Sin[θ]^2 Sin[ϕ]^2) == 0},
{θ, 0, Pi + .01}, {ϕ, 0, 2 Pi + .01}, Frame -> False, Axes -> True,
Ticks -> {{0, Pi}, {0, Pi, 2 Pi}}, AxesLabel -> {"θ", "ϕ"},
ContourStyle -> Directive[Black, Thick],
AxesStyle -> Directive[Black, Bold, Thick, 12]]


Other values of x and y give rise to discrete points in θ and ϕ, located at the intersections of curves in this plot.

Band Structure Plot

With the information just presented, generating band structure plots is straightforward. To do so, evaluate s

seval = f[{2, 0}] + f[{0, 2}]
(* 1/4 Cos[2 x] (1 - 3 Cos[ϕ]^2 Sin[θ]^2) + 1/4 Cos[2 y] (1 - 3 Sin[θ]^2 Sin[ϕ]^2) *)


define

mplt[θ0_, ϕ0_] := Module[{meval =
Det[M] /. {ab -> abeval, ac -> aceval, bc -> bceval,
s -> seval} /. {θ -> θ0, ϕ -> ϕ0} // Simplify},
Plot3D[Evaluate[w /. Solve[meval == 0, w]], {x, -Pi, Pi}, {y, -Pi, Pi},
AxesLabel -> {"x", "y"}, AxesStyle -> Directive[Black, Bold, 12],
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}, Automatic},
ImageSize -> Large, ViewPoint -> {1.4, -3.0, .35}]]


and invoke it with any pair {θ0, ϕ0} from the plot above or from the three equations that the plot represents. For instance

mplt[Pi/2, ArcSin[Sqrt[1/3]]]


or

mplt[Pi/2, ArcCos[Sqrt[1/3]]]


Plots for values corresponding to interior curve intersections from the {θ0, ϕ0} plot are in a sense degenerate.

 mplt[ArcSin[Sqrt[1/3]], Pi]


mplt[0, 0], of course, gives the plot in the Question.

At the request of the OP in a comment below, here is a more detailed argument that M must be diagonal for its three eigenvalues to be equal. To begin, its Determinant must be proportional to [-(w - w0)^3, where w0 is the three-fold-repeated eigenvalue.

eq0 = Expand[-(w - w0)^3]
(* -w^3 + 3 w^2 w0 - 3 w w0^2 + w0^3 *)


The Determinant of M actually is

eq = Collect[Det[M], w, Simplify]
(* 2 ab ac bc - ab^2 s - s (ac^2 + bc^2 - s^2) + (ab^2 + ac^2 + bc^2 - 3 s^2) w +
3 s w^2 - w^3 *)


Equating the two yields

Collect[eq - eq0, w, Simplify]
(* 2 ab ac bc - ab^2 s - ac^2 s - bc^2 s + s^3 + 3 w^2 (s - w0) - w0^3 +
w (ab^2 + ac^2 + bc^2 - 3 s^2 + 3 w0^2) *)


and the coefficient of every power of w must vanish. Thus, w0 must be equal to s. With that substitution,

Collect[% /. w0 -> s, {w, s}, Simplify]
(* 2 ab ac bc + (-ab^2 - ac^2 - bc^2) s + (ab^2 + ac^2 + bc^2) w *)


The coefficient of w in the last expression vanishes if and only if all three of ab, ac, and bc vanish, in other word, if M is diagonal.

• Thank you! Interestingly enough I had actually already generated your 1st plot for another reason. As a first step I plotted all the angles at which the band structure touches for the high symmetry points in the Brillouin zone. However I had no feel for where (or if there would be) any gapless solutions at other points in momentum space. Thanks to your insight it appears that actually these (approximately) linear touching points only occur specifically at those points of high symmetry. For me it does mean unfortunately these points are quite unstable with regards to (theta,phi) but... – Tom Jul 6 '15 at 9:43
• ...interesting nonetheless :) Could I ask you (for the sake of future readers and myself) just to edit and show proof (or link to proof) of your initial starting point. I understand that a diagonal matrix's eigenvalues are the entries on the diagonal. However are we sure there are no other situations in which the three eigenvalues are equal? I'm pretty sure but want to be sure sure. – Tom Jul 6 '15 at 9:46
• @bbcgodfrey can this be easily extended to find the values of $\theta$ and $\phi$ where 2 or more eigenvalues are degenerate? I realise the initial question called for all 3 eigenvalues to be degenerate, but now I have interest in the case where two or more are degenerate. – Tom Jul 10 '15 at 15:10
• @Tom I believe so, although the process would be somewhat more complicated. I suggest you ask a new question. By the way, if my answer above effectively answers your original question, please accept it. – bbgodfrey Jul 10 '15 at 15:46
• Sure I can do that. I had every intention of reaccepting, I just thought that accepted answers got less attention and responses. – Tom Jul 10 '15 at 19:11

This is an incomplete answer, but we will be able to show that there is no solution for most values of θ and ϕ. We will also be able to draw a plot of the regions of interest that you should check further to find solutions, should they exist.

M = {{s - w, ab, ac}, {ab, s - w, bc}, {ac, bc, s - w}};
{d, c, b, a} = CoefficientList[Det[M], w];
disc = Discriminant[Det@M, w] // FullSimplify;
disc0 = b^2 - 3 a c;


Now, if we want three Real (equal) solutions both disc and disc0 must vanish simultaneously. As both discriminants are always positive, we can explore when the minimum of their sum is zero:

f[e_] := 2 Norm[e]^-3 (1 - 3 Sin[theta]^2 Cos[phi - 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}];

Plot3D[FindMinValue[disc + disc0 /. {theta -> t, phi -> p},
{{x, 0, Pi}, {y, 0, Pi}}],
{t, 0, Pi/2}, {p, 0, Pi/2}, PlotRange -> {0, .1}, ClippingStyle -> None]


So there you have the few regions to explore (they are lines, BTW)

Starting with bbgodfrey's excellent suggestion to solve ab == ac == bc == 0, we can obtain a fairly compact list of all of the solutions. If we Reduce the equations with conditions on the variables we get a complicated result, so it's easier to Reduce first and apply conditions after:

Reduce the equations and throw out some obviously inconsistent results:

sols = List @@ (
LogicalExpand@FullSimplify@Reduce[Thread[Simplify@{ab, ac, bc} == 0]] /.
{___ && Cos[x_] == 0 && ___ && Sin[x_] == 0 && ___ -> False,
___ && Sin[x_] == 0 && ___ && Cos[x_] == 0 && ___ -> False}
)


The last four solutions are

sols[[-4 ;;]]


We can show that these are not self-consistent by eliminating θ and ϕ like so:

(Or @@ FullSimplify[# /. Solve[#[[;; 2]], {θ, ϕ}], C[_] ∈ Integers] &) /@ %

{False, False, False, False}


The remainder of the solutions are

(sols = Drop[sols, -4]) // TableForm


We can Solve to find some particular solutions and throw away duplicates with

sols1 = Union@Simplify@Flatten[Solve /@ (sols /. _Unequal -> Sequence[]) /. C[_] -> 0, 1];


Then pick out solutions where the fixed parameters are between 0 and π:

sols2 = Pick[sols1, 0 == # & /@ Count[v_ /; ! 0 <= v < π] /@ ({x, y, θ, ϕ} /. sols1)]


Plot the third and fourth solutions for {x -> 0, y -> π/2} in the θ and ϕ plane:

 ContourPlot[Evaluate[Equal @@@ sols2[[3 ;; 4, 3]]], {θ, 0, π}, {ϕ, -π/2, π/2}, MaxRecursion -> 3]


And plot the energies as a function of θ and ϕ:

Plot3D[Evaluate[Re[w /. {x -> 0, y -> π/2}]], {θ, 0, π}, {ϕ, -π/2, π/2}, MaxRecursion -> 5]