# 1D Integral always evaluates to 0 and Replace All

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.

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].

• Thanks for this. Concerning your remark on normalization (last bit of code in your answer), I think that its the ordering of the vectors that's causing the sign change, as the structure of the matrix $o$ ( $o = u^{\dagger} \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} u$) implies that its diagonal elements are invariant under the change in sign of any of the eigenvectors $v_1,v_2$ making up $u = (v_1,v_2)$. Also, now I understand why the plots of the integrand are different but I am still not sure why the integral f[1,1] is zero since a plot of the integrand suggests otherwise. Commented Sep 7, 2021 at 8:48
• Update: so making all the functions from hMf down numerical with ?NumericQ gives a non-zero integral, however now I get error messages like so: "Part specification o[k,2,2][[1,1]] is longer than depth of object" which I guess arise from Part being evaluated prematurely. I am a bit confused by this message because since we are using NIntegrate, shouldnt k always have a numerical value? Commented Sep 7, 2021 at 9:26
• @RuslanMushkaev Try Indexed[o[k.2.2], {1, 1}] instead of Part. NIntegrate inspects the integrand symbolically to determine which method to choose. Good thing you didn’t have o[k, 2, 2][[1]], because that evaluates to k, which would give you strange results and no error. Commented Sep 7, 2021 at 14:08