Numerical derivative (ND) not working for orthonormalized eigenvectors of a matrix

Question

I have problem making ND work. Define relevant needs and vectors

Needs["NumericalCalculus`"]
d2 = 1./2. {-Sqrt[3.], 1.};
d3 = {0., -1.};
a1 = Sqrt[3.] {1., 0.};
a2 = Sqrt[3.]/2. {1., Sqrt[3.]};
a3 = -a1 + a2;
k = {kx, ky};

Define matrices

H[kx_, ky_, \[Alpha]_, \[Beta]_, \[Gamma]_] = {{0., 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3]), \[Alpha]*(1 + 
E^(I*k.a1) + E^(I*k.a2)), 0.}, {2*\[Gamma]*(Sin[k.a1] + 
Sin[k.a2] + Sin[k.a3]), 0., 
0., \[Beta]*(1 + E^(I*k.a1) + E^(I*k.a2))}, {\[Alpha]*(1 + 
   E^(-I*k.a1) + E^(-I*k.a2)), 0., 0., 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3])}, {0., \[Beta]* 
(1 +E^(-I*k.a1) + E^(-I*k.a2)), 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3]), 0.}};
Htbg[kx_, ky_] = 
2*(H[kx, ky, 1., -0.2, 0.2] + 0.3456 IdentityMatrix[4]);

Solve for normalized eigenstates. Transpose makes the column elements the eigenvalues. Sorting is based on the first column (the eigenvalues).

\[Psi][kx_?NumericQ, ky_?NumericQ] := \[Psi][kx, ky] = 
 Orthogonalize[Sort[Transpose[Eigensystem[Htbg[kx, ky]]]] 
[[All, 2]],Conjugate[#1].#2 &];

This gives wildly fluctuating phase values in k space which is not suitable for taking derivatives. Impose smooth gauge by making the 1st element of each eigenvectors real. That is, by multiplying a phase.

\[Psi]sg[kx_, ky_] := Table[Exp[-I*
  ArcTan[Re[\[Psi][kx, ky][[i, 1]]], 
   Im[\[Psi][kx, ky][[i, 1]]]]] \[Psi][kx, ky][[i]], {i, 1, 
4}];

As example, consider the [[1,1]] element. This gives a smooth plot

Plot[Re[\[Psi]sg[x, 1.][[1, 1]]], {x, -1, 1}]

But numerical derivative fails

ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, 0.]

I do not why this is happening and would appreciate if someone can help. Thank you.

P.S. I tried to composed a simpler example with stripped-down irrelevant details, but the resulting simple example worked just fine. I do not know why the above does not work.

Answer 1

Declaring:

[Psi]sg[kx_?NumericQ, ky_?NumericQ] := ...

will reduce the number of error messages you get, although not completely eliminate them.

As you mentioned, his works:

Plot[Re[\[Psi]sg[x, 1.][[1, 1]]], {x, -1, 1}]

If you now say:

ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, 0.]

You get a bunch of error messages of the form:

Part::partd: Part specification \[Psi]sg[kx,1.][[1,1]] is longer than depth of object.

Part::partw: Part 3 of \[Psi][kx,1.] does not exist.

This is because MMA evaluates arguments before feeding them to functions. In your case it evaluates: "[Psi]sg[kx, 1.][[1, 1]]", feeding "[Psi]sg[kx, 1.]" a symbolic argument, what makes it to fail.

If you define "[Psi]sg[kx_?NumericQ, ky_?NumericQ]=.." this will not evaluate and "[Psi]sg[kx,1.][[1,1]]" gives you only one error message, that it can not take part "[[1,1]]". It nevertheless then feeds the argument to "Re" and finally "ND" and gives a result.

This will not happen inside "Plot", because "Plot" is special, it does not evaluate its argument before giving numerical values to the variables. Therefore, we may plot the derivative without problems provided you define "[Psi]sg[kx_?NumericQ, ky_?NumericQ]=..":

Plot[ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, kx0], {kx0, -1, 1}]