Matsubara frequency

General Formalism

The trick to evaluate Matsubara frequency summation is to use a Matsubara weighting function h_η(z) that has simple poles located exactly at z = i ω . The weighting functions in the boson case η=+1 and fermion case η=-1 differ. The choice of weighting function will be discussed later. With the weighting function, the summation can be replaced by a contour integral in the complex plane.

S η = 1 β ∑ i ω g ( i ω ) = 1 2 π i β ∮ g ( z ) h η ( z ) d z .

As in Fig. 1, the weighting function generates poles (red crosses) on the imaginary axis. The contour integral picks up the residue of these poles, which is equivalent to the summation.

By deformation of the contour lines to enclose the poles of g(z) (the green cross in Fig. 2), the summation can be formally accomplished by summing the residue of g(z)h_η(z) over all poles of g(z),

S η = − 1 β ∑ z 0 ∈ g ( z )  poles Res g ( z 0 ) h η ( z 0 ) .

Note that a minus sign is produced, because the contour is deformed to enclose the poles in the clockwise direction, resulting in the negative residue.

Choice of Matsubara Weighting Function

To produce simple poles on boson frequencies z = i ω n , either of the following two types of Matsubara weighting functions can be chosen

h B ( 1 ) ( z ) = β 1 − e − β z = − β n B ( − z ) = β ( 1 + n B ( z ) ) , h B ( 2 ) ( z ) = − β 1 − e β z = β n B ( z ) ,

depending on which half plane the convergence is to be controlled in. h B ( 1 ) ( z ) controls the convergence in the left half plane (Re z<0), while h B ( 2 ) ( z ) controls the convergence in the right half plane (Re z>0). Here n B ( z ) = ( e β z − 1 ) − 1 is the Bose–Einstein distribution function.

The case is similar for fermion frequencies. There are also two types of Matsubara weighting functions that produce simple poles at z = i ω m

h F ( 1 ) ( z ) = β 1 + e − β z = β n F ( − z ) = β ( 1 − n F ( z ) ) , h F ( 2 ) ( z ) = − β 1 + e β z = − β n F ( z ) .

h F ( 1 ) ( z ) controls the convergence in the left half plane (Re z<0), while h F ( 2 ) ( z ) controls the convergence in the right half plane (Re z>0). Here n F ( z ) = ( e β z + 1 ) − 1 is the Fermi–Dirac distribution function.

In the application to Green's function calculation, g(z) always have the structure

g ( z ) = G ( z ) e − z τ ,

which diverges in the left half plane given 0<τ<β. So as to control the convergence, the weighting function of the first type is always chosen h η ( z ) = h η ( 1 ) ( z ) . However there is no need to control the convergence if the Matsubara summation does not diverge, in that case, any choice of the Matsubara weighting function will lead to identical results.

Table of Matsubara Frequency Summations

The following table concludes the Matsubara frequency summations for some simple rational functions g(z).

S η = 1 β ∑ i ω g ( i ω ) .

η=±1 marks the statistical sign.

[1] Since the summation does not converge, the result may differ by a constant upon different choice of the Matsubara weighting function.

[2] (1↔2) denotes the same expression as the before but with index 1 and 2 interchanged.

Zero Temperature Limit

In this limit β → ∞ , the Matsubara frequency summation is equivalent to the integration of imaginary frequency over imaginary axis.

1 β ∑ i ω = ∫ − i ∞ i ∞ d ( i ω ) 2 π .

Some of the integrals do not converge. They should be regularized by introducing the frequency cutoff Ω , and then subtracting the divergent part ( Ω -dependent) from the integral before taking the limit of Ω → ∞ . For example, the free energy is obtained by the integral of logarithm,

η lim Ω → ∞ [ ∫ − i Ω i Ω d ( i ω ) 2 π ( ln ⁡ ( − i ω + ξ ) − π ξ 2 Ω ) − Ω π ( ln ⁡ Ω − 1 ) ] = { 0 ξ ≥ 0 − η ξ ξ < 0 ,

meaning that at zero temperature, the free energy simply relates to the internal energy below the chemical potential. Also the distribution function is obtained by the following integral

η lim Ω → ∞ ∫ − i Ω i Ω d ( i ω ) 2 π ( 1 − i ω + ξ − π 2 Ω ) = { 0 ξ ≥ 0 − η ξ < 0 ,

which shows step function behavior at zero temperature.

Time Domain

Consider a function G(τ) defined on the imaginary time interval (0,β). It can be given in terms of Fourier series,

G ( τ ) = 1 β ∑ i ω G ( i ω ) e − i ω τ ,

where the frequency only takes discrete values spaced by 2π/β.

The particular choice of frequency depends on the boundary condition of the function G(τ). In physics, G(τ) stands for the imaginary time representation of Green's function

G ( τ ) = − ⟨ T τ ψ ( τ ) ψ ∗ ( 0 ) ⟩ .

It satisfies periodic boundary condition G(τ+β)=G(τ) for boson field. While for fermion field the boundary condition is anti-periodic G(τ+β)=-G(τ).

Given the Green's function G(iω) in the frequency domain, its imaginary time representation G(τ) can be evaluated by Matsubara frequency summation. Depending on the boson or fermion frequencies that is to be summed over, the resulting G(τ) can be different. To distinguish, define

G η ( τ ) = { G B ( τ ) , if  η = + 1 G F ( τ ) , if  η = − 1 ,

with

G B ( τ ) = 1 β ∑ i ω n G ( i ω n ) e − i ω n τ , G F ( τ ) = 1 β ∑ i ω m G ( i ω m ) e − i ω m τ .

Note that τ is restricted in the principal interval (0,β). The boundary condition can be used to extend G(τ) out of the principal interval. Some frequently used results are concluded in the following table.

Operator Switching Effect

The small imaginary time plays a critical role here. The order of the operators will change if the small imaginary time changes sign.

⟨ ψ ψ ∗ ⟩ = ⟨ T τ ψ ( τ = 0 + ) ψ ∗ ( 0 ) ⟩ = − G η ( τ = 0 + ) = − 1 β ∑ i ω G ( i ω ) e − i ω 0 + . ⟨ ψ ∗ ψ ⟩ = η ⟨ T τ ψ ( τ = 0 − ) ψ ∗ ( 0 ) ⟩ = − η G η ( τ = 0 − ) = − η β ∑ i ω G ( i ω ) e i ω 0 +

Distribution Function

The evaluation of distribution function becomes tricky because of the discontinuity of Green's function G(τ) at τ=0. To evaluate the summation

G ( 0 ) = ∑ i ω ( i ω − ξ ) − 1 ,

both choices of the weighting function are acceptable, but the results are different. This can be understood if we push G(τ) away from τ=0 a little bit, then to control the convergence, we must take h η ( 1 ) ( z ) as the weighting function for G ( τ = 0 + ) , and h η ( 2 ) ( z ) for G ( τ = 0 − ) .

Bosons

G B ( τ = 0 − ) = 1 β ∑ i ω n e i ω n 0 + i ω n − ξ = − n B ( ξ ) , G B ( τ = 0 + ) = 1 β ∑ i ω n e − i ω n 0 + i ω n − ξ = − ( n B ( ξ ) + 1 ) .

Fermions

G F ( τ = 0 − ) = 1 β ∑ i ω m e i ω m 0 + i ω m − ξ = n F ( ξ ) , G F ( τ = 0 + ) = 1 β ∑ i ω m e − i ω m 0 + i ω m − ξ = − ( 1 − n F ( ξ ) ) .

Free Energy

Bosons

1 β ∑ i ω n ln ⁡ ( β ( − i ω n + ξ ) ) = 1 β ln ⁡ ( 1 − e − β ξ ) ,

Fermions

− 1 β ∑ i ω m ln ⁡ ( β ( − i ω m + ξ ) ) = − 1 β ln ⁡ ( 1 + e − β ξ ) .

Diagrams Evaluation

Frequently encountered diagrams are evaluated here with the single mode setting. Multiple modes problem can be approached by spectral function integral.

Fermion Self Energy

Σ ( i ω m ) = − 1 β ∑ i ω n 1 i ω m + i ω n − ϵ 1 i ω n − Ω = n F ( ϵ ) − n F ( Ω ) i ω m − ϵ + Ω .

Particle-Hole Bubble

Π ( i ω n ) = 1 β ∑ i ω m 1 i ω m + i ω n − ϵ 1 i ω m − ϵ ′ = − n F ( ϵ ) − n F ( ϵ ′ ) i ω n − ϵ + ϵ ′ .

Particle-Particle Bubble

Π ( i ω n ) = − 1 β ∑ i ω m 1 i ω m + i ω n − ϵ 1 − i ω m − ϵ ′ = 1 − n F ( ϵ ) − n F ( ϵ ′ ) i ω n − ϵ − ϵ ′ .

Distribution Functions

The general notation n η stands for either Bose (η=+1) or Fermi (η=-1) distribution function

n η ( ξ ) = 1 e β ξ − η .

If necessary, the specific notations n_B and n_F are used to indicate Bose and Fermi distribution functions respectively

n η ( ξ ) = { n B ( ξ ) , if  η = + 1 n F ( ξ ) , if  η = − 1 .

Relation to hyperbolic functions

The Bose distribution function is related to hyperbolic cotangent function by

n B ( ξ ) = 1 2 ( c o t h β ξ 2 − 1 ) .

The Fermi distribution function is related to hyperbolic tangent function by

n F ( ξ ) = 1 2 ( 1 − t a n h β ξ 2 ) .

Parity

Both distribution functions do not have definite parity,

n η ( − ξ ) = − η − n η ( ξ ) .

Another formula is in terms of the c η function

n η ( − ξ ) = n η ( ξ ) + 2 ξ c η ( 0 , ξ ) .

However their derivatives have definite parity.

Bose-Fermi Transmutation

Bose and Fermi distribution functions transmute under a shift of the variable by the fermionic frequency,

n η ( i ω m + ξ ) = − n − η ( ξ ) .

However shifting by bosonic frequencies does not make any difference.

First order

n B ′ ( ξ ) = − β 4 c s c h 2 β ξ 2 , n F ′ ( ξ ) = − β 4 s e c h 2 β ξ 2 .

In terms of product:

n η ′ ( ξ ) = − β n η ( ξ ) ( 1 + η n η ( ξ ) ) .

In the zero temperature limit:

n η ′ ( ξ ) = η δ ( ξ ) as β → ∞ .

Second order

n B ′ ′ ( ξ ) = β 2 4 c s c h 2 β ξ 2 c o t h β ξ 2 , n F ′ ′ ( ξ ) = β 2 4 s e c h 2 β ξ 2 t a n h β ξ 2 .

Formula of difference

n η ( a + b ) − n η ( a − b ) = − s i n h β b c o s h β a − η c o s h β b .

Case a=0

n B ( b ) − n B ( − b ) = c o t h β b 2 , n F ( b ) − n F ( − b ) = − t a n h β b 2 .

Case a→0

n B ( a + b ) − n B ( a − b ) = c o t h β b 2 + n B ′ ′ ( b ) a 2 + ⋯ , n F ( a + b ) − n F ( a − b ) = − t a n h β b 2 + n F ′ ′ ( b ) a 2 + ⋯ .

Case b→0

n B ( a + b ) − n B ( a − b ) = 2 n B ′ ( a ) b + ⋯ , n F ( a + b ) − n F ( a − b ) = 2 n F ′ ( a ) b + ⋯ .

The function cη

Definition:

c η ( a , b ) ≡ − n η ( a + b ) − n η ( a − b ) 2 b .

For Bose and Fermi type:

c B ( a , b ) ≡ c + ( a , b ) , c F ( a , b ) ≡ c − ( a , b ) .

Relation to hyperbolic functions

c η ( a , b ) = s i n h β b 2 b ( c o s h β a − η c o s h β b ) .

It is obvious that c F ( a , b ) is positive definite.

To avoid overflow in the numerical calculation, the tanh and coth functions are used

c B ( a , b ) = 1 4 b ( c o t h β ( a − b ) 2 − c o t h β ( a + b ) 2 ) , c F ( a , b ) = 1 4 b ( t a n h β ( a + b ) 2 − t a n h β ( a − b ) 2 ) .

Case a=0

c B ( 0 , b ) = − 1 2 b c o t h β b 2 , c F ( 0 , b ) = 1 2 b t a n h β b 2 .

Case b=0

c B ( a , 0 ) = β 4 c s c h 2 β a 2 , c F ( a , 0 ) = β 4 s e c h 2 β a 2 .

Low temperature limit

For a=0: c F ( 0 , b ) = 1 2 | b | .

For b=0: c F ( a , 0 ) = δ ( a ) .

In general, c F ( a , b ) = { 1 2 | b | , if  | a | < | b | 0 , if  | a | > | b |