Constrained Hilbert spaces and conserved quantum numbers

The basis states of a Hilbert space can be restricted to a subset of states and organized into sectors with different quantum numbers.

You can apply a constraint either when constructing a space, hilbert_space(f, labels, constraint), when doing tensor products tensor_product((H1, H2, ...); constraint), or afterwards with constrain_space(H, constraint).

We define two common types of constraints for conserved quantum numbers:

• NumberConservation — Conserves a (possibly weighted) particle number
• ParityConservation — Conserves parity

As well as some types to help with constructing custom constraints:

• FilterConstraint — Flexible filtering using a boolean test on states (iterates over all states)
• SectorConstraint — Groups states into sectors using a custom rule (iterates over all states)
• AdditiveConstraint — Generalization of number conservation, which can be used for more complicated additive quantum numbers (can be used by iterating over all states, or by branch pruning)
• BranchConstraint — Efficient generation of states using branch pruning

One can also combine constraints by taking products.

NumberConservation and ParityConservation

• NumberConservation(total=missing, subspaces=missing, weights=missing): Conserves a (possibly weighted) particle number.
  • total can be a single integer or a collection of allowed values.
  • subspaces can restrict the count to selected subspaces (or modes).
  • weights lets each selected subspace contribute with a different integer weight.
  Examples: NumberConservation(), NumberConservation(2), NumberConservation(0:1, Hup), NumberConservation(-1:1, [Hl, Hr], [1, -1]).

using FermionicHilbertSpaces
@fermions f
Hf = hilbert_space(f, 1:3, NumberConservation(2)) # Single sector with 2 particles

3-dimensional SectorHilbertSpace
Parent: (f[1, 2, 3])
Sectors: [2 (3-dim)]

@boson b
Hb = hilbert_space(b, 10)
H = tensor_product(Hf, Hb; constraint = NumberConservation(0:1, [Hf,Hb] , [1,-1])) # the number of fermions in Hf minus the number of bosons in Hb must be 0 or 1
basisstates(H)

6-element Vector{FermionicHilbertSpaces.ProductState{Tuple{FockNumber{UInt64}, FermionicHilbertSpaces.BosonicState}}}:
 |11⟩⨯|2⟩
 |101⟩⨯|2⟩
 |011⟩⨯|2⟩
 |11⟩⨯|1⟩
 |101⟩⨯|1⟩
 |011⟩⨯|1⟩

• ParityConservation(parities=[-1, 1], subspaces=missing): Conserves fermionic parity.
  • Allowed parities are -1 (odd) and 1 (even).
  • ParityConservation() keeps both sectors (odd first, then even).
  • ParityConservation([1]) or ParityConservation(1) keeps only even parity.
  • subspaces lets you enforce parity on a specific part of a tensor-product space.

constrain_space(H, ParityConservation(1, [Hb])) # keep only even parity of the bosonic part

3-dimensional SectorHilbertSpace
Parent: SectorHilbertSpace(ConstrainedSpace(ProductSpace(80-dim, 2 factors), 30-dim), 6-dim)
Sectors: [1 (3-dim)]

FilterConstraint and SectorConstraint

• FilterConstraint(...): A flexible filtering constraint.

  • FilterConstraint(reducer) applies reducer(state) and keeps states where it returns true.
  • FilterConstraint(subspaces, subspace_functions, reducer) first applies the subspace_functions to the states in each subspace, then passes them to the reducer function.

  Use this for custom selection rules that are easiest to express as a boolean test on complete states.

• SectorConstraint(...): Like FilterConstraint, but for grouping states into sectors.

  • The reducer should return a sector label/key (or missing to discard a state).
  • States with the same returned key are collected into the same sector.

  Use this when you want an explicit block structure from a custom rule.

@bosons b
H = hilbert_space(b, 1:2, 3)
using FermionicHilbertSpaces: FilterConstraint, SectorConstraint, particle_number
constraint = FilterConstraint(issorted, particle_number, factors(H))
basisstates(constrain_space(H, constraint)) # keep only states with sorted particle numbers

6-element Vector{FermionicHilbertSpaces.ProductState{Tuple{FermionicHilbertSpaces.BosonicState, FermionicHilbertSpaces.BosonicState}}}:
 |0⟩⨯|0⟩
 |0⟩⨯|1⟩
 |1⟩⨯|1⟩
 |0⟩⨯|2⟩
 |1⟩⨯|2⟩
 |2⟩⨯|2⟩

constraint = SectorConstraint(numbers -> issorted(numbers) ? first(numbers) : missing,  particle_number, factors(H)) #  keep only states with sorted particle numbers, organize them into sectors according to the number of particles in the first mode
constrain_space(H, constraint).qn_to_states

OrderedCollections.OrderedDict{Int64, Vector{FermionicHilbertSpaces.ProductState{Tuple{FermionicHilbertSpaces.BosonicState, FermionicHilbertSpaces.BosonicState}}}} with 3 entries:
  0 => [|0⟩⨯|0⟩, |0⟩⨯|1⟩, |0⟩⨯|2⟩]
  1 => [|1⟩⨯|1⟩, |1⟩⨯|2⟩]
  2 => [|2⟩⨯|2⟩]

Product

Constraints can be combined by multiplying them:

Hs = [hilbert_space(f,2k-1:2k) for k in 1:3]
constraint = prod(NumberConservation(0:1, H) for H in Hs) * NumberConservation(2) # each mode can be occupied by at most one fermion, and the total number of fermions must be 2
H = tensor_product(Hs; constraint)

12-dimensional SectorHilbertSpace
Parent: (f[1, 2, 3, 4, 5, 6])
Sectors: [[0, 1, 1, 2] (4-dim), [1, 0, 1, 2] (4-dim), [1, 1, 0, 2] (4-dim)]

BranchConstraint

• BranchConstraint(f): For efficient constrained generation using branch pruning.
  • f(partial_state, depth, spaces) -> Bool is evaluated during state construction.
  • Returning false prunes the branch immediately.
  When branches can be pruned early, this avoids the exponentially large full Hilbert space.

Examples

Spin

This package does not know anything about spin, but one can treat spin just as an extra label as follows:

using FermionicHilbertSpaces
@fermions f
Hup = hilbert_space(f, [(i, :↑) for i in 1:4])
Hdn = hilbert_space(f, [(i, :↓) for i in 1:4])
H = tensor_product(Hup, Hdn)

256-dimensional FermionicSpace
Modes: (f[(1, :↑), (2, :↑), (3, :↑), (4, :↑), (1, :↓), (2, :↓), (3, :↓), (4, :↓)])

If spin is conserved, one can use

Hsectors = constrain_space(H, NumberConservation(Hup)*NumberConservation(Hdn))

256-dimensional SectorHilbertSpace
Parent: (f[(1, :↑), (2, :↑), ..., (3, :↓), (4, :↓)])
Sectors: 25 total [[0, 0] (1-dim), ..., [4, 4] (1-dim)]

to sort states according to the number of fermions with spin up and down. However, this package can't help to sort states into sectors with different total angular momentum, because that requires taking superpositions of different fock states.

To pick out the sector with 2 fermions with spin up and 0 fermions with spin down, one can extract it from the hilbert space defined above using sector as

sector([2,0], Hsectors)

6-dimensional ConstrainedSpace
Parent: (f[(1, :↑), (2, :↑), ..., (3, :↓), (4, :↓)])

or more efficiently by only constructing that specific sector in the first place

constrain_space(H, NumberConservation(2, Hup)*NumberConservation(0, Hdn))

6-dimensional SectorHilbertSpace
Parent: (f[(1, :↑), (2, :↑), ..., (3, :↓), (4, :↓)])
Sectors: [[2, 0] (6-dim)]

Hubbard model

For N fermions, the full hilbert space is exponentially large in N. However, due to conservation laws, we may be interested in only a small subspace. If you use a quantum number which consists of only products of number conservations, this package attempts to find the subspaces without enumerating the full hilbert space.

As an example, consider the Hubbard model on N sites

\[H = -t \sum_{i,\sigma} (c_{i,\sigma}^\dagger c_{i+1,\sigma} + \mathrm{h.c}) + U \sum_i n_{i,↑} n_{i,↓}.\]

which conserves the number of spin up and spin down fermions separately. Let's define a function to get the hamiltonian with symbolic fermions

using FermionicHilbertSpaces
function hubbard_hamiltonian(c, N, t, U)
    spins = (:↑,:↓)
    sum(-t * c[i,σ]' * c[i+1,σ] + hc for σ in spins for i in 1:N-1) + sum(U * c[i,:↑]'c[i,:↑] * c[i,:↓]'c[i,:↓] for i in 1:N)
end

hubbard_hamiltonian (generic function with 1 method)

Let's find the matrix representation of the hamiltonian in the sector with N_up spin up fermions and N_down spin down fermions. To find this subspace we do

@fermions f
N = 20
Nup = 2
Ndn = 1
Hup = hilbert_space(f, [(i, :↑) for i in 1:N], NumberConservation(Nup))
Hdn = hilbert_space(f, [(i, :↓) for i in 1:N], NumberConservation(Ndn))
H = tensor_product(Hup, Hdn)

3800-dimensional SectorHilbertSpace
Parent: ConstrainedSpace((f[(1, :↑), (2, :↑), ..., (19, :↓), (20, :↓)]), 3800-dim)
Sectors: [[2, 1] (3800-dim)]

The full hilbert space is of size 4^20 ≈ 10^12, but the sector with 2 spin up and 1 spin down fermion is only of size 3800 and is generated without constructing the full hilbert space. Finally, we can get the matrix representation of the hamiltonian in this sector as

symham = hubbard_hamiltonian(f, N, 1.0, 4.0)
ham = matrix_representation(symham, H)

3800×3800 SparseArrays.SparseMatrixCSC{Float64, Int64} with 21280 stored entries:
⎡⠻⣦⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⠻⣦⎦

No double occupation

When the onsite Coulomb interaction is very strong, there is a large energy penalty for double occupation of a site. In that case, we can restrict the Hilbert space to not allow double occupation of any site. Consider the site k, which has two labels (k, :↑) and (k, :↓). We can use number_conservation(0:1, label -> label[1] == k) which says that the sum of occupation numbers of all labels where the first element of the label equals k is contained in the set 0:1. To impose this for all sites, we take the product over all sites.

no_double_occ = prod(NumberConservation(0:1, [f[(k, :↑)], f[(k, :↓)]]) for k in 1:N)
H_ndo = constrain_space(H, no_double_occ)

3420-dimensional SectorHilbertSpace
Parent: SectorHilbertSpace(ConstrainedSpace((f[(1, :↑), (2, :↑), ..., (19, :↓), (20, :↓)]), 3800-dim), 3800-dim)
Sectors: 1140 total [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1] (3-dim), ..., [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] (3-dim)]

This quantum number is a product of number conservations, so the sector is constructed without enumerating the full Hilbert space.

The matrix representation of the hamiltonian in this sector can be constructed as before, but now we need to specify projection = true as the symbolic hamiltonian maps states in the subspace to states outside the subspace. The keyword projection = true says to ignore those terms.

symham = hubbard_hamiltonian(f, N, 1, 0)
ham_ndo = matrix_representation(symham, H_ndo; projection = true)

3420×3420 SparseArrays.SparseMatrixCSC{Int64, Int64} with 17442 stored entries:
⎡⢻⣶⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠙⢿⣷⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠻⡿⣯⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢮⠻⣦⡙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠳⣌⠻⣦⡉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⡈⠻⣦⡈⠓⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠈⠻⣦⡀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠈⠻⣦⡀⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠈⠻⣦⡀⠈⠳⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠈⠻⣦⡀⠀⠙⢦⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡄⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢦⡀⠀⠉⠻⣦⡀⠀⠈⠓⢦⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠈⠻⣦⡀⠀⠀⠙⢦⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠈⢿⣷⡀⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠈⠻⣦⣄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠙⢿⣷⎦

Fractionalized hilbert space with BranchConstraint

The $t-J_z$ model is

\[H = -t \sum_{i,\sigma} (c_{i,\sigma}^\dagger c_{i+1,\sigma} + \mathrm{h.c}) + J_z \sum_{i} S^z_i S^z_{i+1},\]

where $S^z_i = n_{i,↑} - n_{i,↓}$ and double occupation is forbidden. This model features a fractionalized hilbert space where the Hilbert space splits into exponentially many dynamically disconnected sectors, see [1910.06341]. We implement the hamiltonian as

function tjz(c, N,t,Jz)
    spins = (:↑,:↓)
    Sz(i) = c[i,:↑]'c[i,:↑] - c[i,:↓]'c[i,:↓]
    -t*sum(c[i,σ]'c[i+1,σ] + hc for σ in spins for i in 1:N-1) + Jz*sum(Sz(i)Sz(i+1) for i in 1:N-1)
end

tjz (generic function with 1 method)

To construct the hilbert space in a specific sector, we need to go beyond simple particle number constraints. Each sector is defined by a spin ordering, e.g. '[:↑, :↑, :↓, :↑, :↓]'. Each state in this sector has three spin up electrons occupied, two spin down electrons occupied, and they come in that specific order with possible holes between them.

This package can help construct these small sectors while avoiding the exponentially large hilbert space, by

function valid_state(partials, depth::Int, spaces, expected_order)
    # partials is a list of 2N focknumbers, but only the ones up to 'depth' are assigned.
    # spaces is a list of 2N fermionic modes
    exc_count = 0
    prev_site = 0
    for n in 1:depth
        if count_ones(partials[n]) == 1
            exc_count += 1
            site, spin = only(spaces[n].modes).label
            exc_count > length(expected_order) && return false # too many excitations
            site == prev_site && return false # two excitations on the same site (this uses the fact that we ordered the spaces by site)
            prev_site = site
            expected_order[exc_count] == spin || return false # check if the spin of the excitation matches the expected order
        end
    end
    depth == length(spaces) && return exc_count == length(expected_order)
    return true
end
using FermionicHilbertSpaces: BranchConstraint
spin_order_constraint(order) = BranchConstraint((p,d,s) -> valid_state(p,d,s,order))

spin_order_constraint (generic function with 1 method)

N = 16
labels = [(i, s) for i in 1:N for s in (:↑,:↓)]
Hfull = hilbert_space(f,  labels)

4294967296-dimensional FermionicSpace
Modes: 32 total (f[(1, :↑), (1, :↓), (2, :↑), ..., (15, :↓), (16, :↑), (16, :↓)])

This space is too large to deal with directly, but we can constrain it to the sector with spin order '[:↑, :↑, :↓, :↑, :↓]' as

Hfrac = hilbert_space(f, labels, spin_order_constraint([:↑, :↑, :↓, :↑, :↓]))

4368-dimensional ConstrainedSpace
Parent: (f[(1, :↑), (1, :↓), ..., (16, :↑), (16, :↓)])

which has only dimension 4368. We can then construct hamiltonian in this sector by

symham = tjz(f, N, 1, 1/4)
ham = matrix_representation(symham, Hfrac; projection = true)

4368×4368 SparseArrays.SparseMatrixCSC{Float64, Int64} with 33606 stored entries:
⎡⢻⣶⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠻⣿⣿⡶⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠘⢯⡻⣮⣅⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⢥⡙⠿⣧⡈⠳⠄⢠⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠙⢦⡈⢿⣷⣤⡀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⣁⠀⠻⠿⣧⡀⠀⠙⢦⡀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⢻⣶⣄⡀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠹⣿⣿⣦⡀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠈⠻⡿⣯⡀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠈⠻⣦⣄⠀⠀⠀⠉⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠙⣿⣿⣦⡀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠈⠻⣿⣿⡲⣄⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠃⠀⠀⠀⠘⢮⡻⣮⡓⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣄⠀⠀⠀⠀⠀⠙⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠙⠆⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠈⢿⣷⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠻⣿⣿⣦⣀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠈⢻⡻⣮⡳⣄⣀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠙⢮⠻⣦⡈⠳⣄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠁⠀⠀⠀⠀⠀⠘⢦⡈⠻⣦⡈⠁⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠆⠈⢿⣷⎦

We can use subregion to find the hilbert space of a subsystem, taking into account the constraint. The full hilbert space of the left half of the system is

Hsub_full = hilbert_space(f,  [(i, s) for i in 1:div(N,2) for s in (:↑,:↓)])

65536-dimensional FermionicSpace
Modes: 16 total (f[(1, :↑), (1, :↓), (2, :↑), ..., (7, :↓), (8, :↑), (8, :↓)])

but taking into account the constraint, we have

Hsub_frac = subregion(Hsub_full, Hfrac)

219-dimensional ConstrainedSpace
Parent: (f[(1, :↑), (1, :↓), ..., (8, :↑), (8, :↓)])

and now one can do partial traces partial_trace(mat, Hfrac => Hsub_frac).