Lattice

H2DChernhex(Lx, Ly, t, t2, ms; uniformhaldane=false, uniformsemenoff=false, ...) -> MPO

Haldane-type Chern insulator on a hexagonal lattice.

• NN hopping t (intra-row hex + vertical inter-row)
• Complex NNN hopping +/-i*T2(x,y) (checkerboard alternation) for next-nearest
• Semenoff mass Ms(x,y) on-site term
• Domain wall in t2 and Ms along x = Nx/2 by default (uniformhaldane=true and uniformsemenoff=true override to uniform fields)

HAAH(L, V, phi, t; b=(1+-)/2, tol_quantics=1e-8, maxbonddim_quantics=50) -> MPO

Aubry-Andre-Harper quasicrystal: H = t * sum c{i+1}ci + V * cos(2pi b i + phi) * n_i

HChern8(Lx, Ly, V, t; a=5/64*2^Lx, t2=0.2t, tol_quantics=1e-8, ...) -> MPO

8-fold "Chern mosaic" Hamiltonian: uniform intra/inter-row hoppings modulated by a spatially varying 8-fold pattern using 4 rotated k-vectors.

HQC2Dsquare(Lx, Ly, t; tol_quantics=1e-8, maxbonddim_quantics=100, cutoff=1e-10) -> MPO

Quasicrystal-modulated square lattice. The hopping amplitude at each bond is evaluated at the bond midpoint using an 8-fold modulation with two competing wavevectors b1 = 5- a/2 and b2 = - Nx a/16.

HSSH(L, t, d; tol_quantics=1e-8, maxbonddim_quantics=10, nn=1) -> MPO

SSH (Su-Schrieffer-Heeger) Hamiltonian: dimerized hopping t+/-d on alternating bonds.

HUniform(L, t; v=1e-6, tol_quantics=1e-8, maxbonddim_quantics=10, nn=1) -> MPO

Uniform-hopping tight-binding chain on 2^L sites with an optional uniform onsite potential v. A small nonzero v is required to avoid TCI failure on constant functions.

HUniform2Dhex(Lx, Ly, t; tol_quantics=1e-8, maxbonddim_quantics=10, cutoff=1e-10) -> MPO

Uniform tight-binding Hamiltonian on a 2^Lx x 2^Ly hexagonal lattice. Intra-row uses kineticintra2DNNhex (checkerboard mask); inter-row uses kineticNNN(- Nx).

HUniform2Dsquare(Lx, Ly, t; tol_quantics=1e-8, maxbonddim_quantics=10, cutoff=1e-10) -> MPO

Uniform tight-binding Hamiltonian on a 2^Lx x 2^Ly square lattice (row-major encoding). Intra-row: kineticintra2DNNN(- nn=1). Inter-row: kineticNNN(- nn=Nx).

HUniform2Dtri(Lx, Ly, t; tol_quantics=1e-8, maxbonddim_quantics=10, cutoff=1e-10) -> MPO

Uniform tight-binding Hamiltonian on a 2^Lx x 2^Ly triangular lattice. Three kinetic terms:

• kineticintra2DNNN(- 1) -intra-row NN
• kineticinterNNNtriSWNE(- Nx+1) -SW->NE diagonal
• kineticinterNNNtriSENW(- Nx-1) -SE->NW diagonal

HUniform2Dtri_bravais(Lx, Ly, t; tol_quantics=1e-8, maxbonddim_quantics=10, cutoff=1e-10) -> MPO

Uniform tight-binding Hamiltonian on a 2^Lx x 2^Ly triangular lattice with proper Bravais vectors a1=(1,0), a2=(1/2,-/2). Exactly three bond types per unit cell:

• (dix=+1, diy= 0): intra-row x via kineticintra2DNNN
• (dix= 0, diy=+1): y-hop via kineticNNN(- Nx)
• (dix=+1, diy=-1): Bravais diag via kineticinterNNNtri_bravais_diag

_col_select_mpo(Lx, Ly, sites; keep=:even) -> MPO

Diagonal mask that retains only even or odd columns of a 2^Lx x 2^Ly grid.

• keep = :even ->1 where ix % 2 == 1 (0-based; LSB of ix = 1)
• keep = :odd ->1 where ix % 2 == 0 (0-based; LSB of ix = 0)

_parse_param_string(s) -> Dict{Symbol,Any}

Parse "key1=val1, key2=val2, - into a Dict. Values are auto-typed as Bool, Int, Float64, or String.

_row_break_mpo(Lx, Ly, sites; which) -> MPO

Diagonal mask that zeroes wrap-around couplings at row boundaries of a 2^Lx x 2^Ly grid (row-major encoding).

• which = :xplus ->0 where ix == 2^Lx -1 (end of each row)
• which = :xplain ->0 where ix == 0 (start of each row)

Multiply a kinetic MPO by this mask on the appropriate side to suppress the bond that crosses a row boundary.

_row_checker_mpo(Lx, Ly, sites) -> MPO

Diagonal checkerboard mask: 1 where (ix + iy) is even, 0 otherwise. Equivalent to projecting onto LSB(ix) == LSB(iy), i.e. both qubits agree: proj{iy-LSB=0, ix-LSB=0} + proj{iy-LSB=1, ix-LSB=1}

_row_select_mpo(_, Ly, sites; keep=:even) -> MPO

Diagonal mask that retains only even or odd rows of a 2^Lx x 2^Ly grid.

• keep = :even ->1 where iy % 2 == 1 (0-based; LSB of iy = 1)
• keep = :odd ->1 where iy % 2 == 0 (0-based; LSB of iy = 0)

build_hamiltonian(model, L; mparams="", mparam_dict=Dict()) -> MPO          (1D)
build_hamiltonian(model, Lx, Ly; mparams="", mparam_dict=Dict()) -> MPO     (2D)

Build a Hamiltonian MPO by model name using the MODEL_REGISTRY.

H = build_hamiltonian("aah",  8; mparams="V=2.0, phi=0.0, t=1.0")
H = build_hamiltonian("square_2d", 4, 4; mparams="t=1.0")
H = build_hamiltonian("chernhex",  4, 4; mparam_dict=Dict(:t=>1.0, :t2=>0.3, :ms=>0.0))

Known models: aah, chern8, chernhex, hex2d, qc2dsquare, square2d, ssh, triangular2d, triangularbravais, uniform

dice_hamiltonian(Lx, Ly[, t]; t_AB, t_AC, cutoff, maxdim) -> TBHamiltonian

Build a dice (T3) tight-binding Hamiltonian as a TBHamiltonian.

Encoding (L+1 sites, L = Lx+Ly):

• Sites 1..L : position qubits for 2^L unit cells on a triangular Bravais lattice
• Site L+1 : dim-3 "Dice" sublattice index A=1 (hub), B=2 (rim), C=3 (rim)

Bond amplitudes

The hub A has coordination 6 (three B neighbors, three C neighbors). Each kwarg controls all bonds of that type (intra- and inter-cell):

Both default to t (uniform dice).

H = dice_hamiltonian(Lx, Ly; t_AB=1.0, t_AC=0.7)  # hub-to-B ~=hub-to-C

Spectrum: doubly degenerate flat band at E=0; dispersive bands reaching +/-3t. Real-space coordinates: dice_positions(Lx, Ly). H.sublattice_s stores the dim-3 index; H.aux_side = :post.

dice_positions(Lx, Ly) -> Matrix{Float64}

Return the (3*2^L x 2) real-space atom-position matrix for a dice (T3) lattice of 2^Lx x 2^Ly unit cells (L = Lx+Ly), consistent with the MPO site ordering.

For total 1-indexed site i: ncell = div(i-1, 3) (0-indexed unit cell, row-major) s = (i-1) % 3 + 1 (sublattice: A=1 hub, B=2 rim, C=3 rim) ix = ncell % Nx, iy = n_cell div Nx

Atom positions (triangular Bravais vectors a(1,0), a(1/2,-/2)): A: (ix + iy/2, iy-/2 ) at 0(aa/3 B: (ix + iy/2 + 1/2, iy-/2 + -/6) at 1(aa/3 C: (ix + iy/2 + 1, iy-/2 + -/3) at 2(aa/3

even_skeleton(L_chain, num_site, sites) -> MPO

Diagonal mask: 1 where iy is odd (0-based, LSB of iy = 1). Selects the odd rows for the lower honeycomb inter-row bond.

even_template(L_chain, num_site, sites) -> MPO

Diagonal mask: 1 where ix is odd AND ix ~=1. Selects B-sublattice columns (odd ix, excluding the boundary ix = 1 column which wraps into the next row). Computed as proj{odd ix} -proj{ix=1}.

generate_kin_d(sites, num_site) -> MPO

Binary-decrement MPO: |n>->|n-1>(mod 2^L). Hermitian conjugate of generate_kin_u; each term handles one borrow level.

generate_kin_u(sites, num_site) -> MPO

Binary-increment MPO: |n>->|n+1>(mod 2^L) on L = log2(numsite) qubits. Each term i handles one carry level: sigmaplus at bit i, sigma_minus on all lower bits (the bits that were 1 and get reset by the carry).

honeycomb_nnn_hamiltonian(Lx, Ly[, t[, t2]]; cutoff, maxdim) -> TBHamiltonian

Build a honeycomb tight-binding Hamiltonian with both nearest-neighbor (NN) and next-nearest-neighbor (NNN) hopping, as a TBHamiltonian.

Encoding is identical to honeycomb_sublattice_hamiltonian: L+1 sites, with the last site being the dim-2 sublattice index (A=1, B=2, postpended).

Hopping structure

NN (amplitude t): same three bonds as honeycomb_sublattice_hamiltonian (intra-cell A->B, x-shift B->A, y-shift B->A).

NNN (amplitude t2): connects same-sublattice atoms along the three triangular Bravais directions. The sublattice operator is the 2x2 identity (both A->A and B->B hop with the same amplitude t2):

• x-direction (shift +/-1): A(n) ->A(n+/-1), B(n) ->B(n+/-1)
• y-direction (shift +/-Nx): A(n) ->A(n+/-Nx), B(n) ->B(n+/-Nx)
• diagonal (shift +/-(1-Nx)): A(n) ->A(n+/-(1-Nx)), same for B

t2 may be complex; conj(t2) is used for the backward hop so that the Hamiltonian is Hermitian. For Haldane-type NNN (sublattice-dependent phases) construct the NN and NNN terms manually.

honeycomb_sublattice_hamiltonian(Lx, Ly[, t]; cutoff, maxdim) -> TBHamiltonian

Build a uniform honeycomb tight-binding Hamiltonian with an explicit 2-component sublattice index, as a TBHamiltonian.

Encoding (L+1 sites total, L = Lx+Ly):

• Sites 1..L : L position qubits for 2^L unit cells on a triangular Bravais lattice (row-major: n = ix + iy*2^Lx)
• Site L+1 : dim-2 "Honeycomb" sublattice index (A=1, B=2), postpended

Hopping structure (uniform amplitude t):

Intra-cell -one bond per unit cell: A-B (same unit cell)

Inter-cell: x (shift +1 ): B(n) ->A(n+1) -break at ix=Nx-1 y (shift +Nx): B(n) ->A(n+Nx) -no x-break needed (pure y step)

The spectrum has two Dirac cones touching at E=0 (gapless for uniform t). Use honeycomb_sublattice_positions(Lx, Ly) for real-space atom coordinates. The sublattice index is stored in H.sublattice_s; H.aux_side = :post.

honeycomb_sublattice_positions(Lx, Ly) -> Matrix{Float64}

Return the (2*2^L x 2) real-space atom-position matrix for a honeycomb lattice of 2^Lx x 2^Ly unit cells (L = Lx+Ly), consistent with the MPO site ordering.

For total 1-indexed site i: ncell = (i-1) div 2 (0-indexed unit cell, row-major) s = (i-1) % 2 + 1 (sublattice: A=1, B=2) ix = ncell % Nx, iy = n_cell div Nx

Atom positions (triangular Bravais vectors a(1,0), a(1/2,-/2)): A: (ix + iy/2, iy-/2 ) B: (ix + iy/2 + 1/2, iy-/2 + -/6 ) displaced along the intra-cell bond

interchain_hopping_honeycomb(L_chain, num_site, sites) -> MPO

Inter-row hopping for a honeycomb lattice with L_chain sites per row. Two inequivalent inter-row bonds, each with a forward and backward term:

• Upper bond (shift +Lchain+1): connectivity mask = oddtemplate * odd_skeleton
• Lower bond (shift +Lchain-1): connectivity mask = eventemplate * even_skeleton

interchain_hopping_square(L_chain, num_site, sites; hopping=Id, t=1) -> MPO

NN hopping along columns (y-direction) of a square lattice. One column step = linear shift by Lchain sites = ku composed Lchain times.

interchain_hopping_square_2nd_minus(L_chain, num_site, sites; hopping=Id, t2=1) -> MPO

NNN hopping in the (+x,-y) diagonal direction (linear shift -(L_chain-1)).. The single x-step is masked to prevent row wrap-around.

interchain_hopping_square_2nd_plus(L_chain, num_site, sites; hopping=Id, t2=1) -> MPO

NNN hopping in the (+x,+y) diagonal direction (linear shift +L_chain+1). The single x-step is masked to prevent row wrap-around.

interchain_hopping_triangle(L_chain, num_site, sites) -> MPO

Inter-row hopping for a triangular lattice. Two diagonal bonds:

• SW->NE (shift +L_chain): pure row hop, no masking needed.
• SE->NW (shift +L_chain-1): row hop + one step back in x; skeleton suppresses the spurious ix=1 wrap-around entry.

intrachain_hopping(L_chain, num_site, sites; hopping=Id, t=1) -> MPO

NN hopping along rows (x-direction) of a 2D lattice with L_chain sites per row. Hops that would wrap ix = Nx-1 ->0 are suppressed by _row_break_mpo.

kagome_hamiltonian(Lx, Ly[, t]; t_AB, t_AC, t_BC, cutoff, maxdim) -> TBHamiltonian

Build a kagome tight-binding Hamiltonian as a TBHamiltonian.

Encoding (L+1 sites, L = Lx+Ly):

• Sites 1..L : position qubits for 2^L unit cells (row-major: n = ix + iy*2^Lx)
• Site L+1 : dim-3 "Kagome" sublattice index A=1, B=2, C=3 (postpended)

Bond amplitudes

Each bond type controls both the intra-cell matrix entry and the matching inter-cell hopping term along the corresponding lattice direction:

All three default to t (uniform kagome). For anisotropic / breathing kagome pass individual values:

H = kagome_hamiltonian(Lx, Ly; t_AB=1.0, t_AC=0.8, t_BC=0.6)

Flat band: at E = -t (uniform case); dispersive bands reach up to +4t. Boundary wrapping is suppressed. Real-space coordinates: kagome_positions(Lx, Ly). H.sublattice_s stores the dim-3 sublattice index; H.aux_side = :post.

kagome_positions(Lx, Ly) -> Matrix{Float64}

Return the (3*2^L x 2) real-space atom-position matrix for a kagome lattice of 2^Lx x 2^Ly unit cells (L = Lx+Ly), consistent with the MPO site ordering.

For total 1-indexed site i: ncell = div(i-1, 3) (0-indexed unit cell, row-major) s = (i-1) % 3 + 1 (sublattice: A=1, B=2, C=3) ix = ncell % Nx, iy = n_cell div Nx

Atom positions (lattice vectors a(1,0), a(1/2,-/2)): A: (ix + iy/2, iy-/2 ) B: (ix + iy/2 + 1/2, iy-/2 ) C: (ix + iy/2 + 1/4, iy*-/2 + -/4)

kineticinterNNNSENW(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

Long-range inter-row hopping along the SE->NW diagonal. Row start-wrap suppressed by _row_break_mpo(:xplain).

kineticinterNNNSWNE(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

Long-range inter-row hopping along the SW->NE diagonal of a 2^Lx x 2^Ly square lattice. Row end-wrap suppressed by _row_break_mpo(:xplus).

kineticinterNNNtriSENW(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

SE->NW diagonal inter-row hopping for a triangular lattice. Applies _row_break_mpo(:xplain) and _row_select_mpo(:odd).

kineticinterNNNtriSWNE(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

SW->NE diagonal inter-row hopping for a triangular lattice. Applies _row_break_mpo(:xplus) and _row_select_mpo(:even) to restrict hops to the correct sublattice rows.

kineticinterNNNtri_bravais_diag(Lx, Ly, sites, hopping; apply_kwargs) -> MPO

Bravais triangular-lattice third-bond hopping: (dix=+1, diy=-1), linear shift 1-Nx. Mirrors kineticinterNNNSWNE with kd/ku swapped. Row x-wrap at ix=Nx- is suppressed by _row_break_mpo(:xplus).

kineticintra2DNNN(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

Long-range intra-row hopping on a 2^Lx x 2^Ly square lattice (nn bonds along x). Row wrap-around at ix = Nx-1 is suppressed by _row_break_mpo(:xplus).

kineticintra2DNNhex(Lx, Ly, sites, hopping, nn; apply_kwargs) -> MPO

Intra-row hopping for a honeycomb lattice. Applies _row_break_mpo(:xplus) and _row_checker_mpo to implement the alternating A/B sublattice pattern.

lieb_hamiltonian(Lx, Ly[, t]; t_AB, t_AC, cutoff, maxdim) -> TBHamiltonian

Build a Lieb tight-binding Hamiltonian as a TBHamiltonian.

Encoding (L+1 sites, L = Lx+Ly):

• Sites 1..L : position qubits for 2^L unit cells on a square Bravais lattice
• Site L+1 : dim-3 "Lieb" sublattice index A=1 (corner), B=2 (x-edge), C=3 (y-edge)

Bond amplitudes

No B-C bond exists (corner connects to edges only). Both default to t.

H = lieb_hamiltonian(Lx, Ly; t_AB=1.0, t_AC=0.5)  # anisotropic Lieb

Flat band at E=0; dispersive bands at +/-2-tAB^2+tAC^2)/- (approx +/-2t uniform). Real-space coordinates: lieb_positions(Lx, Ly). H.sublattice_s stores the dim-3 index; H.aux_side = :post.

lieb_positions(Lx, Ly) -> Matrix{Float64}

Return the (3*2^L x 2) real-space atom-position matrix for a Lieb lattice of 2^Lx x 2^Ly unit cells (L = Lx+Ly), consistent with the MPO site ordering.

For total 1-indexed site i: ncell = div(i-1, 3) (0-indexed unit cell, row-major) s = (i-1) % 3 + 1 (sublattice: A=1, B=2, C=3) ix = ncell % Nx, iy = n_cell div Nx

Atom positions (lattice vectors a(1,0), a(0,1)): A: (ix, iy ) corner B: (ix + 0.5, iy ) x-edge center C: (ix, iy + 0.5) y-edge center

odd_skeleton(L_chain, num_site, sites) -> MPO

Diagonal mask: 1 where iy is even (0-based, LSB of iy = 0). Selects the even rows for the upper honeycomb inter-row bond.

odd_template(_, num_site, sites) -> MPO

Diagonal mask: 1 where ix is odd (LSB of ix = 1). Selects the A-sublattice columns within each row for the honeycomb inter-row bonds.

skeleton(L_chain, num_site, sites) -> MPO

Diagonal mask: 0 where ix == 1 (second column, 0-indexed), 1 elsewhere. Used in the triangular lattice to exclude the wrap-around bond that enters at ix = 1 when shifting by L_chain -1 sites.

ssh_sublattice_hamiltonian(L[, t[, d]]; cutoff, maxdim) -> TBHamiltonian

Build an SSH (Su-Schrieffer-Heeger) tight-binding Hamiltonian with an explicit 2-component sublattice index, as a TBHamiltonian.

Encoding (L+1 sites total):

• Sites 1..L : L position qubits for 2^L unit cells
• Site L+1 : dim-2 "SSH" sublattice index (A=1, B=2), postpended

Hopping structure:

• Intra-cell (amplitude t+d): A->B within each unit cell
• Inter-cell (amplitude t-d): B(n) ->A(n+1)

Geometry (unit cell width = 1, 1-indexed site i over 2*2^L atoms):

• A atom in unit cell n = (i-1)div2: position [n]
• B atom in unit cell n: position [n + 0.5]

geometry_uc returns [n] for every atom in unit cell n (same for A and B).

The chain has periodic boundary conditions (B(N-1) ->A(0) inter-cell bond from the binary-increment wrap-around), consistent with all other QTT Hamiltonians.

_nth_shell_disps(H, n, Lx) -> Vector{NTuple{4,Int}}

Return the canonical set of displacement tuples (dx, dy, from_s, to_s) for the n-th nearest-neighbor shell of H.

Lattice vectors a1, a2 and sublattice offsets are extracted from three anchor cells and used to build a (2n+3)^2 reference patch analytically. This avoids large atom-index queries and only assumes Lx, Ly >= 1.

Only one member of each Hermitian-conjugate pair is returned, filtered by:

dx > 0  OR  (dx=0 AND dy>0)  OR  (dx=dy=0 AND from_s < to_s)

_shift_mpo(dx, dy, ku, kd, Id, brk_xp, Nx; apkw) -> MPO

Kinetic MPO that shifts the quantics unit-cell index by dx + dy*Nx.

• Builds the full linear displacement directly with build_shift_mpo.
• A source-cell mask removes bonds whose destination would leave the finite 2^Lx × 2^Ly grid, including row wrap-arounds.

ku, kd, Id, brk_xp must be pre-built from the same pos_sites for backward compatibility with the old call signature; the shift-MPO path only uses Id to recover the site indices.

add_hopping_2D!(H, f; Lx, Ly, nn=1, layer=nothing,
                lattice=nothing, geometry=nothing,
                maxdim=50, tol=1e-8) -> H

Add the nn-th nearest-neighbor hopping to a 2D TBHamiltonian.

Arguments

• f: hopping amplitude. Six forms are accepted:

  • Number: uniform amplitude for every bond in the shell.
  • f(dx, dy, from_s, to_s): direction-dependent scalar; called once per canonical displacement. The Hermitian conjugate uses conj(f(...)). Call get_shell_disps(H, nn; Lx, Ly) to see which tuples will be passed.
  • f(n): spatially varying amplitude as a function of the 0-based unit-cell index n = ix + iy*Nx. The modulation MPO is built once via QTCI.
  • f(ix, iy): same but with explicit 0-based 2D coordinates.
  • f(n, dx, dy, from_s, to_s): direction- and position-dependent; a separate modulation MPO is built per displacement via QTCI.
  • f(ix, iy, dx, dy, from_s, to_s): same but with explicit 2D coordinates.

  For spatially varying forms the modulation MPO V (built via get_diagonal_mpo) is applied at the destination of the forward hop: forward term uses apply(V, K_fwd), the Hermitian conjugate uses apply(K_bwd, dag(V)).

• Lx, Ly: qubit factorisation of H.L into x-bits and y-bits (required).

• layer: for layered Hamiltonians, nothing applies the hopping to every layer, an integer targets one layer, and a collection targets those layers.

• lattice / geometry: optional geometry hints for layered Hamiltonians whose H.geometry is not set. lattice may be :square, :triangular, or :honeycomb; geometry may be either a function i -> r_i or a coordinate matrix with one row per position site, or per atom when H.sublattice_s is set.

• nn: neighbor shell index (1 = nearest, 2 = next-nearest, etc.).

Compatibility Works for all 2D geometries registered in get_Hamiltonian:

• No explicit sublattice (square_2d, triangular_2d): hopping added directly as a position-space kinetic MPO.
• Explicit sublattice (honeycomb, honeycomb_nnn, kagome, lieb, dice): each displacement is augmented with the sublattice transition matrix via postpend_op.

Complexity Shell detection uses a reference patch of (2nn+3)² unit cells — O(nn²) work, independent of system size. MPO construction uses a direct build_shift_mpo per displacement, plus a simple diagonal source mask to remove finite-grid wrap-around bonds. QTCI-based spatial modulation adds O(2^L) function evaluations.

H = get_Hamiltonian("honeycomb_nnn", (t=1.0, t2=0.0); L=8, Lx=4, Ly=4)
add_hopping_2D!(H, 0.1; Lx=4, Ly=4, nn=3)

# Spatially modulated (Gaussian envelope on unit-cell index):
Nx = 2^4
add_hopping_2D!(H, n -> 0.1 * exp(-((n % Nx - Nx/2)^2)/50); Lx=4, Ly=4, nn=1)

get_shell_disps(H, nn; Lx, Ly) -> Vector{NTuple{4,Int}}

Print and return the canonical displacement tuples (dx, dy, from_s, to_s) for the nn-th nearest-neighbor shell of H.

Use this before writing a direction-dependent amplitude function for add_hopping_2D! to see exactly which tuples f(dx, dy, from_s, to_s) will be called with. The Hermitian-conjugate bond (-dx, -dy, to_s, from_s) is handled automatically and is not listed here.

H = get_Hamiltonian("honeycomb", (t=1.0,); L=4, Lx=2, Ly=2)
get_shell_disps(H, 1; Lx=2, Ly=2)
# nn=1 shell: 3 canonical bond type(s)  [Nx=4, Ny=4]
#   #   dx  dy  from_s ->to_s
#   1     0    0       1 ->2
#   2     1    0       2 ->1
#   3     0    1       2 ->1

add_hopping_2D!(H, (dx, dy, fs, ts) -> dx == 0 && dy == 0 ? 1.2 : 0.9;
                Lx=2, Ly=2, nn=1)

_aa_interlayer_mpo(sites; t_inter=1.0) -> MPO

Build the interlayer coupling MPO for AA stacking: each site in layer 1 couples on-site to the same site in layer 2. The coupling operator is simply t_inter * Identity.

_bernal_interlayer_mpo(L, sites; t_inter=1.0, cutoff=1e-8) -> MPO

Build the interlayer coupling MPO for Bernal (AB) stacking on a honeycomb lattice. In the quantics site ordering, sublattice-A sites have 1-based odd indices and sublattice-B sites have 1-based even indices.

Bernal stacking places each A site in layer 1 directly above a B site in layer 2 (the B site of the same unit cell, index A+1). The interlayer operator in position space is therefore

V = t_inter · (K_u D_A + D_A K_d)

where DA is the A-sublattice projector and Ku/Kd are the ±1 shift operators. This is symmetric (Hermitian for real tinter).

bilayer_hamiltonian(lattice, Lx, Ly;
    stacking=:AA, t_intra=1.0, t_inter=0.3,
    cutoff=1e-8, maxdim=200) -> (MPO, Vector{<:Index})

Build a bilayer tight-binding Hamiltonian as an MPO.

Site encoding (L+1 sites total, L = Lx + Ly):

• Site 1 : layer index (dim = 2)
• Sites 2…L+1 : L position qubits (quantics binary, row-major)

Arguments

• lattice : :square, :triangular, or :honeycomb
• Lx, Ly : each layer has 2^Lx × 2^Ly sites

Keyword arguments

• stacking : :AA (on-site) or :Bernal (A₁↔B₂, honeycomb only)
• t_intra : intra-layer NN hopping amplitude
• t_inter : interlayer hopping amplitude
• cutoff : MPO truncation cutoff
• maxdim : maximum bond dimension of the final MPO

The assembled Hamiltonian is

H = Σₖ Pₖ ⊗ H_mono  +  (|1⟩⟨2| + |2⟩⟨1|) ⊗ V

where V is the exact interlayer MPO for the chosen stacking.

Returns (H_total, ext_sites) where ext_sites[1] is the layer index and ext_sites[2:end] are the L position qubits.

interlayer_mpo(lattice, stacking, Lx, Ly, sites;
               t_inter=1.0, cutoff=1e-8) -> MPO

Build the position-space interlayer coupling MPO for the given stacking. The returned operator V satisfies

H_inter = t_inter · (|k⟩⟨l| ⊗ V + |l⟩⟨k| ⊗ V)

for each pair of adjacent layers k, l.

Supported stackings

• :AA — on-site (identity in position space); any lattice
• :Bernal — A₁↔B₂ coupling within each unit cell; :honeycomb only

For general (non-commensurate) interlayer functions, pass a function f(i,j) to hopping2MPO directly and use prepend_layer_hopping.

multilayer_hamiltonian(lattice, Lx, Ly, n_layers;
    stacking=:AA, t_intra=1.0, t_inter=0.3,
    cutoff=1e-8, maxdim=200) -> (MPO, Vector{<:Index})

Generalisation of bilayer_hamiltonian to n_layers layers. The same stacking and t_inter are used for every adjacent pair.

Returns (H_total, ext_sites) with the same site encoding as bilayer_hamiltonian, extended to a dim = n_layers layer index.

lattice_positions(lattice, Lx, Ly; angle_deg=0.0) -> Matrix{Float64}

Return an N×2 matrix of real-space positions for a 2^Lx × 2^Ly patch of lattice ∈ {:square, :triangular, :honeycomb}, optionally rotated by angle_deg degrees about the geometric centroid of the lattice.

Site ordering matches the quantics row-major encoding: n = ix + iy·2^Lx (0-indexed).

monolayer_hamiltonian(lattice, Lx, Ly, sites; t=1.0, cutoff=1e-8) -> MPO

NN tight-binding Hamiltonian MPO on L = Lx+Ly qubit sites for lattice ∈ {:square, :triangular, :honeycomb} with uniform hopping t. Delegates to build_hamiltonian and replaces internal site indices with sites.

postpend_layer_hopping(H_mpo, layer_s, k, l) -> MPO

Append the off-diagonal operator |k⟩⟨l| on layer_s (1-based) to the end of H_mpo.

postpend_layer_projector(H_mpo, layer_s, k) -> MPO

Append the diagonal projector |k⟩⟨k| on layer_s (1-based) to the end of H_mpo.

prepend_layer_hopping(H_mpo, layer_s, k, l) -> MPO

Prepend the off-diagonal operator |k⟩⟨l| on layer_s (1-based). Equivalent to prepend_op(H_mpo, layer_s, k, l).

prepend_layer_projector(H_mpo, layer_s, k) -> MPO

Prepend the diagonal projector |k⟩⟨k| on layer_s (1-based). Equivalent to prepend_op(H_mpo, layer_s, k).

twisted_bilayer_hamiltonian(lattice, Lx, Ly, θ_deg;
    t_intra=1.0, t_inter=0.3, α_decay=1/16.0,
    tol=1e-6, cutoff=1e-8, maxdim=200) -> (MPO, Vector{<:Index})

Convenience wrapper: two layers, layer 1 at 0° and layer 2 at θ_deg. Equivalent to twisted_multilayer_hamiltonian(lattice, Lx, Ly, [0.0, θ_deg]; …).

twisted_multilayer_hamiltonian(lattice, Lx, Ly, angles_deg;
    t_intra=1.0, t_inter=0.3, α_decay=1/16.0,
    tol=1e-6, cutoff=1e-8, maxdim=200) -> (MPO, Vector{<:Index})

Build a twisted multilayer tight-binding Hamiltonian as an MPO.

Site encoding (L+1 sites, L = Lx+Ly):

• Site 1 : layer index (dim = n_layers) ← replaces the bilayer qubit
• Sites 2…L+1 : L position qubits (quantics binary, row-major)

Each layer k is rigidly rotated to angle angles_deg[k] (degrees). The interlayer coupling is exponentially decaying: V_{kl}(i,j) = t_inter * exp(−α_decay * |r_k[i] − r_l[j]|).

The assembled Hamiltonian is: H = Σ_k P_k ⊗ H_mono + Σ_{k<l} (|k⟩⟨l| ⊗ V_{kl} + |l⟩⟨k| ⊗ V_{lk})

which is Hermitian for real V{kl} since V{lk}(i,j) = V{kl}(j,i) = V{kl}^T.

Arguments

• lattice : :square, :triangular, or :honeycomb
• Lx, Ly : 2^Lx × 2^Ly sites per layer
• angles_deg : twist angles (°) for each layer; length = n_layers

Keyword arguments

• t_intra : intra-layer NN hopping amplitude (same for all layers)
• t_inter : interlayer hopping amplitude
• α_decay : exponential decay constant for interlayer hopping
• tol : QTCI tolerance passed to hopping2MPO
• cutoff : MPO truncation cutoff used throughout assembly
• maxdim : maximum bond dimension of the final MPO

Returns (H_total, ext_sites) where ext_sites[1] is the layer index and ext_sites[2:end] are the L position qubits.

mask_hamiltonian(H, sdf; sigma=0.3, tol=1e-8, maxdim=200, cutoff=1e-8)
    -> TBHamiltonian

Restrict H to a flake geometry defined by sdf by applying the smooth diagonal mask

M_i = σ( sdf(rᵢ) / sigma ),   H_flake = M · H · M

where rᵢ = H.geometry(i) is the real-space position of site i (1-indexed), σ is the logistic sigmoid, and sigma controls boundary sharpness in the same coordinate units as H.geometry.

The mask MPS is learned via QTCI, so even sub-lattice-spacing smoothing (sigma ≪ 1) is captured accurately without explicit enumeration.

Keyword arguments

• sigma : sigmoid half-width (lattice units). Default 0.3.
• tol : QTCI tolerance for the mask MPS. Default 1e-8.
• maxdim : max bond dim during the two M·H·M products. Default 200.
• cutoff : SVD truncation cutoff. Default 1e-8.

Restrictions

• Requires H.geometry to be set (all preset geometries provide this).
• Must be called before add_spin!, add_superconductivity!, or bilayer construction.

Examples

H = get_Hamiltonian("square_2d", 1.0; L=8, Lx=16)

# Disk flake
H_disk = mask_hamiltonian(H, sdf_disk(8.0, 8.0, 6.0); sigma=0.3)

# Triangular flake (intersection of three half-planes)
tri = sdf_intersect(
    sdf_intersect(sdf_halfplane(0, 1, 2.0),
                  sdf_halfplane(-sqrt(3)/2, -0.5, -8.0)),
                  sdf_halfplane( sqrt(3)/2, -0.5, -8.0))
H_tri = mask_hamiltonian(H, tri; sigma=0.2)

# Ring (disk minus inner disk)
H_ring = mask_hamiltonian(H, sdf_annulus(8.0, 8.0, 3.0, 7.0); sigma=0.3)

sdf_annulus(cx, cy, r_in, r_out) -> (x, y) -> Float64

Signed distance function for a circular annulus with inner radius r_in and outer radius r_out, centred at (cx, cy).

sdf_convex_polygon(vertices) -> (x, y) -> Float64

Signed distance function for the convex polygon with the given vertices (a vector of (x, y) tuples listed in counter-clockwise order). Positive inside.

Implemented as the minimum signed distance to each edge's inward half-plane.

sdf_disk(cx, cy, r) -> (x, y) -> Float64

Signed distance function for a disk of radius r centred at (cx, cy). Positive inside, negative outside.

sdf_halfplane(nx, ny, d) -> (x, y) -> Float64

Signed distance function for the half-plane nx·x + ny·y ≥ d. (nx, ny) need not be unit-length — normalisation is applied internally. Positive on the side where the inequality holds.

sdf_intersect(f, g) -> SDF

Boolean intersection: positive where both f and g are positive (min(f, g)).

sdf_interval(x_lo, x_hi) -> x -> Float64

Signed distance function for the interval [x_lo, x_hi]. Positive inside, negative outside.

sdf_rect(cx, cy, w, h) -> (x, y) -> Float64

Signed distance function for an axis-aligned rectangle with full width w and full height h, centred at (cx, cy). Positive inside.

sdf_subtract(f, g) -> SDF

Boolean subtraction: positive inside f and outside g (min(f, -g)).

sdf_union(f, g) -> SDF

Boolean union: positive where either f or g is positive (max(f, g)).

add_tjunction!(H, t_j; junction_site=0, coupling=nothing,
               cutoff=1e-8, maxdim=200) -> H

Extend a 1D TBHamiltonian to a T/Y-junction by postpending a dim-3 "TJunction" index. The three states |1⟩≡|-1⟩, |2⟩≡|0⟩, |3⟩≡|+1⟩ label three identical branches; all share the same intra-chain Hamiltonian.

A junction coupling is added at junction_site (0-indexed):

H_TJ = I₃ ⊗ H_chain  +  J(t_j) ⊗ |junction_site⟩⟨junction_site|

The default 3×3 coupling matrix is the symmetric Y-junction:

J = t_j * [0 1 1 ; 1 0 1 ; 1 1 0]

Pass a custom 3×3 matrix via coupling to override (e.g., for a linear T-junction where only two of the three branches are mutually coupled).

Requirements: H must have no existing auxiliary indices (spin, Nambu, layer, sublattice).

Geometry: H.geometry is set to a 2D Y-junction layout with the three branches radiating at 0°, 120°, 240° from the origin (junction site). See tjunction_positions for the atom index convention.

After the call, H.sublattice_s stores the TJunction index and H.aux_side = :post. All caches are invalidated.

Examples

H = get_Hamiltonian("chain_1d", 1.0; L=8)

# Symmetric Y-junction at the left end (site 0)
add_tjunction!(H, 0.5)

# Junction at the right end (site N-1)
H2 = get_Hamiltonian("chain_1d", 1.0; L=8)
add_tjunction!(H2, 0.5; junction_site=H2.N - 1)

# Linear T-junction: branch 1 ↔ 2 and 2 ↔ 3, but not 1 ↔ 3
H3 = get_Hamiltonian("chain_1d", 1.0; L=8)
add_tjunction!(H3, 0.5; coupling=[0 1 0; 1 0 1; 0 1 0])

tjunction_hamiltonian(L, t, t_j; junction_site=0, coupling=nothing,
                      boundary=:open, cutoff=1e-8, maxdim=200)
    -> TBHamiltonian

Build a T/Y-junction Hamiltonian: three identical 1D nearest-neighbour chains with hopping amplitude t meeting at a common junction site.

Encoding (L+1 sites total):

• Sites 1…L : L position qubits (chain length = 2^L per branch)
• Site L+1 : dim-3 "TJunction" index (postpended), stored in H.sublattice_s

Chain-label states (in H.sublattice_s): 1 ↔ branch |-1⟩, 2 ↔ branch |0⟩, 3 ↔ branch |+1⟩

Geometry (H.geometry): 2D Y-junction with branches at 0°, 120°, 240° from the junction site (placed at the origin). See tjunction_positions.

Arguments

• L : number of position qubits (chain length = 2^L per branch)
• t : intra-chain NN hopping amplitude
• t_j : junction coupling amplitude

Keyword arguments

• junction_site : 0-indexed site where the chains meet (default 0, left end)
• coupling : custom 3×3 coupling matrix (default: symmetric Y-junction t_j * [0 1 1; 1 0 1; 1 1 0])
• boundary : :open (default) or :periodic for the intra-chain hoppings
• cutoff : MPO truncation cutoff
• maxdim : maximum MPO bond dimension

Examples

# Symmetric Y-junction, junction at left end (site 0)
H = tjunction_hamiltonian(8, 1.0, 0.5)

# Junction at the right end (site 2^L - 1)
H = tjunction_hamiltonian(8, 1.0, 0.5; junction_site=255)

# Linear T-junction (branches 1↔2 and 2↔3 only)
H = tjunction_hamiltonian(8, 1.0, 0.5; coupling=[0 1 0; 1 0 1; 0 1 0])

tjunction_index() -> Index

Create a dim-3 Index tagged "TJunction" for a Y/T-junction of three chains. States: 1 = branch |-1⟩, 2 = branch |0⟩, 3 = branch |+1⟩.

tjunction_lattice_hamiltonian(Lx, Ly, L, t, t_j, t_inter;
                               junction_site=0, coupling=nothing,
                               boundary=:open, cutoff=1e-8, maxdim=200)
    -> TBHamiltonian

Build a triangular lattice of T/Y-junctions with 2^Lx × 2^Ly unit cells, where each unit cell is a T-junction with three arms of 2^L sites.

Site ordering (Lx + Ly + L + 1 sites total):

• Sites 1..Ly : iy bits of the triangular lattice (MSB first)
• Sites Ly+1..Ly+Lx : ix bits of the triangular lattice (MSB first)
• Sites Ly+Lx+1..Ly+Lx+L : chain position qubits within each arm (MSB first)
• Site Ly+Lx+L+1 : dim-3 TJunction sublattice index (branch 1/2/3)

Hamiltonian structure:

1. Intra-arm: nearest-neighbor hopping t along each chain.
2. Junction coupling: t_j * [0 1 1; 1 0 1; 1 1 0] at junction_site (acts on the branch index at the shared chain site).
3. Inter-cell: hopping t_inter between the far end (2^L - 1) of one arm and the far end of the appropriate arm in the adjacent unit cell:
   • x-shift (-1) : arm 1 (330°) ↔ arm 3 (210°)
   • y-shift (+Nx) : arm 3 (210°) ↔ arm 2 (90°)
   • diagonal (+Nx-1) : arm 2 (90°) ↔ arm 1 (330°)

Arguments

• Lx, Ly : log₂ of unit cell count along x and y
• L : log₂ of sites per arm
• t : intra-arm NN hopping amplitude
• t_j : junction coupling amplitude
• t_inter : inter-unit-cell hopping amplitude at arm ends

Keyword arguments

• junction_site : 0-indexed junction position within each arm (default 0)
• coupling : custom 3×3 coupling matrix (overrides t_j default)
• boundary : :open or :periodic for intra-arm hopping
• cutoff, maxdim : MPO truncation parameters

H.sublattice_s stores the dim-3 TJunction index; H.aux_side = :post. H.L = Lx + Ly + L, H.N = 2^(Lx+Ly) * 2^L (unitcells × armlength).

Examples

# 4×4 triangular lattice of T-junctions, arms of length 8
H = tjunction_lattice_hamiltonian(2, 2, 3, 1.0, 0.5, 0.8)

# Weaker inter-cell coupling
H = tjunction_lattice_hamiltonian(2, 2, 3, 1.0, 0.5, 0.3; junction_site=0)

tjunction_positions(N, junction_site) -> Matrix{Float64}

Return a (3N × 2) real-space position matrix for a Y-junction of three chains, each of length N, meeting at junction_site (0-indexed).

Atom index convention (branch-fast, matching kagome): atom i (1-indexed): n = (i-1) ÷ 3 — 0-indexed chain position (0…N-1) s = (i-1) % 3 + 1 — branch (1=|-1⟩, 2=|0⟩, 3=|+1⟩)

Branch directions (branches radiate symmetrically from the junction): branch 1: angle 0° branch 2: angle 120° branch 3: angle 240°

The junction site is placed at the origin; position along each branch is the signed distance n − junction_site (positive = away from junction).