diff --git a/src/interaction.jl b/src/interaction.jl
index a45d56c..fa8f406 100644
--- a/src/interaction.jl
+++ b/src/interaction.jl
@@ -525,4 +525,201 @@ function RPA_total(q, n, param; pifunc=Polarization0_ZeroTemp, Vinv_Bare=coulomb
     return Ws, Wa
 end
 
+"""
+    function oldTmatrix(q, n, param; pifunc = Polarization0_ZeroTemp, landaufunc = landauParameterTakada, Vinv_Bare = coulombinv, regular = false, kwargs...)
+
+Dynamic part of the T-matrix. Returns the spin symmetric part and asymmetric part separately.
+
+#Arguments:
+ - q: momentum
+ - n: matsubara frequency given in integer s.t. ωn=2πTn (either one frequency or array of frequencies)
+ - param: other system parameters
+ - pifunc: caller to the polarization function 
+ - landaufunc: caller to the Landau parameter (exchange-correlation kernel)
+ - Vinv_Bare: caller to the bare Coulomb interaction
+ - regular: regularized RPA or not
+
+# Return:
+```math
+    \\frac{1} {\\frac{m}{4\\pi a_s} +\\Gamma(q, n)} - \\frac{4\\pi a_s}{m}.
+```
+"""
+function oldTmatrix(q, n, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...)
+    as = 4π * param.as / param.me
+    if param.as < 1e-6
+        return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...))
+    else
+        return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...))
+    end
+end
+
+"""
+    function oldTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; int_type=:ko,
+        pifunc=Polarization0_ZeroTemp, landaufunc=landauParameterTakada, Vinv_Bare=coulombinv, kwargs...) where {SGT}
+
+Return dynamic part and instant part of KO-interaction Green's function separately. Each part is a MeshArray with inner state 
+(1: spin symmetric part, 2: asymmetric part), and ImFreq and q-grid mesh.
+
+#Arguments:
+ - Euv: Euv of DLRGrid
+ - rtol: rtol of DLRGrid
+ - sgrid: momentum grid
+ - param: other system parameters
+ - pifunc: caller to the polarization function
+ - landaufunc: caller to the Landau parameter (exchange-correlation kernel)
+ - Vinv_Bare: caller to the bare Coulomb interaction
+"""
+function oldTmatrix_wrapped(Euv, rtol, sgrid::SGT, param;
+    ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT}
+    # TODO: innerstate should be in the outermost layer of the loop. Hence, the functions such as KO and Vinv_Bare should be fixed with inner state as argument.
+    @unpack β = param
+    as = 4π * param.as / param.me
+    wn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph)
+    green_dyn_r = GreenFunc.MeshArray(wn_mesh_ph, sgrid; dtype=ComplexF64)
+    green_tail_r = GreenFunc.MeshArray(wn_mesh_ph, sgrid; dtype=ComplexF64)
+
+    wn_mesh_pha = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:pha)
+    green_dyn_i = GreenFunc.MeshArray(wn_mesh_pha, sgrid; dtype=ComplexF64)
+    green_tail_i = GreenFunc.MeshArray(wn_mesh_pha, sgrid; dtype=ComplexF64)
+
+    for (ki, k) in enumerate(sgrid)
+        for (ni, n) in enumerate(wn_mesh_ph.grid)
+            green_tail_r[ni, ki] = real.(Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...))
+            green_dyn_r[ni, ki] = real.(Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...)) - green_tail_r[ni, ki]
+        end
+        for (ni, n) in enumerate(wn_mesh_pha.grid)
+            green_tail_i[ni, ki] = 1im .* imag.(Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...))
+            green_dyn_i[ni, ki] = 1im .* imag.(Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...)) - green_tail_i[ni, ki]
+        end
+    end
+
+    return green_dyn_r, green_dyn_i, green_tail_r, green_tail_i
+end
+
+
+"""
+    function T0(q::Float64, n::Int, param; ; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...)
+
+Vacuum two-particle T-matrix in matsubara frequency. 
+#Arguments:
+ - q: total incoming momentum
+ - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies)
+ - param: other system parameters
+"""
+
+
+function T0matrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...)
+    as = (4π * param.as) / param.me
+    if param.as < 1e-6
+        return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...))
+    else
+        return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...) )
+    end
+end
+
+"""
+    function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.LadderVacuum0_FiniteTemp, kwargs...)
+
+Many-body T-matrix in matsubara frequency. 
+#Arguments:
+ - q: total incoming momentum
+ - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies)
+ - param: other system parameters
+ - ladderfunc: Particle-particle bubble (ladder-diagram of G0G0)
+"""
+
+function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...)
+    as = (4π * param.as) / param.me
+    if param.as < 1e-6
+        return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...))
+    else
+        return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...) )
+    end
+end
+
+"""
+    function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT}
+
+Difference between Many-body T-matrix and vacuum two-particle T-matrix in imaginary time. 
+#Arguments:
+- Euv: Euv of DLRGrid
+- rtol: rtol of DLRGrid
+- sgrid: momentum grid
+- param: other system parameters
+- ladderfunc: Particle-Particle bubble (ladder-diagram of G0G0)
+"""
+
+function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT}
+    β= param.β
+    Ωn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph)
+    δTmatrix_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64)
+
+    for (ki, k) in enumerate(sgrid)
+        for (ni, n) in enumerate(Ωn_mesh_ph.grid)
+            δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc) - T0matrix(k, n, param))
+        end
+    end
+    δTmatrix_tau = δTmatrix_n |> to_dlr |> to_imtime
+    return δTmatrix_tau
+end
+
+"""
+    function T0maxtrix_imtime(tau::Float64, q::Float64, param)
+
+Difference between Many-body T-matrix and vacuum two-particle T-matrix in imaginary time. 
+#Arguments:
+- tau: imaginary time
+- q: total incoming momentum
+- param: other system parameters
+"""
+
+function T0matrix_imtime(tau::Float64, q::Float64,  param)
+    """ Two particle scattering in a vacuum T-matrix in imaginary time """
+    β, me, as = param.β, param.me, param.as
+    B = me^(3/2)/(4π)
+    m = me * as^2
+
+    # Pole Term contribution
+    poletermnumerator = exp(tau/ m)
+    poletermdenominator = 1- exp(β * (1. / (me*as^2) - q^2/ (4 * me)))
+    poleterm = - poletermnumerator/ poletermdenominator
+    # Integral Terms
+
+    function integrand1(vars, tau, param)
+        me, as = param.me, param.as
+        m = me * as^2
+        u = vars
+        y = u/(1-u)
+        fraction = y^2/(y^2 + tau/ m)
+        return 2 * exp(-y^2) * fraction
+    end
+    
+    function integrand2(vars, q, tau, param)
+        β, me = param.β, param.me
+        u = vars
+        y = u/(1-u)
+        factor = 1/(exp(β * (y^2/ tau + q^2/ (4 * me)))-1)
+        return factor * integrand1(vars, tau, param)
+    end
+
+    tgrid = CompositeGrid.LogDensedGrid(
+        :gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb
+        [0.0, 1.0],# The grid is defined on [0.0, 1.0] s.t. the integral in the original variable is from [0.0 , infinity]
+        [0.0],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter.
+        5,# N of log grid
+        0.005, # minimum interval length of log grid
+        5 # N of bottom layer
+    )
+    data = [(integrand1(t, tau, param)) for (ti, t) in  enumerate(tgrid.grid)]
+    # data = [0.0 for (ti,t) in  enumerate(tgrid.grid)]
+    int_result = Interp.integrate1D(data, tgrid)
+    integralfactor = 1/(B * pi * sqrt(tau))
+
+    # Propagator factor
+    propfactor = exp(-q^2/ (4 * me) * tau)
+
+    # return propfactor*(poleterm + int_result * integralfactor)
+    return propfactor*(poleterm + int_result * integralfactor)
+end
+
 end
diff --git a/src/legendreinteraction.jl b/src/legendreinteraction.jl
index d06cc1a..e5c11cc 100644
--- a/src/legendreinteraction.jl
+++ b/src/legendreinteraction.jl
@@ -169,7 +169,7 @@ function helper_function_grid(ygrid, intgrid, n::Int, W, param)
 end
 
 
-struct DCKernel{C}
+struct DCKernel{C,T}
     int_type::Symbol
     spin_state::Symbol
     channel::Int
@@ -180,13 +180,16 @@ struct DCKernel{C}
     qgrids::Vector{CompositeGrid.Composite}
     dlrGrid::DLRGrid
 
-    kernel_bare::Array{Float64,2}
-    kernel::Array{Float64,3}
+    kernel_bare::Array{T,2}
+    kernel::Array{T,3}
 
     function DCKernel(int_type, spin_state, channel, param, kgrid::C, qgrids, dlrGrid, kernel_bare, kernel) where {C}
-        return new{C}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel)
+        return new{C,Float64}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel)
     end
 
+    function DCKernel{T}(int_type, spin_state, channel, param, kgrid::C, qgrids, dlrGrid, kernel_bare, kernel) where {C,T}
+        return new{C,T}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel)
+    end
 end
 
 function DCKernel_2d(param, Euv, rtol, Nk, maxK, minK, order, int_type, channel, spin_state=:auto;
@@ -382,4 +385,91 @@ function DCKernel0(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state
     return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel)
 end
 
+
+function DCKernel_Ladder_r(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...)
+    return DCKernel_Ladder_r(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...)
+end
+
+function DCKernel_Ladder_r(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma;
+    kgrid=CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order),
+    kwargs...)
+    # use helper function
+    @unpack kF, β = param
+    channel = 0
+
+    if spin_state == :sigma
+        # for self-energy, always use ℓ=0
+        channel = 0
+    elseif spin_state == :auto
+        error("not implemented!")
+    end
+
+    bdlr = DLRGrid(Euv, β, rtol, false, :ph)
+    qgrids = [CompositeGrid.LogDensedGrid(:gauss, [0.0, maxK], [k, kF], Nk, minK, order) for k in kgrid.grid]
+    qgridmax = maximum([qg.size for qg in qgrids])
+    #println(qgridmax)
+
+    kernel_bare = zeros(ComplexF64, (length(kgrid.grid), (qgridmax)))
+    kernel = zeros(ComplexF64, (length(kgrid.grid), (qgridmax), length(bdlr.n)))
+
+    helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1 * maxK], [0.0, kF], 2Nk, 0.01minK, 2order)
+    intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0, kF], 2Nk, 0.01minK, 2order)
+
+    # dynamic
+    for (ni, n) in enumerate(bdlr.n)
+        helper = helper_function_grid(helper_grid, intgrid, 1, u -> real(Interaction.Tmatrix(u, n, param; kwargs...)), param)
+        for (ki, k) in enumerate(kgrid.grid)
+            for (pi, p) in enumerate(qgrids[ki].grid)
+                Hp, Hm = Interp.interp1D(helper, helper_grid, k + p), Interp.interp1D(helper, helper_grid, abs(k - p))
+                kernel[ki, pi, ni] = (Hp - Hm)
+            end
+        end
+    end
+
+    return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel)
+end
+
+function DCKernel_Ladder_i(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...)
+    return DCKernel_Ladder_i(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...)
+end
+
+function DCKernel_Ladder_i(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma;
+    kgrid=CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order),
+    kwargs...)
+    # use helper function
+    @unpack kF, β = param
+    channel = 0
+
+    if spin_state == :sigma
+        # for self-energy, always use ℓ=0
+        channel = 0
+    elseif spin_state == :auto
+        error("not implemented!")
+    end
+
+    bdlr = DLRGrid(Euv, β, rtol, false, :ph)
+    qgrids = [CompositeGrid.LogDensedGrid(:gauss, [0.0, maxK], [k, kF], Nk, minK, order) for k in kgrid.grid]
+    qgridmax = maximum([qg.size for qg in qgrids])
+    #println(qgridmax)
+
+    kernel_bare = zeros(ComplexF64, (length(kgrid.grid), (qgridmax)))
+    kernel = zeros(ComplexF64, (length(kgrid.grid), (qgridmax), length(bdlr.n)))
+
+    helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1 * maxK], [0.0, kF], 2Nk, 0.01minK, 2order)
+    intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0, kF], 2Nk, 0.01minK, 2order)
+
+    # dynamic
+    for (ni, n) in enumerate(bdlr.n)
+        helper = helper_function_grid(helper_grid, intgrid, 1, u -> imag(Interaction.Tmatrix(u, n, param; kwargs...)), param)
+        for (ki, k) in enumerate(kgrid.grid)
+            for (pi, p) in enumerate(qgrids[ki].grid)
+                Hp, Hm = Interp.interp1D(helper, helper_grid, k + p), Interp.interp1D(helper, helper_grid, abs(k - p))
+                kernel[ki, pi, ni] = (Hp - Hm)
+            end
+        end
+    end
+
+    return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel)
+end
+
 end
diff --git a/src/parameter.jl b/src/parameter.jl
index b431d3f..37f5416 100644
--- a/src/parameter.jl
+++ b/src/parameter.jl
@@ -24,6 +24,7 @@ using ..Parameters
     EF::Float64 = 1.0     #kF^2 / (2me)
     β::Float64 = 200 # bare inverse temperature
     μ::Float64 = 1.0
+    as::Float64 = 2.5   # s-wave scattering length (local contact interaction)
 
     # artificial parameters
     Λs::Float64 = 0.0   # Yukawa-type spin-symmetric interaction  ~1/(q^2+Λs)
@@ -130,7 +131,7 @@ end
 # end
 
 """
-    function fullUnit(ϵ0, e0, me, EF, β)
+    function fullUnit(ϵ0, e0, me, EF, β, dim=3, spin=2; kwargs...)
 
 generate Para with a complete set of parameters, no value presumed.
 
@@ -140,6 +141,8 @@ generate Para with a complete set of parameters, no value presumed.
  - me: electron mass
  - EF: Fermi energy
  - β: inverse temperature
+- dim: dimension
+- spin: spin
 """
 @inline function fullUnit(ϵ0, e0, me, EF, β, dim=3, spin=2; kwargs...)
     # μ = try
diff --git a/src/polarization.jl b/src/polarization.jl
index 81e6c87..bc6dca5 100644
--- a/src/polarization.jl
+++ b/src/polarization.jl
@@ -5,7 +5,7 @@ using ..Parameter
 using ..Parameters
 using ..Convention
 
-export Polarization0_ZeroTemp, Polarization0_FiniteTemp
+export Polarization0_ZeroTemp, Polarization0_FiniteTemp, Ladder0_FiniteTemp
 
 # Analytical calculated integrand of Π0 in 2D.
 @inline function _ΠT2d_integrand(k, q, ω, param)
@@ -50,14 +50,21 @@ end
     @unpack me, β, μ = param
     # ω only appears as ω^2 so no need to check sign of ω
     k = x / (1 - x)
-    nk = 1.0 / (exp(β * (k^2 / 2 / me - μ)) + 1)
+    ek = β * (k^2 / 2 / me - μ)
+    nk = ek < 0.0 ? 1.0 / (exp(ek) + 1) : exp(-ek) / (1 + exp(-ek))
 
     if q > 1e-6 * param.kF
         a = ω * 1im - k^2 / me - q^2 / 2 / me + 2μ
         b = q * k / me
         return me / (2π)^2 * (k / q * (log(a + b) - log(a - b)) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
     else
-        return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+        if ω != 0.0 || abs(ek) > 1e-4
+            me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+        else # ω==0.0 && abs(ek)<1e-4
+            #expansion of 1/ek*(2/(exp(ek)+1)-1)
+            f = -0.5 + ek^2 / 24 - ek^4 / 240 + 17 * ek^6 / 40320
+            return me / (2π)^2 * (2k^2 / me * (f / 2 * β) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+        end
     end
 
     # if q is too small, use safe form
@@ -74,6 +81,43 @@ end
     # end
 end
 
+# Analytical calculated integrand of the vacuum ladder in 3D.
+@inline function _LadderTVacuum3d_integrand(x, q, ω, param)
+    me, β= param.me, param.β
+    # μ = 0.0 for vacuum
+    μ = 0.0
+    # ω only appears as ω^2 so no need to check sign of ω
+    k = x / (1 - x)
+    ek = β * (k^2 / 2 / me - μ)
+    nk = ek < 0.0 ? 1.0 / (exp(ek) + 1) : exp(-ek) / (1 + exp(-ek))
+
+    if q > 1e-6 * param.kF
+        a = ω * 1im - k^2 / me - q^2 / 2 / me + 2μ
+        b = q * k / me
+        return me / (2π)^2 * (k / q * (log(a + b) - log(a - b)) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+    else
+        if ω != 0.0 || abs(ek) > 1e-4
+            return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+        else # ω==0.0 && abs(ek)<1e-4
+            #expansion of 1/ek*(2/(exp(ek)+1)-1)
+            f = -0.5 + ek^2 / 24 - ek^4 / 240 + 17 * ek^6 / 40320
+            return me / (2π)^2 * (2k^2 / me * (f / 2 * β) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian
+        end
+    end
+
+    # if q is too small, use safe form
+    # if q < 1e-16 && ω ≈ 0
+    #     if abs(q - 2 * k)^2 < 1e-16
+    #         return 0.0
+    #     else
+    #         return -k * me / (4 * π^2) * nk * ((8 * k) / ((q - 2 * k)^2))
+    #     end
+    # elseif q < 1e-16 && !(ω ≈ 0)
+    #     return -k * me / (4 * π^2) * nk * ((8 * k * q^2) / (4 * me^2 * ω^2 + (q^2 - 2 * k * q)^2))
+    # else
+    #     return -k * me / (4 * π^2 * q) * nk * log1p((8 * k * q^3) / (4 * me^2 * ω^2 + (q^2 - 2 * k * q)^2))
+    # end
+end
 
 function finitetemp_kgrid(q::Float64, kF::Float64, maxk=20, scaleN=20, minterval=1e-6, gaussN=10)
     mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF)
@@ -90,14 +134,15 @@ function finitetemp_kgrid_ladder(μ::Float64, m::Float64, q::Float64, kF::Float6
     mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF)
     mixX = mink / (1 + mink)
     if q >= sqrt(8 * μ * m)
-        kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0], [0.0, kF], scaleN, mink, gaussN)
+        x = kF / (1 + kF)
+        kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0 - 1e-8], [0.0, x], scaleN, mink, gaussN)
     else
         k1 = abs(q - sqrt(8 * μ * m - q^2)) / 2
         k2 = (q + sqrt(8 * μ * m - q^2)) / 2
         x1 = k1 / (1 + k1)
         x2 = k2 / (1 + k2)
         x3 = kF / (1 + kF)
-        kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0], sort([x1, x2, x3]), scaleN, mink, gaussN)
+        kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0 - 1e-8], sort([x1, x2, x3]), scaleN, mink, gaussN)
     end
     return kgrid
 end
@@ -130,6 +175,72 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6
         error("No support for finite-temperature polarization in $dim dimension!")
     end
 
+    # check sign of q, use -q if negative
+    if q < 0
+        q = -q
+    end
+    Ωn = 2π * n / β
+    # mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF)
+    # kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, maxk * kF], [0.5 * q, kF], scaleN, mink, gaussN)
+    kgrid = finitetemp_kgrid_ladder(μ, me, q, kF, scaleN, minterval, gaussN)
+    integrand = zeros(ComplexF64, kgrid.size)
+    if dim == 3
+        for (ki, k) in enumerate(kgrid.grid)
+            #integrand[ki] = _LadderT3d_integrand(k, q, Ωn, param)
+            if Ωn <100
+                integrand[ki] = _LadderT3d_integrand(k, q, Ωn, param)
+                @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n"
+            else
+                integrand[ki] = 0 # Do not calculate integrand for high frequency tail
+            end
+        end
+    end
+    result = Interp.integrate1D(integrand, kgrid)
+    if Ωn> 100
+        result = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωn - 2μ)) # Replace high frequency tail with two-particle particle-particle term
+    end
+    return result
+end
+
+function Ladder0_FreeElectron(q::Float64, n::Int, param)
+    @unpack dim, kF, β, me, μ = param
+    if dim != 3
+   
+        error("No support for finite-temperature polarization in $dim dimension!")
+    end
+    return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β - 2μ)
+end
+
+"""
+    function Ladder0Vacuum_FiniteTemp(q::Float64, n::Int, param, scaleN=20, minterval=1e-6, gaussN=10)
+
+Finite temperature vacuum ladder function for matsubara frequency and momentum. Analytically sum over total incoming frequency and angular
+dependence of momentum, and numerically calculate integration of magnitude of momentum.
+Assume G_vacuum^{-1} = iω_n - (k^2/(2m))
+
+The ladder function is defined as 
+```math
+\\int \\frac{d^3 \\vec{p}}{\\left(2π^3\\right)} T \\sum_{i ω_n} \\frac{1}{iω_n+iΩ_n-\\frac{(\\vec{k}+\\vec{p})^2}{2 m}} \\frac{1}{-iω_n-\\frac{p^2}{2 m}} - \\frac{m}{2π^2}Λ
+```
+where we subtract the UV divergent term that is a constant proportional to the UV cutoff Λ.
+
+#Arguments:
+ - q: momentum
+ - n: matsubara frequency given in integer s.t. Ωn=2πTn
+ - param: other system parameters
+ - maxk: optional, upper limit of integral -> maxk*kF
+ - scaleN: optional, N of Log grid in LogDensedGrid, check CompositeGrids for more detail
+ - minterval: optional, actual minterval of grid is this value times min(q,kF)
+ - gaussN: optional, N of GaussLegendre grid in LogDensedGrid.
+"""
+function Ladder0Vacuum_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6, gaussN=10)
+    dim, kF, β, me = param.dim, param.kF, param.β, param.me
+    # μ = 0.0 for vacuum
+    μ = 0.0
+    if dim != 3
+        error("No support for finite-temperature polarization in $dim dimension!")
+    end
+
     # check sign of q, use -q if negative
     if q < 0
         q = -q
@@ -141,8 +252,8 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6
     integrand = zeros(ComplexF64, kgrid.size)
     if dim == 3
         for (ki, k) in enumerate(kgrid.grid)
-            integrand[ki] = _LadderT3d_integrand(k, q, 2π * n / β, param)
-            @assert !isnan(integrand[ki]) "nan at k=$k, q=$q"
+            integrand[ki] = _LadderTVacuum3d_integrand(k, q, 2π * n / β, param)
+            @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n"
         end
     end
 
@@ -243,13 +354,38 @@ function Ladder0_FiniteTemp(q::Float64, n::AbstractVector, param; scaleN=20, min
     if dim == 3
         for (ki, k) in enumerate(kgrid.grid)
             for (mi, m) in enumerate(n)
-                integrand[ki, mi] = _LadderT3d_integrand(k, q, 2π * m / β, param)
-                @assert !isnan(integrand[ki, mi]) "nan at k=$k, q=$q"
+              Ωm = 2π * m / β
+              #integrand[ki, mi] = _LadderT3d_integrand(k, q, Ωm, param)
+                if Ωm < 1000
+                    integrand[ki, mi] = _LadderT3d_integrand(k, q, Ωm, param)
+                    @assert !isnan(integrand[ki, mi]) "nan at k=$k, q=$q"
+                else 
+                    integrand[ki, mi] = 0 # Don't integrate for high frequencies, we use abrupt two-particle expression first
+                end
             end
         end
     end
+    result = Interp.integrate1D(integrand, kgrid; axis=1)
+    for (mi, m) in enumerate(n)
+       Ωm = 2π * m / β
+        if Ωm > 100
+            result[mi] = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωm - 2μ)) # Replace with two-particle particle-particle high frequency tail
+        end
+    end
 
-    return Interp.integrate1D(integrand, kgrid; axis=1)
+    return result
+end
+
+function Ladder0_FreeElectron(q::Float64, n::AbstractVector, param)
+    @unpack dim, kF, β, me, μ = param
+    if dim != 3
+        error("No support for finite-temperature polarization in $dim dimension!")
+    end
+    result = zeros(ComplexF64, length(n))
+    for (mi, m) in enumerate(n)
+        result[mi] = -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * m / β - 2μ)
+    end
+    return result
 end
 
 
diff --git a/src/selfenergy.jl b/src/selfenergy.jl
index 64b7511..01c1768 100644
--- a/src/selfenergy.jl
+++ b/src/selfenergy.jl
@@ -255,6 +255,47 @@ function G0W0(param, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing;
     return G0W0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid; kwargs...)
 end
 
+function G0Γ0(param, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; Euv=100 * param.EF, rtol=1e-14, Nk=12, maxK=6 * param.kF, minK=1e-8 * param.kF, order=8, int_type=:none,
+    kwargs...)
+    return G0Γ0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid; kwargs...)
+end
+
+function G0Γ0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; kwargs...)
+    @unpack dim = param
+    # kernel = SelfEnergy.LegendreInteraction.DCKernel_old(param;
+    # Euv = Euv, rtol = rtol, Nk = Nk, maxK = maxK, minK = minK, order = order, int_type = int_type, spin_state = :sigma)
+    # kernel = SelfEnergy.LegendreInteraction.DCKernel0(param;
+    #     Euv = Euv, rtol = rtol, Nk = Nk, maxK = maxK, minK = minK, order = order, int_type = int_type, spin_state = :sigma)
+    # G0 = G0wrapped(Euv, rtol, kernel.kgrid, param)
+
+    kGgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order)
+    if dim == 2
+        error("No support for G0Γ0 in 2d dimension!")
+    elseif dim == 3
+        if isnothing(kgrid)
+            kernel_r = SelfEnergy.LegendreInteraction.DCKernel_Ladder_r(param;
+                Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kwargs...)
+            kernel_i = SelfEnergy.LegendreInteraction.DCKernel_Ladder_i(param;
+                Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kwargs...)
+        else
+            if (kgrid isa AbstractVector)
+                kgrid = SimpleG.Arbitrary{eltype(kgrid)}(kgrid)
+            end
+            kernel_r = SelfEnergy.LegendreInteraction.DCKernel_Ladder_r(param;
+                Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kgrid=kgrid, kwargs...)
+            kernel_i = SelfEnergy.LegendreInteraction.DCKernel_Ladder_i(param;
+                Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kgrid=kgrid, kwargs...)
+        end
+        G0 = G0wrapped(Euv, rtol, kGgrid, param)
+        Σr, Σ_ins = calcΣ_3d(G0, kernel_r)
+        Σi, Σ_ins = calcΣ_3d(G0, kernel_i)
+    else
+        error("No support for G0W0 in $dim dimension!")
+    end
+
+    return Σr, Σi
+end
+
 """
     function zfactor(param, Σ::GreenFunc.MeshArray; kamp=param.kF, ngrid=[0, 1])
     
diff --git a/test/galileaninvariance.jl b/test/galileaninvariance.jl
new file mode 100644
index 0000000..126f777
--- /dev/null
+++ b/test/galileaninvariance.jl
@@ -0,0 +1,111 @@
+# using ElectronGas
+using Parameters
+using GreenFunc
+using CompositeGrids
+# print(dirname(pwd()))
+include("./ElectronGas.jl/src/convention.jl")
+include("./ElectronGas.jl/src/parameter.jl")
+include("./ElectronGas.jl/src/polarization.jl")
+include("./ElectronGas.jl/src/interaction.jl")
+include("./ElectronGas.jl/src/twopoint.jl")
+using Plots
+
+para = Parameter.rydbergUnit(0.01, 1, 3)
+Ωn_mesh_ph = GreenFunc.ImFreq(para.β, BOSON; Euv=100.0, symmetry=:ph)
+sgrid = 0.0
+
+δladderarray_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64)
+δTarray_n2 = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64)
+
+for (ki, k) in enumerate(sgrid)
+    for (ni, n) in enumerate(Ωn_mesh_ph.grid)
+        δTarray_n2[ni, ki] = Interaction.Tmatrix(k, n, para) - Interaction.T0matrix(k, n, para)
+    end
+end
+
+
+
+
+# Many Body Particle-Particle
+mbladderarray = Polarization.Ladder0_FiniteTemp(0.0, Ωn_mesh_ph.grid, para)
+
+# Vacuum Particle-Particle
+vacladderarray = Polarization.Ladder0_FreeElectron(0.0,  Ωn_mesh_ph.grid, para)
+
+for (ki, k) in enumerate(sgrid)
+    for (ni, n) in enumerate(Ωn_mesh_ph.grid)
+        δladderarray_n[ni, ki] = mbladderarray[ni,ki] - vacladderarray[ni,ki] 
+    end
+end
+
+# for (ki, k) in enumerate(sgrid)
+#     for (ni, n) in enumerate(Ωn_mesh_ph.grid)
+#         δTarray_n[ni, ki] = Interaction.Tmatrix(k, n, para) - Interaction.T0matrix(k, n, para)
+#     end
+# end
+
+δladderarray_τ = δladderarray_n |> to_dlr |> to_imtime
+# δTarray_τ = δTarray_n |> to_dlr |> to_imtime
+
+# Analytic expression for T0(k, τ) #
+
+function T0_τ(q::Float64, τ::Float64, param)
+
+    me, as, μ, β  = param.me, param.as, param.μ, param.β
+    B = me^(3/2) / (4 * π)
+    # Generating a log densed composite grid with LogDensedGrid()
+    ygrid =  CompositeGrid.LogDensedGrid(:gauss,[0.0,50.],[0.0],5,0.005,5)
+
+    # Functions to be integrated 
+    f(y) = ( exp(- y^2 ) * y^2 ) / (τ / (me * as) + y^2)
+    g(y) = ( exp(- y^2 ) * y^2 ) / (τ / (me * as) + y^2) * 1. / ( exp( β * ( y^2/ τ + (q^2 / (4me) - 2μ)) ) - 1)
+
+    # Generate data for integration
+    fdata = [f(y) for (yi, y) in enumerate(ygrid.grid)]
+    gdata = [g(y) for (yi, y) in enumerate(ygrid.grid)]
+
+    fint_result = Interp.integrate1D(fdata, ygrid)
+    gint_result = Interp.integrate1D(gdata, ygrid)
+
+    return exp(-(q^2 / (4 * me) - 2μ) * τ) * (2/( B* π * sqrt(τ) ) * ( fint_result + gint_result ) - exp( τ / (me*as^2) ) / ( 1 - exp( β* ( 1 / (me*as^2) - (q^2 / (4me) - 2μ)) ) ))
+    # return  - exp( τ / (me*as^2) ) / ( 1 - exp( β* ( 1 / (me*as^2) - (q^2 / (4me) - 2μ)) ) )
+
+end
+
+function Σ_τ(q::Float64, param)
+    # Output self energy for the mesh imaginary values
+
+    kgrid = Polarization.finitetemp_kgrid_ladder(param.μ, param.me, q, param.kF, 20, 1e-6, 10)
+    δTarray_n = GreenFunc.MeshArray(Ωn_mesh_ph, kgrid; dtype=ComplexF64)
+
+    for (ki, k) in enumerate(kgrid)
+        for (ni, n) in enumerate(Ωn_mesh_ph.grid)
+            δTarray_n[ni, ki] = Interaction.Tmatrix(k+q, n, param) - Interaction.T0matrix(k+q, n, param)
+        end
+    end
+
+    δTarray_τ = δTarray_n |> to_dlr |> to_imtime
+    integrand = GreenFunc.MeshArray(δTarray_τ, kgrid; dtype=ComplexF64)
+    result = GreenFunc.MeshArray(δTarray_τ.mesh[1]; dtype=ComplexF64)
+
+    for (ti, t) in enumerate(δTarray_τ.mesh[1][:])
+        for (ki, k) in enumerate(kgrid.grid)
+            integrand[ti, ki] = (T0_τ(k + q, t, param) + δTarray_τ[ti]) * exp(t * ( k^2 / (2*param.me)  - param.μ)) * 1 / (exp( param.β * (k^2 / (2* param.me) - param.μ)) + 1)
+        end
+        result[ti] = Interp.integrate1D(integrand[ti, :], kgrid)
+    end
+    return result
+end
+
+qarray = LinRange(0,2*para.kF, 5)
+#test_τ = #[Σ_τ(qi, para) for qi in qarray]
+#test_ω = test_τ |> to_dlr |> to_imfreq
+#plot(test_ω.mesh[1][:], imag(test_ω))
+
+#G0_τ = GreenFunc.MeshArray(test_τ.mesh[1]; dtype=ComplexF64)
+
+#for (ti, t) in enumerate(G0_τ.mesh[1][:])
+#    G0_τ[ti] = T0_τ(0.0, t, para) #exp(t * ( 0.0^2 / (2*para.me)  - para.μ)) * 1 / (exp( para.β * (0.0^2 / (2* para.me) - para.μ)) + 1)
+
+#end
+#G0_ω = G0_τ |> to_dlr |> to_imfreq
\ No newline at end of file
diff --git a/test/ladder.jl b/test/ladder.jl
index cc985cc..70aa7f1 100644
--- a/test/ladder.jl
+++ b/test/ladder.jl
@@ -1,9 +1,9 @@
 using MCIntegration
 using ElectronGas
 
-const para = Parameter.rydbergUnit(1.0 / 10, 4.0, 3)
+const para = Parameter.rydbergUnit(1.0 / 100, 1.0, 3)
 const n = 0
-const k = 0.1 * para.kF
+const k = 0.0 * para.kF
 
 function integrand(vars, config)
     n, k, para = config.userdata