# Solve computes discontinuous eigenvalues of parameter-dependent matrix

So I have a family of unitary matrices $$m(x,y)$$, which depend on two parameters $$x,y \in [0, 2 \pi)$$. Its eigenvalues should be continuous in $$(x,y)$$. Since $$m(x,y)$$ is a unitary matrix, its eigenvalues are of norm $$1$$, which is why I want to use Mathematica to plot their argument over the domain $$[0, 2 \pi)^2$$. For this I basically use Solve to get the roots of the characteristic polynomial of $$m(x,y)$$. The matrix $$m(x,y)$$ is a product of four $$3 \times 3$$ unitary matrices and is generated by the following code:

S12[x_, y_] = {{Exp[I y], 0., 0.}, {0., Exp[I (y - x)], 0.}, {0., 0., 1.}};
S21[x_, y_] = {{Exp[I (x - y)], 0., 0.}, {0., 1., 0.}, {0., 0., Exp[-I y]}};
CDFT = 1/(3^(1/2)) {{1, 1, 1}, {1, Exp[I 2 Pi/3], Exp[I 4 Pi/3]}, {1,Exp[I 4 Pi/3],
Exp[I 8 Pi/3]}};
m[x_, y_] = S12[x, y] . CDFT . S21[x, y] . CDFT;


I use the following code to compute its eigenvalues depending on $$x$$ and $$y$$:

charpoly[x_, y_] = CharacteristicPolynomial[m[x, y], k];
zeros[x_, y_] = k /. Solve[charpoly[x, y] == 0, k];
argzeros[x_, y_] = Arg[zeros[x, y]];
Plot3D[{argzeros[x, y][[1]], argzeros[x, y][[2]], argzeros[x, y][[3]]},
{x, 0, 2 Pi}, {y, 0, 2 Pi}, PlotLegends -> Automatic]


$m(x,y)$" />

As you can see, the bands have a discontinuity exactly when the phase of an eigenvalue has to transition from $$2 \pi$$ to $$0$$. Technically the two small blue parts should be green, the small green parts should be orange and the small orange parts should be blue. Is there any way to avoid this discontinuity, so that Solve "labels" the eigenvalues correctly and the colors in the plot align?

• Welcome to the Mathematica StackExchange! Please share a complete example (including $m(x,y)$) so that people can experiment and propose solutions. Apr 7, 2023 at 17:41
• The Arg Function is a non-continuous function. Lines of negative Eigenvalues produce steps by 2 pi in the graphics. m[x,y] is not given. The form of the jump lines as parabolas suggest a 2x2 matrix Apr 7, 2023 at 18:16

After a bit of manual fiddling with transformations, it works out: Starting with exact expressions (removing the machine-precision 0. and 1. symbols in the original expression),

S12[x_, y_] = {{Exp[I y], 0, 0}, {0, Exp[I (y - x)], 0}, {0, 0, 1}};
S21[x_, y_] = {{Exp[I (x - y)], 0, 0}, {0, 1, 0}, {0, 0, Exp[-I y]}};
CDFT = 1/(3^(1/2)) {{1, 1, 1}, {1, Exp[I 2 Pi/3], Exp[I 4 Pi/3]},
{1, Exp[I 4 Pi/3], Exp[I 8 Pi/3]}};
m[x_, y_] = S12[x, y] . CDFT . S21[x, y] . CDFT // FullSimplify;


Eigenvalues: using complex logarithms instead of the argument function, so that we can simplify more easily later on,

s[x_, y_] = -I*Log[Eigenvalues[m[x, y]]] // FullSimplify
(*    {-I Log[1/3 E^(-2 I (x + y)) Root[..., 1]],
-I Log[1/3 E^(-2 I (x + y)) Root[..., 2]],
-I Log[1/3 E^(-2 I (x + y)) Root[..., 3]]}    *)


Transform the Root objects manually to avoid branch cuts: Note that the transformation function $$q(x,y)$$ captures the prefactor of the root object inside of the logarithm and we end up with Log[Root[...]] expressions. This is the goal: we want to avoid branch cuts in both the Root and Log functions!

q[x_, y_] = 3 I E^(2 I (x + y));
t[x_, y_] = s[x, y] /.
Root[f_, i_] :> q[x, y] Root[Evaluate[f[q[x, y] #]] &, i] /.
Log[a_] :> Log[-I a] - Log[-I]
(*    {-I (I*π/2 + Log[Root[..., 1]]),
-I (I*π/2 + Log[Root[..., 2]]),
-I (I*π/2 + Log[Root[..., 3]])}    *)

Plot3D[Evaluate@t[x, y], {x, 0, 2π}, {y, 0, 2π}, PlotLegends -> Automatic, PlotPoints -> 100]


## some more explanations

We are looking for three values $$s_k(x,y)$$ such that the $$e^{i s_k(x,y)}$$ are eigenvalues of $$m(x,y)$$. We impose two continuity conditions:

1. No permutations hoppings: the ordering of the $$s_k(x,y)$$ shall stay consistent and not jump (permute) as a function of $$(x,y)$$, and
2. No branch cuts: the $$s_k(x,y)$$ shall not do any $$\pm2\pi$$ jumps, even though such jumps have no effect on the values of $$e^{i s_k(x,y)}$$.

For this purpose we first calculate the form of the $$s_k(x,y)$$ in the variable s[x,y]; but this result has both permutation hops and branch cuts.

To eliminate the permutation hops (point 1), we need to scale the Root objects. We note that for any root $$z$$ satisfying $$f(z)=0$$, we can define a scaled root $$\zeta=z/q$$ that satisfies $$f(q\zeta)=0$$. We thus define a scaled function $$g(\zeta)=f(q\zeta)$$, find the root $$\zeta$$ of $$g$$, and then compute $$z=q\zeta$$. Seems trivial, yes; but it helps greatly to reduce or even eliminate permutation jumps in the Root objects. With $$q(x,y)=3i e^{2i(x+y)}$$ we thus eliminate the permutation jumps.

But we still have $$\pm2\pi$$ branch cuts (point 2). To eliminate these, we shift the complex logarithm by an amount so that the branch cut happens in a place that none of the solutions ever reaches. In this case, it suffices to shift the argument by $$\pi/2$$ by multiplying all roots by $$-i$$ and then subtracting the shift: Log[a_] :> Log[-I a] - Log[-I].

• Thanks for the help, that is perfect! Could you please elaborate a bit on the purpose of $q$ and $t$? I am new to Mathematica, so I am having quite some trouble understanding what especially $t$ does, and why you use it. Apr 9, 2023 at 12:09
• @Andreas132 I've added some more explanations. This is more math than Mathematica. Apr 10, 2023 at 14:46