Defining a function that outputs a matrix, and later finding its eigenvalues

Question

I am trying to do the following:

I have a simple 2x2 matrix that depends on three parameters (physically -- momentum coordinates kx, ky, kz). Then I want to replace each of these parameters by a matrix block of size "cut x cut" (physically -- Peirels substitutaion with a cut-off). The results is a "2*cut x 2*cut" matrix. I want to study eigenvalues of this matrix.

Let me show how I implement this, step by step, in my code. The problem is that in the end I get wrong results, and I don't understand why.

First, I define the small matrix that depends on three parameters. For example (although in general I also have powers of kx,ky,kz):

H[kx_,ky_,kz_]:=kx*PauliMatrix[1]+ky*PauliMatrix[2]+kz*PauliMatrix[3];

Next, I know I will replace parameters kx, ky, kz (and potentially also their powers) by some matrix blocks. I define these blocks as follows:

a[cut_] := Module[ 
   {c = cut, t},
   t = ConstantArray[0, {c, c}];
   Do[t[[i, i + 1]] = Sqrt[i], {i, 1, c - 1}]; 
   t
   ];

aDag[cut_] := ConjugateTranspose[a[cut]];

eye[cut_] := IdentityMatrix[cut];

Now I substitute the parameters kx,ky,kz by the blocks defined above:

HLLz[kz_, cut_] := 
  ArrayFlatten[
   H[kx, ky, kz] /. {kx -> (a[cut] + aDag[cut])/Sqrt[2], 
     ky -> I*(a[cut] - aDag[cut])/Sqrt[2], 
     kz -> kz*IdentityMatrix[cut]}];

For example,

HLLz[kz, 6] // MatrixForm
Eigenvalues[HLLz[kz, 3]] // MatrixForm

correctly output

These eigenvalues can be plotted as a function of kz

Plot[Sort[{-kz, kz, -Sqrt[2 + kz^2], Sqrt[2 + kz^2], -Sqrt[4 + kz^2], 
   Sqrt[4 + kz^2]}, Greater], {kz, -2, 2}]

Leading again to a correct output

However, when I calculate the eigenvalues and plot them in one step (for example when I want to apply this calculations to more complicated models which do not have analytic solution for eigenvalues), I get something wrong. This piece of code:

Plot[Sort[Eigenvalues[HLLz[kz, 3]], Greater], {kz, -2, 2}]

outputs the following different and wrong plot of eigenvalues:

Could anybody explain to me why the two pieces of code yield different plots? How do I correct the latter piece of code where I calculate the eigenvalues and plot them in a single cell/step?

Answer 1

It's hard to trace out exactly what is happening (for instance by using Trace: there is just too much there to sift through), but problem is in the definition of HLLz and the difference between calling it with a symbolic vs a numerical kz. Note that

Eigenvalues[HLLz[kz, 3]] /. kz -> 0.1
(* {-0.1, 0.1, -1.41774, 1.41774, -2.0025, 2.0025} *)

is correct whereas

Eigenvalues[HLLz[0.1, 3]]
(* {-2.01133, 2.00884, 1.41814, -1.41714, 0.1, -0.0985211} *)

is not. The issue is due to calling H[kx, ky, kz] inside the definition of HLLz and then making the the replacement

kz -> kz*IdentityMatrix[cut]

Since kz is the argument for the original function, the numerical value of kz gets put in first before the replacement is made, and so the matrix and hence the eigenvalues will be different.

To fix this in the simplest way, make the following modifications: In the definition of HLLz[kz_, cut_], use HLLz[kzp_, cut_] instead and modify kz -> kz*IdentityMatrix[cut] to kz -> kzp*IdentityMatrix[cut], i.e. define

HLLz[kzp_, cut_] := ArrayFlatten[
  H[kx, ky, kz] /. {
    kx -> (a[cut] + aDag[cut])/Sqrt[2], 
    ky -> I*(a[cut] - aDag[cut])/Sqrt[2], 
    kz -> kzp*IdentityMatrix[cut]
   }
 ];