5
$\begingroup$

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:

Integrand of s1

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:

enter image description here

Leading to:

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:

Integrand of s2

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

Improve this question
$\endgroup$
11
  • $\begingroup$ What's a typical value for $\beta$ like? $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 10, 2015 at 16:19
  • $\begingroup$ 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. $\endgroup$
    – zakk
    Oct 10, 2015 at 16:22
  • $\begingroup$ I just tried $\beta=10$ and $\beta=20$, the convergence problems are unchanged. They disappear in the low-temperature limit, $\beta=100$, for example. $\endgroup$
    – zakk
    Oct 10, 2015 at 16:23
  • $\begingroup$ 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... $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 10, 2015 at 16:47
  • 1
    $\begingroup$ g is undefined in the question. I presume it should be g[z_, bval_] := bval/2 Tanh[z bval/2] $\endgroup$
    – bbgodfrey
    Oct 10, 2015 at 16:55

1 Answer 1

Reset to default
6
$\begingroup$

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[%]

enter image description here

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.

Addendum

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]]
(* -beta ((-0.0772029 + 0.0405577 I) Tan[(0.259456 + 0.333305 I) beta] - 
          (0.0772029 - 0.172992 I) Tan[(0.509456 + 0.301283 I) beta] - 
          (0.172992 - 0.0772029 I) Tanh[(0.301283 + 0.509456 I) beta] - 
          (0.0405577 - 0.0772029 I) Tanh[(0.333305 + 0.259456 I) beta]) *)
(% /. beta -> .2) // Chop
(* 0.0100333 *)

as expected.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.