Physics

cdw_hartree_builder(U; background=0.5, maxdim=100, cutoff=1e-8)

Return a simple spinless CDW Hartree builder:

V_H[n] = U * (rho[n] - background)

This is intentionally local and minimal; pass a custom hartree_builder to meanfield for long-range interactions or other channels.

dense_hartree_builder(V, L, sites; background=nothing, kwargs...)

Return a Hartree builder backed by a dense interaction kernel. V can be an existing interaction MPO or a function V(i, j) on 0-indexed coordinates.

density_profile_from_dm(density_mpo, sites; mode=:direct) -> MPS

Extract a density profile MPS from the diagonal of a density-matrix MPO. mode=:complement returns 1 - diag(D).

fock_mpo_from_density(density_mpo, interaction_op, sites;
                          sign=-1, maxdim=100, cutoff=1e-8) -> MPO

Build a Fock/exchange MPO from a density matrix and a two-body interaction kernel by element-wise multiplication:

F[i, j] = sign * V[i, j] * rho[i, j]

For the usual exchange term use the default sign=-1.

get_scf(H0, U, channel; kwargs...) -> NamedTuple

Convenience front-end for the self-consistent mean-field solvers.

Examples:

res = get_scf(H0, U, :CDW; tol=1e-5, maxiters=30, mixing=0.4)
Hmf = res.ham

Supported channels are :CDW, :magnetic, :swave, and :pwave. The wrapper keeps the lower-level SCF routines available while providing a single get_* style entry point for notebooks.

get_scf(H0::TBHamiltonian, channel::Symbol; interaction=:dense, kwargs...) -> NamedTuple

Like get_scf(H0, U, channel) but reads the interaction from H0.interaction_mpo (stored via add_interaction!). Supported channels: :cdw, :magnetic.

For :cdw the default interaction=:dense passes the stored MPO directly to dense_hartree_builder. For :magnetic the stored MPO is handled by magnetic_hubbard which dispatches on U::Union{Number,MPO}.

Calling with :swave or :pwave raises an informative error — those channels require an explicit coupling constant.

hartree_mpo_from_density(density_mps, interaction_op, sites;
                             background=nothing, maxdim=100, cutoff=1e-8)

Apply an interaction kernel MPO to a density profile and convert the resulting coefficient MPS into a diagonal Hartree MPO.

initial_guess_Neel_dn(Lx, Ly, sites) -> (MPO, MPS)

Checkerboard spin-down density seed for 2D Hubbard: occupation 1 where (ix+iy) is odd.

initial_guess_Neel_up(Lx, Ly, sites) -> (MPO, MPS)

Checkerboard spin-up density seed for 2D Hubbard: occupation 1 where (ix+iy) is even.

initial_guess_trivial_down_1D(L, sites) -> (MPO, MPS)

Diagonal density MPO with occupation (x+1) % 2 on site x (spin-down Néel seed for 1D).

initial_guess_trivial_up_1D(L, sites) -> (MPO, MPS)

Diagonal density MPO with occupation x % 2 on site x (spin-up Néel seed for 1D).

pair_distance_hartree_builder(L, sites, distance, weight=1; kwargs...)
pair_distance_hartree_builder(L, sites, terms; kwargs...)

Return a Hartree builder for sparse pair interactions generated from non-cyclic shift MPOs.

pair_distance_interaction_mpo(L, sites, distance, weight=1; kwargs...)
pair_distance_interaction_mpo(L, sites, terms; kwargs...)

Build an interaction MPO for arbitrary fixed-distance site pairs. For a bond weight w(i) on the pair (i, i+d), the operator contains row=i+d,col=i with weight w(i) and row=i,col=i+d with weight w(i). The shifts are non-cyclic, so boundary-wrapping bonds are absent.

pwave_equalspin_anomalous_profiles(density_mpo, Hbdg; distance=1)
    -> (F_up, F_dn)

Extract nearest-neighbor equal-spin triplet anomalous bond profiles from a spinful BdG density matrix. The profile value at i corresponds to the bond (i, i + distance) in the same non-cyclic convention as build_shift_mpo.

scf_magnetic_hubbard(H0, U; kwargs...) -> NamedTuple

Two-channel collinear magnetic mean-field loop for the on-site Hubbard model. If H0 is spinful, the spin-up and spin-down one-body blocks are obtained by projecting out the spin core, matching the magnetic RPA convention. If H0 is spinless, the previous two-copy behavior is retained.

H_up = H0_up + U * diag(n_down - background)
H_dn = H0_dn + U * diag(n_up   - background)

meanfield(H0, hartree_builder; kwargs...) -> NamedTuple

Generic self-consistent mean-field loop.

Workflow per iteration:

1. Build H_MF = H0 + H_Hartree.
2. Compute the density matrix by density_method = :kpm | :mcweeny | :sp2.
3. Extract a density profile MPS.
4. Compute RMS error, mix density, rebuild Hartree term.

Returns a named tuple with density, Hartree term, final Hamiltonian, and history.

scf_pwave_equalspin(H0, V; kwargs...) -> NamedTuple

Minimal spinful equal-spin p-wave triplet SCF. The pairing fields live on nearest-neighbor bonds by default:

Delta_up(i) = pairing_sign * V * F_up(i, i + distance)
Delta_dn(i) = pairing_sign * V * F_dn(i, i + distance)

eta_down=-1 gives a helical initial seed, while eta_down=1 starts the two equal-spin channels with the same sign.

scf_swave_hubbard(H0, U; kwargs...) -> NamedTuple

Attractive on-site Hubbard s-wave SCF. This is a convenience wrapper around swave_superconducting with the Hartree channel enabled. Positive U means attraction by default:

Delta_i = -U F_i
V_up    = -U (rho_dn - background)
V_dn    = -U (rho_up - background)

scf_swave_superconducting(H0, g; kwargs...) -> NamedTuple

Spin-singlet s-wave pairing SCF. By default this keeps the normal Hamiltonian fixed and updates only the pairing channel. Set include_hartree=true to also self-consistently update spin-resolved onsite Hartree potentials,

Delta_new(i) = pairing_sign * g * F_singlet(i)
V_up(i)      = hartree_sign * hartree_coupling * (rho_dn(i) - background)
V_dn(i)      = hartree_sign * hartree_coupling * (rho_up(i) - background)

where F_singlet is extracted from the Nambu off-diagonal block of the BdG density matrix. If H0 is spinless, a spin degree of freedom is added to a working copy before building BdG Hamiltonians.

staggered_magnetic_initial(H; amplitude=0.05, background=0.5)
    -> (rho_up, rho_dn)

Build a simple antiferromagnetic initial guess for a Hubbard mean-field loop. If H is spinful, the spin core is first projected out so the returned profiles live on the same spinless position/sublattice sites as each spin block.

swave_anomalous_profile(density_mpo, Hbdg; maxdim=100, cutoff=1e-8) -> MPS

Extract the local singlet anomalous profile from a spinful s-wave BdG density matrix. The returned MPS is the diagonal of the Nambu particle-hole block projected onto the spin-singlet tensor (i sigma_y).

swave_normal_profiles(density_mpo, Hbdg; mode=:particle,
                           maxdim=100, cutoff=1e-8)
    -> (rho_up, rho_dn)

Extract spin-resolved local normal densities from a spinful BdG density matrix. mode=:particle uses the particle-particle block directly; mode=:hole_complement uses 1 - diag(hole-hole block).

_project_spin_sector(H, sector) -> TBHamiltonian

Project a spinful TBHamiltonian onto spin sector sector (1 = ↑, 2 = ↓) by contracting the spin site tensor with the projector |sector⟩⟨sector|.

The spin index is identified by its "Spin" tag, so the function is robust to whether spin is prepended or postpended. The contracted tensor is absorbed into its neighbour, leaving a valid L-qubit MPO.

Returns a new TBHamiltonian with spin_s = nothing and fresh (empty) caches; scale and center are reset to 0.0 so _ensure_scale! will re-estimate them on the first KPM call.

apply_interleave_swaps(W, sites; cutoff, maxdim, verbose) -> MPO

Re-order the sites of W from [1…N, N+1…2N] to the interleaved order [1, N+1, 2, N+2, …] by composing a sequence of adjacent SWAP gates. Truncates after each swap to control bond dimension growth.

build_bubble_mpo(ω; Tn_list1, Tn_list2, Tn_listeff,
                    sites, sites2, finalsites, finalfinalsites,
                    chi, ϵF, a, maxdim) -> MPO

Compute the L-site polarization bubble Π₀(ω) from pre-computed Chebyshev lists.

chebyshev2d_gf_coeffs(ω, scale1, center1, scale2, center2, η, N) -> Matrix{ComplexF64}

Compute 2D Chebyshev expansion coefficients c_{mn} for the scalar Green's function

f(x, y) = 1 / (ω + iη − (scale2·y + center2 − scale1·x − center1))

on the domain [-1,1]×[-1,1], using an N×N Chebyshev-Gauss grid and 2D DCT-II.

C[m+1, n+1] = c{mn} in the expansion f(x,y) ≈ Σ{m,n} c{mn} Tm(x) T_n(y).

N should equal length(Tn_list) from KPM_Tn (i.e. Ncheb + 1).

collapse_mpo_pairs(mpo2L, out_sites) -> MPO

Merge each consecutive pair of sites (2n-1, 2n) in a 2L-site MPO into a single site of an L-site MPO by contracting and tying the shared bra and ket indices to out_sites[n].

convert_mpo(old_mps, new_sites) -> MPO

Convert a 2N-site MPS (typically from QTCI) into an N-site MPO by merging pairs and remapping physical indices to new_sites.

eval_haydock_cf(a, b, z) -> ComplexF64

Evaluate the Haydock continued fraction ⟨seed|(z−H)⁻¹|seed⟩ via backward recursion. b[1] must be norm0 = ||seed||F (as returned by `haydockcf`).

G(z) = b[1]² / (z − a[1] − b[2]²/(z − a[2] − b[3]²/…))

Calling with truncated arrays a[1:N], b[1:N] gives the N-th CF convergent, whose sequence over N is suitable for Wynn ε-acceleration.

extract_diagonal_to_mps(M) -> MPS

Extract the diagonal of an MPO M as an MPS by projecting each site tensor onto the subspace where bra and ket indices are equal.

get_Tnlists(H, H2, sites, sites2, N; a, maxdim)
    -> (Tn_list1, Tn_list2, Tn_listeff)

Build the three Chebyshev moment lists needed for the bubble:

• Tn_list1 : moments for H₁ (system 1)
• Tn_list2 : moments for H₂ (system 2)
• Tn_listeff : moments for H_eff = I⊗H₂ − H₁⊗I on the combined 2L-site space

get_bubble_diag_cheb2d(H1, H2, ωlist; Ncheb, maxdim, cutoff,
                        ϵF, P_method, purify_*, η, coeff_tol,
                        qft_tol, qft_maxdim, verbose) -> Vector{MPS}

Diagonal-only variant of get_bubble_mpo_cheb2d.

Returns the k-space diagonal of the non-interacting polarization bubble, diag_Π₀(k, ω) = ⟨k| Π₀(ω) |k⟩, as a Vector{MPS} (one MPS per ω in ωlist) ready for direct plotting.

Compared to get_bubble_mpo_cheb2d, this function:

• QFTs each D_mn once (inside the (m,n) loop) and extracts its k-space diagonal as an MPS.
• Accumulates c_mn(ω) · diag(D_mn) as MPS sums instead of MPO sums.
• Skips the per-ω conjugate_by_qft + extract_diagonal_to_mps steps.

This follows the same paradigm as get_bands: the expensive MPO-level work (Hadamard products, QFT conjugation) is done once per (m,n) pair and shared across all frequencies; per-ω cost is a cheap scalar-weighted MPS addition.

Use get_bubble_mpo_cheb2d when you need the full off-diagonal MPO (e.g. for RPA resummation). Use this function when only χ₀(k,ω) is needed.

Keyword arguments — identical to get_bubble_mpo_cheb2d, plus:

• qft_tol : truncation tolerance inside conjugate_by_qft. Default 1e-9.
• qft_maxdim : max bond dimension inside conjugate_by_qft. Default 100.

get_bubble_diag_cheb2d_svd(H1, H2, ωlist; ..., svd_tol, svd_maxrank) -> Vector{MPS}

Per-ω SVD-accelerated variant of get_bubble_diag_cheb2d.

For each frequency ω the coefficient matrix C[m,n](ω) is rank-truncated via its own SVD:

`C[m,n](ω) = Σ_s  S_s(ω) · U[m,s](ω) · conj(V[n,s](ω))`   (exact up to truncation)

The per-ω rank r(ω) (typically 2–5 for smooth Lorentzian kernels) is usually much smaller than the Tucker/joint-SVD rank, which must span all frequencies simultaneously. For each (ω, s) one Hadamard product and one QFT are performed, giving

`diag_Π[ω] = Σ_s S_s(ω) · diag(QFT( A_s ⊙ B_s − C_s ⊙ E_s ))`

where A_s = Σ_m U[m,s]·TP1[m], B_s = Σ_n conj(V[n,s])·Tn2[n], etc.

Total Hadamard+QFT operations: Σω r(ω) — compared to N² for the plain variant or rm·rn for Tucker. When r(ω) ≪ rTucker the per-ω SVD is both faster and more accurate, because it uses the optimal basis for each individual C(ω).

Additional keyword arguments (beyond get_bubble_diag_cheb2d):

• svd_tol : relative singular-value cutoff per ω (fraction of σ_max(ω)). Default 1e-6.
• svd_maxrank: hard cap on the per-ω rank. Default 20.
• kernel : Chebyshev damping kernel applied before SVD. :jackson (default) suppresses Gibbs oscillations and dramatically reduces per-ω rank (typically 2–5 instead of ~26). Use :none for exact Chebyshev coefficients.

get_bubble_diag_cheb2d_tucker(H1, H2, ωlist; ..., tucker_tol, tucker_maxrank, kernel) -> Vector{MPS}

Tucker (HOSVD) variant of get_bubble_diag_cheb2d.

Finds a global low-rank basis in the (m, n) indices shared across all frequencies by stacking the coefficient matrices and performing two mode-SVDs:

• Mode-1 SVD of [C(ω₁) | C(ω₂) | … | C(ωₙ)] (shape N × N·nω) → basis Um (N × rm)
• Mode-2 SVD of [C(ω₁)ᵀ | … | C(ωₙ)ᵀ] (shape N × N·nω) → basis Vn (N × rn)

Then for each (s1, s2) ∈ {1…rm} × {1…rn}, one Hadamard product and one QFT are computed (both ω-independent), giving r_m × r_n such operations in total. Per-ω cost reduces to cheap scalar-weighted MPS additions over the core tensor Ã[s1, s2, ω] = Umᵀ C(ω) Vn.

Scaling comparison (assuming rank saturation with Jackson kernel):

• Plain diagonal: N² Hadamard+QFT
• Per-ω SVD: nω × rperω Hadamard+QFT
• Tucker: rm × rn Hadamard+QFT ← frequency-independent

Tucker wins when rm × rn < nω × rperω, which holds for large nω or when the global (m,n) structure is very low-dimensional (as it is after Jackson damping).

Additional keyword arguments (beyond get_bubble_diag_cheb2d):

• tucker_tol : relative singular-value cutoff for both mode SVDs. Default 1e-3.
• tucker_maxrank: hard cap on rm and rn. Default 20.
• kernel : :jackson (default) or :none. Jackson damping is essential here — without it the Tucker rank is high and results are inaccurate.
• coeff_tol : skip core-tensor entries |Ã[s1,s2,ω]| below this threshold. Default 1e-12.
• hooi_iters : number of Higher-Order Orthogonal Iteration refinement steps after the initial HOSVD. Each iteration re-optimises Um and Vn jointly, giving the globally optimal Tucker bases for the given rank (vs. the independent mode-SVDs of plain HOSVD). Default 3; set to 0 for HOSVD only.

get_bubble_mpo(H1::TBHamiltonian, H2::TBHamiltonian, ω; ...) -> MPO

Compute the non-interacting polarization bubble Π₀(ω) on H1.sites.

Keyword arguments

• ϵF : Fermi energy (physical units). Default 0.0.
• P_method : :purification (default) or :kpm — how to compute density matrices. With :purification, H._density_cache is reused if present.
• GF_method : :kpm (default) or :krylov — how to compute G_eff(ω).
• Ncheb : Chebyshev order (KPM methods only). Default 150.
• maxdim : Maximum bond dimension. Default 200.
• cutoff : SVD truncation cutoff. Default 1e-8.
• purify_method : :mcweeny (default) or :sp2.
• purify_maxdim : Max bond dim during purification. Default 40.
• purify_maxiters: Max purification iterations. Default 30.
• purify_tol : Purification convergence tolerance. Default 1e-5.
• η : Lorentzian broadening for the GF. Default 1e-3.
• krylov_nsweeps : DMRG sweeps for Krylov solver. Default 12.
• krylov_maxdim : Max bond dim for Krylov solver. Default 100.
• krylov_cutoff : SVD cutoff for Krylov solver. Default 1e-8.
• verbose : Print progress. Default false.

get_bubble_mpo_cheb2d(H1, H2, ωlist; Ncheb, maxdim, cutoff,
                       ϵF, P_method, purify_*, η, verbose) -> Vector{MPO}

Compute the non-interacting polarization bubble Π₀(ω) for each ω in ωlist using the double Chebyshev decomposition.

Instead of building the 2L-site effective Hamiltonian Heff = I⊗H₂ − H₁⊗I and running KPM on it (where bond dimension grows at each Chebyshev step due to entanglement between subsystems), this routine decomposes G_eff as

G_eff(ω) ≈ Σ_{mn} c_{mn}(ω) · T_m(H̃₁) ⊗ T_n(H̃₂)

where Tm, Tn are Chebyshev polynomials of the L-site rescaled Hamiltonians and c_{mn}(ω) are scalar 2D Chebyshev coefficients (cheap, via DCT-II).

The bubble on L-site MPOs is assembled as

Π₀(ω) = Σ_{mn} c_{mn}(ω) · D_{mn}

where D{mn} = (Tm(H̃₁)·P₁) ⊙ Tn(H̃₂) − Tm(H̃₁) ⊙ (Tn(H̃₂)·P₂) and ⊙ is the site-wise Hadamard product (`hadamardmpo`).

Online multi-ω sweep: All coefficient matrices C[m,n](ω) are precomputed at once (cheap DCT scalars). The (m,n) double loop runs once; each D{mn} is computed once and accumulated into every Π(ω) simultaneously using the scalar c{mn}(ω). This matches the KPM "online" paradigm: the expensive MPO work (Hadamard products) is done once and shared across all frequencies.

Keyword arguments

• Ncheb : Chebyshev expansion order. Default 50.
• maxdim : Max bond dimension throughout. Default 200.
• cutoff : SVD truncation cutoff. Default 1e-8.
• ϵF : Fermi energy. Default 0.0.
• P_method : :purification (default) or :kpm.
• purify_method : :mcweeny (default) or :sp2.
• purify_maxdim, purify_maxiters, purify_tol : purification controls.
• η : Lorentzian broadening. Default 1e-3.
• coeff_tol : Skip (m,n) pairs where |C[m,n]| < coeff_tol, and entire m rows where the row maximum is below coeff_tol. For smooth integrands (large η) this prunes most of the N² terms at negligible accuracy cost. Default 1e-12.
• verbose : Print progress. Default false.

get_bubble_mpo_cheb2d_tucker(H1, H2, ωlist; Ncheb, maxdim, cutoff,
                              ϵF, P_method, purify_*, η, coeff_tol,
                              tucker_tol, tucker_maxrank, kernel,
                              hooi_iters, verbose) -> Vector{MPO}

Tucker-accelerated variant of get_bubble_mpo_cheb2d.

Returns the full non-interacting polarization bubble Π₀(ω) as an MPO at each frequency in ωlist, suitable for RPA resummation and Wynn acceleration on the Dyson geometric series.

The Tucker-2 decomposition of the stacked coefficient tensor finds global bases Um (N×rm) and Vn (N×rn) satisfying

`C[m,n](ω) ≈ Σ_{s₁,s₂} G[s₁,s₂,ω] · (U_m)_{ms₁} · conj((V_n)_{ns₂})`

with C ≈ Um G(ω) Vn†. The rm + rn frequency-independent weighted MPO sums and the rm × rn Hadamard products are computed once; per-ω cost is only rm × rn cheap scalar-weighted MPO additions.

Speedup over get_bubble_mpo_cheb2d: N²→rm·rn Hadamard products.

Additional keyword arguments (beyond get_bubble_mpo_cheb2d):

• tucker_tol : relative singular-value cutoff for both mode SVDs. Default 1e-3.
• tucker_maxrank: hard cap on rm and rn. Default 20.
• kernel : :jackson (default) or :none. Jackson damping reduces Tucker rank.
• hooi_iters : HOOI refinement iterations after HOSVD initialisation. Default 3.

get_bubble_mpo_haydock(H1, H2, ωlist; N_steps, η, maxdim, cutoff,
                        ϵF, P_method, purify_method, purify_maxdim,
                        purify_maxiters, purify_tol, Ncheb, verbose)
    -> Vector{MPO}

Compute the bare polarization bubble Π₀(ω) as an L-site MPO for each frequency in ωlist using Haydock recursion on H_eff = I⊗H₂ − H₁⊗I.

The Krylov basis is built once from the seed Φ₀ = I⊗P₂ − P₁⊗I. For each ω, the resolvent (ω+iη−H_eff)⁻¹Φ₀ is recovered by solving the N×N Lanczos tridiagonal system and forming a linear combination of the stored basis MPOs.

Returns a Vector{MPO} compatible with rpa_wynn_from_bubbles.

Keyword arguments

• N_steps : Haydock recursion depth. Default 30.
• η : Lorentzian broadening. Default 1e-2.
• maxdim : Maximum bond dimension. Default 200.
• cutoff : SVD truncation cutoff. Default 1e-8.
• ϵF : Fermi energy. Default 0.0.
• P_method : :purification (default) or :kpm.
• purify_method : :mcweeny (default) or :sp2.
• purify_maxdim : Max bond dim during purification. Default 40.
• purify_maxiters: Max purification iterations. Default 30.
• purify_tol : Purification convergence tolerance. Default 1e-5.
• Ncheb : Chebyshev order for :kpm P_method. Default 150.
• verbose : Print progress. Default false.

get_bublle_expanded_from_Tn(Tn_list1, Tn_list2, Tn_listeff,
                             sites1, sites2, N, ω, ϵF;
                             a, maxdim) -> MPO

Compute the non-interacting polarization bubble Π₀(ω) as a 2L-site MPO using the Lehmann representation in terms of the Chebyshev moments.

The bubble is: Π₀(ω) = (P₁⊗I₂ − I₁⊗P₂) · Geff(ω) where P₁, P₂ are density matrices (filled-band projectors) and Geff is the retarded Green's function of H_eff = I⊗H₂ − H₁⊗I.

get_magnon_bubble(H, ω; ...) -> MPO

Non-interacting transverse spin polarization bubble Π₀^{+−}(ω) for a spinful TBHamiltonian.

The spin degree of freedom is projected out analytically: the spin-↑ and spin-↓ blocks of H become two independent L-qubit Hamiltonians H_↑ and H_↓, which are passed to get_bubble_mpo(H_↑, H_↓, ω). This corresponds to the Kubo S⁺S⁻ bubble

Π₀^{+−}(ω) = ∑_k (f_{k↓} − f_{k↑}) / (ω − (ε_{k↑} − ε_{k↓}) + iη)

Errors if H.spin_s === nothing. All get_bubble_mpo keyword arguments are accepted and forwarded.

get_magnon_susceptibility(H, MPOV, ω; ...) -> MPS

RPA transverse spin susceptibility χ^{+−}_RPA(ω) for a spinful TBHamiltonian.

Builds Π₀^{+−}(ω) via get_magnon_bubble, then solves the Dyson equation (I − Π₀ V) χ = Π₀. For a Hubbard-like interaction the interaction MPO is MPOV = U · Id on the orbital sites.

Keyword arguments

• rpa_nsweeps, rpa_maxdim, rpa_cutoff : Dyson linsolve parameters.
• All get_bubble_mpo keywords forwarded via kwargs....

get_magnon_susceptibility_wynn(H, MPOV, ωlist; K_max, ...) -> (chi_partial, chi_wynn)

Wynn ε-accelerated transverse spin RPA susceptibility over a frequency list.

The spin-↑/↓ sector projections are performed once outside the ω loop; the Neumann-series / Wynn logic is identical to get_rpa_susceptibility_wynn.

Returns (chi_partial, chi_wynn) with the same layout as get_rpa_susceptibility_wynn.

Keyword arguments

• K_max, maxdim_apply, cutoff_apply, verbose : series / Wynn control.
• All get_bubble_mpo keywords forwarded via kwargs....

get_rpa_susceptibility(H::TBHamiltonian, MPOV, ω; mode, ...) -> MPS

Compute the RPA susceptibility χ(ω) for a system described by H with interaction MPO MPOV. Returns a 2L-site MPS encoding the diagonal χ_{iijj}^RPA.

mode keyword

• :charge (default) — density–density bubble χ^{ρρ}: calls get_bubble_mpo(H, H, ω).
• :magnetic — transverse spin bubble χ^{+−}: projects H onto its spin-↑ and spin-↓ blocks and calls get_magnon_bubble. Requires H.spin_s !== nothing.

Internally solves the Dyson equation (I − Π₀V) χ = Π₀. All get_bubble_mpo keyword arguments are accepted and forwarded.

Additional keywords (Dyson solve)

• rpa_nsweeps : sweeps for the RPA linsolve. Default 20.
• rpa_maxdim : max bond dim for the RPA linsolve. Default 400.
• rpa_cutoff : cutoff for the RPA linsolve. Default 1e-8.

get_rpa_susceptibility_wynn(H, MPOV, ωlist; mode, K_max, maxdim_apply,
                             cutoff_apply, verbose, <bubble kwargs>)
                             -> (chi_partial, chi_wynn)

Compute the RPA susceptibility χ_RPA(q,ω) for all frequencies in ωlist using the Wynn ε-algorithm for Padé acceleration of the geometric (bubble) series.

mode keyword

• :charge (default) — density–density channel; uses get_bubble_mpo(H, H, ω).
• :magnetic — transverse spin channel S⁺S⁻; uses get_magnon_bubble(H, ω). The spin-↑/↓ projections are performed once before the ω loop. Requires H.spin_s !== nothing.

Key idea: instead of inverting (I − Π₀V), build the Neumann series T₀ = Π₀, Tₙ = Tₙ₋₁·V·Π₀ (so Σ Tₙ → χRPA as K→∞), extract scalars sₙ(q,ω) = −Im⟨q|Tₙ(ω)|q⟩/π via `getspect_k`, and apply Wynn ε to the partial-sum sequence per (q,ω) for fast convergence.

Returns

• chi_partial[k+1, i_ω, q] : partial sum Σₙ₌₀ᵏ sₙ(q,ω)
• chi_wynn[m, i_ω, q] : Wynn ε_{2m}(0) estimate (uses 2m+1 terms)

Keyword arguments

• mode : :charge (default) or :magnetic.
• K_max : highest order in the series (total K_max+1 terms). Default 6.
• maxdim_apply : bond dim for the Tₙ·V·Π₀ products. Default 200.
• cutoff_apply : truncation cutoff for those products. Default 1e-8.
• verbose : print per-ω progress. Default false.
• All get_bubble_mpo keywords (ϵF, P_method, GF_method, Ncheb, maxdim, cutoff, purify_*, η, krylov_*) are accepted and forwarded.

get_spect_k(W; tol, maxdim) -> Vector{ComplexF64}

Extract the k-space diagonal of MPO W as a dense vector of 2^L values.

Conjugates W by the QFT (giving W̃ = QFT·W·QFT†), extracts the diagonal as an MPS, then evaluates each ⟨k|W̃|k⟩ for the 2^L quantics k-indices. Uses LSB-first quantics convention (site 1 = least significant bit).

haydock_cf(H_mpo, seed, N_steps; maxdim, cutoff, verbose)
    -> (a, b, basis, norm0)

Haydock (Lanczos) recursion with Hmpo acting on MPO vectors from the left. Starting from seed, builds an orthogonal Krylov basis under Hmpo using the Frobenius (Hilbert-Schmidt) inner product (A, B) = Tr[A† B].

Three-term recurrence (Φ₀ = seed / β₀, β₀ = ||seed||_F):

Φₙ₊₁ = H·Φₙ − aₙ·Φₙ − bₙ·Φₙ₋₁    (b₁ = 0)

Returns:

• a : diagonal coefficients a[1..N]
• b : b[1] = norm0 = ||seed||_F; b[2..N] = off-diagonal βₙ
• basis: normalized Krylov MPOs {Φ₀, …, Φₙ₋₁}
• norm0: sqrt(inner(seed, seed))

The scalar projected GF ⟨seed|(z−H)⁻¹|seed⟩ is recovered via eval_haydock_cf(a, b, z). The full resolvent MPO (z−H)⁻¹|seed⟩ is recovered via haydock_resolve_mpo(a, b, basis, z).

haydock_resolve_mpo(a, b, basis, z; maxdim, cutoff) -> MPO

Reconstruct (z−H)⁻¹|seed⟩ as an MPO by solving the N×N Lanczos tridiagonal system and forming a linear combination of the Krylov basis MPOs:

(z·I − T) c = b[1]·e₁,   Π₀(z) = Σₙ c[n]·basis[n]

where T has diagonal a and off-diagonal b[2:], and b[1] = norm0.

interleave_mpo(target_mpo, phys_sites, n) -> MPO

Embed an L-site MPO into a 2L-site space by interleaving it with identity operators. phys_sites must have length 2L.

• n = 0 : operator sits at even positions (2, 4, 6, …), identities at odd
• n = 1 : operator sits at odd positions (1, 3, 5, …), identities at even

Note: phys_sites must be interleaved as [A[1], B[1], A[2], B[2], …] so that each operator site lands on an index with the correct dimension. For heterogeneous site spaces (layer, sublattice, …), use interleave_mpo_tb which builds the interleaved site list automatically.

interleave_mpo_tb(op, sites_A, sites_B, which) -> MPO

Generalization of interleave_mpo for heterogeneous site spaces (Layer, Qubit, Honeycomb, etc.). Embeds an N-site MPO op into the 2N-site product space whose physical indices are ordered as

[sites_A[1], sites_B[1], sites_A[2], sites_B[2], …, sites_A[N], sites_B[N]]

so every consecutive pair (sites_A[i], sites_B[i]) shares the same dimension regardless of site type.

• which = :A : op acts on sites_A (odd positions), identity on sites_B
• which = :B : op acts on sites_B (even positions), identity on sites_A

merge_mps_to_mpo(mps) -> MPO

Contract each consecutive pair (2i-1, 2i) of an MPS into a single MPO tensor. The resulting MPO has length(mps) ÷ 2 sites.

mpo_kron(A, B) -> MPO

Concatenate two MPOs into a single MPO on the combined site space, joined by a bond-dimension-1 link.

mps_kron(A, B) -> MPS

Concatenate two MPS into a single MPS on the combined site space, joined by a bond-dimension-1 link.

rpa_from_bubble_diag(Π, MPOV, finalsites, finalfinalsites;
                     nsweeps, maxdim, cutoff) -> MPS

Solve the RPA Dyson equation (I − Π₀V) χ = Π₀ for the interacting susceptibility χ using DMRG-style linear solve.

Returns a 2L-site MPS encoding the diagonal χ_{iijj}^RPA.

rpa_wynn_from_bubbles(Π0_list, MPOV; K_max, maxdim_apply, cutoff_apply, verbose)
-> (chi_partial, chi_wynn)

Wynn ε-accelerated RPA susceptibility from a pre-computed list of bubble MPOs.

Accepts the output of get_bubble_mpo_cheb2d_tucker (or any Vector{MPO}) directly, skipping the internal bubble calculation. All other logic is identical to get_rpa_susceptibility_wynn: Neumann series T₀ = Π₀, Tₙ = Tₙ₋₁·V·Π₀, followed by per-k-point Wynn ε-acceleration of the partial-sum sequence.

Returns

• chi_partial[k+1, iω, q] : partial sum Σₙ₌₀ᵏ (−Im⟨q|Tₙ(ω)|q⟩)
• chi_wynn[m, iω, q] : Wynn ε_{2m}(0) estimate

Keyword arguments

• K_max : series order (K_max+1 terms). Default 6.
• maxdim_apply : bond dim for Tₙ·V·Π₀ products. Default 200.
• cutoff_apply : truncation cutoff for those products. Default 1e-8.
• verbose : print per-ω progress. Default false.

swap_every_other_legs(MPOin, newsites) -> MPO

Replace site indices and additionally swap bra↔ket on every even-numbered site. Used to convert the 2L-site bubble MPO from the interleaved ordering into the form expected by collapse_mpo_pairs.

wynn_epsilon(s) -> Vector{ComplexF64}

Wynn ε-algorithm applied to a scalar sequence s = [s₀, s₁, ..., sK]. Returns the even-column first-row Padé estimates [ε₂(0), ε₄(0), ...]. Uses only additions and reciprocals — no matrix operations.

_get_projector(H; method, fermi, Nchebychev, maxdim, cutoff, Nel) -> MPO

Compute or retrieve the ground-state projector P for H.

• method=:KPM (default): uses the cached H._tn_cache if present; otherwise runs KPM_Tn(H, Nchebychev). The Fermi level fermi (in physical units) is rescaled internally.
• method=:mcweeny: returns H._density_cache if set; otherwise runs McWeeny purification. fermi is ignored.
• method=:sp2: same but uses SP2 purification. Nel sets the target electron count (default: H.N ÷ 2).
• maxdim, cutoff: bond dimension and truncation threshold forwarded to the underlying method.

get_C(H::TBHamiltonian, xfunc=nothing, yfunc=nothing;
      method=:KPM, fermi=0.0, l=nothing, Λ=10,
      Nchebychev=300, maxdim=500, cutoff=1e-8,
      Nel=nothing, quenched=true) -> Function

High-level wrapper: compute the ground-state projector via method and return the Chern marker closure from get_C_op_MPO_from_P.

xfunc(i, L_chain) and yfunc(i, L_chain) accept a 0-indexed physical site number i and return raw x/y coordinates. Both default to nothing, in which case they are auto-derived:

• If H.geometry_uc is set (sublattice models: honeycomb, kagome, lieb, dice, sshsublattice): uses `geometryuc(i+1)[1/2]`, which returns the same Bravais unit-cell position for all sublattice atoms in the same UC.
• Otherwise falls back to H.geometry(i+1)[1/2].

When xfunc/yfunc are auto-derived for a sublattice model, α in the returned closure ranges over 1 … n_sub·2^L (all physical sites).

Reuses H._tn_cache or H._density_cache when available. maxdim and cutoff are forwarded uniformly to the projector computation and to all MPO multiplications in the Chern marker assembly.

See get_C_op_MPO_from_P for full documentation of the remaining arguments.

Returns

calculate_chern_number(uc::Int) -> ComplexF64 where uc is a 1-indexed unit cell number; the closure sums the marker over all n_sub sublattice atoms in that UC. Take real(·) for the density.

get_C_op_MPO_from_P(P, L, sites, xfunc, yfunc;
                    l=nothing, Λ=10, maxdim=500, cutoff=1e-8,
                    quenched=true) -> Function

Build the real-space Chern marker and return a closure calculate_chern_number(α) that evaluates it at any physical site α (1-indexed).

Coordinate functions

xfunc(i, L_chain) and yfunc(i, L_chain) accept a 0-indexed physical site number i and return the raw coordinate (not yet quenched).

• Plain models (length(sites) == L): i ∈ 0…2^L−1. Examples: xfunc(i, Lchain) = Float64(mod(i, Lchain)) # square x yfunc(i, Lchain) = Float64(div(i, Lchain)) # square y
• Sublattice models (length(sites) == L+1, last index has dim n_sub): i ∈ 0…n_sub·2^L−1. The function should return the same coordinate for all n_sub atoms within the same unit cell — typically the Bravais position. The auto-derived functions from H.geometry_uc satisfy this automatically. Internally, position MPOs are built on the L position qubits only and extended to the full chain via postpend_op(⋅, sub_s, I).

Quenched mode (quenched=true, default)

Position operators are quenched: sin(xfunc/Λ), cos(xfunc/Λ), and similarly for y. The Chern marker uses a 4-term trig decomposition that pre-computes 4 α-independent MPOs (C1–C4) and combines them per site in the closure:

C(α) = 2πi Λ² [ cos_xα cos_yα ⟨α|C1|α⟩ + sin_xα sin_yα ⟨α|C2|α⟩
               − cos_xα sin_yα ⟨α|C3|α⟩ − sin_xα cos_yα ⟨α|C4|α⟩ ]

Exploits sin(θ_r − θ_α) = sinθ_r cosθ_α − cosθ_r sinθ_α to reduce cost from O(N) MPO products to 4 products computed once. Λ² restores physical units.

Flat mode (quenched=false)

Position operators use xfunc/yfunc directly:

C_op = 2πi · (Q x P y Q − P x Q y P)

Most accurate for OBC systems or bulk-averaged quantities (no per-site centring).

Arguments

• P : ground-state projector MPO (over sites)
• L : number of position qubits; system has 2^L unit cells
• sites : full ITensor site index list (length == L or L+1 for sublattice)
• xfunc, yfunc : coordinate functions; i is 0-indexed over all physical sites
• l : qubits per spatial direction; inferred as L ÷ 2 if nothing
• Λ : quenching period (quenching angle = coord / Λ)
• maxdim : MPO bond dimension during all multiplications
• cutoff : truncation threshold during all multiplications and subtractions
• quenched: true = 4-term sin/cos decomposition; false = flat operators

Returns

calculate_chern_number(uc::Int) -> ComplexF64 where uc is a 1-indexed unit cell number (1 … 2^L). For sublattice models the value is the sum of the Chern marker over all n_sub atoms within that unit cell. Take real(·) for the Chern number density.

Example — quenched square lattice

L_chain = 2^(L ÷ 2)
xfunc(i, _) = Float64(mod(i, L_chain))
yfunc(i, _) = Float64(div(i, L_chain))
C_at  = get_C_op_MPO_from_P(P, L, sites, xfunc, yfunc; Λ=L_chain, maxdim=100)
chern = real(sum(C_at(α) for α in 1:2^L)) / L_chain^2

Example — honeycomb via get_C (auto-derived geometry)

C_at  = get_C(H)   # xfunc/yfunc from H.geometry_uc automatically
N_sub = 2 * H.N
chern = real(sum(C_at(α) for α in 1:N_sub)) / (2^(H.L ÷ 2))^2

get_W(H::TBHamiltonian, xfunc=nothing;
      method=:KPM, fermi=0.0, Nchebychev=300,
      maxdim=15, cutoff=1e-8, Nel=nothing,
      quenched=true, l=nothing, Λ=10) -> Function

Compute the real-space winding number and return a closure calculate_winding(α::Int) -> ComplexF64 that evaluates the local winding number density at any site α (1-indexed).

The winding number operator is

W_op = σ_z · (P x̂ Q + Q x̂ P)

where P is the occupied-band projector, Q = I − P, x̂ is the position operator, and σ_z is the sublattice chirality (A → +1, B → -1).

H must have a 2-component sublattice index (H.sublattice_s with dim 2). σ_z is built automatically as the diagonal operator diag(+1, −1) on that index tensored with identity on the position qubits.

Coordinate function

xfunc(i, L_chain) -> Float64 accepts a 0-indexed physical site number and returns the raw x-coordinate. Defaults to nothing, in which case it is auto-derived from H.geometry_uc (preferred) or H.geometry: xfunc(i, _) = geom(i+1)[1]. Because geometry_uc returns the same Bravais position for both sublattice atoms in a unit cell, this correctly assigns the same x-coordinate to both A and B sites of each UC.

Quenched vs flat mode

• quenched=true (default): pre-computes two α-independent MPOs W1 and W2. W(α) = Λ · ⟨α| cos(xα/Λ) W1 − sin(xα/Λ) W2 |α⟩
• quenched=false: builds a global position operator and returns a closure over the resulting W_op MPO.

Arguments

• method : :KPM, :mcweeny, or :sp2 (see _get_projector).
• fermi : Fermi level in physical energy units (KPM only).
• Nchebychev: Chebyshev order when method=:KPM and no cache is present.
• maxdim : MPO bond dimension during all multiplications.
• cutoff : truncation threshold during all multiplications.
• Nel : target electron count for SP2 (default H.N ÷ 2).
• l : qubits per direction; inferred as H.L ÷ 2 if nothing.
• Λ : quenching period (angle = xfunc/Λ); sets the Λ prefactor.

Returns

calculate_winding(uc::Int) -> ComplexF64 where uc is a 1-indexed unit cell number (1 … 2^L). The value is the sum of the winding marker over both sublattice atoms (A and B) within that unit cell.

get_cosx_op(L, sites, L_chain, Λ, xfunc) -> MPO

Diagonal MPO for cos(xfunc(i, L_chain) / Λ). See get_sinx_op.

get_cosy_op(L, sites, L_chain, Λ, yfunc) -> MPO

Diagonal MPO for cos(yfunc(i, L_chain) / Λ). See get_sinx_op.

get_pump_xop(L, sites, xfunc; quenched=false, Λ=nothing) -> MPO

Diagonal position operator MPO for the Thouless pump formula.

xfunc(i, N) accepts a 0-indexed site index i ∈ 0…N-1 and the chain length N = 2^L, and returns the raw coordinate. For a 1-indexed chain: xfunc(i, N) = Float64(i + 1).

• quenched=false (default): diagonal entries are xfunc(i, N) directly.
• quenched=true: entries are Λ * sin(xfunc(i, N) / Λ), which smooths the discontinuity at PBC at the cost of a Λ prefactor. Λ defaults to N (one full period), giving sin(x/N) * N ≈ x for x ≪ N.

get_sinx_op(L, sites, L_chain, Λ, xfunc) -> MPO

Diagonal MPO for sin(xfunc(i, L_chain) / Λ) over the L-qubit position chain given by sites. xfunc(i, L_chain) receives a 0-indexed site/UC number and returns the raw x-coordinate.

get_siny_op(L, sites, L_chain, Λ, yfunc) -> MPO

Diagonal MPO for sin(yfunc(i, L_chain) / Λ). See get_sinx_op.

get_thouless_pump(H_of_t, Nt, T, xfunc;
                  P_method=:mcweeny, Nchebychev=200,
                  maxdim=100, cutoff=1e-8,
                  quenched=false, Λ=-1.0,
                  r_center=nothing, Nel=nothing, verbose=false) -> Float64

High-level Thouless pump: build P(t_k) for k = 0…Nt-1 via P_method, then compute the M1Q invariant C = M1Q(T) − M1Q(0).

Arguments

• H_of_t : t -> TBHamiltonian — all calls must share the same site indices (pass ref_sites to get_Hamiltonian inside the factory).
• Nt : number of time steps.
• T : period of the pump cycle.
• xfunc : coordinate function (i, N) -> Float64, 0-indexed.
• P_method : :mcweeny, :sp2, or :KPM.
• r_center : 0-indexed bulk site for M1Q evaluation; defaults to N ÷ 2.
• quenched : false = flat x̂; true = sin-quenched (removes PBC discontinuity).
• verbose : print progress.

get_valley_C(H, xfunc=nothing, yfunc=nothing;
             valley=:K, method=:mcweeny, fermi=0.0, l=nothing, Λ=10,
             Nchebychev=300, maxdim=500, cutoff=1e-8,
             Nel=nothing, quenched=true) -> Function

Compute the valley-resolved Chern marker and return a closure calculate_valley_chern(uc::Int) -> ComplexF64.

The valley operator V is built automatically via get_valley_operator(H). Its sign S = sign(V) is computed via sign_mpo (McWeeny purification), giving eigenvalues in {−1, +1}. The occupied K-valley projector is then

PK = P · (I ± S) / 2 · P

where P is the ground-state projector and the sign is chosen by valley (:K → +, :K_prime → −). The PK is symmetrized and passed as pk_mpo to get_C_op_MPO_from_P, which evaluates the Chern marker as ⟨α|PK C_op PK|α⟩.

Arguments

• valley : :K or :K_prime.
• use_sign : if true (default), sharpen V to eigenvalues {−1,+1} via sign_mpo before forming the valley projector. Set to false to use V directly (cheaper but less accurate when the valley operator spectrum is not already close to ±1).

All remaining arguments are forwarded to _get_projector and get_C_op_MPO_from_P; see those functions for documentation.

get_valley_operator(H::TBHamiltonian; maxdim=500, cutoff=1e-8) -> MPO

Build the valley operator V for a 2D honeycomb Hamiltonian.

V is constructed as the Haldane NNN Hamiltonian (φ = π/2, NN term zeroed) multiplied by an additional sublattice sign η_i: +1 on sublattice A (index 1), −1 on sublattice B (index 2). The global prefactor is −i/(3√3):

t(dx,dy,fs,ts) = (−i / 3√3) · ν_{dx,dy,fs,ts} · η_{fs}

where ν ∈ {±1} is the Haldane chirality (counterclockwise = +1).

H must be a 2D honeycomb model with H.Lx set and a 2-component sublattice index. The returned MPO shares the same site indices as H.mpo.

get_valley_projectors(V_mpo, sites; maxdim=500, cutoff=1e-8) -> (PK, PK_prime)

Return the two valley projectors from the valley operator V_mpo:

PK       = (I + V) / 2   (K  valley)
PK_prime = (I − V) / 2   (K′ valley)

thouless_pump(P_array, dt, x_op, sites; r_center, maxdim, cutoff,
              verbose, return_trajectory) -> Float64 or (Float64, Vector{Float64})

Compute the Thouless pump invariant (Chern number) using the local M1Q marker:

M1Q(r, t) = ⟨r| P(t) U†(t) x̂ U(t) P(t) |r⟩

where U(t) is the adiabatic evolution operator generated by h(t) = [∂P/∂t, P(t)], propagated with a second-order Taylor step:

U(k) = (I + h(k−1) dt + h(k−1)² dt²/2) U(k−1),  U(0) = I

The invariant is C = M1Q(r, T) − M1Q(r, 0), evaluated at site r_center. Finite differences for ∂P/∂t use central differences (one-sided at endpoints).

Arguments

• P_array : Vector{MPO} of length Nt, one instantaneous projector per step.
• dt : time step (T / Nt).
• x_op : position operator MPO from get_pump_xop.
• sites : physical site indices (for identity MPO and matrix_checker).
• r_center : 0-indexed bulk site at which to evaluate M1Q.
• maxdim : max bond dimension for all MPO operations.
• cutoff : SVD truncation threshold.
• verbose : print M1Q(0), M1Q(T), and bond dimension at each step.
• return_trajectory : if true, return (C, M1Q_traj) where M1Q_traj is a Vector{Float64} of length Nt+1 with M1Q at each step (index 1 = t=0, index k+1 = t=k·dt). Default false.

_idempotency_error(ρ, ρ2) -> Float64

Compute ‖ρ² - ρ‖ / ‖ρ‖ as a measure of how far ρ is from a projection. Returns 0 for an exact density matrix.

_mpo_sq(ρ, sites; maxdim, cutoff) -> MPO

Compute ρ² via apply and truncate immediately. The intermediate bond dimension of apply is controlled by maxdim.

get_density(H::TBHamiltonian; method, ϵF, Ncheb, kernel, lambda,
            maxdim, cutoff, Nel, maxiters, tol, verbose) -> MPO

Compute and cache the zero-temperature density matrix P = θ(ϵF − H).

If H._density_cache is already populated it is returned immediately. Set H._density_cache = nothing to force a fresh computation.

method

• :mcweeny (default) — McWeeny purification P{n+1} = 3Pn² − 2P_n³
• :sp2 — SP2 purification (1 MPO product/step), requires Nel
• :kpm — KPM Chebyshev expansion of the Fermi step function

Keyword arguments

• ϵF : Fermi energy in physical units (:kpm only). Default 0.0.
• Ncheb : Chebyshev order (:kpm only). Default 150. Reuses H._tn_cache if already built at order ≥ Ncheb; otherwise calls KPM_Tn to build and cache it.
• kernel : KPM kernel — :jackson (default). HODC is not meaningful for the step function so only convolution kernels are supported.
• lambda : Jackson kernel damping parameter. Default 4.0.
• maxdim : Maximum bond dimension. Default 40.
• cutoff : SVD truncation cutoff. Default 1e-8.
• Nel : Target electron count (:sp2 only). Default H.N ÷ 2.
• maxiters: Maximum purification iterations. Default 30.
• tol : Idempotency convergence tolerance (purification). Default 1e-5.
• verbose : Print iteration progress. Default false.

get_dos_drho(H::TBHamiltonian, ω; dmu=0.05, maxdim=40, cutoff=1e-8,
             maxiters=30, tol=1e-5, verbose=false) -> Float64

get_dos_drho(H::TBHamiltonian, ωs::AbstractVector; ...) -> Vector{Float64}

Compute the total density of states at energy ω as the finite-difference derivative of Tr[ρ(μ)] with respect to μ:

DOS(ω) ≈ (Tr[ρ(ω + δμ)] − Tr[ρ(ω − δμ)]) / (2δμ)

This is the cheapest variant: only the scalar trace of each purified density matrix is needed, so no MPO or MPS subtraction is performed. All other arguments are identical to get_ldos_drho.

get_ldos_drho(H::TBHamiltonian, ω; mode=:mpo, dmu=0.05, maxdim=40,
              cutoff=1e-8, maxiters=30, tol=1e-5, verbose=false) -> MPO or MPS

get_ldos_drho(H::TBHamiltonian, ωs::AbstractVector; ...) -> Vector{MPO} or Vector{MPS}

Compute the local density-of-states at energy ω as the finite-difference derivative of the McWeeny density matrix with respect to μ:

A(ω) ≈ [ρ(ω + δμ) − ρ(ω − δμ)] / (2δμ)

where ρ(μ) = θ(μ − H) is the purified density matrix at Fermi level μ.

Modes

• mode=:mpo (default): returns the full MPO A(ω). Its diagonal elements give the LDoS: LDoS(r, ω) = ⟨r|A(ω)|r⟩. Subtraction and truncation are performed on the full MPO bond structure.

• mode=:mps: extracts the diagonal of each ρ as an MPS via extract_diagonal_to_mps before taking the difference. The returned MPS encodes diag(ρ(ω+δμ)) − diag(ρ(ω−δμ)) scaled by 1/(2δμ). This is cheaper than :mpo since the subtraction and all subsequent operations stay in the MPS bond space (no MPO–MPO off-diagonal contributions).

Arguments

• ω : energy (or ωs: vector) at which to evaluate LDoS
• mode : :mpo (full operator) or :mps (diagonal only, cheaper)
• dmu : finite-difference step / broadening
• maxdim : bond dimension for purification and the final combination
• cutoff : SVD truncation threshold throughout
• maxiters: maximum McWeeny iterations per purification call
• tol : idempotency convergence threshold ‖ρ²−ρ‖/‖ρ‖
• verbose : print McWeeny residuals

mcweeny_purify(ρ0; maxiters=30, maxdim=40, cutoff=1e-8,
              tol=1e-5, verbose=false) -> MPO

Iterate the McWeeny map ρ{n+1} = 3ρn² - 2ρ_n³ until the idempotency residual ‖ρ² - ρ‖/‖ρ‖ < tol or maxiters is reached.

Arguments

• ρ0 : initial density matrix MPO (need not be idempotent)
• maxiters : maximum number of iterations
• maxdim : maximum MPO bond dimension during multiplication
• cutoff : cutoff during truncation
• tol : convergence threshold on ‖ρ²−ρ‖/‖ρ‖
• verbose : print residual at each iteration

Returns

Purified density matrix MPO.

mcweeny_purify(H::TBHamiltonian; ϵF, maxiters, maxdim, cutoff, tol, verbose) -> MPO

High-level overload: builds the initial guess from H, runs McWeeny purification, caches the result in H._density_cache, and returns the purified density matrix.

ϵF shifts the Fermi level of the initial guess ρ₀ = (I − (H − ϵF·I)/scale) / 2, allowing purification to target a band other than half-filling. Default 0.0.

purification_initial_guess(H_mpo, scale, sites; maxdim=40, cutoff=1e-8) -> MPO
purification_initial_guess(H::TBHamiltonian; maxdim=40, cutoff=1e-8) -> MPO

Construct the simplest valid initial guess for purification:

ρ₀ = (I - (H − center·I)/scale) / 2

This maps the rescaled spectrum ∈ [-1, 1] to ρ₀ eigenvalues ∈ [0, 1], the required input range for both mcweeny_purify and sp2_purify.

The TBHamiltonian overload calls _ensure_scale! automatically and accounts for a non-zero spectral center.

sign_mpo(A::MPO, sites; scale=1.0, maxdim=500, cutoff=1e-8,
         maxiters=30, tol=1e-5, verbose=false) -> MPO

Compute sign(A) for a Hermitian MPO A via two McWeeny purifications:

sign(A) = θ(A) − θ(−A)

where θ(A) is the projector onto the positive-eigenvalue subspace of A. The two initial guesses

ρ₀₊ = (I + A/scale) / 2   →  McWeeny  →  θ(A)
ρ₀₋ = (I − A/scale) / 2   →  McWeeny  →  θ(−A)

map eigenvalues of ±A from [−scale, +scale] to [0, 1] before purification. scale must be ≥ the spectral radius of A; it defaults to 1.0, which is correct when A has already been constructed with normalised eigenvalues (e.g. from get_valley_operator whose spectrum lies in (−1, +1)).

The returned MPO has eigenvalues in {−1, +1}.

sp2_purify(ρ0, Nel; maxiters=40, maxdim=40, cutoff=1e-8,
           tol=1e-5, verbose=false) -> MPO

Iterate the SP2 map until convergence:

if Tr(ρ_n²) ≥ N_el:   ρ_{n+1} = ρ_n²          (contract)
else:                  ρ_{n+1} = 2ρ_n - ρ_n²   (expand)

Each step costs one MPO-MPO product. The direction rule drives Tr(ρ) toward Nel and simultaneously pushes eigenvalues to 0 or 1. Convergence is quadratic.

The spectrum of ρ0 must lie in [0, 1]; normalise with ρ0 = (Id - H/scale) / 2 if starting from scratch.

Arguments

• ρ0 : initial density matrix MPO with spectrum ⊆ [0,1]
• Nel : target electron number (Tr of the converged projector)
• remaining kwargs: same as mcweeny_purify

sp2_purify(H::TBHamiltonian; Nel, maxiters, maxdim, cutoff, tol, verbose) -> MPO

High-level overload: builds the initial guess from H, runs SP2 purification, caches the result in H._density_cache, and returns the purified density matrix. Nel defaults to half-filling (H.N ÷ 2).

Exciton_Hamiltonian(H_c, H_v, Ufunc; on_site, tol_quantics, maxbonddim_quantics,
                    tol, cutoff, maxdim) -> MPO

Build the exciton Hamiltonian on the interleaved 2L-site electron-hole space:

H_exc = (H_c ⊗ I_h − I_e ⊗ H_v) + U

where U is the contact interaction diagonal MPO built from Ufunc.

H_c and H_v are TBHamiltonian objects for the electron and hole single-particle sectors (any geometry: "chain_1d", "square_2d", etc.). Both must have the same L and distinct site indices. Ufunc(x) gives the interaction strength at site x ∈ {1, …, 2^L} (1-indexed).

on_site keyword (optional): a function V(x) representing the conduction band edge modulation, compressed via QTCI. Applied as +V to the electron sector and −V to the valence sector, so that the hole also feels +V after the H_c − H_v subtraction (type-I semiconductor convention: both carriers confined). tol_quantics and maxbonddim_quantics control the QTCI compression of V.

Examples

# 1D, uniform hopping, contact interaction, Gaussian confinement
H_c = get_Hamiltonian("chain_1d", t; L=L)
H_v = get_Hamiltonian("chain_1d", t; L=L)
H_exc = Exciton_Hamiltonian(H_c, H_v, x -> -U;
                             on_site = x -> -V0 * exp(-((x - N/2)^2) / (2σ^2)))

# 2D square lattice
H_c = get_Hamiltonian("square_2d", t; L=Lx+Ly, Lx=Lx, Ly=Ly)
H_v = get_Hamiltonian("square_2d", t; L=Lx+Ly, Lx=Lx, Ly=Ly)
H_exc = Exciton_Hamiltonian(H_c, H_v, Ufunc; on_site = x -> dot_potential(x))

build_interaction_op_exciton(L, sites, Ufunc) -> MPO

Build the electron-hole interaction MPO on the 2L-site interleaved space. Ufunc(x) gives the interaction strength when electron and hole are both at position x (contact interaction). sites must be the 2L-site interleaved site index vector.

exciton_hamiltonian(geometry, params, Ufunc; L, [on_site, scale, tol_quantics,
    maxbonddim_quantics, tol, cutoff, maxdim, kwargs...]) -> TBHamiltonian

exciton_hamiltonian(H_c, H_v, Ufunc; [on_site, scale, tol_quantics,
    maxbonddim_quantics, tol, cutoff, maxdim]) -> TBHamiltonian

Build an exciton Hamiltonian and wrap it in a TBHamiltonian for use with TensorBinding's KPM, DMRG, and spectral tools.

Site encoding (2L sites total, L = H_c.L):

• Odd sites (1, 3, …) : electron position qubits
• Even sites (2, 4, …) : hole position qubits (interleaved)

TBHamiltonian.L = L counts position qubits per sector; TBHamiltonian.sites holds all 2L interleaved MPO sites.

Two calling modes

1. Geometry string — both bands built from the same geometry and parameters via get_Hamiltonian; any extra kwargs are forwarded to it:

   H = exciton_hamiltonian("square_2d", t, Ufunc; L=Lx+Ly, Lx=Lx, Ly=Ly,
                            on_site = x -> V(x))

2. Pre-built sectors — pass explicit TBHamiltonian objects for the electron (H_c) and hole (H_v) bands when the two sectors differ (different hopping, disorder, external fields):

   H_c = get_Hamiltonian("chain_1d", t_c; L=L)
   H_v = get_Hamiltonian("chain_1d", t_v; L=L)
   H   = exciton_hamiltonian(H_c, H_v, Ufunc)

Keyword arguments

• on_site : conduction band edge modulation V(x), 1-indexed. Applied as +V to electron and −V to valence sector (type-I confinement). Compressed via QTCI.
• scale : exciton spectral half-bandwidth (0.0 → lazy DMRG estimate).
• tol_quantics : QTCI tolerance for Ufunc and on_site. Default 1e-8.
• maxbonddim_quantics : QTCI max bond dimension. Default 100.
• tol : MPO assembly truncation tolerance. Default 1e-8.
• cutoff : SVD cutoff for MPO arithmetic. Default 1e-8.
• maxdim : max bond dimension of the final MPO. Default 200.

hopping_mpo_exciton(H_c, H_v) -> MPO

Embed the electron (H_c) and hole (H_v) single-particle Hamiltonians into the interleaved 2L-site electron-hole space:

H_kin = H_c ⊗ I_h  -  I_e ⊗ H_v

H_c sits at odd sites (1, 3, 5, …) and H_v at even sites (2, 4, 6, …). Both TBHamiltonian objects must have the same L and distinct site indices.

mpsexcitonKQ(k, Q, sites) -> MPS

Product momentum-basis electron-hole continuum state at fixed total momentum:

|k, Q - k>

Both k and Q are 1-indexed (1:2^LPhys) in the public API. Internally they are shifted to 0-indexed modular arithmetic, and the hole momentum is computed as h = Q - k (mod N).

This is the incoherent continuum probe complement to mpsexcitonQ. mpsexcitonQ(Q, sites) builds the coherent pair state (1 / sqrt(N)) * sum_k |k, Q-k>, while mpsexcitonKQ(k, Q, sites) builds one rank-1 product basis vector in that fixed-Q sector. Stochastic continuum traces sample many such k values and average the resulting MPS-KPM spectra.

The site convention matches the QFT momentum convention: site pair 1 carries bit 0 (LSB-first), with electron bits on odd sites and hole bits on even sites.

mpsexcitonQ(Q, sites) -> MPS

Normalized fixed-total-momentum exciton state on the interleaved electron-hole chain:

|Q> = (1 / sqrt(N)) * sum_k |k, -k + Q>

Q is 1-indexed (Q in 1:2^LPhys), matching mpsexciton and the rest of the public exciton sampling API. Internally the modular momentum arithmetic is 0-indexed and uses the QFT momentum convention: site pair 1 carries bit 0 (LSB-first).

The construction encodes the constraint k + h = Q (mod N) as a small finite-state MPS rather than explicitly summing over all N = 2^LPhys momenta. For each electron-hole bit pair (a_i, b_i) it keeps a binary carry on the MPS bond and allows only local configurations satisfying

a_i + b_i + carry_in = Q_i + 2 * carry_out

where Q_i is the ith bit of Q - 1. The first bit has carry_in = 0, and the final carry is not fixed, giving addition modulo 2^LPhys.

Scaling: the probe has 2LPhys physical sites, bond dimension at most 2 between bit pairs (the carry), and a local intermediate link of dimension 4 inside each electron-hole pair. Building the probe is therefore O(LPhys) in storage and time, not O(2^LPhys). For very large systems the expensive part is the subsequent MPS-KPM recursion with H_QFT; this compact probe construction should not be the bottleneck.

mpsexcitonQTrace(Q, sites; rng=Random.default_rng()) -> MPS

Random-phase fixed-total-momentum trace probe for the electron-hole continuum:

|r_Q> = (1 / sqrt(N)) * sum_k eta_k |k, Q-k>

where the phases are Rademacher signs eta_k = +/-1. The signs are generated as a product over the bits of k, so the state still has the same compact finite-state structure as mpsexcitonQ: it enforces k + h = Q (mod N) with a carry bond of dimension at most 2, while adding random signs to the electron-bit tensors.

This is the MPS probe used for stochastic continuum traces. For an operator restricted to the fixed-Q sector,

E[<r_Q|A|r_Q>] = (1 / N) * sum_k <k,Q-k|A|k,Q-k>

because the random signs satisfy E[eta_k eta_l] = delta_kl. Building one probe is O(LPhys) and avoids launching one KPM recursion per explicit relative momentum k.

NonHermitianHamiltonian

Wrapper for a non-Hermitian TBHamiltonian together with its hermitized block Hamiltonian.

Fields

• parent : original TBHamiltonian, not modified in-place
• z : complex reference point used in zI - H
• block_s : dim-2 auxiliary Index tagged "NHBlock"
• hermitized : Hermitian TBHamiltonian on [parent.sites...; block_s]

The hermitized Hamiltonian can be passed to existing MPO/KPM routines. Avoid calling tight-binding mutation helpers like add_onsite! on hermitized; mutate parent first and call hermitize again.

_nh_diag_online(NH, n; scale, maxdim, cutoff) -> (A_mps, dos)

Online NH KPM diagonal spectral function: run the partial Chebyshev recursion and accumulate the site-resolved diagonal MPS A(r, z) in a single pass, keeping only two partial MPOs in memory at a time.

Compared with nh_kpm_partials + nh_reconstruct_spectral_mps:

• Memory: O(2 χP²) instead of O(N χP²).
• Diagonal extractions: N/2 (only even Julia-index partials contribute).

_nh_kpm_mps_ldos(NH, n, probe_site; scale, maxdim, cutoff) -> Real

Online MPS NH KPM: compute the site-resolved spectral weight A(probe_site, z) using the dual-chain MPS partial recursion, keeping only 4 MPS in memory at a time.

probe_site is a 0-indexed physical site. The probes are localized basis states: ketprobe = |block=1⟩ ⊗ |probesite⟩ braprobe = |block=2⟩ ⊗ |probesite⟩

so inner(braprobe, pk) = ⟨2, probesite | Pk | 1, probesite⟩, which is the diagonal element of block{2,1}(Pk) at site probesite — the correct LDOS contribution at that site.

Two chains are propagated on the hermitized block space: |tk⟩ = Tk(A)|ketprobe⟩ (Chebyshev, A = Hherm / scale) |pk⟩ = Pk|ket_probe⟩ (NH partial sum)

with partial recurrence: |p0⟩ = 0, |p1⟩ = S|ketprobe⟩ |pk⟩ = 2S|t{k-1}⟩ + 2A|p{k-1}⟩ − |p_{k-2}⟩

Cost per z-point: O(Ncheb × χH × χψ) instead of O(Ncheb × χ_P²) for MPO mode.

_nh_kpm_probe_mps(sites, block_s, block_state, site_r) -> MPS

Product-state MPS on the hermitized block space. The block site carries block_state (1-indexed); position sites are set to the big-endian binary encoding of the 0-indexed physical site site_r. Bond dimension 1.

Works for both :pre (sites = [block_s; pos_sites...]) and :post (sites = [pos_sites...; block_s]) layouts — placement is auto-detected.

_nh_scalar_online(NH, n; scale, maxdim, cutoff) -> ComplexF64

Online NH KPM scalar DOS: run the partial Chebyshev recursion and accumulate Tr[block{2,1}(Pk)] contributions in a single pass, keeping only two partial MPOs in memory at a time.

Avoids the O(N·χP²) memory cost of `nhkpmpartials` and skips building the intermediate Amps entirely — each contributing step adds only one scalar to the accumulator.

_nh_stochastic_online(NH, n; scale, n_random=10, maxdim, cutoff) -> Real

Stochastic trace NH KPM DOS: estimate Tr[block{2,1}(Pk)] via Monte Carlo averaging over n_random random product-state probes on the position sites.

Each realization draws |φ⟩ = ⊗i (random local state) and runs the dual-chain MPS recursion with ket = |1⟩block ⊗ |φ⟩, bra = |2⟩block ⊗ |φ⟩. The same position state is shared between bra and ket so the estimator is unbiased: E[⟨2,φ|Pk|1,φ⟩] = Tr[P_k] / D, D = 2^L

Cost: O(nrandom × Ncheb × χH × χ_ψ) — no MPO×MPO products.

add_loss!(H, f; coefficient=-1im, space=:full, ...)

Add a loss/gain term coefficient * diag(f) to the original Hamiltonian MPO. This only modifies H.mpo; it does not create the hermitized NH block.

add_nh_imag_onsite!(H, f; prefactor=1im, space=:full, ...)

Backward-compatible name for adding prefactor * diag(f) to the original Hamiltonian MPO. This is not hermitization.

add_nh_loss!(H, f; prefactor=-1im, space=:full, ...)

Convenience wrapper for onsite loss. By default this adds -im * diag(f). Use prefactor=1im for gain or for the convention i*f(x).

add_nh_nonreciprocal_hopping!(H, t_forward, t_backward; nn=1, ...)

Add asymmetric hopping directly to H. This is the skin-effect helper: choose, for example, t_forward=t*exp(g) and t_backward=t*exp(-g).

add_nh_onsite!(H, v; Lx=nothing, tol=1e-8, maxdim=200, type=ComplexF64)

Add a possibly complex onsite potential to a position-only TBHamiltonian.

v may be:

• a number, e.g. 1im * gamma
• v(n) with n = 0, ..., H.N-1
• v(ix, iy) with 0-indexed coordinates, requiring Lx=...

This is the non-Hermitian counterpart of add_onsite!; unlike the generic version it preserves complex values by default.

add_nh_skin_hopping!(H, t, g; nn=1, convention=:exp)

Convenience wrapper for non-reciprocal skin hopping.

• convention=:exp uses t_R = t * exp(g), t_L = t * exp(-g).
• convention=:linear uses t_R = t + g, t_L = t - g.

contract_nh_block(W, block_s; row=2, col=1) -> MPO

Extract the (row, col) block of an MPO whose first or last site is block_s, returning an MPO on the remaining sites. The block position is auto-detected.

hermitize(NH; z=NH.z, cutoff=1e-8, maxdim=200, scale=0.0)

Rebuild the hermitized block Hamiltonian from NH.parent, optionally at a new reference point z.

hermitize(H; z=0, cutoff=1e-8, maxdim=200, scale=0.0)
    -> NonHermitianHamiltonian

Build a NonHermitianHamiltonian wrapper without modifying H.

hermitized_hamiltonian(H; z=0, block_s=nh_block_index(), cutoff=1e-8,
                       maxdim=200, scale=0.0, block_placement=:post) -> TBHamiltonian

Return the Hermitian block Hamiltonian

[ 0        zI - H ;
  (zI-H)'  0      ]

as a TBHamiltonian. block_placement controls where the auxiliary block site lives:

• :post (default) — site order [H.sites...; block_s]; position qubits occupy 1:L directly
• :pre — site order [block_s; H.sites...]; original layout before postpend change

scale=0.0 keeps the usual lazy KPM spectral-bound estimation.

loss_profile_mpo(H, f; Lx=nothing, type=Float64, space=:full)

Build the real diagonal profile MPO diag(f) used for loss/gain terms. This function does not multiply by im and does not hermitize anything.

By default space=:full, so the diagonal is built on the full Hamiltonian site space H.sites, with basis coordinate n = 0, ..., 2^length(H.sites)-1. Use space=:position only when the profile should live on position qubits before any auxiliary spaces are attached.

nh_block_index() -> Index

Create the dim-2 auxiliary block index used for non-Hermitian hermitization. State 1 is the upper block, state 2 is the lower block. The index keeps a "Qubit" tag so generic MPO(sites, "Id") construction remains compatible with the existing KPM code paths.

nh_block_source(NH; row=2, col=1) -> MPO
nh_block_source(Hh, block_s; row=2, col=1) -> MPO

Build the off-diagonal block source |row><col| x I. The default row=2, col=1 matches the old I_ldn source used in the non-Hermitian KPM recursion.

nh_imag_onsite_mpo(H, f; prefactor=1im, Lx=nothing, type=Float64,
                   space=:full)

Build prefactor * diag(f). The real profile MPO is constructed first by loss_profile_mpo, then the imaginary prefactor is multiplied afterward. This returns a term on the original Hilbert space, not an NH hermitized block.

nh_kpm_partials(Hh, n; source, scale, maxdim=100, cutoff=1e-8)
    -> Vector{MPO}

Compute the auxiliary "partial" Chebyshev recursion used by the old non-Hermitian spectral algorithm. If A = Hh / scale and S is the block source, the recurrence is

P_0 = 0
P_1 = S
P_k = 2 S T_{k-1}(A) + 2 A P_{k-1} - P_{k-2}

while T_k(A) is advanced in parallel. The returned vector has length 2n and stores P_0, P_1, ..., P_{2n-1}.

nh_kpm_scale(H, z_points; scale=nothing, padding=1.05, maxdim=200,
             cutoff=1e-8, dmrg_nsweeps=5, dmrg_maxdim=[10,20,40],
             dmrg_linkdim=4, printinfo=false)

Return one universal, zero-centered Chebyshev scale for NH KPM over all complex points in z_points.

For scale=nothing or scale=0, the parent operator norm is estimated from the Hermitian zero-shift block Hamiltonian

[ 0  -H ; -H'  0 ]

whose spectral radius is ||H||_2. The scale then uses the triangle bound

||zI - H||_2 <= |z| + ||H||_2

and returns padding * (maximum(abs(z_points)) + ||H||_2). This is deliberately conservative and uses the same scale for every point, keeping values comparable across a grid. Passing a positive numeric scale bypasses the estimator.

nh_nonreciprocal_hopping_mpo(H, t_forward, t_backward; nn=1, ...)

Build the position-space MPO

t_forward  * K_+^nn + t_backward * K_-^nn

without imposing Hermiticity. t_forward and t_backward may be numbers or site-dependent one-argument functions t(n) with 0-indexed n.

nh_preprocess_partials(partials, block_s; row=2, col=1) -> Vector{MPS}

For each partial MPO, extract the requested NH block and then extract its diagonal as an MPS. This is the old pre_process step without global state.

nh_reconstruct_spectral_mpo(partials, n, NH; maxdim=400, cutoff=1e-8)
    -> (ldos_mps, dos, rotated_mpo)

High-bond-dimension reconstruction of the non-Hermitian spectral object.

This mirrors the older all-site LDOS workflow: first reconstruct the full partial MPO on the hermitized block space, then left-multiply by |1><2| x I and take its trace / diagonal. Compared with nh_reconstruct_spectral_mps, this keeps the full MPO until the end, so it can represent all sites at once but usually needs a much larger maxdim.

nh_reconstruct_spectral_mps(partials, n, block_s; maxdim=100,
                            row=2, col=1) -> (A_mps, dos)

Apply the old Jackson reconstruction to the even partial terms:

A = 2/(pi^2 (2n+1)) * sum_l (-1)^(l/2-1) g_l diag(P_l)

where l = 2, 4, ..., 2n in one-based Julia indexing of the partial list. Returns the diagonal spectral MPS and its summed value.

nh_spectral_function(NH, n; scale, maxdim=100, cutoff=1e-8,
                     source_row=2, source_col=1, block_row=2, block_col=1)
    -> (A_mps, dos, partials)

Convenience wrapper for the full non-Hermitian KPM spectral calculation at the reference point stored in NH.z.

nh_spectral_function_allsite_mpo(NH, n; scale, maxdim=400, cutoff=1e-8)
    -> (ldos_mps, dos, rotated_mpo, partials)

Alternative non-default NH KPM path that reconstructs the full spectral MPO before extracting the LDOS. It can compute all sites in one object, but needs larger bond dimensions for accuracy.

nh_spectrum_grid(H, xlims, nx, ylims, ny, n; scale=nothing,
                 nh_scale_padding=1.05, convention=:z_minus_H,
                 mode=:scalar, probe_site=0, n_random=10,
                 maxdim=100, cutoff=1e-8, verbose=false)

Evaluate the NH KPM spectral weight on a rectangular complex energy grid.

Modes

• :scalar — full NH partial MPO recursion; total DOS. O(Ncheb × χ_P²).
• :mps — dual-chain MPS at a single site (probe_site, 0-indexed). LDOS at that site. O(χH × χψ) per step.
• :diag — same as :scalar but also extracts site-resolved diagonal MPS A(r,z). Extra return Z_spatial has shape (H.N, ny, nx).
• :stochastic — Monte Carlo trace: average over n_random random product-state probes. Total DOS estimate. O(nrandom × Ncheb × χH × χ_ψ). No MPO×MPO products.

Set verbose=true to print one progress line per Re(z) column. If scale is omitted or zero, a single conservative scale is estimated from nh_kpm_scale and reused over the whole grid.

_impurity_mpo(x0, L, sites, V) -> MPO

Build the on-site impurity potential V · |x₀−1⟩⟨x₀−1| as a bond-dimension-1 MPO using OpSum — exact, no QTCI required.

x0 is 1-indexed. The quantics encoding is big-endian (matching binary_to_MPS): site i carries bit (L−i) of the 0-indexed address x0−1, so site 1 holds the MSB.

Uses "projUp" (|0⟩⟨0|) and "projDn" (|1⟩⟨1|) operators defined for SiteType"Qubit".

get_qpi(H, Ncheb, ω_phys_vals;
        impurity_site=nothing, V=1.0,
        window_fraction=nothing, window_sigma=1.5,
        kernel=:jackson, lambda=4.0,
        maxdim=100, cutoff=1e-8, verbose=false)
    -> Matrix{Float64}

Compute the Quasiparticle Interference (QPI) pattern from a single on-site impurity.

Algorithm

1. Place impurity V · |x₀−1⟩⟨x₀−1| at site x0 (default: geometric center via central_index). The impurity MPO is built analytically from OpSum (no QTCI).
2. Run KPM Chebyshev expansion for the clean (H) and perturbed (H_imp) systems and extract the site-resolved diagonal LDOS as Vector{MPS} via get_ldos_spectrum.
3. For each energy ω:
   • Subtract LDOS MPS: δA_mps = ldos_imp[ω] − ldos_clean[ω].
   • If window_fraction is set, multiply δA_mps element-wise by a smooth disk mask (built via QTCI using the same SDF machinery as mask_hamiltonian). The mask equals 1 deep in the bulk and rolls off to 0 at the boundary over window_sigma lattice units, suppressing Gibbs-like ringing from open edges.
   • Apply the forward QFT MPO (position → momentum) directly to δA_mps. This works because we already have the diagonal of the spectral function as an MPS; no conjugate_by_qft (which operates on full MPOs) is needed.
   • Read out QPI(k, ω) = |⟨k|δÃ_mps⟩|² using LSB-first bit encoding (site 1 = bit 0), consistent with the TCI.reverse QFT convention.

Arguments

• H : TBHamiltonian with no auxiliary DOFs (spinless, no Nambu/layer/sublattice).
• Ncheb : Number of Chebyshev moments for KPM.
• ω_phys_vals : Physical energies to evaluate.

Keyword arguments

• impurity_site : 1-indexed site index for the impurity (default: central_index(H)).
• V : Impurity potential strength.
• impurity_mode : :delta (default) — exact rank-1 projector V·|x₀⟩⟨x₀| via OpSum; :gaussian — smooth Gaussian potential V·exp(-|r−r₀|²/2σ²) built via add_onsite! (lower bond dimension for the KPM recursion). Requires H.geometry to be set when :gaussian.
• sigma : Gaussian half-width in physical distance units (only used when impurity_mode=:gaussian). Default: 1.5.
• window_fraction : If set, build a circular (2D) or interval (1D) apodization mask with radius window_fraction × min_half_extent of the bounding box. Requires H.geometry to be set. Default: nothing (no windowing).
• window_sigma : Sigmoid roll-off half-width in lattice units. Default: 1.5.
• kernel : KPM kernel (:jackson or :lorentz).
• lambda : Lorentz kernel width.
• maxdim : Maximum bond dimension for MPS/MPO operations.
• cutoff : SVD truncation cutoff.
• verbose : Print progress.

Returns Matrix{Float64} of shape (Nω, N) where N = H.N. Column index k+1 (1-based) corresponds to 0-indexed momentum k.

Example

H     = get_Hamiltonian("square_2d", 1.0; L=8, Lx=4)
ωlist = [0.0]
qpi   = get_qpi(H, 80, ωlist; V=2.0, window_fraction=0.8)
# Reshape to 2D: reshape(qpi[1,:], 2^Lx, 2^Ly)'

_autoenable_proj(H, nambu_proj, spin_proj, layer_proj, sublat_proj)
    -> (nambu_proj, spin_proj, layer_proj, sublat_proj)

Enable projection flags for any auxiliary DOF detected on H, printing one info line per auto-enabled flag. Called at the top of every TBHamiltonian spectral method before any aux-index logic runs.

_embed_displacement_in_full_sites(H, mpo_pos) -> MPO

Like _embed_in_full_sites but pads auxiliary sites with all-ones matrices instead of identity. Required for current-operator construction: the displacement (xᵣ − xᵣ′) depends only on position, so it must be broadcast uniformly across all auxiliary (sublattice, layer, spin, Nambu) index pairs, including off-diagonal ones where physical hoppings exist.

Using identity at an aux site with off-diagonal hoppings (e.g. sublattice A↔B in honeycomb, or inter-layer tunneling in bilayers) would set those current matrix elements to zero and give σ = 0.

_embed_in_full_sites(H, mpo_pos) -> MPO

Embed mpo_pos (which lives on _pos_sites(H)) into the full H.sites space by prepending/appending dim-1 identity tensors at each auxiliary site.

_eval_diag_mps(A, x) -> Float64

Evaluate the diagonal MPS A at the 0-indexed position x using a LSB-first bit encoding (site 1 = bit 0 of x). Equivalent to inner(binary_MPS(x), A) but avoids constructing the full basis MPS.

_kpm_weight_matrix(Ncheb, ω_vals; kernel=:jackson, lambda=4.0) -> Matrix{Float64}

Precompute the full KPM weight matrix W[n, iω] for fast in-loop accumulation.

W[n, iω] = c_n · g_n · cos((n-1) · arccos(ω_iω))

• c_n = 1 for n=1, c_n = 2 otherwise (Chebyshev expansion factor)
• g_n = kernel damping: Jackson (default, finite-size ringing suppressed) or Lorentz (controlled width lambda, smoother tails)
• Entries for |ω| ≥ 1 are set to zero (outside the spectral support)

Pre-computing W avoids recomputing cos((n-1)·arccos(ω)) inside the inner loop, which is called Ncheb × Nω times.

aux_site(H, which) -> (Index, Symbol)

Return the auxiliary Index and its position side (:pre or :post) for the named auxiliary degree of freedom in H.

which ∈ :spin, :sublattice, :nambu, :layer.

Useful for passing the correct arguments to project_aux without manually inspecting H.sites.

s, side = aux_site(H, :sublattice)
W_A = project_aux(W, s, 1; side=side)   # sublattice-A channel

conjugate_by_qft(W; tol=1e-9, maxdim=100) -> MPO

Return U · W · U† where U is the Quantum Fourier Transform MPO built from QuanticsTCI.quanticsfouriermpo (normalised, with TCI.reverse applied).

TCI.reverse places the LSB at site 1, matching the quantics encoding used throughout this codebase. The resulting k-space MPO has the same site structure as W but with momenta as the diagonal degree of freedom.

Calling this on the Chebyshev spectral MPO T_n and then extracting the diagonal gives the k-resolved contribution to A(k,ω).

conjugate_by_qft(H::TBHamiltonian, W::MPO; tol, maxdim) -> MPO

TBHamiltonian-aware version of conjugate_by_qft. Applies U·W·U† where U is the QFT acting only on the position (Qubit) sites of H, with identity operators at all auxiliary sites (Layer, spin, sublattice, Nambu).

Use this overload whenever W lives on the full H.sites space (including aux indices), as is the case in the bubble pipeline after replace_sites.

conjugate_by_qft_exciton(H::TBHamiltonian, W; tol=1e-9, maxdim=100) -> MPO

Two-particle (electron–hole) analogue of conjugate_by_qft for the interleaved 2L-site exciton encoding produced by Exciton_Hamiltonian / exciton_hamiltonian, where one carrier occupies the odd sites (1, 3, …) and the other the even sites (2, 4, …).

Returns U · W · U† where U = U_e ⊗ U_h is the product of two independent Quantum Fourier Transforms, one acting on each carrier's R = H.L position qubits. Conjugating the exciton spectral MPO T_n by U rotates both carriers into momentum space; extracting the diagonal then resolves the joint (kₑ, k_h) weight (apply a relative-/centre-of-mass-momentum shift first to read a fixed total momentum q).

The site list comes from _pos_sites(H) — the same canonical accessor the single-particle conjugate_by_qft(H, W) overload uses — rather than being reverse-engineered from siteinds(W).

The construction mirrors the single-particle conjugate_by_qft but builds the full two-particle transform first:

1. Two single-particle QFT cores from QuanticsTCI.quanticsfouriermpo (TCI.reverse, normalize=true) — sign=-1 and sign=+1, identical to conjugate_by_qft so the per-carrier convention matches.
2. The same core is interleaved onto both carrier registers — the even sites (interleave_mpo(…, 0)) and the odd sites (interleave_mpo(…, 1)) — and the two halves are multiplied into the full transforms U and U† (= U_e ⊗ U_h). The product contracts away the interleaved identity tensors, leaving dense cores.
3. fix_sites re-canonicalises the combined MPOs onto the physical sites, then swapprime undoes the bra↔ket swap interleave_mpo introduces (needed because the reversed QFT is not symmetric, so the swap is a genuine transpose).
4. The conjugation U · W · U† is applied exactly as in the single-particle routine (swapprime left-multiply); the result stays on sites.

Verified to machine precision against the dense U·W·U† (with U taken from the single-particle conjugate_by_qft) for R = 2, 3, 4. As with the 2L-site MPOs involved, raise maxdim (and/or tighten tol) for larger R; the default maxdim=100 starts truncating around R ≳ 4.

See also conjugate_by_qft, Exciton_Hamiltonian, interleave_mpo, fix_sites.

get_bands(H_mpo, scale, center, sites, Ncheb, D, ω_vals; kwargs...) -> Matrix{Float64}

Memory-efficient band structure via online Chebyshev KPM accumulation. See the section 4 block comment above for the full four-step projection pipeline.

Arguments

• H_mpo : unscaled Hamiltonian MPO on all sites (position + any aux).
• scale, center: energy rescaling so that H̃ = (H−center)/scale ∈ (−1, 1).
• sites : the full site list of H_mpo including any aux indices. Position-only site count is inferred as L_pos = L − n_aux.
• Ncheb : number of Chebyshev moments.
• D : spatial dimension (1 or 2).
• ω_vals : rescaled energies ∈ (−1, 1) at which to evaluate A(k,ω).

Projection keyword arguments

Each projection flag is independent; any combination is valid.

Nambu (BdG particle/hole) projection — Step 0:

• nambu_proj : project each T_n onto Nambu sectors (default false).
• proj_nambu : 1 = particle only, 2 = hole only, nothing = sum both.
• nambu_s : the Nambu Index (auto-detected from H.nambu_s via the TBHamiltonian overload).
• nambu_side : :pre (default) or :post — position of the Nambu site.

Spin projection — Step 1:

• spin_proj : project each T_n onto spin channels (default false).
• proj_s : 1 = ↑ only, 2 = ↓ only, nothing = sum both.
• spin_s_aux : explicit spin Index; when nothing, falls back to sites[1]. Set automatically by the TBHamiltonian overload so that spin is correctly identified even when Nambu is also prepended at site 1.

Layer projection — Step 1c (bilayer / multilayer):

• layer_proj : project each T_n onto individual layers (default false).
• proj_layer : k = layer k only, nothing = sum all layers.
• layer_s : the layer Index (auto-detected from H.layer_s).
• layer_side : :pre (default) — layer is always prepended.

Sublattice auxiliary projection — Step 1b (kagome, Lieb, honeycomb):

• sublat_proj : project each T_n onto sublattice aux sectors (default false).
• proj_sl : k = sublattice k only, nothing = sum all. Shared with the legacy sublattice flag (Step 2).
• sublat_s : the sublattice Index (auto-detected from H.sublattice_s).
• sublat_side : :post (default) or :pre — position of the sublattice site.

Legacy sublattice mask projection — Step 2 (2-sublattice models without aux index):

• sublattice : apply a mask sandwich mask · T_n · mask (default false).
• proj_sl : 1 = mask A only, 2 = mask B only, nothing = both.

k-point sampling keyword arguments

• xmin, xmax, num_x : grid in x (1D) or kx (2D). Default: full range, 10 pts.
• ymin, ymax, num_y : grid in ky (2D only).
• num_avg : number of offset points averaged around each center (default 1).

Truncation and performance

• kernel : KPM broadening kernel (:jackson or :lorentz).
• lambda : Lorentz kernel width (ignored for Jackson).
• tol, maxdim, cutoff : MPO truncation parameters passed to apply and truncate!.
• printinfo : print maxlinkdim every 10 Chebyshev steps (default false).

Returns

Matrix{Float64} of shape (Nω, num_x).

get_bands(H, Ncheb, ω_phys_vals; kwargs...)

Convenience overload that infers the spatial dimension D from H.geometry (via length(H.geometry(1))), so callers do not need to pass D explicitly. All keyword arguments are forwarded unchanged to the 4-argument form.

Errors if H.geometry is nothing (custom or geometry-free Hamiltonians must still pass D explicitly via the 4-argument form).

get_bands(H, Ncheb, D, ω_phys_vals; kwargs...)
    -> Matrix{Float64}  or  NamedTuple(Ak, ticks, labels)

High-level overload of get_bands for a TBHamiltonian.

Physical energies ω_phys_vals are rescaled via H.scale and H.center.

Auto-detection: All auxiliary site Indices (Nambu, spin, layer, sublattice) and their positions (:pre/:post) are read from the struct fields and excluded from position k-space automatically — no manual index passing required.

Projection kwargs (forwarded verbatim to the low-level MPO method): spin_proj, proj_s, nambu_proj, proj_nambu, layer_proj, proj_layer, sublat_proj, proj_sl, sublattice.

High-symmetry k-path shortcut — replaces the manual hsk_* + kpath_2d

• k_groups_override boilerplate with a single call:

res = get_bands(H, Ncheb, 2, omega;
                kpath=[:G, :M, :Kp, :G], kpath_lattice=:honeycomb, num_x=30)

• kpath : symbol vector defining the path. Use Latin aliases: G=Γ, M, K, Kp=K', X — no special characters needed.
• kpath_lattice : :honeycomb, :square, or :triangular.
• kpath_Lx : Lx for the 2D grid; defaults to H.L ÷ 2.
• num_x : reused as npts_per_segment when kpath is given.

When kpath is set the return value is a NamedTuple: (Ak = Matrix{Float64}(Nω×Nk), ticks = Vector{Int}, labels = Vector{String})

julia heatmap(1:size(res.Ak,2), omega, res.Ak; xticks=(res.ticks, res.labels)) vline!(p, res.ticks; ls=:dash, color=:white)

Otherwise returns Matrix{Float64} as usual (backward-compatible).

get_exciton_bands(H_QFT, H, Ncheb, omega_phys_vals; Q_list, Q_groups,
                  num_q, num_avg, q_start, q_end, kernel,
                  lambda, eta, m_order, maxdim, cutoff,
                  verbose, printinfo) -> Matrix{Float64}

CPU momentum-space exciton bands from an already-QFT-conjugated exciton MPO. This is the MPS-KPM analogue of get_exciton_ldos_spatial: for each total momentum label Q it runs an online MPS Chebyshev recursion from mpsexcitonQ(Q, H.sites) = (1 / sqrt(N)) * sum_k |k, -k + Q> and accumulates all requested energies in one pass. The exciton must not use the MPO-KPM get_bands pipeline.

H_QFT supplies the MPO used in the recurrence. The original real-space H::TBHamiltonian supplies metadata (sites, N, L, scale, center) and sets the basis convention for mpsexcitonQ.

Rows are energies, columns are total momenta/groups. Q_list selects momenta directly. Q_groups (or alias q_groups) averages several momentum probes into one output column. K_list, K_groups, k_groups, num_k, k_start, and k_end are accepted as backward-compatible aliases.

get_exciton_continuum(H_QFT, H, Ncheb, omega_phys_vals;
                      Q_list, num_q, q_start, q_end,
                      N_sample, k_list, seed, normalize,
                      kernel, lambda, eta, m_order, maxdim, cutoff,
                      verbose, printinfo) -> Matrix{Float64}

Stochastic MPS-KPM trace over the electron-hole continuum at fixed total momentum. For each selected total momentum Q, this estimates

A_cont(Q, omega) = (1 / N) * sum_k <k, Q-k| delta(omega - H_QFT) |k, Q-k>

with random-phase trace probes built by mpsexcitonQTrace(Q, H.sites). This is an incoherent trace over relative momentum, unlike get_exciton_bands, which probes the coherent pair state mpsexcitonQ(Q, H.sites).

H_QFT supplies the already-QFT-conjugated MPO used in the MPS recursion. The original H::TBHamiltonian supplies metadata (sites, N, L, scale, center) and the exciton site convention. The exciton continuum remains an MPS-KPM calculation: no exciton MPO-KPM / QFT_tk.get_bands path is used.

Rows are energies, columns are total momenta. If k_list is provided, those 1-indexed relative momenta are used deterministically for every Q; otherwise N_sample compact random-phase trace probes are drawn for each Q. Each trace probe is a randomized superposition over all |k,Q-k> states, so this avoids a full KPM recursion for every explicit relative momentum.

With normalize=true (default), the result estimates the per-relative-momentum average (1/N) * Tr_Q. With normalize=false, it is multiplied by H.N and estimates the total fixed-Q trace Tr_Q.

hsk_honeycomb(Lx, Ly) -> NamedTuple

High-symmetry k-points for the honeycomb sublattice Hamiltonian (honeycomb_sublattice_hamiltonian) in quantics integer units (kxidx ∈ [0, 2^Lx−1], kyidx ∈ [0, 2^Ly−1]).

Derivation. The Bloch off-diagonal element is h(k) = t ( 1 + e^{2πi kx/Nx} + e^{2πi ky/Ny} ) Dirac points h=0 require θx = 2π/3 AND θy = 4π/3 (or their conjugates): K : (kxidx, kyidx) = (Nx/3, 2Ny/3) → Cartesian (2π/3, 2π/√3) K' : (kxidx, kyidx) = (2Nx/3, Ny/3) → Cartesian (4π/3, 0)

Because Nx = 2^Lx is never divisible by 3, the K/K' indices are rounded to the nearest integer. Use a large Lx (≥4) for a good approximation.

M is the edge midpoint adjacent to K' along the kx axis: M : (Nx/2, Ny/4) → Cartesian (π, 0) — |h|=1, saddle point

Standard G–M–Kp–G path (along the kx direction, K' corner at (4π/3,0)):

hs = hsk_honeycomb(Lx, Ly)
kg, ticks = kpath_2d([hs.G, hs.M, hs.Kp, hs.G], Lx; npts_per_segment=30)
Ak = get_bands(H, Ncheb, 2, omega; k_groups_override=kg)
# or using the kpath shortcut:
res = get_bands(H, Ncheb, 2, omega; kpath=[:G, :M, :Kp, :G],
                kpath_lattice=:honeycomb, num_x=30)

hsk_square(Lx, Ly) -> NamedTuple

High-symmetry k-points for a 2D square lattice in quantics integer units.

Standard path: [:G, :X, :M, :G]

hsk_triangular(Lx, Ly) -> NamedTuple

Approximate high-symmetry k-points for a 2D triangular lattice in quantics integer units.

Standard path: [:G, :M, :K, :G]

ilinspace(xmin, xmax, num_x) -> Vector{Int}

Return num_x as almost evenly spaced integers in [xmin, xmax], inclusive, with a preference for the endpoints. Used to build the k-point center grid for band-structure sampling.

kpath_2d(hs_points, Lx; npts_per_segment=20) -> (k_groups, tick_positions)

Build k-groups for get_bands(…; D=2, k_groups_override=k_groups) sampling along a high-symmetry path.

hs_points is an ordered list of (kx_idx, ky_idx) integer tuples defining the path vertices (0-based, kxidx ∈ [0, 2^Lx−1]). `nptsper_segment` points are linearly interpolated between each consecutive pair of vertices.

Returns:

• k_groups : Vector{Vector{Int}} — single-element groups in the (ky << Lx) | kx linear-index format expected by get_bands. Pass as k_groups_override.
• tick_positions : 1-based indices into k_groups at each vertex; use as x-axis tick positions when plotting.

# Low-level: explicit tuples
hs = hsk_honeycomb(Lx, Ly)
kg, ticks = kpath_2d([hs.G, hs.M, hs.Kp, hs.G], Lx; npts_per_segment=30)
Ak = get_bands(H, Ncheb, 2, omega; k_groups_override=kg)

# High-level shortcut (preferred):
res = get_bands(H, Ncheb, 2, omega;
                kpath=[:G, :M, :Kp, :G], kpath_lattice=:honeycomb, num_x=30)
# heatmap(1:size(res.Ak,2), omega, res.Ak; xticks=(res.ticks, res.labels))

kpath_setup(lattice, Lx, Ly, path_syms; npts_per_segment=20)
    -> (k_groups, ticks, labels)

Build k-path inputs for get_bands from a list of high-symmetry symbols. path_syms is a vector of symbols such as [:G, :M, :Kp, :G] (use G for Γ — the Latin alias avoids non-ASCII input). npts_per_segment is the number of points between each consecutive pair.

Returns (k_groups, ticks, labels) ready to pass to get_bands as k_groups_override, and to heatmap as xticks=(ticks, labels).

project_aux(W, aux_s, σ; side=:pre) -> MPO

Remove an auxiliary site from MPO W by projecting onto state σ.

• side=:pre — aux site is at position 1 (prepended, e.g. spin).
• side=:post — aux site is at the last position (postpended, e.g. sublattice).

Contracts the projector |σ⟩⟨σ| on both bra and ket physical indices of the aux tensor; the resulting dim-1 link is absorbed into the adjacent position site. Returns an (L−1)-site MPO suitable for conjugate_by_qft.

project_spin(W, spin_s, σ) -> MPO

Convenience wrapper for project_aux when the auxiliary site is prepended (spin at site 1). Equivalent to project_aux(W, spin_s, σ; side=:pre).

projop_1DSL(O, sites, Lx, SL) -> MPO

Project MPO O (on a 2^Lx chain) onto sublattice SL (1 = even sites, 2 = odd sites) by sandwiching with the corresponding diagonal mask: mask · O · mask.

projop_2DSL(O, sites, Lx, Ly, SL) -> MPO

Project MPO O (on a 2^Lx × 2^Ly lattice) onto sublattice SL (1 = even checkerboard, 2 = odd checkerboard) by sandwiching with the corresponding diagonal mask: mask · O · mask.

sample_diag(Tn_k, ikstart, ikend) -> Vector{Float32}

Extract the diagonal of MPO Tn_k as an MPS and evaluate it at every integer index in ikstart:ikend. Convenience wrapper around _eval_diag_mps; used when all k-points in a contiguous range are needed.

bdg_hamiltonian(H_kin, H_pair, nambu_s; cutoff=1e-8) -> MPO

Build a spinless Bogoliubov–de Gennes Hamiltonian on [nambu_s; pos_sites…]:

H_BdG = τ_z ⊗ H_kin  +  τ_+ ⊗ H_pair  +  τ_- ⊗ dag(H_pair)

H_kin is the single-particle kinetic/hopping MPO measured from the chemical potential (H_kin = H_tb − μ·I). H_pair encodes the pairing amplitude Δ(i,j).

Note: for spinless fermions c_i† c_j† = −c_j† c_i†, so the pairing matrix must satisfy Δ(i,j) = −Δ(j,i). On-site (s-wave) pairing is therefore forbidden; the minimal allowed symmetry is p-wave (nearest-neighbour, antisymmetric). For the Kitaev chain with uniform p-wave amplitude Δ:

H_pair = hopping2MPO(L, [(i, i+1, Δ) for i in 0:N-2], sites)   # Δ(i,i+1) = +Δ

The h.c. term τ_- ⊗ H_pair† (which carries Δ(j,i) = −Δ) is built automatically.

The result is Hermitian for any H_kin = H_kin† and any complex H_pair.

bdg_spin_hamiltonian(H_kin_up, H_kin_down, H_pair, spin_s, nambu_s;
                     H_soc=nothing, cutoff=1e-8) -> MPO

Build a spin-½ singlet BdG Hamiltonian on [nambu_s, spin_s; pos_sites…].

Nambu–spin convention: Ψ = (c_↑, c_↓, c†_↓, −c†_↑)ᵀ (standard BCS).

H_BdG = τ_z⊗P_↑ ⊗ H_kin_up  +  τ_z⊗P_↓ ⊗ H_kin_down
      + τ_+⊗(i·σ_y) ⊗ H_pair  +  τ_-⊗(−i·σ_y) ⊗ H_pair†
      [+ τ_z⊗S_z ⊗ H_soc]

The first two lines are the kinetic energy (allowing spin-dependent fields, e.g. Zeeman: pass H_kin_up = H_tb − (μ+h)·I, H_kin_down = H_tb − (μ−h)·I).

The pairing lines implement singlet Cooper-pair creation via the antisymmetric spin factor i·σ_y = [[0,1],[−1,0]]. For on-site s-wave pairing, build H_pair as a diagonal MPO with Δ on the diagonal.

H_soc (optional) adds an Ising-type spin-orbit coupling τ_z⊗S_z⊗H_soc.

The result is Hermitian for real or complex H_pair and any H_kin_up/dn.

nambu_index() -> Index

Create a dim-2 Index tagged "Nambu" (state 1 = particle, state 2 = hole). Pass the result as nambu_s to all prepend_nambu calls.

pairing2MPO(f, N, sites; tol=1e-8, type=ComplexF64,
            initial_positions=[], unfoldingscheme=:interleaved) -> MPO

Antisymmetric analogue of hopping2MPO for superconducting pairing.

Compresses a general N×N pairing matrix Δ[i,j] = f(i,j) into an MPO via Quantics TCI. f is passed directly to TCI — the caller is responsible for providing an antisymmetric function (f(i,j) = −f(j,i)).

For nearest-neighbour pairing prefer pairingNNN — it is exact and avoids TCI. Use pairing2MPO for longer-range or d-wave / p±ip patterns.

H_pair is passed directly to bdg_hamiltonian or bdg_spin_hamiltonian.

pairingNNN(L, sites, hopping, nn; apply_kwargs=NamedTuple()) -> MPO

Antisymmetric analogue of kineticNNN for superconducting pairing.

Returns

H_pair = hopping · Kf^nn  −  Kb^nn · dag(hopping)

where Kf / Kb are the quantics forward / backward unit-shift operators. The sign flip relative to kineticNNN (+ → −) produces an antisymmetric matrix Δ(i,j) = −Δ(j,i), satisfying the Fermi constraint c_i†c_j† = −c_j†c_i†.

For uniform pairing hopping = Δ·I:

H_pair = Δ · (Kf^nn − Kb^nn)

H_pair is passed directly to bdg_hamiltonian or bdg_spin_hamiltonian; the τ_- ⊗ H_pair† (hole-particle) term is built automatically there.

postpend_nambu(H_mpo, nambu_s, op) -> MPO

Append a Nambu operator on nambu_s to the end of H_mpo. op is a Symbol (same table as prepend_nambu) or an explicit 2×2 matrix. Equivalent to postpend_op(H_mpo, nambu_s, op).

postpend_spin(H_mpo, spin_s, op) -> MPO

Append a spin-½ operator on spin_s to the end of H_mpo. op is a Symbol (same table as prepend_spin) or an explicit 2×2 matrix. Equivalent to postpend_op(H_mpo, spin_s, op).

prepend_nambu(H_mpo, nambu_s, op) -> MPO

Prepend a Nambu (particle–hole) operator on nambu_s (created with nambu_index()). op is a Symbol from the table below or an explicit 2×2 matrix. Equivalent to prepend_op(H_mpo, nambu_s, op).

Basis: state 1 = particle, state 2 = hole.

prepend_op(H_mpo, s, op::Symbol) -> MPO
postpend_op(H_mpo, s, op::Symbol) -> MPO

Symbol dispatch for indices tagged "Spin" or "Nambu". Looks up op in the appropriate operator dictionary and calls the matrix-form prepend_op / postpend_op.

Spin ops (index tagged "Spin"): :Id, :Pup, :Pdn, :Sz, :Sp, :Sm, :Sx, :Sy, :iSy, :miSy

Nambu ops (index tagged "Nambu"): :Id, :Pp, :Ph, :tz, :tx, :ty, :tp, :tm

prepend_spin(H_mpo, spin_s, op) -> MPO

Prepend a spin-½ operator on spin_s (created with spin_index()). op is a Symbol from the table below or an explicit 2×2 matrix. Equivalent to prepend_op(H_mpo, spin_s, op).

Basis: state 1 = ↑, state 2 = ↓.

spin_hamiltonian(H_up, H_down, spin_s;
                 H_Zeeman=nothing, cutoff=1e-8) -> MPO

Build a spin-resolved Hamiltonian on [spin_s; pos_sites…]:

H = P_↑ ⊗ H_up  +  P_↓ ⊗ H_down  [+  S_z ⊗ H_Zeeman]

H_up, H_down are MPOs on the same position sites (they may differ for spin-orbit coupling or magnetic exchange). H_Zeeman is an optional position-MPO encoding a local magnetic field h(x); it enters as S_z ⊗ H_Zeeman so spin-↑ gains +½ h(x) and spin-↓ gains −½ h(x).

spin_index() -> Index

Create a dim-2 Index tagged "Spin" (state 1 = ↑, state 2 = ↓). Pass the result as spin_s to all prepend_spin calls.