diff --git a/src/Quantica.jl b/src/Quantica.jl
index 70794573..b1b384b6 100644
--- a/src/Quantica.jl
+++ b/src/Quantica.jl
@@ -28,7 +28,7 @@ export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian,
        plusadjoint, neighbors, siteselector, hopselector, diagonal, sitepairs,
        unflat, torus, transform, translate, combine,
        spectrum, energies, states, bands, subdiv,
-       greenfunction, selfenergy, attach,
+       greenfunction, selfenergy, attach, matsubara,
        plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults,
        conductance, josephson, ldos, current, transmission, densitymatrix,
        OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes, siteindexdict,
diff --git a/src/docstrings.jl b/src/docstrings.jl
index c8d5d115..e1927853 100644
--- a/src/docstrings.jl
+++ b/src/docstrings.jl
@@ -2054,6 +2054,13 @@ are passed to the `QuadGK.quadgk` integration routine. See below for additiona `
 
 As above with `ωmin = -ωmax`.
 
+    densitymatrix(gs::GreenSlice, matsubara(approximate_half_bandwidth::Real); opts...)
+
+As above, but using a vertical integration path from `µ + 0im` to `µ + Inf*im`. At the
+moment, this method does not support finite temperatures. The integral is typically more
+efficient if `approximate_half_bandwidth` is approximately the half-bandwidth of the system,
+but the result should not depend on this number.
+
 ## Full evaluation
 
     ρ(μ = 0, kBT = 0; params...)   # where ρ::DensityMatrix
@@ -2215,6 +2222,18 @@ Return the vertices of the polygonal integration path used to compute `O(args...
 """
 path
 
+"""
+    matsubara(half_bandwdith::Real)
+
+Represent an integration method along a vertical path in the complex-ω plane, from the real
+axis to imaginary infinity. This is intended to be used with `densitymatrix`, and currently
+only supports zero temperatures.
+
+# See also
+    `densitymatrix`
+"""
+matsubara
+
 
 """
     OrbitalSliceArray <: AbstractArray
diff --git a/src/observables.jl b/src/observables.jl
index 2728782a..bc5f1ba2 100644
--- a/src/observables.jl
+++ b/src/observables.jl
@@ -4,8 +4,8 @@
 
 abstract type IndexableObservable end  # any type that can be sliced into (current and ldos)
 
-fermi(ω::C, β = Inf; atol = sqrt(eps(real(C)))) where {C} =
-    isinf(β) ? ifelse(abs(ω) < atol, C(0.5), ifelse(real(ω) <= 0, C(1), C(0))) : C(1/(exp(β * ω) + 1))
+fermi(ω::C, β = Inf) where {C} =
+    isinf(β) ? ifelse(real(ω) <= 0, C(1), C(0)) : C(1/(exp(β * ω) + 1))
 
 normal_size(h::AbstractHamiltonian) = normal_size(blockstructure(h))
 
@@ -79,11 +79,9 @@ end
 
 #region ## Constructor ##
 
-function Integrator(result, f, pts; imshift = missing, post = identity, slope = 0, callback = Returns(nothing), quadgk_opts...)
-    imshift´ = imshift === missing ?
-        sqrt(eps(promote_type(typeof.(float.(real.(pts)))...))) : float(imshift)
+function Integrator(result, f, pts; applyshifts = true, imshift = missing, post = identity, slope = 0, callback = Returns(nothing), quadgk_opts...)
     sanitize_integration_points(pts)
-    pts´ = apply_complex_shifts(pts, imshift´, slope)
+    pts´ = maybe_apply_complex_shifts(pts, imshift, slope)
     quadgk_opts´ = NamedTuple(quadgk_opts)
     return Integrator(f, pts´, result, quadgk_opts´, callback, post)
 end
@@ -104,14 +102,14 @@ function sanitize_integration_points(pts)
     return pts
 end
 
-# If all points are real, apply sawtooth with slope
-apply_complex_shifts(pts::NTuple{<:Any,Real}, imshift, slope) =
-    iszero(slope) ? pts .+ (im * imshift) : triangular_sawtooth(imshift, slope, pts)
-apply_complex_shifts(pts::AbstractVector{<:Real}, imshift, slope) =
-    iszero(slope) ? (pts .+ (im * imshift)) : triangular_sawtooth(imshift, slope, pts)
+maybe_apply_complex_shifts(pts, imshift, slope) = pts
 
-# Otherwise do not modify
-apply_complex_shifts(pts, imshift, slope) = pts
+function maybe_apply_complex_shifts(pts::Union{NTuple{<:Any,Real},AbstractVector{<:Real}}, imshift, slope::Number)
+    imshift´ = imshift === missing ?
+        sqrt(eps(promote_type(typeof.(float.(real.(pts)))...))) : float(imshift)
+    pts´ = iszero(slope) ? (pts .+ (im * imshift´)) : triangular_sawtooth(imshift´, slope, pts)
+    return pts´
+end
 
 triangular_sawtooth(is, sl, ωs::Tuple) = _triangular_sawtooth(is, sl, (), ωs...)
 _triangular_sawtooth(is, sl, ::Tuple{}, ω1, ωs...) = _triangular_sawtooth(is, sl, (ω1 + im * is,), ωs...)
@@ -469,11 +467,13 @@ end
 #region
 
 # this produces gs(ω; params...) * f(ω-mu). Use post = gf_to_rho! after integration
-struct DensityMatrixIntegrand{G<:GreenSlice,T,O}
+struct DensityMatrixIntegrand{G<:GreenSlice,T,O,F1,F2}
     gs::G
     mu::T
     kBT::T
-    omegamap::O
+    omegamap::O                 # to change system parameters for each ω
+    omega_transform::F1         # to apply to ω at each point
+    integrand_transform!::F2    # to apply to integrand at each point
 end
 
 # Default solver (integration in complex plane)
@@ -488,6 +488,14 @@ end
 
 #region ## Constructors ##
 
+# default omega_transform, integrand_transform! and omegamap
+DensityMatrixIntegrand(gs, mu, kBT, omegamap = identity, omega_transform = identity) =
+    DensityMatrixIntegrand(gs, mu, kBT, omegamap, omega_transform, identity)
+
+#endregion
+
+#region ## Call ##
+
 # generic fallback (for other solvers)
 (ρ::DensityMatrix)(mu = 0, kBT = 0; params...) =
     ρ.solver(mu, kBT; params...)
@@ -509,6 +517,49 @@ densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) =
 # default integrator solver
 densitymatrix(gs::GreenSlice, ωmax::Number; opts...) = densitymatrix(gs, (-ωmax, ωmax); opts...)
 
+#endregion
+
+#region ## API ##
+
+(gf::DensityMatrixIntegrand)(ω; params...) = copy(call!(gf, ω; params...))
+
+function call!(gf::DensityMatrixIntegrand, ω; params...)
+    ω´ = gf.omega_transform(ω)
+    output = call!(gf.gs, ω´; gf.omegamap(ω´)..., params...)
+    serialize(output) .*= fermi(ω´ - gf.mu, inv(gf.kBT))
+    gf.integrand_transform!(output)
+    return output
+end
+
+function gf_to_rho!(x::AbstractMatrix)
+    x .= x .- x'
+    x .*= -1/(2π*im)
+    return x
+end
+
+# For diagonal indexing
+function gf_to_rho!(x::AbstractVector)
+    x .= x .- conj.(x)
+    x .*= -1/(2π*im)
+    return x
+end
+
+integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = integrand(ρ.solver(mu, kBT))
+
+path(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = path(ρ.solver(mu, kBT))
+
+temperature(D::DensityMatrixIntegrand) = D.kBT
+
+chemicalpotential(D::DensityMatrixIntegrand) = D.mu
+
+Base.parent(ρ::DensityMatrix) = ρ.gs
+
+call!_output(ρ::DensityMatrix) = call!_output(ρ.gs)
+
+#endregion
+
+#region ## Default integration solver
+
 function densitymatrix(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), imshift = missing, atol = 1e-7, opts...) where {T}
     # check_nodiag_axes(gs)
     result = copy(call!_output(gs))
@@ -561,42 +612,47 @@ override_path!(override_path, ptsvec, pts) =
 
 #endregion
 
-#region ## API ##
+#region ## Alternative integration solvers - Matsubara
 
-(gf::DensityMatrixIntegrand)(ω; params...) = copy(call!(gf, ω; params...))
-
-function call!(gf::DensityMatrixIntegrand, ω; params...)
-    output = call!(gf.gs, ω; gf.omegamap(ω)..., params...)
-    serialize(output) .*= fermi(ω - gf.mu, inv(gf.kBT))
-    return output
+struct Matsubara{T<:AbstractFloat}
+    half_bandwidth::T
 end
 
-function gf_to_rho!(x::AbstractMatrix)
-    x .= x .- x'
-    x .*= -1/(2π*im)
-    return x
-end
+matsubara(h::Real) = Matsubara(float(h))
+matsubara(h...) = argerror("Expected matsubara(half_bandwidth::Real)")
 
-# For diagonal indexing
-function gf_to_rho!(x::AbstractVector)
-    x .= x .- conj.(x)
-    x .*= -1/(2π*im)
-    return x
+function densitymatrix(gs::GreenSlice{T}, c::Matsubara; omegamap = Returns((;)), imshift = missing, atol = 1e-7, opts...) where {T}
+    result = copy(call!_output(gs))
+    opts´ = (; imshift, slope = nothing, post = cauchy_post_transform, atol, opts...)
+    hW = c.half_bandwidth
+    hW > 0 || argerror("The estimated half-bandwidth for Matsubara integration should be positive, got $hW")
+    integrand_transform! = cauchy_integrand_tranform(hW)
+    function ifunc(mu, kBT, override_path)
+        iszero(kBT) || argerror("The Matsubara integrator currently requires zero temperature, got $kBT")
+        override_path === missing || argerror("The Matsubara integrator doesn't allow path overrides")
+        omega_transform = cauchy_omega_tranform(mu, hW)
+        ρd = DensityMatrixIntegrand(gs, T(mu), zero(T), omegamap, omega_transform, integrand_transform!)
+        pts = (zero(T), T(Inf))
+        return Integrator(result, ρd, pts; opts´...)
+    end
+    return DensityMatrix(DensityMatrixIntegratorSolver(ifunc), gs)
 end
 
-integrand(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = integrand(ρ.solver(mu, kBT))
-
-path(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver}, mu = 0.0, kBT = 0.0) = path(ρ.solver(mu, kBT))
-
-temperature(D::DensityMatrixIntegrand) = D.kBT
-
-chemicalpotential(D::DensityMatrixIntegrand) = D.mu
+cauchy_omega_tranform(mu, hW) = x -> mu + hW * x * im
 
-Base.parent(ρ::DensityMatrix) = ρ.gs
+cauchy_integrand_tranform(hW) = x -> begin
+    serialize(x) .*= -im*hW
+    gf_to_rho!(x)
+    return x
+end
 
-call!_output(ρ::DensityMatrix) = call!_output(ρ.gs)
+function cauchy_post_transform(x)
+    return x + (I*0.5)
+end
 
 #endregion
+#endregion
+
 
 ############################################################################################
 # josephson