Solvers

KPM_Tn(H_mpo, N, sites; scale=nothing, center=0.0, maxdim=40,
       dmrg_nsweeps, dmrg_maxdim, dmrg_linkdim) -> (Tn_list, scale, center)

Build the list of Chebyshev MPOs T_n((H−center·I)/scale) for n = 0…N.

Scale / center arguments

• If scale=nothing (default): spectral bounds estimated automatically via _estimate_spectral_bounds (two short DMRG runs).
• If scale is provided: used directly; center defaults to 0.0 but can be set explicitly for non-symmetric spectra.

High-level overload

Pass a TBHamiltonian as the first argument to skip manual rescaling entirely: Tn, scale, center = KPM_Tn(H, Ncheb; maxdim=100) H.scale and H.center are computed lazily on the first call and cached.

Return value

Returns (Tn_list, scale, center). To convert a physical energy ω: ω_r = (ω − center) / scale ∈ (−1, 1)

KPM_Tn(H::TBHamiltonian, Ncheb::Int, X::Int; maxdim, cutoff, verbose, ...)
    -> (Tn_mps_list, scale, center)

MPS-based Chebyshev expansion for the exciton state |X,X⟩ on the exciton TBHamiltonian. X ∈ {1, …, 2^L} is the 1-indexed exciton site.

Handles rescaling automatically from H.scale/H.center, caches the result in H._tn_mps_cache and H._tn_Ncheb.

No MPO-mode is provided for exciton systems: the 2L-site MPO chain makes the MPO Chebyshev variant prohibitively expensive.

KPM_Tn(H::TBHamiltonian, Ncheb; mode=:mpo, psi0=nothing,
       maxdim=40, cutoff=1e-8, dmrg_nsweeps, dmrg_maxdim, dmrg_linkdim)
    -> (Tn_list, scale, center)

High-level Chebyshev expansion for a TBHamiltonian.

Lazily determines H.scale and H.center via DMRG if not already set, builds the rescaled Chebyshev list, caches the result on H, and returns (Tn_list, H.scale, H.center).

mode keyword

mode=:mps requires psi0 (a reference MPS). The MPS pathway is more memory-efficient when a single reference state is sufficient.

Overview of the three KPM pathways

Pathway 1 — MPO × MPO cache  [legacy / rarely used]
  KPM_Tn(H, Ncheb; mode=:mpo)           # build and cache {T_n(H̃)} MPOs
  → get_ldos_spectrum(H, ωlist)          # all ω in one pass → Vector{MPS}
  → get_ldos(H, ω; mode=:diag)           # single ω → diagonal site-LDOS MPS
  → get_ldos(H, ω; mode=:mpo)            # single ω → full off-diagonal MPO

  ⚠ Bond dimension χ_T grows at each MPO × MPO step; memory scales as
  O(Ncheb × χ_T²).  Prefer Pathways 3 or 4 unless the cache is reused for
  multiple downstream calls.  Kept mainly for legacy compatibility.

Pathway 2 — MPS cache  [legacy / fixed reference state]
  KPM_Tn(H, Ncheb; mode=:mps, psi0=ψ₀)  # cache {T_n(H̃)|ψ₀⟩} MPS on H
  → get_ldos(H, ω; mode=:mps, psi0=ψ₀)  # μₙ = ⟨ψ₀|T_n(H̃)|ψ₀⟩ → scalar

  Propagates a single reference MPS and stores the full trajectory
  {|φ_n⟩ = T_n(H̃)|ψ₀⟩} for repeated re-use across many energy queries on the
  same state.  Kept for advanced workflows; the public LDOS helpers below avoid
  storing this cache.

Pathway 3 — Online MPO × MPO  [k-space and spatial spectral functions]
  get_bands(H, Ncheb, D, ωlist; …)       # k-resolved A(k,ω), QFT-conjugated
  get_ldos_spatial(H, Ncheb, ωlist; …)   # real-space LDOS heatmap (default mode)

  Runs the MPO × MPO Chebyshev recursion online with no prior KPM_Tn call;
  only 3 MPOs alive at a time (truncated after each step).  Preferred for
  computing band structures and spatial LDOS over many positions simultaneously.

Pathway 4 — Online MPO × MPS  [single-particle default, most memory-efficient]
  get_ldos_online(H, Ncheb, X, ωlist; …) # LDOS at one site, all ω
  get_ldos_spatial(H, Ncheb, ωlist; mode=:mps; …)  # per-position MPS recursion
  get_exciton_ldos_spatial(H, Ncheb, ωlist; …)      # bound-pair exciton LDOS
  get_exciton_ldos(H, X, ωlist; …)                  # one-position wrapper
  get_dos_stochastic(H, Ncheb, ωlist; …) # stochastic trace DOS

  Propagates MPS states rather than full MPOs: only 3 MPS alive per sample/site.
  For single-particle problems this is almost always the best choice —
  memory cost is O(χ_H × χ_ψ) instead of O(χ_T²).  get_dos_stochastic
  applies this with random initial states for a stochastic trace estimate.

All pathways share the same KPM kernel and normalization. Aux-DOF projections (spin_proj, nambu_proj, sublat_proj, …) are supported in Pathways 1, 3, and 4. Pathway 2 operates on a fixed reference state and does not expose per-DOF projection keywords.

KPM_Tn_mps(H_mpo, N, psi0, sites; scale=nothing, center=0.0, maxdim=40,
           dmrg_nsweeps, dmrg_maxdim, dmrg_linkdim, cutoff, verbose)
-> (Tn_mps_list, scale, center)

MPS-based Chebyshev expansion. Instead of storing Chebyshev MPOs Tn(H) (as `KPMTn` does), this builds the projected MPS states

|φ_n⟩ = T_n((H − center·I)/scale) |ψ₀⟩,   n = 0 … N

via the three-term recurrence

|φ₀⟩ = |ψ₀⟩,   |φ₁⟩ = H̃|ψ₀⟩,   |φ_k⟩ = 2H̃|φ_{k-1}⟩ − |φ_{k-2}⟩.

This is more memory-efficient than the full MPO version when only a single reference state is needed (e.g. site-resolved LDoS). Moments and spectral weights are then obtained as inner(ref_mps, Tn_mps_list[n+1]).

psi0 is normalised internally. scale/center follow the same convention as KPM_Tn: if scale=nothing the spectral bounds are estimated via DMRG. Returns (Tn_mps_list, scale, center) where Tn_mps_list[n+1] = |φ_n⟩.

_aux_setup(H, nambu_proj, proj_nambu, spin_proj, proj_s,
           layer_proj, proj_layer, sublat_proj, proj_sl) -> NamedTuple

Detect all auxiliary DOF indices from H and compute sector iteration ranges. Returns a NamedTuple with fields: nambu_s_det, nambu_side_det, spin_s_det, layer_s_det, layer_side_det, sublat_s_det, sublat_side_det, nambu_range, spin_range, layer_range, sl_range, any_aux_proj.

_dos_weight_matrix(Ncheb, ω_vals; kernel, lambda, eta, m_order)
    -> (W::Matrix, denom::Vector)

Stochastic-DOS reconstruction weights W[n, iω] and per-ω normalisation denom[iω] for a given KPM kernel. The DOS is recovered from the (sample- averaged) Chebyshev moments μ_n as Σ_n W[n,iω] μ_n / denom[iω].

• Convolution kernels (:jackson, :lorentz, :fejer, :dirichlet): W follows _kpm_weight_matrix and denom = π²·Ncheb·√(1−ω²), matching get_ldos_from_mun.
• :hodc: the contour weights νₙ(ω) from get_hodc_weights already carry the full normalisation (denom = 1), matching get_ldos_hodc_from_mun. eta=0 falls back to 1/(Ncheb+1).

Entries with |ω| ≥ 1 are zeroed in W (outside the rescaled spectral support).

_ensure_scale!(H::TBHamiltonian; dmrg_nsweeps, dmrg_maxdim, dmrg_linkdim)

If H.scale == 0 (sentinel meaning "not yet determined"), run _estimate_spectral_bounds and store the results in H.scale and H.center. No-op if H.scale > 0 (analytic estimate already set at construction or a previous KPM call already ran DMRG).

_estimate_spectral_bounds(H_mpo, sites; dmrg_nsweeps, dmrg_maxdim, dmrg_linkdim)
    -> (scale, center)

Run two short DMRG sweeps (minimising H and −H) to find the spectral edges Emin and Emax, then return: center = (Emax + Emin) / 2 scale = (Emax − Emin) / 2 × 1.1 (10 % buffer)

_get_exciton_ldos_cached(H::TBHamiltonian, X::Int, omega; kwargs...) -> Real

Internal legacy cache-backed scalar exciton LDOS helper. Public exciton LDOS now routes through get_exciton_ldos_spatial, which performs online MPS recursion without storing a Chebyshev cache on H.

_ldos_make_psi0(H, x, σ_n, σ_s, σ_l, σ_sl) -> MPS

Product-state MPS over all H.sites for KPM evaluation in LDOS :mps mode.

• Position sites encode x-1 in big-endian binary (first position site = MSB).
• Auxiliary sites are set to 1-based sector indices: σ_n (nambu), σ_s (spin), σ_l (layer), σ_sl (sublattice). Indices for absent aux dofs are ignored.

_run_kpm_mps!(Ham_n, psi0, Ncheb, W, valid, accum;
              weight=1.0, cutoff=1e-8, maxdim=100,
              verbose=false, label="") -> Int

Online MPS Chebyshev KPM recursion. Computes μ_n = ⟨psi0|T_n(Ham_n)|psi0⟩ for n = 1…Ncheb and accumulates W[n,iω] × μ_n × weight into accum[iω] for each valid energy index. Returns the maxlinkdim of the final state.

compute_dos_ldos_hodc(N, tn_lis, en_num; eta, m_order, maxdim)
    -> (dos_vec, ldos_mpo_list)

Compute DOS and energy-resolved LDOS MPOs from a list of Chebyshev moment MPOs tn_lis using the HODC kernel.

exciton_radius2(rho, d_list) -> Matrix{Float64}

Energy-resolved exciton radius

xi^2(R, E) = sum_d d^2 * P(d | R, E),   P(d | R, E) = rho(d,R,E) / sum_d' rho(d',R,E)

from a rho array as returned by get_exciton_ldos_separation (size(rho) == (Nomega, length(d_list), n_R)). Returns an (Nomega, n_R) matrix. Entries where any rho[:, :, R] along d is NaN (e.g. R+d outside the chain under boundary=:open), or where the normalisation sum_d rho(d,R,E) is zero, are returned as NaN.

get_dos_stochastic(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
                   N_sample, N_bound, seed, normalize, dos_weighting,
                   kernel, lambda, eta, m_order, maxdim, cutoff, verbose,
                   nambu_proj, proj_nambu, spin_proj, proj_s,
                   layer_proj, proj_layer, sublat_proj, proj_sl)
    -> Vector{Float64}

Stochastic full DOS via random trace estimation (MPS Chebyshev, 3 MPS per sample).

Normalization

• dos_weighting=:trace (default): returns the trace DOS. With normalize=false this is the total spectral weight Tr[δ(ω-H)]; with normalize=true it is divided by the traced Hilbert-space dimension.
• dos_weighting=:sample: returns the unweighted sample signal. For exciton stratified runs this is avg_full + avg_bound (when N_bound > 0), with no phase-space factor multiplying the continuum. This is intended for visualising the bound peak; normalize is ignored in this mode.

Auxiliary DOF projections

Unlike get_ldos_spatial, projections are not auto-enabled here. The default is full Hilbert-space sampling over all D states — always correct and cheapest. Projections must be requested explicitly:

• layer_proj=true, proj_layer=k — DOS on layer k only.
• sublat_proj=true, proj_sl=k — DOS on sublattice k only.
• Combining multiple *_proj=true flags is supported.

When any flag is set the function samples from position-basis states with fixed aux sectors (N_phys effective states), which is n_sectors× slower than the default. Only use projections when you actually need a sector-resolved DOS.

Exciton stratification (no aux projections)

For exciton Hamiltonians (length(H.sites) == 2*H.L), stratified sampling dedicates N_bound samples to the bound sector |x,x⟩ and N_sample to the full Hilbert space, combining with proper weights: DOS = N_phys × avg_bound + (D − N_phys) × avg_scatter. N_bound = 0 (default) = uniform sampling over all D states. Set dos_weighting=:sample to inspect the sampled spectral signal before these sector-size weights are applied.

Reconstruction kernel

kernel=:hodc uses the Higher-Order Delta Chebyshev contour reconstruction (eta, m_order control it; eta=0 → 1/(Ncheb+1)), whose weights already carry the full KPM normalisation. Other values are convolution kernels (:jackson default, :lorentz with lambda, :fejer, :dirichlet).

Examples

# Per-state DOS, size-independent
dos = get_dos_stochastic(H, 100, ωlist; N_sample=50, normalize=true)

# Sublattice-resolved per-state DOS (kagome)
dos_A = get_dos_stochastic(H_kg, 100, ωlist; sublat_proj=true, proj_sl=1,
                            N_sample=50, normalize=true)

# BdG: particle + hole combined per-state DOS
dos   = get_dos_stochastic(H_bdg, 100, ωlist; nambu_proj=true,
                            N_sample=50, normalize=true)

get_exciton_ldos_separation(H, Ncheb, omega_phys_vals; d_list, R_list=1:H.N,
                            boundary=:open, kernel=:jackson, lambda=4.0, eta=0.0,
                            m_order=4, maxdim=100, cutoff=1e-8, verbose=false,
                            printinfo=false) -> Array{Float64,3}

Relative-separation-resolved exciton LDOS

rho(d, R, E) = <R+d, R| delta(E - H) |R+d, R>

generalising get_exciton_ldos_spatial (the d = 0 slice) to electron-hole probes with electron at R+d and hole at R (both 1-indexed in 1:H.N).

For each (d, R) pair the probe |R+d, R> is a single product state (mpsexciton(R+d, R, H.sites)) — no superposition / carry construction is needed, unlike the momentum-space probes mpsexcitonQ/mpsexcitonQTrace. An online MPS Chebyshev recursion accumulates all energies in one pass, exactly as in get_exciton_ldos_spatial.

boundary=:open (default): pairs with R+d outside 1:H.N are left as NaN. boundary=:periodic: R+d is wrapped modulo H.N.

Returns an (Nomega, length(d_list), length(R_list)) array. Cost is length(d_list) * length(R_list) Chebyshev recursions, so keep these (and Ncheb) small for a first pass.

get_exciton_ldos_spatial(H, Ncheb, omega_phys_vals; X_list, X_groups,
                         num_x, num_avg, x_start, x_end, kernel,
                         lambda, eta, m_order, maxdim, cutoff,
                         verbose, printinfo) -> Matrix{Float64}

CPU spatial exciton LDOS. For each bound exciton position X (electron = hole = X, 1-indexed in 1:H.N) this runs an online MPS Chebyshev recursion from |X,X> and accumulates all requested energies in one pass. No Chebyshev cache is stored on H.

Rows are energies, columns are positions/groups. X_list selects positions directly. X_groups (or alias x_groups) averages several bound-pair probes into one output column. If no explicit positions are provided, num_x coarse groups are generated over x_start:x_end, with num_avg subpositions per group.

kernel=:hodc uses the HODC reconstruction (eta, m_order); otherwise the standard KPM kernels are available (:jackson, :lorentz, :fejer, :dirichlet).

get_ldos(H::TBHamiltonian, ω_phys; mode, psi0, kernel, lambda, eta, m_order,
         maxdim, cutoff, zl, wl)

Compute the local density of states at physical energy ω_phys using the Chebyshev expansion cached in H by a prior KPM_Tn or KPM_Tn_mps call.

Modes

• :mpo (legacy) — calls get_ldos_w_from_Tn and returns a full spectral-weight MPO at ω_phys. Requires KPM_Tn(H, N; mode=:mpo). Retains off-diagonal information; use when spatial correlations are needed.

• :diag — calls get_ldos_diag_from_Tn and returns an MPS encoding only the diagonal A(r, ω_phys). Requires KPM_Tn(H, N; mode=:mpo). Much cheaper than :mpo for typical LDOS use-cases; mirrors the get_bands momentum-space pattern. Use get_ldos_spectrum to compute all energies in a single pass. Use get_ldos_online to avoid storing the Tn cache entirely.

• :mps — computes moments μₙ = ⟨ψ₀|φₙ⟩ from the MPS Chebyshev cache and calls get_ldos_from_mun, returning a Real. Requires KPM_Tn(H, N; mode=:mps, psi0=...).

Physical energies are converted via E = (ω_phys − H.center) / H.scale.

get_ldos_diag_from_Tn(Tn_list, N, ω_vals; kernel, lambda, maxdim, cutoff)
    -> Vector{Union{Nothing, MPS}}

Compute site-resolved LDOS at every energy in ω_vals from a stored Chebyshev MPO list Tn_list, using the same online-accumulation + diagonal-extraction pattern as get_bands in momentum space.

At each Chebyshev step n, the diagonal of Tn_list[n] is extracted as an MPS via extract_diagonal_to_mps and its KPM-weighted contribution is accumulated into the LDOS for every energy point simultaneously — avoiding construction and storage of full weighted MPOs.

Returns a Vector of length length(ω_vals). Each entry is either:

• an MPS encoding the site-resolved LDOS A(r, ω) at that (rescaled) energy, or
• nothing for energies with |ω| ≥ 1 (outside the rescaled spectral support).

ω_vals must be rescaled energies in (−1, 1) — convert from physical units with E = (ω_phys − center) / scale. The legacy per-energy full-MPO path is still available via get_ldos_w_from_Tn.

get_ldos_from_mun(mun_list, N, E; kernel=:jackson, lambda=4.0) -> Real

Reconstruct the local spectral weight at rescaled energy E ∈ (−1, 1) from a list of Chebyshev moments μ_n = ⟨ψ₀|T_n(H̃)|ψ₀⟩ produced by KPM_Tn_mps.

Equivalent to computing ⟨ψ₀|δ(E − H̃)|ψ₀⟩ via the KPM expansion:

A(E) ≈ [g₀μ₀ + 2 Σ_{n≥1} gₙ Tₙ(E) μₙ] / (π √(1−E²))

where gₙ are the kernel damping weights (Jackson by default). Supported kernel values: :jackson, :lorentz (requires lambda), :fejer, :dirichlet. Returns 0 for |E| ≥ 1.

To convert a physical energy ω: E = (ω − center) / scale. To obtain the density of states per site, sum over all sites and divide by N.

get_ldos_hodc_from_mun(mun_list, N, E; eta=0.02, m_order=6) -> Real

HODC (Higher-Order Delta Chebyshev) variant of get_ldos_from_mun. Uses a contour-based kernel that gives sharper spectral features than the Jackson kernel, at the cost of m_order extra parameters.

The HODC weights νₖ (from compute_hodc_params / get_hodc_weights) already carry the full KPM normalisation, so no extra denominator is needed:

A_hodc(E) ≈ ν₁μ₁ + 2 Σ_{n≥2} νₙ μₙ

Returns 0 for |E| ≥ 1.

get_ldos_online(H::TBHamiltonian, Ncheb::Int, X::Int, ω_phys_vals;
                kernel, lambda, maxdim, cutoff, verbose,
                nambu_proj, proj_nambu, spin_proj, proj_s,
                layer_proj, proj_layer, sublat_proj, proj_sl)
    -> Vector{Float64}

Online real-space LDOS at unit-cell position X for all physical energies in ω_phys_vals. Never stores more than 3 MPS simultaneously (no Chebyshev cache).

Algorithm: MPS Chebyshev recursion |φ_k⟩ = T_k(H̃)|X⟩ with moment accumulation μ_k = ⟨X|φ_k⟩. Auxiliary DOF sectors are summed by running the recursion once per requested sector.

X ∈ {1, …, H.N} is the 1-indexed unit-cell position.

Auxiliary DOF projections (same interface as get_bands and get_ldos_spatial):

• spin_proj, nambu_proj, layer_proj, sublat_proj — enable projection of the corresponding auxiliary DOF auto-detected from H.
• proj_s, proj_nambu, proj_layer, proj_sl — sector selector: nothing sums all sectors of that DOF; an integer selects a single sector (1-based).
• Contributions from all requested sectors are accumulated into a single result vector.

Returns Vector{Float64} of length Nω with 0.0 outside the spectral support.

Examples

ωlist = range(-3.0, 3.0; length=300)
ldos  = get_ldos_online(H, 200, 2^(H.L-1), ωlist)          # no aux DOF

# Spin-summed LDOS at site 16
ldos_tot = get_ldos_online(H_spin, 200, 16, ωlist; spin_proj=true)

# Spin-↑ LDOS only
ldos_up  = get_ldos_online(H_spin, 200, 16, ωlist; spin_proj=true, proj_s=1)

get_ldos_spatial(H, Ncheb, ω_phys_vals;
                 num_x, num_avg, mode, x_start, x_end, x_groups,
                 kernel, lambda, maxdim, cutoff, verbose,
                 nambu_proj, proj_nambu, spin_proj, proj_s,
                 layer_proj, proj_layer, sublat_proj, proj_sl)
    -> Matrix{Float64}

Spatially-resolved LDOS, real-space analogue of get_bands.

Sampling procedures (reduce) — full detail in spatial_sampling_plan.

• :point (default) — read the LDOS at num_x[×num_y] sample cells; with box_half > 0 each pixel is the mean over a (2·box_half+1)² box. Cheap, but a grid coarser than a feature's width aliases it (thin in-gap edge / domain-wall channels can fall between pixels and be missed).
• :block — partition the system into num_x × num_y blocks (powers of two) and report the integral over each block, computed by tracing out the within-block position bits (a partial contraction, cost independent of block size). Gap-free: every cell belongs to one block, so a thin feature on a gapped background cannot be missed. Use it for large-scale maps of edge networks. Output columns are row-major over coarse pixels (col = ixp + iyp·num_x + 1); block centres come from the plan. :mpo mode only.

Return shape (sublattice geometry-awareness)

For a multi-atom unit cell (H.sublattice_s set) the layout depends on the sampling scale, decided by spatial_sampling_plan from the local stride (see sublattice below):

• resolved (atomic scale — every unit cell probed, or proj_sl=k): (Nω × ng×n_sub), columns interleaved in atom order matching the *_positions functions — [A₀, B₀, A₁, B₁, …] (2-sublattice), [A₀, B₀, C₀, …] (3-sublattice). proj_sl=k fills only sublattice k.
• averaged (large scale — the grid skips unit cells): (Nω × ng), the sublattice is traced out into one value per unit cell (mean over the n_sub atoms).

With no sublattice DOF the shape is always (Nω × ng), ng = num_x.

sublattice (resolve vs average)

• :auto (default) — resolve when sampling at full unit-cell resolution (stride == 1, e.g. zooming a small window and probing every cell); average whenever the grid is coarser or box_half > 0.
• :resolve — always emit per-atom columns. :average — always trace the sublattice to one value per cell. proj_sl=k always resolves that one atom.

Sampling parameters

• num_x/num_y : sample counts (per axis for grid=true; num_x is the total for the default 1D linear sweep). Default H.N = full resolution.
• num_avg : sub-samples per coarse block for local averaging (1D, default 1).
• x_start, x_end : 1-indexed linear position range (1D layout).
• grid : true lays centers on a 2D num_x × num_y unit-cell grid.
• xwin, ywin : 0-indexed unit-cell (lo, hi) windows for grid=true (e.g. zoom into a patch of a large system).
• x_groups : explicit Vector{Vector{Int}} override (treated as atomic).
• box_half : 2D neighbourhood half-width (averages, forces sublattice averaging).
• reduce : :point (sample/box) or :block (block-integrate; see above).

Modes

• :mpo (default) — single Chebyshev pass; evaluates all positions simultaneously. Cost ∝ Ncheb × (MPO×MPO), independent of num_x or n_sub.
• :mps — independent MPS recursion per (position, sector) combination.

Other auxiliary DOF projections (same interface as get_bands): nambu_proj/proj_nambu, spin_proj/proj_s, layer_proj/proj_layer.

Examples

# Standard 1D chain — shape (Nω × 8)
ldos = get_ldos_spatial(H, 200, ωlist; num_x=8)

# Honeycomb, large-scale map — sublattice averaged, shape (Nω × 64)
ldos_uc  = get_ldos_spatial(H_hc, 200, ωlist; num_x=64)

# Honeycomb, atomic zoom into a 50×50 patch — sublattice resolved (Nω × 50*50*2)
ldos_zoom = get_ldos_spatial(H_hc, 200, ωlist; grid=true,
                             xwin=(1000, 1049), ywin=(1000, 1049))

# Large 2^14×2^14 system: 128×128 block-integrated map — catches thin in-gap
# edge channels that point sampling would alias away. Shape (Nω × 128*128).
ldos_blk = get_ldos_spatial(H_big, 200, ωlist; reduce=:block, num_x=128, num_y=128)

# Kagome: sublattice A only — only A columns filled, B/C columns zero
ldos_A   = get_ldos_spatial(H_kg, 200, ωlist; proj_sl=1, num_x=H_kg.N)

get_ldos_spectrum(H::TBHamiltonian, ω_phys_vals; kernel, lambda, maxdim, cutoff)
    -> Vector{Union{Nothing, MPS}}

Compute the site-resolved LDOS at all physical energies in ω_phys_vals in a single pass over the cached Chebyshev MPO list — the real-space equivalent of get_bands.

At each Chebyshev step the diagonal of T_n is extracted once and its weighted contribution is accumulated into every energy slot simultaneously, so the cost scales as O(Ncheb) MPO operations regardless of how many energy points are requested.

Returns Vector{Union{Nothing, MPS}} of length length(ω_phys_vals). Each MPS encodes the site-resolved LDOS A(r, ω) at the corresponding physical energy; entries are nothing for energies outside the spectral support.

Requires KPM_Tn(H, Ncheb; mode=:mpo) to have been called first.

Example

KPM_Tn(H, 200; mode=:mpo, maxdim=100)

ωlist    = range(-4.0, 4.0; length=200)
ldos_vec = get_ldos_spectrum(H, ωlist)

# Evaluate LDOS at site x=16 for each energy:
ldos_at_16 = [l === nothing ? 0.0 : _eval_diag_mps(l, 15) for l in ldos_vec]

ldos_exc_KPM_Tn(H, N, X; cutoff, maxdim) -> Vector

Low-level: Chebyshev moment list ⟨X|Tn(H)|X⟩ for the exciton state |X,X⟩. X ∈ {1, …, 2^LPhys} (1-indexed). H must already be rescaled so its spectrum lies in (−1, 1). Prefer `KPMTn(H::TBHamiltonian, Ncheb, X)` for the high-level interface which handles rescaling and caching automatically.

_vec_mps_from_mpo(G_mpo, sites2; cutoff, maxdim) -> MPS

Convert an L-site MPO into the 2L-site interleaved vectorized MPS used as the initial guess for get_green_krylov.

Each MPO tensor at site k is split by SVD into two MPS tensors at positions (2k-1, 2k): the bra (primed) physical index maps to the odd row site and the ket (unprimed) index maps to the even column site, matching the encoding produced by custom_mpo and _identity_vec_mps.

get_green_krylov(H_mpo, sites, ω_phys; η, nsweeps, maxdim, cutoff,
                 x0_mpo, ishermitian, tol, maxiter, krylovdim, verbose) -> MPO

Low-level: compute the retarded Green's function G(ω) = (ω + iη − H)⁻¹ for a raw MPO H_mpo defined on sites, via the vectorized linear system

[(ω + iη − H) ⊗ I] |G⟩⟩ = |I⟩⟩

See the TBHamiltonian overload for the full keyword-argument reference.

get_green_krylov(H::TBHamiltonian, ω_phys; η, nsweeps, maxdim, cutoff,
                 x0_mpo, ishermitian, tol, maxiter, krylovdim, verbose) -> MPO

Compute the retarded Green's function

G(ω) = (ω + iη − H)⁻¹

as an MPO by solving the vectorized linear system

[(ω + iη − H) ⊗ I] |G⟩⟩ = |I⟩⟩

using ITensorMPS.linsolve (DMRG-like Krylov solver). The Hamiltonian is used unscaled — no KPM Chebyshev expansion required.

Keyword arguments

• η : Lorentzian broadening. Default 1e-2.
• nsweeps : Number of DMRG sweeps for the linear solver. Default 12.
• maxdim : Maximum bond dimension of the solution MPS. Default 100.
• cutoff : SVD truncation cutoff. Default 1e-8.
• x0_mpo : Optional MPO initial guess for G(ω). When provided it is vectorized via _vec_mps_from_mpo and passed as x0 to linsolve, replacing the default identity-matrix guess. Typical use: pass a low-accuracy KPM Green's function to warm-start the Krylov iteration. Default nothing.
• ishermitian : Set true only when the shifted operator is Hermitian (requires purely imaginary η = 0, not physical for GF). Default false.
• tol : Krylov solver convergence tolerance. Default 1e-10.
• maxiter : Maximum Krylov iterations per site. Default 600.
• krylovdim : Krylov subspace dimension. Default 30.
• verbose : Print progress messages. Default false.

Usage

G = get_green_krylov(H, ω; η=0.05, nsweeps=20, maxdim=200)
dos  = -imag(tr(G)) / π                              # total DoS
ldos = real(inner(psi_i, apply(G, psi_i)))           # LDoS at site i
gij  = inner(psi_i, apply(G, psi_j))                 # off-diagonal element

# Warm-start from a cheap KPM estimate
TensorBinding.KPM_Tn(H, 15; maxdim=50)
G_kpm = TensorBinding.get_Green_retarded_from_Tn(H._tn_cache, 15, ω; η=η, maxdim=50)
G_ws  = get_green_krylov(H, ω; x0_mpo=G_kpm, nsweeps=6, maxdim=200)

build_K(H_mpo, sites, ω, η; maxdim_K, cutoff_K) -> MPO

Build the resolvent-squared MPO

K(ω,η) = (H − ω I)² + η² I

The ground state energy of K satisfies E_K ≥ η², with E_K → η² when ω coincides with an eigenvalue of H. The ground state concentrates on the eigenstate of H nearest to ω.

maxdim_K and cutoff_K control truncation of the intermediate MPO product; larger maxdim_K gives a more accurate K at the cost of DMRG wall time.

dmrg_gs(H_mpo, sites; nsweeps, maxdim, cutoff, noise, linkdim_init) -> (E, ψ)

Find the ground state and energy of H_mpo using DMRG.

Keyword arguments

• linkdim_init : bond dimension of the random initial MPS
• nsweeps : total number of DMRG sweeps
• maxdim : max bond dimension per sweep (scalar or vector)
• cutoff : SVD truncation cutoff
• noise : perturbative noise per sweep (scalar or vector; aids convergence)

Returns (E, ψ).

dmrg_spectral(H_mpo, sites, ω, η; ...) -> (E_K, ψ)

Find the ground state of K(ω,η) = (H−ωI)² + η²I using DMRG.

The ground state |ψ(ω)⟩ concentrates on the eigenstate of H closest to ω. Use local_weight(ψ, i, L, sites) to read off the LDoS proxy |⟨i|ψ(ω)⟩|².

Keyword arguments

• ψ0 : warm-start MPS (random if nothing)
• maxdim_K, cutoff_K : truncation for building K (see build_K)
• linkdim_init : bond dimension of the random initial MPS (if ψ0=nothing)
• nsweeps, maxdim, cutoff, noise : DMRG sweep parameters

Returns (E_K, ψ) where E_K ≈ η² at spectral peaks.

local_weight(ψ, i, L, sites) -> Float64

Return |⟨i|ψ⟩|² where |i⟩ is the position-basis state for 0-based integer i (big-endian quantics encoding over L qubit sites).

When ψ = ψ(ω) is the DMRG ground state of K(ω,η), this gives the LDoS proxy at site i and energy ω:

ρ(ω, i) ≈ |⟨i|ψ(ω)⟩|²

which peaks at eigenvalues of H that have support on site i.

For a BdG or spin-extended system pass the full ext_sites and set L = length(ext_sites).

apply_mpo_to_mps(U_mpo, psi; cutoff, maxdim, normalize) -> MPS

Apply a propagator MPO U_mpo to the MPS psi with optional truncation and normalisation. Used to advance a state by one time step when U_mpo was prebuilt by build_tdvp_propagator_mpo.

basis_amplitude(psi, n, L, sites) -> ComplexF64

Return the amplitude ⟨n|ψ⟩ where |n⟩ is the n-th computational basis state (0-indexed big-endian quantics encoding over L qubit sites).

bond_current_x(ρ, j, tx, L, sites) -> ComplexF64

Compute the x-direction bond current

Jⱼˣ = i tₓ (ρⱼ,ⱼ₊₁ − ρⱼ₊₁,ⱼ)

for a single bond at 0-indexed site j in density-matrix MPO ρ. tx is the hopping amplitude.

bond_current_x_trajectory(states, j, tx, L, sites; dt) -> Vector{ComplexF64}

Compute the x-direction bond current Jⱼˣ(t) along a trajectory of density-matrix MPOs.

tx may be a scalar or a callable t -> hopping amplitude for a time-dependent drive. States are assumed spaced by dt: t_n = (n-1)*dt.

build_tdvp_propagator_mpo(H, dt, L, sites; maxdim, cutoff, reverse_step,
                          outputlevel, nsite, cross_tol, initial_positions,
                          use_diagonal_pivots, interpolation_type) -> MPO

Build an MPO approximation of the short-time propagator U(dt) = e^{-iH dt} by sampling matrix elements ⟨i|U(dt)|j⟩ via TDVP and compressing with TCI.

H must already be multiplied by -im for Schrödinger evolution. The diagonal is dominant for small dt; by default TCI is seeded with all diagonal pivots (use_diagonal_pivots=true) so the near-identity structure is captured first.

Keyword arguments

• maxdim, cutoff : TDVP truncation parameters.
• cross_tol : TCI interpolation tolerance.
• use_diagonal_pivots : Seed TCI with the N diagonal positions. Default true.
• interpolation_type : Element type for TCI sampling. Default ComplexF64.

A TBHamiltonian overload applies -im internally: build_tdvp_propagator_mpo(H::TBHamiltonian, dt; ...).

central_x_bond(L; Nx) -> Int

Return the 0-indexed site j of the central x-direction bond (j → j+1).

• 1D (Nx=nothing): central bond at j = 2^L ÷ 2 − 1.
• 2D row-major (Nx = 2^Lx): bond at the centre column of the centre row, j = (Ny÷2)*Nx + Nx÷2 − 1 where Ny = 2^L ÷ Nx.

check_tdvp_vs_U_mpo(H, U_mpo, dt, L, sites; test_states, ...) -> (max_overlap_error, max_phase_error)

Validate that U_mpo agrees with direct TDVP on a set of computational basis states. Prints per-state overlap errors and phase-aligned distances, then returns the maxima.

A TBHamiltonian overload is available.

compare_propagator_and_tdvp_heatmaps(U_mpo, H, psi0, L, sites, nsteps; ...)

Full comparison of MPO-propagator and TDVP trajectories: evolves psi0 with both methods for nsteps steps, renders heatmaps of |⟨x|ψ(t)⟩| and |⟨x|ψ(t)⟩|², and returns all trajectory data and agreement metrics as a named tuple.

A TBHamiltonian overload is available.

compute_basis_overlaps(states, L, sites)
    -> NamedTuple(overlaps, abs_overlaps, probabilities, norms)

For each MPS in states, compute overlaps with all 2^L computational basis states.

Returns a named tuple with fields:

• overlaps : (nsteps+1) × 2^L matrix of ComplexF64 amplitudes ⟨j|ψ(t)⟩
• abs_overlaps : element-wise absolute values
• probabilities : |⟨j|ψ(t)⟩|²
• norms : ⟨ψ(t)|ψ(t)⟩ at each step

dm_expect(O, ρ) -> Float64

Compute Tr(O ρ) for a Hermitian operator MPO O and density-matrix MPO ρ, using inner(O, ρ) = Tr(O† ρ) = Tr(O ρ).

evolve_rk4_dm_nh(Hoft, ρ0, nsteps, dt; maxdim, cutoff,
                 truncate_intermediates, verbose) -> Vector{MPO}

RK4 evolution of a density-matrix MPO under the non-Hermitian von Neumann equation

dρ/dt = -i(H(t) ρ − ρ H(t)†)

Hoft is a callable t::Float64 -> MPO. For H = H₀ - iΓ the anti-Hermitian part causes Tr(ρ) to decay whenever Γ > 0, modelling lossy open systems.

Returns [ρ(0), ρ(dt), ..., ρ(nsteps*dt)] as a Vector{MPO}. See evolve_rk4_dm_timedep for keyword-argument descriptions.

evolve_rk4_dm_timedep(Hoft, ρ0, nsteps, dt; maxdim, cutoff,
                      truncate_intermediates, verbose) -> Vector{MPO}

Evolve a density-matrix MPO ρ0 under dρ/dt = -i[H(t), ρ] for nsteps steps of size dt using RK4.

Hoft is a callable t::Float64 -> MPO returning the physical Hamiltonian at time t. Returns the full trajectory [ρ(0), ρ(dt), ..., ρ(nsteps*dt)] as a Vector{MPO}.

Keyword arguments

• maxdim, cutoff : Truncation for intermediate MPO sums and products.
• truncate_intermediates : Truncate after each RK4 sub-step to control bond growth.
• verbose : Print step/bond-dim progress.

evolve_with_propagator(U_mpo, psi0, nsteps; normalize_each_step,
                       cutoff, maxdim) -> Vector{MPS}

Apply the fixed MPO propagator U_mpo repeatedly for nsteps steps, returning the full trajectory [psi(0), psi(1*dt), ..., psi(nsteps*dt)].

Useful when the same short-time propagator is reused at every step (time-independent H). For efficiency the MPO is built once via build_tdvp_propagator_mpo; this function then applies it nsteps times.

evolve_with_tdvp(H, psi0, nsteps, dt; normalize_each_step, maxdim,
                 cutoff, reverse_step, outputlevel, nsite) -> Vector{MPS}

Run a TDVP loop for nsteps steps of size dt under a fixed Hamiltonian H, returning [psi(0), psi(dt), ..., psi(nsteps*dt)].

H must carry the -im prefactor for Schrödinger evolution.

A TBHamiltonian overload applies -im internally: evolve_with_tdvp(H::TBHamiltonian, psi0, nsteps, dt; ...).

evolve_with_tdvp_timedep(Hoft, psi0, nsteps, dt; normalize_each_step,
                         maxdim, cutoff, reverse_step, outputlevel,
                         nsite, krylovdim, tol) -> Vector{MPS}

TDVP loop for a time-dependent Hamiltonian H(t).

Hoft is a callable t::Float64 -> MPO. On each interval [t, t+dt] the Hamiltonian is frozen at the midpoint t + dt/2 (midpoint rule). Hoft must return the physical Hamiltonian; the -im prefactor is applied internally.

Returns [psi(0), psi(dt), ..., psi(nsteps*dt)].

observables_trajectory(ops, states) -> Dict

Measure a named collection of Hermitian operator MPOs along a density-matrix trajectory.

ops is a NamedTuple or Dict mapping labels to MPOs, e.g. (J=J_mpo, N=N_mpo). Returns a Dict mapping each label to a Vector{Float64} of expectation values.

phase_aligned_distance(psi_a, psi_b) -> Float64

Phase-insensitive distance between two (unnormalised) MPS states:

d = min_{φ} ‖â − e^{iφ} b̂‖ = √(2 − 2|⟨â|b̂⟩|)

where â = psi_a/‖psi_a‖. Returns Inf when either state has zero norm. Useful for comparing TDVP and MPO propagator trajectories independent of any global phase accumulated during time evolution.

purity(ρ) -> Float64

Return the purity Tr(ρ²) = inner(ρ, ρ). Equals 1 for a pure state and decreases as the system becomes mixed.

purity_trajectory(states) -> Vector{Float64}

Return purity(ρ) for each MPO in states.

rk4_step_dm_nh(Hoft, ρ, t, dt; maxdim, cutoff, truncate_intermediates) -> MPO

Single RK4 step for the non-Hermitian von Neumann equation

dρ/dt = -i(H(t) ρ − ρ H(t)†)

Hoft is a callable t -> MPO; pass (_ -> H.mpo) for a static NH Hamiltonian.

rk4_step_dm_timedep(Hoft, ρ, t, dt; maxdim, cutoff,
                    truncate_intermediates) -> MPO

Single RK4 step for dρ/dt = -i[H(t), ρ] with a time-dependent Hamiltonian MPO. H is evaluated at t, t+dt/2, and t+dt per the classical RK4 tableau.

tdvp_evolve(H, psi, dt; maxdim, cutoff, normalize, reverse_step,
            outputlevel, nsite) -> MPS

Apply one TDVP step to psi under Hamiltonian H for time dt.

H must already carry the -im prefactor for Schrödinger evolution. When normalize=true the output is renormalised after the step.

A TBHamiltonian overload applies -im internally: tdvp_evolve(H::TBHamiltonian, psi, dt; ...).

timedep_observable_trajectory(Oft, states, dt) -> Vector{Float64}

Measure a time-dependent Hermitian operator O(t) along a density-matrix trajectory.

Oft is a callable t -> MPO. State n is assigned time t_n = (n-1)*dt. Covers e.g. ⟨H(t)⟩ = Tr(H(t) ρ(t)) under a driven Hamiltonian.