diff --git a/docs/src/tutorial/greenfunctions.md b/docs/src/tutorial/greenfunctions.md
index d5454b90..aca5362c 100644
--- a/docs/src/tutorial/greenfunctions.md
+++ b/docs/src/tutorial/greenfunctions.md
@@ -67,18 +67,21 @@ The currently implemented `GreenSolver`s (abbreviated as `GS`) are the following
 
   Uses a sparse `LU` factorization to compute the inverse of `⟨i|ω - H - Σ(ω)|j⟩`, where `Σ(ω)` is the self-energy from the contacts.
 
+
 - `GS.Spectrum(; spectrum_kw...)`
 
-  For bounded (`L=0`) AbstractHamiltonians.
+  For bounded (`L=0`) Hamiltonians. This solver does not accept ParametricHamiltonians. Convert to Hamiltonian with `h(; params...)` first.
 
   Uses a diagonalization of `H`, obtained with `spectrum(H; spectrum_kw...)`, to compute the `G⁰ᵢⱼ` using the Lehmann representation `∑ₖ⟨i|ϕₖ⟩⟨ϕₖ|j⟩/(ω - ϵₖ)`. Any eigensolver supported by `spectrum` can be used here. If there are contacts, it dresses `G⁰` using a T-matrix approach, `G = G⁰ + G⁰TG⁰`.
 
+
 - `GS.KPM(order = 100, bandrange = missing, kernel = I)`
 
   For bounded (`L=0`) Hamiltonians, and restricted to sites belonging to contacts (see the section on Contacts).
 
   It precomputes the Chebyshev momenta, and incorporates the contact self energy with a T-matrix approach.
 
+
 - `GS.Schur(boundary = Inf)`
 
   For 1D (`L=1`) AbstractHamiltonians with only nearest-cell coupling. Default for `L=1`.
@@ -86,6 +89,7 @@ The currently implemented `GreenSolver`s (abbreviated as `GS`) are the following
   Uses a deflating Generalized Schur (QZ) factorization of the generalized eigenvalue problem to compute the unit-cell self energies.
   The Dyson equation then yields the Green function between arbitrary unit cells, which is further dressed using a T-matrix approach if the lead has any attached self-energy.
 
+
 - `GS.Bands(bandsargs...; boundary = missing, bandskw...)`
 
   For unbounded (`L>0`) Hamiltonians.
@@ -94,6 +98,7 @@ The currently implemented `GreenSolver`s (abbreviated as `GS`) are the following
 
   To retrieve the bands from a `g::GreenFunction` that used the `GS.Bands` solver, we may use `bands(g)`.
 
+
 ## Attaching Contacts
 
 A self energy `Σ(ω)` acting of a finite set of sites of `h` (i.e. on a `LatticeSlice` of `lat = lattice(h)`) can be incorporated using the `attach` command. This defines a new Contact in `h`. The general syntax is `oh = attach(h, args...; sites...)`, where the `sites` directives define the Contact `LatticeSlice` (`lat[siteselector(; sites...)]`), and `args` can take a number of forms.
diff --git a/docs/src/tutorial/observables.md b/docs/src/tutorial/observables.md
index 59a29316..7d9ac13e 100644
--- a/docs/src/tutorial/observables.md
+++ b/docs/src/tutorial/observables.md
@@ -14,23 +14,81 @@ See the corresponding docstrings for full usage instructions. Here we will prese
 
 ## Local density of states (LDOS)
 
-Let us compute the LDOS in a cavity like in the previous section. Instead of computing the Green function between a contact to an arbitrary point, we can construct an object `ρ = ldos(g(ω))` without any contacts. By using a small imaginary part in `ω`, we broaden the discrete spectrum, and obtain a finite LDOS. Then, we can pass `ρ` directly as a site shader to `qplot`
+Let us compute the LDOS in a cavity like in the previous section. Instead of computing the Green function between a contact to an arbitrary point, we can construct an object `d = ldos(g(ω))` without any contacts. By using a small imaginary part in `ω`, we broaden the discrete spectrum, and obtain a finite LDOS. Then, we can pass `d` directly as a site shader to `qplot`
 ```julia
 julia> h = LP.square() |> onsite(4) - hopping(1) |> supercell(region = r -> norm(r) < 40*(1+0.2*cos(5*atan(r[2],r[1]))));
 
-julia> g = h|> greenfunction;
+julia> g = h |> greenfunction;
 
-julia> ρ = ldos(g(0.1 + 0.001im))
+julia> d = ldos(g(0.1 + 0.001im))
 LocalSpectralDensitySolution{Float64} : local density of states at fixed energy and arbitrary location
   kernel   : LinearAlgebra.UniformScaling{Bool}(true)
 
-julia> qplot(h, hide = :hops, sitecolor = ρ, siteradius = ρ, minmaxsiteradius = (0, 2), sitecolormap = :balance)
+julia> qplot(h, hide = :hops, sitecolor = d, siteradius = d, minmaxsiteradius = (0, 2), sitecolormap = :balance)
 ```
 ```@raw html
 <img src="../../assets/star_shape_ldos.png" alt="LDOS" width="400" class="center"/>
 ```
 
-Note that `ρ[sites...]` produces a vector with the LDOS at sites defined by `siteselector(; sites...)` (`ρ[]` is the ldos over all sites). We can also define a `kernel` to be traced over orbitals to obtain the spectral density of site-local observables (see `diagonal` slicing in the preceding section).
+Note that `d[sites...]` produces a vector with the LDOS at sites defined by `siteselector(; sites...)` (`d[]` is the ldos over all sites). We can also define a `kernel` to be traced over orbitals to obtain the spectral density of site-local observables (see `diagonal` slicing in the preceding section).
+
+We can also compute the convolution of the density of states with the Fermi distribution `f(ω)=1/(exp((ω-μ)/kBT) + 1)`, which yields the density matrix in thermal equilibrium, at a given temperature `kBT` and chemical potential `μ`. This is computed with `ρ = densitymatrix(gs, (ωmin, ωmax); kBT = kBT, mu = μ)`. Here `gs = g[sites...]` is a `GreenSlice`, and `(ωmin, ωmax)` are integration bounds (they should span the full bandwidth of the system). Then, `ρ(; params...)` will yield a matrix over the selected `sites` for a set of model `params`.
+```julia
+julia> ρ = densitymatrix(g[region = RP.circle(1)], (-0.1, 8.1), mu = 4)
+Integrator: Complex-plane integrator
+  Integration path    : (-0.1 + 1.4901161193847656e-8im, 1.95 + 2.050000014901161im, 4.0 + 1.4901161193847656e-8im)
+  Integration options : (atol = 1.0e-7,)
+  Integrand           : GFermi{Float64}
+    kBT               : 0.0
+
+julia> @time ρ(4)
+  5.436825 seconds (57.99 k allocations: 5.670 GiB, 1.54% gc time)
+5×5 Matrix{ComplexF64}:
+          0.5+0.0im          -7.34893e-10-3.94035e-15im  0.204478+1.9366e-14im   -7.34889e-10-1.44892e-15im  -5.70089e-10+5.48867e-15im
+ -7.34893e-10+3.94035e-15im           0.5+0.0im          0.200693-2.6646e-14im   -5.70089e-10-1.95251e-15im  -7.34891e-10-2.13804e-15im
+     0.204478-1.9366e-14im       0.200693+2.6646e-14im        0.5+0.0im              0.200693+3.55692e-14im      0.204779-4.27255e-14im
+ -7.34889e-10+1.44892e-15im  -5.70089e-10+1.95251e-15im  0.200693-3.55692e-14im           0.5+0.0im          -7.34885e-10-3.49861e-15im
+ -5.70089e-10-5.48867e-15im  -7.34891e-10+2.13804e-15im  0.204779+4.27255e-14im  -7.34885e-10+3.49861e-15im           0.5+0.0im
+```
+
+Note that the diagonal is `0.5`, indicating half-filling.
+
+The default algorithm used here is slow, as it relies on numerical integration in the complex plane. Some GreenSolvers have more efficient implementations. If they exist, they can be accessed by omitting the `(ωmin, ωmax)` argument. For example, using `GS.Spectrum`:
+```julia
+julia> @time g = h |> greenfunction(GS.Spectrum());
+ 38.024189 seconds (2.81 M allocations: 2.929 GiB, 0.33% gc time, 1.86% compilation time)
+
+julia> ρ = densitymatrix(g[region = RP.circle(1)])
+DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver
+
+julia> @time ρ(4)
+  0.000723 seconds (8 allocations: 430.531 KiB)
+5×5 Matrix{ComplexF64}:
+          0.5+0.0im  -2.21437e-15+0.0im  0.204478+0.0im   2.67668e-15+0.0im   3.49438e-16+0.0im
+ -2.21437e-15+0.0im           0.5+0.0im  0.200693+0.0im  -1.40057e-15+0.0im  -2.92995e-15+0.0im
+     0.204478+0.0im      0.200693+0.0im       0.5+0.0im      0.200693+0.0im      0.204779+0.0im
+  2.67668e-15+0.0im  -1.40057e-15+0.0im  0.200693+0.0im           0.5+0.0im   1.81626e-15+0.0im
+  3.49438e-16+0.0im  -2.92995e-15+0.0im  0.204779+0.0im   1.81626e-15+0.0im           0.5+0.0im
+```
+
+Note, however, that the computation of `g` is much slower in this case, due to the need of a full diagonalization. A better algorithm choice in this case is `GS.KPM`. It requires, however, that we define the region for the density matrix beforehand, as a `nothing` contact.
+```julia
+julia> @time g = h |> attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 10000, bandrange = (0,8)));
+Computing moments: 100%|█████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01
+  1.704793 seconds (31.04 k allocations: 11.737 MiB)
+
+julia> ρ = densitymatrix(g[1])
+DensityMatrix: density matrix on specified sites with solver of type DensityMatrixKPMSolver
+
+julia> @time ρ(4)
+  0.006521 seconds (2 allocations: 992 bytes)
+5×5 Matrix{ComplexF64}:
+         0.5+0.0im  2.15097e-17+0.0im   0.20456+0.0im  2.15097e-17+0.0im   3.9251e-17+0.0im
+ 2.15097e-17+0.0im          0.5+0.0im  0.200631+0.0im  1.05873e-16+0.0im  1.70531e-18+0.0im
+     0.20456+0.0im     0.200631+0.0im       0.5+0.0im     0.200631+0.0im      0.20482+0.0im
+ 2.15097e-17+0.0im  1.05873e-16+0.0im  0.200631+0.0im          0.5+0.0im  1.70531e-18+0.0im
+  3.9251e-17+0.0im  1.70531e-18+0.0im   0.20482+0.0im  1.70531e-18+0.0im          0.5+0.0im
+```
 
 ## Current
 
diff --git a/src/Quantica.jl b/src/Quantica.jl
index 3504c369..3f748d67 100644
--- a/src/Quantica.jl
+++ b/src/Quantica.jl
@@ -29,7 +29,7 @@ export sublat, bravais_matrix, lattice, sites, supercell,
        spectrum, energies, states, bands, subdiv,
        greenfunction, selfenergy, attach, contact, cellsites,
        plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults,
-       conductance, josephson, ldos, current, transmission
+       conductance, josephson, ldos, current, transmission, densitymatrix
 
 export LatticePresets, LP, RegionPresets, RP, HamiltonianPresets, HP
 export EigenSolvers, ES, GreenSolvers, GS
diff --git a/src/docstrings.jl b/src/docstrings.jl
index 54783e75..46d25eb8 100644
--- a/src/docstrings.jl
+++ b/src/docstrings.jl
@@ -1347,10 +1347,11 @@ Add self-energy `Σᵢⱼ(ω)` defined by a `model` composed of parametric terms
 
     attach(h, nothing; sites...)
 
-Add null self-energy `Σᵢⱼ(ω) = 0` on selected sites, which in effect simply amounts to
-defining a contact on said sites, but does not lead to any dressing the Green function. This
-is useful for some `GreenFunction` solvers such as `GS.KPM` (see `greenfunction`), which
-need to know the sites of interest beforehand (the contact sites in this case).
+Add a `nothing` contact with a null self-energy `Σᵢⱼ(ω) = 0` on selected sites, which in
+effect simply amounts to labeling those sites with a contact number, but does not lead to
+any dressing the Green function. This is useful for some `GreenFunction` solvers such as
+`GS.KPM` (see `greenfunction`), which need to know the sites of interest beforehand (the
+contact sites in this case).
 
     attach(h, g1D::GreenFunction; reverse = false, transform = identity, sites...)
 
@@ -1521,6 +1522,7 @@ possible keyword arguments are
 - `GS.SparseLU()` : Direct inversion solver for 0D Hamiltonians using a `SparseArrays.lu(hmat)` factorization
 - `GS.Spectrum(; spectrum_kw...)` : Diagonalization solver for 0D Hamiltonians using `spectrum(h; spectrum_kw...)`
     - `spectrum_kw...` : keyword arguments passed on to `spectrum`
+    - This solver does not accept ParametricHamiltonians. Convert to Hamiltonian with `h(; params...)` first. Contact self-energies that depend on parameters are supported.
 - `GS.Schur(; boundary = Inf)` : Solver for 1D Hamiltonians based on a deflated, generalized Schur factorization
     - `boundary` : 1D cell index of a boundary cell, or `Inf` for no boundaries. Equivalent to removing that specific cell from the lattice when computing the Green function.
 - `GS.KPM(; order = 100, bandrange = missing, kernel = I)` : Kernel polynomial method solver for 0D Hamiltonians
@@ -1531,7 +1533,7 @@ possible keyword arguments are
 - `GS.Bands(bands_arguments; boundary = missing, bands_kw...)`: solver based on the integration of bandstructure simplices
     - `bands_arguments`: positional arguments passed on to `bands`
     - `bands_kw`: keyword arguments passed on to `bands`
-    - `boundary`: either `missing` (no boundary), or `dir => cell_pos`, where `dir::Integer` is the Bravais vector normal to the boundary, and `cell_pos::Integer` the value of cell indices `cells[dir]` that define the boundary (i.e. `cells[dir] <= cell_pos` are vaccum)
+    - `boundary`: either `missing` (no boundary), or `dir => cell_pos` (single boundary), where `dir::Integer` is the Bravais vector normal to the boundary, and `cell_pos::Integer` the value of cell indices `cells[dir]` that define the boundary (i.e. `cells[dir] <= cell_pos` are vaccum)
     - This solver only allows zero or one boundary. WARNING: if a boundary is used, the algorithm may become unstable for very fine band meshes.
 
 """
@@ -1564,7 +1566,7 @@ julia> g(1)[diagonal(:)]                            # g diagonal on all contact
  -0.10919028169083328 - 0.08398577667508969im
    -0.109190281684109 - 0.08398577667508969im
 
-julia> g(1)[diagonal(:, kernel = SA[1 0; 0 -1])]    # σz spin density at ω = 1
+julia> g(1)[diagonal(:, kernel = SA[1 0; 0 -1])]    # σz spin spectral density at ω = 1
 2-element Vector{ComplexF64}:
   6.724287793247186e-12 + 2.7755575615628914e-17im
  -6.724273915459378e-12 + 0.0im
@@ -1638,6 +1640,57 @@ true
 """
 ldos
 
+"""
+
+    densitymatrix(gs::GreenSlice; opts...)
+
+Compute a `ρ::DensityMatrix` at thermal equilibrium on sites encoded in `gs`. The actual
+matrix for given system parameters `params`, and for a given chemical potential `mu` and
+temperature `kBT` is obtained by calling `ρ(mu = 0, kBT = 0; params...)`. The algorithm used is
+specialized for the GreenSolver used, if available. In this case, `opts` are options for
+said algorithm.
+
+    densitymatrix(gs::GreenSlice, (ωmin, ωmax); opts...)
+
+As above, but using a generic algorithm that relies on numerical integration along a contour
+in the complex plane, between points `(ωmin, ωmax)`, which should be chosen so as to
+encompass the full system bandwidth. Keywords `opts...` are passed to the `QuadGK.quadgk`
+integration routine.
+
+    densitymatrix(gs::GreenSlice, ωmax::Number; opts...)
+
+As above with `ωmin = -ωmax`.
+
+## Full evaluation
+
+    ρ(μ = 0, kBT = 0; params...)   # where ρ::DensityMatrix
+
+Evaluate the density matrix at chemical potential `μ` and temperature `kBT` (in the same
+units as the Hamiltonian) for the given `g` parameters `params`, if any.
+
+## Algorithms
+
+Currently, the following GreenSolvers implement dedicated densitymatrix algorithms:
+
+- `GS.Spectrum`: based on summation occupation-weigthed eigenvectors. No `opts`.
+- `GS.KPM`: based on the Chebyshev expansion of the Fermi function. Currently only works for zero temperature and only supports `nothing` contacts (see `attach`).
+
+# Example
+```
+julia> g = HP.graphene(a0 = 1) |> supercell(region = RP.circle(10)) |> greenfunction(GS.Spectrum());
+
+julia> ρ = densitymatrix(g[region = RP.circle(0.5)])
+DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver
+
+julia> ρ()  # with mu = kBT = 0 by default
+2×2 Matrix{ComplexF64}:
+       0.5+0.0im  -0.262865+0.0im
+ -0.262865+0.0im        0.5+0.0im
+```
+
+"""
+densitymatrix
+
 """
     current(h::AbstractHamiltonian; charge = -I, direction = 1)
 
@@ -1804,24 +1857,23 @@ true
 transmission
 
 """
-    josephson(gs::GreenSlice, ωmax; kBT = 0.0, phases = missing, imshift = missing, slope = 1, post = real, atol = 1e-7, quadgk_opts...)
+    josephson(gs::GreenSlice, ωmax; phases = missing, imshift = missing, slope = 1, post = real, atol = 1e-7, quadgk_opts...)
 
 For a `gs = g[i::Integer]` slice of the `g::GreenFunction` of a hybrid junction, build a
-partially evaluated object `J::Integrator` representing the equilibrium (static) Josephson
-current `I_J` flowing into `g` through contact `i`, integrated from `-ωmax` to `ωmax` (or
-from `-ωmax` to `0` at zero temperature `kBT = 0`). The result of `I_J` is given in units of
-`qe/h` (`q` is the dimensionless carrier charge). `I_J` can be written as ``I_J = Re ∫ dω
-f(ω) j(ω)``, where ``j(ω) = (qe/h) × 2Tr[(ΣʳᵢGʳ - GʳΣʳᵢ)τz]``.
+`J::Josephson` object representing the equilibrium (static) Josephson current `I_J` flowing
+into `g` through contact `i`, integrated from `-ωmax` to `ωmax`. The result of `I_J` is
+given in units of `qe/h` (`q` is the dimensionless carrier charge). `I_J` can be written as
+``I_J = Re ∫ dω f(ω) j(ω)``, where ``j(ω) = (qe/h) × 2Tr[(ΣʳᵢGʳ - GʳΣʳᵢ)τz]``.
 
 ## Full evaluation
 
-    J(; params...)
+    J(kBT = 0; params...)   # where J::Josephson
 
-Evaluate the Josephson current `I_J` for the given `g` parameters `params`, if any.
+Evaluate the current `I_J` at temperature `kBT` (in the same units as the Hamiltonian) for
+the given `g` parameters `params`, if any.
 
 ## Keywords
 
-- `kBT` : temperature in same energy units as the Hamiltonian
 - `phases` : collection of superconducting phase biases to apply to the contact, so as to efficiently compute the full current-phase relation `[I_J(ϕ) for ϕ in phases]`. Note that each phase bias `ϕ` is applied by a `[cis(-ϕ/2) 0; 0 cis(ϕ/2)]` rotation to the self energy, which is almost free. If `missing`, a single `I_J` is returned.
 - `imshift`: if `missing` the initial and final integration points `± ωmax` are shifted by `im * sqrt(eps(ωmax))`, to avoid the real axis. Otherwise a shift `im*imshift` is applied (may be zero if `ωmax` is greater than the bandwidth).
 - `slope`: if non-zero, the integration will be performed along a piecewise-linear path in the complex plane `(-ωmax, -ωmax/2 * (1+slope*im), 0, ωmax/2 * (1+slope*im), ωmax)`, taking advantage of the holomorphic integrand `f(ω) j(ω)` and the Cauchy Integral Theorem for faster convergence.
@@ -1837,27 +1889,20 @@ julia> glead = LP.square() |> hamiltonian(onsite(0.0005 * SA[0 1; 1 0]) + hoppin
 julia> g0 = LP.square() |> hamiltonian(hopping(SA[1 0; 0 -1]), orbitals = 2) |> supercell(region = r->-2<r[2]<2 && r[1]≈0) |> attach(glead, reverse = true) |> attach(glead) |> greenfunction;
 
 julia> J = josephson(g0[1], 4; phases = subdiv(0, pi, 10))
-Integrator: Complex-plane integrator
-  Integration path    : (-4.0 + 1.4901161193847656e-8im, -2.0 + 2.000000014901161im, 0.0 + 1.4901161193847656e-8im)
-  Integration options : (atol = 1.0e-7,)
-  Integrand:          :
-  JosephsonDensity{Float64} : Equilibrium (dc) Josephson current observable before integration over energy
-    kBT                     : 0.0
-    Contact                 : 1
-    Number of phase shifts  : 10
-
-julia> J()
+Josephson: equilibrium Josephson current at a specific contact using solver of type JosephsonIntegratorSolver
+
+julia> J(0.0)
 10-element Vector{Float64}:
- -8.240920080846876e-16
-  0.0016315088241543722
-  0.003213820056116566
-  0.004699191781509914
-  0.006042752632291697
-  0.007203835441100842
-  0.008147188939637184
-  0.008844017741697798
-  0.009272686486294741
-  9.310666393245707e-13
+ -4.497874185217862e-16
+  0.0016315088241548356
+  0.003213820056116977
+  0.004699191781510337
+  0.006042752632292157
+  0.007203835441101561
+  0.008147188939637824
+  0.008844017741699536
+  0.009272686486297674
+ -1.3101026101755857e-12
 ```
 
 # See also
diff --git a/src/greenfunction.jl b/src/greenfunction.jl
index 52b5df38..853474a7 100644
--- a/src/greenfunction.jl
+++ b/src/greenfunction.jl
@@ -31,30 +31,35 @@ default_green_solver(::AbstractHamiltonian) = GS.Bands()
 (g::GreenSlice)(; params...) = minimal_callsafe_copy(call!(g; params...))
 (g::GreenSlice)(ω; params...) = copy(call!(g, ω; params...))
 
-function call!(g::GreenFunction{T}, ω::Real; params...) where {T}
-    ω´ = retarded_omega(real_or_complex_typed(T, ω), solver(g))
+call!(g::G, ω; params...) where {T,G<:Union{GreenFunction{T},GreenSlice{T}}} =
+    call!(g, real_or_complex_promote(T, ω); params...)
+
+function call!(g::GreenFunction{T}, ω::T; params...) where {T}
+    ω´ = retarded_omega(ω, solver(g))
     return call!(g, ω´; params...)
 end
 
-function call!(g::GreenFunction{T}, ω::Complex; params...) where {T}
-    ω´ = real_or_complex_typed(T, ω)
+function call!(g::GreenFunction{T}, ω::Complex{T}; params...) where {T}
     h = parent(g)
     contacts´ = contacts(g)
     call!(h; params...)
-    Σblocks = call!(contacts´, ω´; params...)
+    Σblocks = call!(contacts´, ω; params...)
     corbs = contactorbitals(contacts´)
-    slicer = solver(g)(ω´, Σblocks, corbs)
+    slicer = solver(g)(ω, Σblocks, corbs)
     return GreenSolution(g, slicer, Σblocks, corbs)
 end
 
 call!(g::GreenSlice; params...) =
-    GreenSlice(call!(greenfunction(g); params...), slicerows(g), slicecols(g))
+    GreenSlice(call!(greenfunction(g); params...), greenindices(g)...)
+
+call!(g::GreenSlice{T}, ω::T; params...) where {T} =
+    call!(g, retarded_omega(ω, solver(parent(g))); params...)
 
-call!(g::GreenSlice{T}, ω; params...) where {T} =
-    call!(greenfunction(g), real_or_complex_typed(T, ω); params...)[slicerows(g), slicecols(g)]
+call!(g::GreenSlice{T}, ω::Complex{T}; params...) where {T} =
+    call!(greenfunction(g), ω; params...)[orbinds_or_contactinds(g)...]
 
-real_or_complex_typed(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω)
-real_or_complex_typed(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω)
+real_or_complex_promote(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω)
+real_or_complex_promote(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω)
 
 retarded_omega(ω::T, s::AppliedGreenSolver) where {T<:Real} =
     ω + im * sqrt(eps(float(T))) * needs_omega_shift(s)
@@ -85,12 +90,7 @@ call!_output(c::Contacts) = selfenergyblocks(c)
 #   We convert any index down to cellorbs to pass to slicer, except contacts (Int, Colon)
 #region
 
-Base.getindex(g::GreenFunction, i::Union{Integer,Colon,DiagIndices}, j::Union{Integer,Colon,DiagIndices} = i) =
-    GreenSlice(g, i, j)
-Base.getindex(g::GreenFunction, i, j) =
-    GreenSlice(g, sites_to_orbs(i, g), sites_to_orbs(j, g))
-Base.getindex(g::GreenFunction, i) =
-    (i´ = sites_to_orbs(i, g); GreenSlice(g, i´, i´))
+Base.getindex(g::GreenFunction, i, j = i) = GreenSlice(g, i, j)
 Base.getindex(g::GreenFunction; kw...) = g[siteselector(; kw...)]
 Base.getindex(g::GreenFunction, kw::NamedTuple) = g[siteselector(; kw...)]
 
@@ -100,11 +100,9 @@ Base.getindex(g::GreenSolution, kw::NamedTuple) = g[getindex(lattice(g); kw...)]
 Base.view(g::GreenSolution, i::Integer, j::Integer = i) = view(slicer(g), i, j)
 Base.view(g::GreenSolution, i::Colon, j::Colon = i) = view(slicer(g), i, j)
 
-Base.getindex(g::GreenSolution, i::Integer, j::Integer = i) = copy_or_missing(view(g, i, j))
-Base.getindex(g::GreenSolution, ::Colon, ::Colon = :) = copy_or_missing(view(g, :, :))
-
-copy_or_missing(::Missing) = missing
-copy_or_missing(x) = copy(x)
+# fastpath for intra and inter-contact
+Base.getindex(g::GreenSolution, i::Integer, j::Integer = i) = copy(view(g, i, j))
+Base.getindex(g::GreenSolution, ::Colon, ::Colon = :) = copy(view(g, :, :))
 
 # conversion down to CellOrbitals. See sites_to_orbs in slices.jl
 Base.getindex(g::GreenSolution, i, j) = getindex(g, sites_to_orbs(i, g), sites_to_orbs(j, g))
@@ -482,8 +480,9 @@ end
 
 function GreenSolutionCache(gω::GreenSolution{T,<:Any,L}) where {T,L}
     cache = Dict{Tuple{SVector{L,Int},SVector{L,Int},Int},Matrix{Complex{T}}}()
+    # if contacts exists, we preallocate columns for each of their sites
     if ncontacts(gω) > 0
-        gmat = gω[:]    # may be missing if solver does not support (or doesn't have) contacts
+        gmat = gω[:]
         g = parent(gω)
         h = hamiltonian(g)
         bs = blockstructure(h)
diff --git a/src/observables.jl b/src/observables.jl
index f076b05d..a6c5def1 100644
--- a/src/observables.jl
+++ b/src/observables.jl
@@ -35,13 +35,13 @@ mul_tau!(tau::Vector, g) = (g .*= tau)
 tauz_diag(i, normalsize) = ifelse(iseven(fld1(i, normalsize)), -1, 1)
 taue_diag(i, normalsize) = ifelse(iseven(fld1(i, normalsize)), 0, 1)
 
-check_contact_slice(gs) = (slicerows(gs) isa Integer && slicecols(gs) isa Integer) ||
+check_contact_slice(gs) = (rows(gs) isa Integer && cols(gs) isa Integer) ||
     argerror("Please use a Green slice of the form `g[i::Integer, j::Integer]` or `g[i::Integer]`")
 
-check_same_contact_slice(gs) = (slicerows(gs) isa Integer && slicecols(gs) === slicerows(gs)) ||
+check_same_contact_slice(gs) = (rows(gs) isa Integer && cols(gs) === rows(gs)) ||
     argerror("Please use a Green slice of the form `g[i::Integer]`")
 
-check_different_contact_slice(gs) = (slicerows(gs) isa Integer && slicecols(gs) != slicerows(gs)) ||
+check_different_contact_slice(gs) = (rows(gs) isa Integer && cols(gs) != rows(gs)) ||
     argerror("Please use a Green slice of the form `g[i::Integer, j::Integer] with `i ≠ j`")
 
 #endregion
@@ -113,7 +113,7 @@ end
 function call!(I::Integrator; params...)
     fx! = (y, x) -> (y .= call!(I.integrand, x; params...))
     result, err = quadgk!(fx!, I.result, I.points...; I.opts...)
-    result´ = I.post.(result)
+    result´ = I.post(result)  # note: post-processing is not element-wise & can be in-place
     return result´
 end
 
@@ -146,7 +146,7 @@ end
 ldos(gω::GreenSolution; kernel = I) = LocalSpectralDensitySolution(gω, kernel)
 
 function ldos(gs::GreenSlice{T}; kernel = I) where {T}
-    slicerows(gs) === slicecols(gs) ||
+    rows(gs) === cols(gs) ||
         argerror("Cannot take ldos of a GreenSlice with rows !== cols")
     return LocalSpectralDensitySlice(gs, kernel)
 end
@@ -175,10 +175,10 @@ Base.getindex(d::LocalSpectralDensitySolution{T}, i) where {T} =
 
 # fallback through LocalSpectralDensitySolution - overload to allow a more efficient path
 function call!(d::LocalSpectralDensitySlice{T}, ω; params...) where {T}
-    sites = slicerows(d.gs)
+    orbslice_or_contact = first(orbinds_or_contactinds(d.gs))
     gω = call!(parent(d.gs), ω; params...)
     l = ldos(gω; kernel = d.kernel)
-    return l[sites]
+    return l[orbslice_or_contact]
 end
 
 #endregion
@@ -186,7 +186,7 @@ end
 
 ############################################################################################
 # current: current density Jᵢⱼ(ω) as a function of a charge operator
-#   d = current(::GreenSolution[, dir]; charge)      -> d[sites...]::SparseMatrixCSC{SVector{E,T}}
+#   d = current(::GreenSolution[, dir]; charge) -> d[sites...]::SparseMatrixCSC{SVector{E,T}}
 #   d = current(::GreenSlice[, dir]; charge) -> d(ω; params...)::SparseMatrixCSC{SVector{E,T}}
 #   Computes the zero-temperature equilibrium current density Jᵢⱼ from site j to site i
 #       Jᵢⱼ(ω) = (2/h) rᵢⱼ Re Tr[(Hᵢⱼgʳⱼᵢ - gʳᵢⱼHⱼᵢ)Q]
@@ -215,10 +215,10 @@ current(gω::GreenSolution; direction = missing, charge = -I) =
     CurrentDensitySolution(gω, charge, GreenSolutionCache(gω), sanitize_direction(direction, gω))
 
 function current(gs::GreenSlice; direction = missing, charge = -I)
-    slicerows(gs) === slicecols(gs) ||
+    rows(gs) === cols(gs) ||
         argerror("Cannot currently take ldos of a GreenSlice with rows !== cols")
     g = parent(gs)
-    orbslice = sites_to_orbs(slicerows(gs), g)
+    orbslice = orbrows(gs)
     return CurrentDensitySlice(g, charge, orbslice, sanitize_direction(direction, g))
 end
 
@@ -310,8 +310,8 @@ end
 
 function conductance(gs::GreenSlice{T}; nambu = false) where {T}
     check_contact_slice(gs)
-    i = slicerows(gs)
-    j = slicecols(gs)
+    i = rows(gs)
+    j = cols(gs)
     g = parent(gs)
     ni = norbitals(contactorbitals(g), i)
     nj = norbitals(contactorbitals(g), j)
@@ -400,6 +400,73 @@ end
 #endregion
 #endregion
 
+############################################################################################
+# densitymatrix: equilibrium (static) ρ::DensityMatrix
+#   ρ = densitymatrix(g::GreenSlice, (ωmin, ωmax); opts...)
+#   ρ(mu, kBT = 0; params...) gives the DensityMatrix that is solved with an integral over (ωmin, ωmax)
+#       ρ(mu, kBT; params...) = -(1/π) Im ∫dω f(ω) g(ω; params...)
+#   ρ = densitymatrix(g::GreenSlice; opts...) uses a GreenSolver-specific algorithm
+#   Keywords opts are passed to QuadGK.quadgk for the integral or the algorithm used
+#region
+
+struct DensityMatrix{S}
+    solver::S
+end
+
+# Default solver (integration in complex plane)
+struct DensityMatrixIntegratorSolver{I}
+    ifunc::I
+end
+
+#region ## Constructors ##
+
+(ρ::DensityMatrix)(mu = 0, kBT = 0; params...) = ρ.solver(mu, kBT; params...)
+
+(s::DensityMatrixIntegratorSolver)(mu, kBT; params...) = s.ifunc(mu, kBT; params...)
+
+# redirects to specialized method
+densitymatrix(gs::GreenSlice; kw...) = densitymatrix(solver(parent(gs)), gs::GreenSlice; kw...)
+
+# generic fallback
+densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) =
+    argerror("Dedicated `densitymatrix` algorithm not implemented for $(nameof(typeof(s))). Use generic one instead.")
+
+# default integrator solver
+densitymatrix(g::GreenSlice, ωmax::Number; opts...) = densitymatrix(g, (-ωmax, ωmax); opts...)
+
+function densitymatrix(gs::GreenSlice{T}, (ωmin, ωmax)::Tuple; imshift = missing, atol = 1e-7, opts...) where {T}
+    result = similar_Matrix(gs)
+    opts´ = (; imshift, slope = 1, post = gf_to_rho!, atol, opts...)
+    integratorfunc(mu, kBT; params...) = iszero(kBT) ?
+        Integrator(result, GFermi(gs, T(mu), T(kBT)), (ωmin, mu); opts´...)(; params...) :
+        Integrator(result, GFermi(gs, T(mu), T(kBT)), (ωmin, mu, ωmax); opts´...)(; params...)
+    return DensityMatrix(DensityMatrixIntegratorSolver(integratorfunc))
+end
+
+function gf_to_rho!(x)
+    x .= x .- x'
+    x .*= -1/(2π*im)
+    return x
+end
+
+#endregion
+
+#region ## API ##
+
+struct GFermi{G<:GreenSlice,T}
+    gs::G
+    mu::T
+    kBT::T
+end
+
+function call!(gf::GFermi, ω; params...)
+    gω = call!(gf.gs, ω; params...)
+    f = fermi(ω - gf.mu, gf.kBT)
+    gω .*= f
+    return gω
+end
+
+#endregion
 
 ############################################################################################
 # josephson
@@ -429,24 +496,37 @@ struct JosephsonDensity{T<:AbstractFloat,P<:Union{Missing,AbstractArray},G<:Gree
     cisτz::Vector{Complex{T}}   # preallocated workspace
 end
 
+struct Josephson{S}
+    solver::S
+end
+
+# default solver (integration in complex plane)
+struct JosephsonIntegratorSolver{I}
+    ifunc::I
+end
+
 #region ## Constructors ##
 
-function josephson(gs::GreenSlice{T}, ωmax; kBT = 0.0, phases = missing, imshift = missing, atol = 1e-7, opts...) where {T}
+(j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...)
+
+(s::JosephsonIntegratorSolver)(kBT; params...) = s.ifunc(kBT; params...)
+
+function josephson(gs::GreenSlice{T}, ωmax; phases = missing, imshift = missing, atol = 1e-7, opts...) where {T}
     check_same_contact_slice(gs)
-    contact = slicerows(gs)
+    contact = rows(gs)
     g = parent(gs)
-    kBT´ = T(kBT)
     Σfull = similar_contactΣ(g)
     Σ = similar_contactΣ(g, contact)
     normalsize = normal_size(hamiltonian(g))
     tauz = tauz_diag.(axes(Σ, 1), normalsize)
     phases´, traces = sanitize_phases_traces(phases, T)
-    jd = JosephsonDensity(g, kBT´, contact, tauz, phases´,
+    jd(kBT) = JosephsonDensity(g, T(kBT), contact, tauz, phases´,
         traces, Σfull, Σ, similar(Σ), similar(Σ), similar(Σ), similar(tauz, Complex{T}))
-    integrator = iszero(kBT) ?
-        Integrator(traces, jd, (-ωmax, 0); imshift, slope = 1, post = real, atol, opts...) :
-        Integrator(traces, jd, (-ωmax, 0, ωmax); imshift, slope = 1, post = real, atol, opts...)
-    return integrator
+    opts´ = (; imshift, slope = 1, post = real, atol, opts...)
+    ifunc(kBT; params...) = iszero(kBT) ?
+        Integrator(traces, jd(kBT), (-ωmax, 0); opts´...)(; params...) :
+        Integrator(traces, jd(kBT), (-ωmax, 0, ωmax); opts´...)(; params...)
+    return Josephson(JosephsonIntegratorSolver(ifunc))
 end
 
 sanitize_phases_traces(::Missing, ::Type{T}) where {T} = missing, missing
diff --git a/src/show.jl b/src/show.jl
index 01595dff..345eb9ba 100644
--- a/src/show.jl
+++ b/src/show.jl
@@ -376,7 +376,7 @@ function Base.show(io::IO, I::Integrator)
     print(io, i, summary(I), "\n",
 "$i  Integration path    : $(points(I))
 $i  Integration options : $(display_namedtuple(options(I)))
-$i  Integrand:          :\n")
+$i  Integrand           : ")
     ioindent = IOContext(io, :indent => i * "  ")
     show(ioindent, integrand(I))
 end
@@ -389,19 +389,46 @@ display_namedtuple(nt::NamedTuple) = isempty(nt) ? "()" : "$nt"
 #endregion
 
 ############################################################################################
-# Josephson
+# Josephson (integrand)
 #region
 
 function Base.show(io::IO, J::JosephsonDensity)
     i = get(io, :indent, "")
-    print(io, i, summary(J), "\n",
+    print(io, summary(J), "\n",
 "$i  kBT                     : $(temperature(J))
 $i  Contact                 : $(contact(J))
 $i  Number of phase shifts  : $(numphaseshifts(J))")
 end
 
-Base.summary(::JosephsonDensity{T}) where {T} =
-    "JosephsonDensity{$T} : Equilibrium (dc) Josephson current observable before integration over energy"
+Base.summary(::JosephsonDensity{T}) where {T} = "JosephsonDensity{$T}"
+
+#endregion
+
+############################################################################################
+# DensityMatrix
+#region
+
+function Base.show(io::IO, d::DensityMatrix)
+    i = get(io, :indent, "")
+    print(io, summary(d))
+end
+
+Base.summary(::DensityMatrix{S}) where {S} =
+    "DensityMatrix: density matrix on specified sites using solver of type $(nameof(S))"
+
+#endregion
+
+############################################################################################
+# Josephson
+#region
+
+function Base.show(io::IO, j::Josephson)
+    i = get(io, :indent, "")
+    print(io, summary(j))
+end
+
+Base.summary(::Josephson{S}) where {S} =
+    "Josephson: equilibrium Josephson current at a specific contact using solver of type $(nameof(S))"
 
 #endregion
 
diff --git a/src/slices.jl b/src/slices.jl
index 137a0591..3a6e5ada 100644
--- a/src/slices.jl
+++ b/src/slices.jl
@@ -9,7 +9,7 @@ Base.getindex(lat::Lattice, ls::LatticeSlice) = ls
 
 Base.getindex(lat::Lattice, ss::SiteSelector) = lat[apply(ss, lat)]
 
-Base.getindex(lat::Lattice, ::UnboundedSiteSelector) = lat[siteselector(; cells = zerocell(lat))]
+Base.getindex(lat::Lattice, ::SiteSelectorAll) = lat[siteselector(; cells = zerocell(lat))]
 
 function Base.getindex(lat::Lattice, as::AppliedSiteSelector)
     L = latdim(lat)
@@ -287,29 +287,36 @@ end
 #       - AnyOrbitalSlice
 #       - AnyCellOrbitalsDict
 #       - AnyCellOrbitals
-# sites_to_orbs_flat: converts sites to orbitals, without site groups
+# sites_to_orbs_nogroups: converts sites to orbitals, without site groups
 #   conversion rules:
 #       - SiteSlice to an OrbitalSlice
 #region
 
+
+#region ## sites_to_orbs
+
 ## no-ops
 
 sites_to_orbs(s::AnyOrbitalSlice, _) = s
 sites_to_orbs(c::AnyCellOrbitalsDict, _) = c
 sites_to_orbs(c::AnyCellOrbitals, _) = c
 
-## convert SiteSlice -> OrbitalSliceGrouped/OrbitalSlice
+# unused
+# sites_to_orbs_nogroups(cs::CellOrbitals, _) =  cs
+
+## DiagIndices
+
+sites_to_orbs(d::DiagIndices, g) = DiagIndices(sites_to_orbs(parent(d), g), kernel(d))
 
+## convert SiteSlice -> OrbitalSliceGrouped
+
+sites_to_orbs(kw::NamedTuple, g) = sites_to_orbs(siteselector(; kw...), g)
 sites_to_orbs(s::SiteSelector, g) = sites_to_orbs(lattice(g)[s], g)
-sites_to_orbs(kw::NamedTuple, g) = sites_to_orbs(getindex(lattice(g); kw...), g)
-sites_to_orbs(i::Integer, g) = orbslice(selfenergies(contacts(g), i))
+sites_to_orbs(i::Union{Colon,Integer}, g) = orbslice(contacts(g), i)
 sites_to_orbs(l::SiteSlice, g) =
     OrbitalSliceGrouped(lattice(l), sites_to_orbs(cellsdict(l), blockstructure(g)))
 
-sites_to_orbs_flat(l::SiteSlice, g) =
-    OrbitalSlice(lattice(l), sites_to_orbs_flat(cellsdict(l), blockstructure(g)))
-
-## convert CellSitesDict to CellOrbitalsGroupedDict/CellOrbitalsDict
+## convert CellSitesDict to CellOrbitalsGroupedDict
 
 sites_to_orbs(c::CellSitesDict, g) = sites_to_orbs(c, blockstructure(g))
 
@@ -319,14 +326,6 @@ function sites_to_orbs(cellsdict::CellSitesDict{L}, os::OrbitalBlockStructure) w
     return CellOrbitalsGroupedDict(co)
 end
 
-sites_to_orbs_flat(c::CellSitesDict, g) = sites_to_orbs_flat(c, blockstructure(g))
-
-function sites_to_orbs_flat(cellsdict::CellSitesDict{L}, os::OrbitalBlockStructure) where {L}
-    # inference fails if cellsdict is empty, so we need to specify eltype
-    co = CellOrbitals{L,Vector{Int}}[sites_to_orbs_flat(cellsites, os) for cellsites in cellsdict]
-    return CellOrbitalsDict(co)
-end
-
 ## convert CellSites -> CellOrbitalsGrouped
 
 sites_to_orbs(c::CellSites, g) = sites_to_orbs(c, blockstructure(g))
@@ -338,12 +337,42 @@ function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure)
     return CellOrbitalsGrouped(cell(cs), orbinds, Dictionary(groups...))
 end
 
-function sites_to_orbs_flat(cs::CellSites, os::OrbitalBlockStructure)
+#endregion
+
+#region ## sites_to_orbs_nogroups
+
+## convert SiteSlice -> OrbitalSlice
+
+sites_to_orbs_nogroups(l::SiteSlice, g) =
+    OrbitalSlice(lattice(l), sites_to_orbs_nogroups(cellsdict(l), blockstructure(g)))
+
+## convert CellSitesDict to CellOrbitalsDict
+
+sites_to_orbs_nogroups(c::CellSitesDict, g) = sites_to_orbs_nogroups(c, blockstructure(g))
+
+function sites_to_orbs_nogroups(cellsdict::CellSitesDict{L}, os::OrbitalBlockStructure) where {L}
+    # inference fails if cellsdict is empty, so we need to specify eltype
+    co = CellOrbitals{L,Vector{Int}}[sites_to_orbs_nogroups(cellsites, os) for cellsites in cellsdict]
+    return CellOrbitalsDict(co)
+end
+
+## convert CellSites -> CellOrbitals
+
+function sites_to_orbs_nogroups(cs::CellSites, os::OrbitalBlockStructure)
     sites = siteindices(cs)
     orbinds = _orbinds(sites, os)
     return CellOrbitals(cell(cs), orbinds)
 end
 
+## convert CellOrbitalsGrouped -> CellOrbitals
+
+# unused
+# sites_to_orbs_nogroups(cs::CellOrbitalsGrouped, _) = CellOrbitals(cell(cs), orbindices(cs))
+
+#endregion
+
+#region ## CORE FUNCTIONS
+
 _groups(i::Integer, os) = [i], [flatrange(os, i)]
 _groups(::Colon, os) = _groups(siterange(os), os)
 
@@ -382,3 +411,4 @@ function _orbinds(sites, os)
 end
 
 #endregion
+#endregion
diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl
index 6f9cdbec..3e814262 100644
--- a/src/solvers/green/bands.jl
+++ b/src/solvers/green/bands.jl
@@ -633,7 +633,7 @@ function coupled_orbslice(condition, h, seedcell, dir)
     ls = grow(lat[CellSites(seedcell, :)], h)
     sc´ = [projected_cellsites(c, dir) for c in cellsdict(ls) if condition(cell(c)[dir])]
     sslice = SiteSlice(lat, CellSitesDict(sc´))
-    oslice = sites_to_orbs_flat(sslice, h)
+    oslice = sites_to_orbs_nogroups(sslice, h)
     return oslice
 end
 
diff --git a/src/solvers/green/kpm.jl b/src/solvers/green/kpm.jl
index 18f0443c..21c8e919 100644
--- a/src/solvers/green/kpm.jl
+++ b/src/solvers/green/kpm.jl
@@ -127,32 +127,6 @@ end
 
 #endregion
 
-############################################################################################
-# KPM observables
-#   h = rescaled(H) = (H - ω0)/β
-#   G_H(ω) = 1/(ω-H) = 1/(ω - β*h - ω0*I) = β⁻¹/((ω - ω0)/β - h) = β⁻¹ G_h((ω-ω0)/β)
-#region
-
-function KPMgreen(momenta::Vector{<:Matrix}, ω, (ω0, β) = (0, 1))
-    β⁻¹ = 1/β
-    ω´ =  (ω - ω0) * β⁻¹
-    g0 = zero(first(momenta))
-    for (i, μi) in enumerate(momenta)
-        g0n = β⁻¹ * KPMgreen_coefficient(i - 1, ω´)
-        g0 .+= g0n .* μi
-    end
-    return g0
-end
-
-function KPMgreen_coefficient(n, ω)
-    σ = ifelse(imag(ω) < 0, -1, 1)
-    ωc = complex(ω)
-    g0n = -2 * im * σ * cis(-n * σ * acos(ωc)) / (ifelse(iszero(n), 2, 1) * sqrt(1 - ωc^2))
-    return g0n
-end
-
-#endregion
-
 ############################################################################################
 # AppliedKPMGreenSolver
 #region
@@ -212,11 +186,87 @@ end
 
 function (s::AppliedKPMGreenSolver{T})(ω, Σblocks, corbitals) where {T}
     g0contacts = KPMgreen(s.momenta, ω, s.bandCH)
-    # We rely on TMatrixSlicer to incorporate contact self-energies
+    # We rely on TMatrixSlicer to incorporate contact self-energie, and to slice contacts
     gslicer = TMatrixSlicer(g0contacts, Σblocks, corbitals)
     return gslicer
 end
 
 #endregion
+#endregion
+
+############################################################################################
+# KPM Green
+#   h = rescaled(H) = (H - ω0)/Δ, where (ω0, Δ) = bandsCH = (center, halfwidth)
+#   G_H(ω) = 1/(ω-H) = 1/(ω - Δ*h - ω0*I) = Δ⁻¹/((ω - ω0)/Δ - h) = Δ⁻¹ G_h((ω-ω0)/Δ)
+#region
 
+function KPMgreen(momenta::Vector{<:Matrix}, ω, (ω0, Δ) = (0, 1))
+    Δ⁻¹ = 1/Δ
+    ω´ =  (ω - ω0) * Δ⁻¹
+    g0 = zero(first(momenta))
+    for (i, μi) in enumerate(momenta)
+        g0n = Δ⁻¹ * KPMgreen_coefficient(i - 1, ω´)
+        g0 .+= g0n .* μi
+    end
+    return g0
+end
+
+function KPMgreen_coefficient(n, ω)
+    σ = ifelse(imag(ω) < 0, -1, 1)
+    ωc = complex(ω)
+    g0n = -2 * im * σ * cis(-n * σ * acos(ωc)) / (ifelse(iszero(n), 2, 1) * sqrt(1 - ωc^2))
+    return g0n
+end
+
+#endregion
+
+############################################################################################
+# densitymatrix
+#   specialized DensityMatrix method for GS.KPM
+#region
+
+struct DensityMatrixKPMSolver{T,M}
+    momenta::Vector{M}
+    bandCH::Tuple{T,T}
+end
+
+## Constructor
+
+function densitymatrix(s::AppliedKPMGreenSolver, gs::GreenSlice{T}) where {T}
+    has_selfenergy(gs) && argerror("The KPM densitymatrix solver currently support only `nothing` contacts")
+    momenta = slice_momenta(s.momenta, gs)
+    solver = DensityMatrixKPMSolver(momenta, s.bandCH)
+    return DensityMatrix(solver)
+end
+
+function slice_momenta(momenta, gs)
+    r = _maybe_contactinds(gs, rows(gs))
+    c = _maybe_contactinds(gs, cols(gs))
+    momenta´ = [view(m, r, c) for m in momenta]
+    return momenta´
+end
+
+_maybe_contactinds(gs, ::Colon) = Colon()
+_maybe_contactinds(gs, i::Integer) = contactinds(gs, i)
+_maybe_contactinds(gs, _) = throw(argerror("KPM doesn't support generic indexing"))
+
+## call
+
+function (d::DensityMatrixKPMSolver)(mu, kBT; params...)
+    ω0, Δ = d.bandCH
+    kBT´ = kBT / Δ
+    mu´ = (mu - ω0) / Δ
+    if kBT´ ≈ 0
+        ϕ = acos(mu´)
+        ρ = copy(first(d.momenta)) * (1-ϕ/π)
+        for n in 1:length(d.momenta)-1
+            @. ρ += d.momenta[n+1] * 2 * (sin(π*n)-sin(n*ϕ))/(π*n)
+        end
+    else
+        throw(argerror("KPM densitymatrix currently doesn't support finite temperatures"))
+    end
+    return ρ
+end
+
+#endregion
 #endregion
diff --git a/src/solvers/green/spectrum.jl b/src/solvers/green/spectrum.jl
index 18fc52e6..5463b610 100644
--- a/src/solvers/green/spectrum.jl
+++ b/src/solvers/green/spectrum.jl
@@ -28,9 +28,12 @@ minimal_callsafe_copy(s::SpectrumGreenSlicer) = SpectrumGreenSlicer(s.ω, s.solv
 
 #region ## apply ##
 
-apply(s::GS.Spectrum, h::AbstractHamiltonian0D, ::Contacts) =
+apply(s::GS.Spectrum, h::Hamiltonian{<:Any,<:Any,0}, ::Contacts) =
     AppliedSpectrumGreenSolver(spectrum(h; s.spectrumkw...))
 
+apply(s::GS.Spectrum, h::ParametricHamiltonian, ::Contacts) =
+    argerror("Cannot use GS.Spectrum with ParametricHamiltonian. Apply parameters with h(;params...) first.")
+
 apply(::GS.Spectrum, h::AbstractHamiltonian, ::Contacts) =
     argerror("Can only use Spectrum with bounded Hamiltonians")
 
@@ -63,4 +66,41 @@ end
 
 #endregion
 
+############################################################################################
+# densitymatrix
+#   specialized DensityMatrix method for GS.Spectrum
+#region
+
+struct DensityMatrixSpectrumSolver{T,R,C}
+    es::Vector{Complex{T}}
+    psirows::R
+    psicols::C
+end
+
+## Constructor
+
+function densitymatrix(s::AppliedSpectrumGreenSolver, gs::GreenSlice)
+    # SpectrumGreenSlicer is 0D, so there is a single cellorbs in dict.
+    # If rows/cols are contacts, we need their orbrows/orbcols (unlike for gs(ω; params...))
+    i, j = only(cellsdict(orbrows(gs))), only(cellsdict(orbcols(gs)))
+    oi, oj = orbindices(i), orbindices(j)
+    es, psis = spectrum(s)
+    solver = DensityMatrixSpectrumSolver(es, _maybe_view(psis, oi), _maybe_view(psis, oj))
+    return DensityMatrix(solver)
+end
+
+## API
+
+## call
+
+function (d::DensityMatrixSpectrumSolver)(mu, kBT; params...)
+    vi = d.psirows
+    vj´ = d.psicols' .* fermi.(d.es .- mu, kBT)
+    return vi * vj´
+end
+
+# if the orbindices cover all the unit cell, use matrices instead of
+_maybe_view(m, oi) = length(oi) == size(m, 1) ? m : view(m, oi, :)
+
+#endregion
 #endregion
diff --git a/src/solvers/selfenergy/generic.jl b/src/solvers/selfenergy/generic.jl
index 9cc3ed3e..1e2b4c0b 100644
--- a/src/solvers/selfenergy/generic.jl
+++ b/src/solvers/selfenergy/generic.jl
@@ -35,9 +35,9 @@ end
 #region ## API ##
 
 function SelfEnergy(hparent::AbstractHamiltonian, gslice::GreenSlice, model::AbstractModel; sites...)
-    slicerows(gslice) === slicecols(gslice) ||
+    rows(gslice) === cols(gslice) ||
         argerror("To attach a Greenfunction with `attach(h, g[cols, rows], coupling; ...)`, we must have `cols == rows`")
-    lsbath = sites_to_orbs(latslice(parent(gslice), slicerows(gslice)), parent(gslice))
+    lsbath = orbrows(gslice)
     lat0bath = lattice0D(lsbath)
     lsparent = sites_to_orbs(getindex(lattice(hparent); sites...), hparent)
     lat0parent = lattice0D(lsparent)
diff --git a/src/solvers/selfenergy/nothing.jl b/src/solvers/selfenergy/nothing.jl
index 0602987f..b7137d7d 100644
--- a/src/solvers/selfenergy/nothing.jl
+++ b/src/solvers/selfenergy/nothing.jl
@@ -21,6 +21,8 @@ call!(s::SelfEnergyEmptySolver, ω; params...) = s.emptymat
 
 call!_output(s::SelfEnergyEmptySolver) = s.emptymat
 
+has_selfenergy(::SelfEnergyEmptySolver) = false
+
 minimal_callsafe_copy(s::SelfEnergyEmptySolver) = s
 
 #endregion
diff --git a/src/types.jl b/src/types.jl
index e38f8933..437eaf46 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -187,7 +187,7 @@ struct SiteSelector{F,S,C}
     cells::C
 end
 
-const UnboundedSiteSelector = SiteSelector{Missing,Missing,Missing}
+const SiteSelectorAll = SiteSelector{Missing,Missing,Missing}
 
 struct AppliedSiteSelector{T,E,L}
     lat::Lattice{T,E,L}
@@ -1736,6 +1736,10 @@ solver(Σ::SelfEnergy) = Σ.solver
 
 plottables(Σ::SelfEnergy) = Σ.plottables
 
+has_selfenergy(s::SelfEnergy) = has_selfenergy(solver(s))
+has_selfenergy(s::AbstractSelfEnergySolver) = true
+# see nothing.jl for override for the case of SelfEnergyEmptySolver
+
 call!(Σ::SelfEnergy; params...) = SelfEnergy(call!(Σ.solver; params...), Σ.orbslice)
 call!(Σ::SelfEnergy, ω; params...) = call!(Σ.solver, ω; params...)
 
@@ -1850,6 +1854,8 @@ offsets(c::ContactOrbitals) = c.offsets
 selfenergies(c::Contacts) = c.selfenergies
 selfenergies(c::Contacts, i::Integer) = check_contact_index(i, c) && c.selfenergies[i]
 
+has_selfenergy(c::Contacts) = any(has_selfenergy, selfenergies(c))
+
 # c::Union{Contacts,ContactOrbitals} here
 check_contact_index(i, c) = 1 <= i <= ncontacts(c) ||
     argerror("Cannot get contact $i, there are $(ncontacts(c)) contacts")
@@ -1934,23 +1940,36 @@ struct DiagIndices{K,V}   # represents Green indices to return only diagonal ele
     kernel::K
 end
 
+struct GreenIndices{I,R}
+    inds::I
+    orbinds::R  # orbinds = sites_to_orbs(inds)
+end
+
 # Obtained with gs = g[; siteselection...]
 # Alows call!(gs, ω; params...) or gs(ω; params...)
 #   required to do e.g. h |> attach(g´[sites´], couplingmodel; sites...)
-struct GreenSlice{T,E,L,G<:GreenFunction{T,E,L},R,C}
+struct GreenSlice{T,E,L,G<:GreenFunction{T,E,L},R<:GreenIndices,C<:GreenIndices}
     parent::G
     rows::R
     cols::C
-    function GreenSlice{T,E,L,G,R,C}(parent, rows, cols) where {T,E,L,G<:GreenFunction{T,E,L},R,C}
-        if rows isa DiagIndices || cols isa DiagIndices
-            rows === cols || argerror("Diagonal indices should be identical for rows and columns")
-        end
-        return new(parent, rows, cols)
+end
+
+#region ## Constuctors ##
+
+function GreenSlice(parent, rows, cols)
+    if rows isa DiagIndices || cols isa DiagIndices
+        rows === cols || argerror("Diagonal indices should be identical for rows and columns")
     end
+    rows´ = greenindices(rows, parent)
+    cols´ = cols === rows ? rows´ : greenindices(cols, parent)
+    return GreenSlice(parent, rows´, cols´)
 end
 
-GreenSlice(parent::G, rows::R, cols::C) where {T,E,L,G<:GreenFunction{T,E,L},R,C} =
-    GreenSlice{T,E,L,G,R,C}(parent, rows, cols)
+# see sites_to_orbs in slices.jl
+greenindices(inds, g) = GreenIndices(inds, sites_to_orbs(inds, g))
+greenindices(g::GreenSlice) = g.rows, g.cols
+
+#endregion
 
 #region ## API ##
 
@@ -1961,20 +1980,21 @@ Base.parent(i::DiagIndices) = i.inds
 
 kernel(i::DiagIndices) = i.kernel
 
-hamiltonian(g::GreenFunction) = hamiltonian(g.parent)
-hamiltonian(g::GreenSolution) = hamiltonian(g.parent)
+hamiltonian(g::Union{GreenFunction,GreenSolution,GreenSlice}) = hamiltonian(g.parent)
 
-lattice(g::GreenFunction) = lattice(g.parent)
-lattice(g::GreenSolution) = lattice(g.parent)
+lattice(g::Union{GreenFunction,GreenSolution,GreenSlice}) = lattice(g.parent)
 
 latslice(g::GreenFunction, i) = orbslice(g.contacts, i)
-function latslice(g::GreenFunction, ls::LatticeSlice)
-    lattice(g) === lattice(ls) || internalerror("latslice: parent lattice mismatch")
-    return ls
-end
+
+# function latslice(g::GreenFunction, ls::LatticeSlice)
+#     lattice(g) === lattice(ls) || internalerror("latslice: parent lattice mismatch")
+#     return ls
+# end
 # latslice(g::GreenFunction, is::SiteSelector) = lattice(g)[is]
 # latslice(g::GreenFunction; kw...) = latslice(g, siteselector(; kw...))
 
+zerocell(g::Union{GreenFunction,GreenSolution,GreenSlice}) = zerocell(lattice(g))
+
 solver(g::GreenFunction) = g.solver
 
 contacts(g::GreenFunction) = g.contacts
@@ -1987,23 +2007,38 @@ slicer(g::GreenSolution) = g.slicer
 
 selfenergies(g::GreenSolution) = g.contactΣs
 
+has_selfenergy(g::Union{GreenFunction,GreenSlice,GreenSolution}) =
+    has_selfenergy(contacts(g))
+
 contactorbitals(g::GreenFunction) = contactorbitals(g.contacts)
 contactorbitals(g::GreenSolution) = g.contactorbs
+contactorbitals(g::GreenSlice) = contactorbitals(parent(g))
 
 blockstructure(g::GreenFunction) = blockstructure(hamiltonian(g))
 blockstructure(g::GreenSolution) = blockstructure(hamiltonian(g))
+blockstructure(g::GreenSlice) = blockstructure(parent(g))
 
 norbitals(g::GreenFunction) = norbitals(g.parent)
 norbitals(g::GreenSlice) = norbitals(g.parent.parent)
 
 contactinds(g::GreenFunction, i...) = contactinds(contacts(g), i...)
-contactinds(g::GreenSolution, i...) = contactinds(contactorbitals(g), i...)
+contactinds(g::Union{GreenSolution,GreenSlice}, i...) = contactinds(contactorbitals(g), i...)
 
 greenfunction(g::GreenSlice) = g.parent
 
-slicerows(g::GreenSlice) = g.rows
+rows(g::GreenSlice) = g.rows.inds
+
+cols(g::GreenSlice) = g.cols.inds
 
-slicecols(g::GreenSlice) = g.cols
+orbrows(g::GreenSlice) = g.rows.orbinds
+
+orbcols(g::GreenSlice) = g.cols.orbinds
+
+# ifelse(rows && cols are contacts, (rows, cols), (orbrows, orbcols))
+# I.e: if rows, cols are contact indices retrieve them instead of orbslices.
+orbinds_or_contactinds(g) = orbinds_or_contactinds(rows(g), cols(g), orbrows(g), orbcols(g))
+orbinds_or_contactinds(r::Union{Colon,Integer}, c::Union{Colon,Integer}, _, _) = (r, c)
+orbinds_or_contactinds(_, _, or, oc) = (or, oc)
 
 Base.parent(g::GreenFunction) = g.parent
 Base.parent(g::GreenSolution) = g.parent
@@ -2015,6 +2050,12 @@ Base.size(g::GreenSolution, i...) = size(g.parent, i...)
 flatsize(g::GreenFunction, i...) = flatsize(g.parent, i...)
 flatsize(g::GreenSolution, i...) = flatsize(g.parent, i...)
 
+function similar_Matrix(gs::GreenSlice{T}) where {T}
+    m = norbitals(orbrows(gs))
+    n = norbitals(orbcols(gs))
+    return Matrix{Complex{T}}(undef, m, n)
+end
+
 copy_lattice(g::GreenFunction) = GreenFunction(copy_lattice(g.parent), g.solver, g.contacts)
 copy_lattice(g::GreenSolution) = GreenSolution(
     copy_lattice(g.parent), g.slicer, g.contactΣs, g.contactbs)
diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl
index 370eff64..25c8bcc7 100644
--- a/test/test_greenfunction.jl
+++ b/test/test_greenfunction.jl
@@ -11,8 +11,8 @@ function testgreen(h, s; kw...)
     z = zero(SVector{L,Int})
     o = Quantica.unitvector(1, SVector{L,Int})
     conts = ntuple(identity, ncontacts(h))
-    locs = (cellsites(z, :), cellsites(z, 2), cellsites(z, 1:2), cellsites(z, (1,2)),
-            cellsites(o, :), CellOrbitals(o, 1), CellOrbitals(z, 1:2), CellOrbitals(z, SA[2,1]),
+    locs = (cellsites(z, :), cellsites(z, 2), cellsites(o, 1:2), cellsites(o, (1,2)),
+            CellOrbitals(o, 1), CellOrbitals(z, 1:2), CellOrbitals(z, SA[2,1]),
             conts...)
     for loc in locs, loc´ in locs
         gs = g[loc, loc´]
@@ -97,6 +97,9 @@ end
     h = supercell(h)
     @test_throws ArgumentError greenfunction(h, GS.Schur())
     @test_throws ArgumentError greenfunction(h, GS.Bands())
+    h = LP.honeycomb() |> @onsite((; o = 1) -> o*I) |> supercell
+    @test_throws ArgumentError greenfunction(h, GS.Spectrum())
+    @test greenfunction(h(), GS.Spectrum()) isa GreenFunction
 end
 
 @testset "greenfunction KPM" begin
@@ -258,6 +261,7 @@ function testjosephson(g0)
 end
 
 @testset "greenfunction observables" begin
+    # single-orbital vs two-orbital
     g1 = LP.square() |> supercell((1,0), region = r->-2<r[2]<2) |> hamiltonian(@hopping((r, dr; B = 0.1) -> I * cis(B * dr' * SA[r[2],-r[1]])), orbitals = 1) |> greenfunction(GS.Schur(boundary = 0));
     g2 = LP.square() |> supercell((1,0), region = r->-2<r[2]<2) |> hamiltonian(@hopping((r, dr; B = 0.1) -> I * cis(B * dr' * SA[r[2],-r[1]])), orbitals = 2) |> greenfunction(GS.Schur(boundary = 0));
     J1 = current(g1[cells = SA[1]])
@@ -265,6 +269,25 @@ end
     @test size(J1(0.2)) == size(J2(0.2)) == (3, 3)
     @test 2*J1(0.2; B = 0.1) ≈ J2(0.2; B = 0.1)
 
+    ρ = densitymatrix(g1[cells = SA[1]], 5)
+    @test all(≈(0.5), diag(ρ(0, 0; B=0.3))) # half filling
+    ρ = densitymatrix(g1[cells = SA[1]], 7)
+    @test all(<(0.96), real(diag(ρ(4, 1; B=0.1)))) # thermal depletion
+    h = LP.honeycomb() |> hopping(1) |> supercell(region = RP.circle(10))
+    reg = (; region = RP.circle(0.5))
+    gLU = h |> greenfunction(GS.SparseLU());
+    gSpectrum = h |> greenfunction(GS.Spectrum());
+    gKPM = h |> attach(nothing; reg...) |> greenfunction(GS.KPM(order = 100000, bandrange = (-3,3)));
+    ρ1, ρ2, ρ3 = densitymatrix(gLU[reg], (-3,3)), densitymatrix(gSpectrum[reg]), densitymatrix(gKPM[1])
+    @test ρ1() ≈ ρ2() atol = 0.00001
+    @test ρ2() ≈ ρ3() atol = 0.00001
+    gLU´ = h |> attach(nothing; reg...) |> greenfunction(GS.SparseLU());
+    ρ1´ = densitymatrix(gLU´[1], (-3, 3))
+    @test ρ1() ≈ ρ1´()
+    gSpectrum´ = h |> attach(nothing; reg...) |> greenfunction(GS.Spectrum());
+    ρ2´ = densitymatrix(gSpectrum´[1])
+    @test ρ2() ≈ ρ2´()
+
     glead = LP.square() |> hamiltonian(hopping(1)) |> supercell((0,1), region = r -> -1 <= r[1] <= 1) |> attach(nothing; cells = SA[10]) |> greenfunction(GS.Schur(boundary = 0));
     contact1 = r -> r[1] ≈ 5 && -1 <= r[2] <= 1
     contact2 = r -> r[2] ≈ 5 && -1 <= r[1] <= 1
diff --git a/test/test_show.jl b/test/test_show.jl
index 7f166c03..a9886274 100644
--- a/test/test_show.jl
+++ b/test/test_show.jl
@@ -35,4 +35,10 @@
     h = first(hs)
     g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing))
     @test nothing === show(stdout, josephson(g[1], 2))
+    @test nothing === show(stdout, densitymatrix(g[1], 2))
+    h = supercell(h, 3) |> supercell
+    g = greenfunction(supercell(h) |> attach(nothing), GS.KPM())
+    @test nothing === show(stdout, densitymatrix(g[1]))
+    g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing), GS.Spectrum())
+    @test nothing === show(stdout, densitymatrix(g[1]))
 end