# A 1D numerical integral Mathematica cannot compute, from physics

A well know result in theoretical physics is that a sum over Matsubara fermionic frequencies, i.e.:

$$S = \sum_{n=-\infty}^{\infty} h(\omega_n) \hspace{32pt} \omega_n=(2n+1)\frac{\pi}{\beta}$$

can be rewritten as a contour integral as

$$S = \frac{-1}{2 \pi \mathrm{i}} \oint_C \mathrm{d} z g(z) h(-\mathrm{i}z) \hspace{32pt} g(z)=\frac{\beta}{2} \tanh(\frac{\beta z}{2})$$

where $C$ is an infinite contour enclosing counterclockwise the imaginary axis in the complex plane. This result is demonstrated, for instance, in Altand & Simons, "Condensed Matter Field Theory", Cambridge University Press (2010).

If the integrand decays fast enough the contour $C$ can be split in two contours, parallel to the imaginary axis. Enough with the physics! To make a long story short the theory implies that the following functions in Mathematica should give the same result (for any choice of dist):

h[omega_] := omega^2/(omega^4 + omega^3 + 1)
omega[n_] := (2*n + 1)*\[Pi]/b
s1[bval_] := Sum[h[omega[n]], {n, -Infinity, Infinity}] /. {b -> bval}
g[z_, b_] := b/2 Tanh[z b/2]
s2[b_, dist_] := 1/(2*\[Pi])*NIntegrate[(g[z1, b]*h[-I*z1] - g[z2, b]*h[-I*z2])
/. {z1 -> I*x + dist, z2 -> I*x - dist}, {x, -Infinity, Infinity},
Method -> "DoubleExponentialOscillatory", MaxRecursion -> 16]


i.e. s1[b] should evaluate to the same result as s2[b,d] for any choice of b and d. However, while the sum in s1 is evaluated instantly, the integral in s2 fails to converge, even if a plot of the integrand looks good enough:

Note: The peaks are in infinite number, periodically spaced, but their height is rapidly decreasing. Probably this confuses Mathematica integration strategies.

I tried to change the integration method, AccuracyGoal and Precision goal, as long as increasing MaxRecursion. Sometimes the results of the integral s2 are very near to the exact result from s1, but I always get the NIntegrate::ncvb error (and most often also NIntegrate::slwcon).

I also tried to convert the infinite integral to a finite one, using:

s3[b_, dist_] := 1/(2*\[Pi])*NIntegrate[(g[z1, b]*h[-I*z1] - g[z2, b]*h[-I*z2])
*(1 + t^2)/(1 - t^2)^2 /. {z1 -> I*t/(1 - t^2) + dist,z2 -> I*t/(1 - t^2) - dist},
{t, -1, 1}]


And the peaks become now an oscillating behaviour at the borders, which apparently Mathematica cannot handle:

This conversion seems to help somehow, but I never get an errorless, precise answer from Mathematica.

• What's a typical value for $\beta$ like? – J. M. will be back soon Oct 10 '15 at 16:19
• Sorry, I forgot to mention it, for the plots I am using $\beta=0.2$. Physically speaking it is the inverse temperature, so $\beta=0.2$ corresponds to high temperatures. – zakk Oct 10 '15 at 16:22
• I just tried $\beta=10$ and $\beta=20$, the convergence problems are unchanged. They disappear in the low-temperature limit, $\beta=100$, for example. – zakk Oct 10 '15 at 16:23
• Do you really need to use contour integration to evaluate your infinite sum? There are a lot of alternative approaches you could use for infinite summation... – J. M. will be back soon Oct 10 '15 at 16:47
• g is undefined in the question. I presume it should be g[z_, bval_] := bval/2 Tanh[z bval/2] – bbgodfrey Oct 10 '15 at 16:55

To understand this integral, it is helpful to know the location of the poles. The poles of g are given by

poleg = Simplify[Solve[Cosh[beta z/2] == 0, z], C[1] ∈ Integers] /. C[1] -> n;
(z /. #) & /@ poleg
(* {(I (-1 + 4 n) π)/beta, (I (π + 4 n π))/beta} *)


which agrees with the expression in the Question. However, h also has poles.

poleh = Solve[Denominator[h[-I z]] == 0, z];
N[z /. #] & /@ poleh
(* {-0.66661 + 0.518913 I, -0.602565 - 1.01891 I, 0.602565 - 1.01891 I,
0.66661 + 0.518913 I} *)


The contour integral is equal to the sum of residues of all enclosed poles. Consider now the value of the integral for beta = 0.2 as dist is varied. (Quiet is applied to eliminate a plethora of error messages that, nonetheless, do not prevent reasonable results.)

Quiet@Table[{dist, s2[2/10, dist] // Chop}, {dist, .5, .7, .01}];
ListPlot[%]


Only integrals for dist < 0.602565 would be expected to equal s1[0.2], and indeed that is the case.

s1[2/10] // N // Chop
(* 0.0100333 *)


The average value of the integral for dist < 0.602565, as shown in the plot, is about 0.01002, which seems like reasonable agreement.

Because the contour integrals enclosing all poles seem to be equal to zero (approximately -10^-5 in the plot), it follows that the sum of the g residues must equal the negative of the sum of the h residues. Symbolic expressions for the h residues are quite lengthy but can be greatly simplified with N.
-Simplify[Total[N[Residue[h[-I z] g[z, beta], {z, z /. #}]] & /@ poleh]]