1D Integral always evaluates to 0 and Replace All

Question

I am a bit puzzled by Mathematica's behaviour. I am trying to evaluate some 1D integrals which depend on parameters phi1,phi2 numerically (functions f1 and f2), however NIntegrate always outputs 0 whereas plotting the integrand with Plot shows that the integral is non-zero.

Here is the code in question:

t = 1;
lambda = 0.3;
hbarOmega = 4;
q = 0.5;
alpha = 4 lambda^2 t (Sin[q/2])/q;

A[k_] := 2 t (Sin[q/2]/q) Sin[k] (1/(hbarOmega + 4 t Sin[q/2] Sin[k]) + 
 1/(hbarOmega - 4 t Sin[q/2] Sin[k]));

B[k_] := 2 t (Sin[q/2]/q) Sin[k] (1/(hbarOmega + 4 t Sin[q/2] Sin[k]) - 
 1/(hbarOmega - 4 t Sin[q/2] Sin[k]));

e1[k_] := 4 t lambda^2  (Sin[q/2]/q) B[k + q/2] Sin[k + q/2];
e2[k_] := -4 t lambda^2  (Sin[q/2]/q) B[k - q/2] Sin[k - q/2];
en[k_] := -2 t Cos[k] + e1[k - q] + e2[k + q];
F[k_, phi1_, phi2_] := - alpha (A[k] phi1 + Sin[k] phi2);
mu[phi1_, phi2_] := - 2 alpha Re[phi1 Conjugate[phi2]];

hMf[k_, phi1_, phi2_] := {{en[k + q/2],F[k, phi1, phi2]}, {Conjugate[F[k, phi1, phi2]], en[k - 
q/2]}};


u[k_, phi1_, phi2_] := Transpose[Eigenvectors[hMf[k, phi1, phi2]]];
enMf[k_, phi1_, phi2_] := ConjugateTranspose[u[k, phi1, phi2]].hMf[k, phi1, phi2].u[k, phi1, 
phi2];

o[k_, phi1_, phi2_] := ConjugateTranspose[u[k, phi1, phi2]].{{0, 0}, {1, 0}}.u[k, phi1,phi2];


f1[phi1_?NumericQ,phi2_?NumericQ] := (1/(2 Pi)) NIntegrate[
Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
     enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
   o[k, phi1, phi2][[2, 2]] Boole[
     enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]), {k, -Pi, Pi}];

f2[phi1_?NumericQ,phi2_?NumericQ] := (1/(2 Pi)) NIntegrate[A[k] (o[k, phi1, phi2][[1, 1]] 
Boole[enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
   o[k, phi1, phi2][[2, 2]] Boole[
     enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]), {k, -Pi,Pi}];

f1[1, 1]

Plot[Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
   enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
 o[k, phi1, phi2][[2, 2]] Boole[
   enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]) /. {phi1 -> 1, phi2 -> 1}, {k, -Pi, Pi}]

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + o[k, 1, 1][[2, 2]] 
Boole[enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

So for example, when I evaluate the integral given by f1[1,1], I get an output of 0 alongside with a warning that Mathematica couldn't find all eigenvectors of a 2x2 matrix which is bizarre, because when I plot the eigenvalues separately with:

Plot[{enMf[k, 1, 1][[1, 1]], enMf[k, 1, 1][[2, 2]]}, {k, -Pi, Pi}]

the warning is not present.

However when I plot the integrand inside f1 with the same parameters phi1=1,phi2=1 with:

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[
  enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + 
o[k, 1, 1][[2, 2]] Boole[
  enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

it shows that the integrand is always positive and that the integral should be non-zero.

Moreover, there is another oddity: replacing the parameters phi1,phi2 within the integrand with /. for the same values phi1=1,phi2=1 like so:

Plot[Sin[k]*(o[k, phi1, phi2][[1, 1]] Boole[
   enMf[k, phi1, phi2][[1, 1]] <= mu[phi1, phi2]] + 
 o[k, phi1, phi2][[2, 2]] Boole[
   enMf[k, phi1, phi2][[2, 2]] <= mu[phi1, phi2]]) /. {phi1 -> 1,phi2 -> 1}, {k, -Pi, Pi}]

gives a different plot compared to just substituting the values phi1=1,phi2=1 directly:

Plot[Sin[k]*(o[k, 1, 1][[1, 1]] Boole[
  enMf[k, 1, 1][[1, 1]] <= mu[1, 1]] + 
o[k, 1, 1][[2, 2]] Boole[
  enMf[k, 1, 1][[2, 2]] <= mu[1, 1]]), {k, -Pi, Pi}]

which I find very bizarre and assume has to do with delayed evaluation. Any suggestions would be appreciated.

Answer 1

You're right that it's a delayed evaluation problem! I traced the problem back to Eigenvectors. Note that it behaves differently when given symbolic arguments and numeric ones:

Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1}

(* Out: {{-4.064089631, 1.}, {0.2460575654, 1.}} *)

Eigenvectors[hMf[0.01, 1, 1]]

(* Out: {{0.2389308952, 0.9710365736}, {-0.9710365736, 0.2389308952}} *)

The deal is that Mathematica both 1) doesn't normalize eigenvectors when symbolic and 2) might order them differently. Note that

Normalize /@ (Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1})

(* Out: {{-0.9710365736, 0.2389308952}, {0.2389308952, 0.9710365736}} *)

Which is the same as Eigenvectors[hMf[0.01, 1, 1]], but in reverse order.

Sort[Normalize /@ (Eigenvectors[hMf[0.01, phi1, phi2]] /. {phi1 -> 1, phi2 -> 1})]
  == Sort[Eigenvectors[hMf[0.01, 1, 1]]]

(* Out: True *)

So, you'll need to decide whether the vectors should be normalized, and whether the order matters!

But there's also another issue: different normalizations can be chosen. Normalized eigenvectors are unique only up to norm 1 elements of your field, which in this case are $\pm 1$ if we're doing this over the reals. Which normalization do you want? For this reason, even after redefining u[k_, phi1_, phi2_] := Transpose[Normalize /@ Eigenvectors[hMf[k, phi1, phi2]]];, we still have

o[0.01, phi1, phi2][[2, 2]] /. {phi1 -> 1, phi2 -> 1}

(* Out: 0.2320106378 *)

o[0.01, 1, 1][[2, 2]] 
(* Out: -0.2320106378 *)

So that's another choice to make as well.

---

However, there is a way to avoid some of this, since I noticed you seem to be using the matrix of eigenvectors for diagonalization. Check out JordanDecomposition and other *Decomposition functions, which return not only the decomposition, but similarity matrices you can use as well! Also look at Eigensystem, which returns a list {evalues, evectors} for which M.Transpose[evectors] == Transpose[evectors].DiagonalMatrix[evalues]. You could use this to save on some computation in calculating enMf by simply using DiagonalMatrix[evalues].