# Series expand root of quartic polynomial containing many real parameters and FullSimplify return Indeterminate

Using Mathematica I can get the Eigenvalues of the following matrix (using Quartics->True):

{{(DE1^2 + EC*DE1 + t1^2)*(1/DE1), ((2*I*DE1*xi + t2*t1)*DE2 + DE1*t1*t2)*(1/(2*DE1*DE2)), t1^2*(1/DE1), -t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), -J*(1/(2*Exp[I*phi])), 0, 0,  0}, {((-2*I*DE1*xi + t2*t1)*DE2 + DE1*t1*t2)*(1/(2*DE1*DE2)), (DE2^2 + DE2*EC + t2^2)*(1/DE2), t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), -t2^2*(1/DE2), 0, -J*(1/(2*Exp[I*phi])), 0, 0}, {t1^2*(1/DE1), t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), (DE1^2 + EC*DE1 + t1^2)*(1/DE1), ((2*I*DE1*xi - t2*t1)*DE2 - DE1*t1*t2)*(1/(2*DE1*DE2)), 0, 0, -J*(1/(2*Exp[I*phi])), 0}, {-t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), -t2^2*(1/DE2), ((-2*I*DE1*xi - t2*t1)*DE2 - DE1*t1*t2)*(1/(2*DE1*DE2)), (DE2^2 + DE2*EC + t2^2)*(1/DE2), 0, 0, 0, -J*(1/(2*Exp[I*phi]))}, {-(1/2)*(Exp[I*phi]*J), 0, 0, 0, (DE1^2 + EC*DE1 + t1^2)*(1/DE1), ((-2*I*DE1*xi + t2*t1)*DE2 + DE1*t1*t2)*(1/(2*DE1*DE2)), t1^2*(1/DE1), -t1* t2*(DE1 + DE2)*(1/(2*DE1*DE2))}, {0, -(1/2)*(Exp[I*phi]*J), 0, 0, ((2*I*DE1*xi + t2*t1)*DE2 + DE1*t1*t2)*(1/(2*DE1*DE2)), (DE2^2 + DE2*EC + t2^2)*(1/DE2), t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), -t2^2*(1/DE2)}, {0, 0, -(1/2)*(Exp[I*phi]*J), 0, t1^2*(1/DE1), t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), (DE1^2 + EC*DE1 + t1^2)*(1/DE1), ((-2*I*DE1*xi - t2*t1)*DE2 - DE1*t1*t2)*(1/(2*DE1*DE2))}, {0, 0, 0, -(1/2)*(Exp[I*phi]*J), -t1*t2*(DE1 + DE2)*(1/(2*DE1*DE2)), -t2^2*(1/ DE2), ((2*I*DE1*xi - t2*t1)*DE2 - DE1*t1*t2)*(1/(2*DE1*DE2)), (DE2^2 + DE2*EC + t2^2)*(1/DE2)}}


If I take now the first eigenvalue and want to series expand it in t1,t2,xi,J I would normally simply replace all four variables by #*epsilon and series expand in epsilon (and use Normal). However, if I try to FullSimplify the resulting expression it returns Indeterminate for the second order.

I figured since calculating eigenvalues is not an uncommon task maybe someone has encountered the same issue and found a way to handle these things. For a similar matrix I wanted to do the same and ended up doing lots of by hand work. However, I want to study many such matrices so something more automatic or less matrix-specific would be great.

Here are my assumptions on the parameters:

Element[EA1 | EA2 | EC | DE1 | DE2 | t1 | t2 | phi | xi | J, Reals]
&& phi > -Pi/200 && phi < Pi/200


EDIT: I meant FullSimplify, not Simplify, I edited the question.

Edit2: I forgot to mention that I use the Quartics->True option. Here is a MWE with a slightly "faster" matrix, but with the same issues:

H = {{EA - c, 0, t1, t2, I*s, 0, 0, 0}, {0, EA - c, t1, -t2, 0, -I*s,
0, 0}, {t1, t1, EC, 0, 0, 0, 0, I*xi}, {t2, -t2, 0, EC, 0,
0, -I*xi, 0}, {-I*s, 0, 0, 0, EA + c, 0, t1, t2}, {0, I*s, 0, 0,
0, EA + c, t1, -t2}, {0, 0, 0, I*xi, t1, t1, EC, 0}, {0, 0, -I*xi,
0, t2, -t2, 0, EC}} /. {s -> J*Sin[phi]/2, c -> J*Cos[phi]/2}
Assuming[Element[t1 | t2 | EA | EC | xi | J | phi, Reals],
FullSimplify[#] & /@
CoefficientList[
Normal@Series[
Eigenvalues[H, 1, Quartics -> True] /. {t1 -> t*t1, t2 -> t*t2,
xi -> xi*t, J -> J*t}, {t, 0, 4}], t]]


Designating the matrix in the Question as m, I ran the following:

fsteig= Eigenvalues[m, 1];
Normal[Series[fsteig /. {t1 -> a t1, t2 -> a t2, xi -> a xi, J -> a J}, {a, 0, 3}]];
FullSimplify[%, Element[EA1 | EA2 | EC | DE1 | DE2 | t1 | t2 | phi | xi | J, Reals]]


Although I received a warning message,

Root::sbr: Because of branch cuts, the series may represent a different root of ...


I received no error messages. Concatenating

&& phi > -Pi/200 && phi < Pi/200 && EA1 != EA2 && EA1 != EC && EA2 != EC


to the second argument of FullSimplify (Equivalent to your definition of FullSimplify2, which I also tried) does not produce an error message either. Perhaps, some other part of your notebook has corrupted your results.

Update

The addition of Quartics->True triggered the error described by the OP, both in his original problem and the faster running one. I analyzed the latter as follows:

scdeig = Eigenvalues[H, 1, Quartics -> True][[1]];
lst = CoefficientList[Normal@Series[
scdeig /. {t1 -> t*t1, t2 -> t*t2, xi -> xi*t, J -> J*t}, {t, 0, 4}], t];


and then applying

FullSimplify[#, Element[t1 | t2 | EA | EC | xi | J | phi, Reals]] &


To each element of lst in turn. The first element yielded Piecewise[{{EA, EA < EC}}, EC], and the second yielded 0. The third, however, yielded two

FullSimplify::infd: Expression ... simplified to ComplexInfinity. >>


error messages, because it contains denominators that vanish. One would expect that putting all of lst[[3]] over a common denominator before applying FullSimplify would resolve this issue, but it does not.

lst3 = Together[lst[[3]]];


Then, applying FullSimplify to Numerator[lst3] yields

(* 6*Sqrt[3]*(EA - EC)^7*(J^2 - 4*xi^2 - (J^2 + 4*xi^2)/Sign[EA - EC]) *)


whereas applying FullSimplify to Denominator[lst3] yields 0.

Alternative Approach (for "fast" matrix)

The difficulties encountered above probably are associated with the many branch cuts occurring there. Instead, expand the CharacteristicPolynomial of H and its eigenvalues z:

cl = Simplify[CharacteristicPolynomial[H /. {t1 -> t*t1, t2 -> t*t2, xi -> xi*t, J -> J*t}, z[t]]];
cls = CoefficientList[Simplify[Normal[Series[cl, {t, 0, 6}]]], t];


Note that Simplify with no assumptions can be used, because we are dealing only with polynomials. We now solve cls == 0 term by term to obtain the expansion of z.

Solve[cls[[1]] == 0, z[0]]


yields z[0] -> {{z[0] -> EA}, {z[0] -> EC}}, each with four-fold degeneracy. Substituting EA, for instance, in place of z[0]

Simplify[cls /. z[0] -> EA]


next yields

{0, 0, 0, 0, 1/16 (EA - EC)^4 (J^2 - 4 z'[0]^2)^2, ...}


Hence, we have {{z'[0] -> J/2}, {z'[0] -> -J/2}}. Substituting either of these into cls further reduces it

Simplify[% /. z'[0] -> J/2]


to

{0, 0, 0, 0, 0, 0, 1/4 (EA - EC)^2 J^2 (2 t1^4 + 2 t2^4 - 2 (t1^2 - t2^2)^2 Cos[2 phi] -
4 EA t2^2 z''[0] + 4 EC t2^2z''[0] + EA^2 z''[0]^2 - 2 EA EC z''[0]^2 +
EC^2 z''[0]^2 + 4 t1^2 (3 t2^2 + (-EA + EC) z''[0]))}


and, finally,

Simplify[Solve[%[[7]] == 0, z''[0]]]


yields z''[0]. This process can be extended to higher order in a straightforward, although tedious, manner.

Undoubtedly, the same approach can be used for the "slow" matrix.

• Ah, I'm sorry, I used FullSimplify, otherwise the assumptions aren't used. I will edit the question – NOhs Jan 29 '15 at 21:27
• FullSimplify also produced no error with your assumptions. However, there are several ways to apply assumptions. Which are you using? – bbgodfrey Jan 29 '15 at 22:10
• FullSimplify2[x_] := FullSimplify[x, Element[EA1 | EA2 | EC | DE1 | DE2 | t1 | t2 | phi | xi | J, Reals] && phi > -Pi/200 && phi < Pi/200 && EA1 != EA2 && EA1 != EC && EA2 != EC], like this – NOhs Jan 30 '15 at 10:34
• I also forgot to mention, that I use Quartics->True – NOhs Feb 2 '15 at 8:46
• Please include in your Question the actual call to Eigenvalues that you used. For completeness, you might move your definition of FullSimplify2[x_] there as well and show how you called it. – bbgodfrey Feb 2 '15 at 13:36