GitHub

Conserved quantum numbers

The basis states of a Hilbert space can be organized into sectors with different quantum numbers. Only quantum numbers which are diagonal in the fock state basis are supported.

Use hilbert_space(labels, qn), where qn can be e.g.

  • ParityConservation(): Puts all states with odd parity first, then all states with even parity.
  • ParityConservation(p::Int): Only contains states with parity p (-1 for odd and 1 for even).
  • number_conservation(): Sorts basis states by the number of fermions.
  • number_conservation(sectors::Union{Int,Vector{Int}}): Only contains states with the number(s) in the list sectors.
  • number_conservation(sectors, weight_function): 'weight_function' is a function that takes a label and returns an integer weight that represents the contribution of that fermion to the total number. The returned states will have total weighted number of fermions contained in sectors.
  • number_conservation(weight_function): As above, but allowing all sectors.
  • Products of the above quantum numbers, which sorts states according to each factor in the product.

Using number_conservation with small sectors avoids the exponentially large Hilbert space.

Spin

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

using FermionicHilbertSpaces
labels = Base.product(1:4, (:↑,:↓))
H = hilbert_space(labels)
256-dimensional SimpleFockHilbertSpace:
8 fermions: [(1, :↑), (2, :↑), (3, :↑), (4, :↑), (1, :↓), (2, :↓), (3, :↓), (4, :↓)]

If spin is conserved, one can use 

H = hilbert_space(labels, number_conservation(label -> :↑ in label) * number_conservation(label -> :↓ in label))
256-dimensional SymmetricFockHilbertSpace:
8 fermions: [(1, :↑), (2, :↑), (3, :↑), (4, :↑), (1, :↓), (2, :↓), (3, :↓), (4, :↓)]
FockSymmetry(Number conservation for 2 subsets)

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, or construct it directly

FermionicHilbertSpaces.sector((2,0), H)
hilbert_space(labels, number_conservation(2, label -> :↑ in label) * number_conservation(0, label -> :↓ in label))
6-dimensional SymmetricFockHilbertSpace:
8 fermions: [(1, :↑), (2, :↑), (3, :↑), (4, :↑), (1, :↓), (2, :↓), (3, :↓), (4, :↓)]
FockSymmetry(Number conservation for 2 subsets)

Small subspaces in large systems

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(N, t, U)
    @fermions c
    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

N = 20
labels = Base.product(1:N, (:↑,:↓))
N_up = 2
N_down = 1
qn_spin = number_conservation(N_up, label -> :↑ in label) * number_conservation(N_down, label -> :↓ in label)
H = hilbert_space(labels, qn_spin)
3800-dimensional SymmetricFockHilbertSpace:
40 fermions: [(1, :↑), (2, :↑), (3, :↑), (4, :↑), (5, :↑), (6, :↑), (7, :↑), (8, :↑), (9, :↑), (10, :↑)  …  (11, :↓), (12, :↓), (13, :↓), (14, :↓), (15, :↓), (16, :↓), (17, :↓), (18, :↓), (19, :↓), (20, :↓)]
FockSymmetry(Number conservation for 2 subsets)

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(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.

qn_no_double_occ = prod(number_conservation(0:1, label -> label[1] == k) for k in 1:N)
qn = qn_spin * qn_no_double_occ
H_ndo = hilbert_space(labels, qn)
3420-dimensional SymmetricFockHilbertSpace:
40 fermions: [(1, :↑), (2, :↑), (3, :↑), (4, :↑), (5, :↑), (6, :↑), (7, :↑), (8, :↑), (9, :↑), (10, :↑)  …  (11, :↓), (12, :↓), (13, :↓), (14, :↓), (15, :↓), (16, :↓), (17, :↓), (18, :↓), (19, :↓), (20, :↓)]
FockSymmetry(Number conservation for 22 subsets)

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(N, 1, 0)
ham_ndo = matrix_representation(symham, H_ndo; projection = true)
3420×3420 SparseArrays.SparseMatrixCSC{Int64, Int64} with 17442 stored entries:
⎡⢻⣶⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠙⢿⣷⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠻⡿⣯⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢮⡻⣮⡳⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢮⠻⣦⡙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠳⣌⠻⣦⡉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⡈⠻⣦⡈⠓⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠈⠻⣦⡀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠈⠻⣦⡀⠉⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠈⠻⣦⡀⠈⠳⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠈⠻⣦⡀⠀⠙⢦⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡄⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢦⡀⠀⠉⠻⣦⡀⠀⠈⠓⢦⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠈⠻⣦⡀⠀⠀⠙⢦⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠈⢿⣷⡀⠀⠀⠙⢦⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠈⠻⣦⣄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠙⢿⣷⎦

More advanced use of conserved quantities: Fractionalized hilbert space

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(N,t,Jz)
    @fermions c
    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 for this model, we first use the same conservation as above to restrict to no double occupation and conserve spin. This gives a conserved quantum numbers which is a product of number conservations and so it is efficient in generating states.

N = 12
N_up = 4
N_down = 1
labels = Base.product(1:N, (:↑,:↓))
qn_spin = number_conservation(N_up, label -> :↑ in label) * number_conservation(N_down, label -> :↓ in label)
qn = qn_spin * qn_no_double_occ
H = hilbert_space(labels, qn)
symham = tjz(N, 1, 1/4)
ham = matrix_representation(symham, H; projection = true)
3960×3960 SparseArrays.SparseMatrixCSC{Float64, Int64} with 26780 stored entries:
⎡⢻⣶⡄⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⢭⡻⣮⣳⠀⢠⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠙⠚⢿⣷⡀⠙⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠲⣄⠈⠻⣦⣄⠙⠂⠀⠐⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠢⣄⠙⠿⣧⡳⣄⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠀⠙⢮⡻⣮⣄⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠰⣄⠀⠀⠀⠙⠻⣦⣀⠀⠀⠀⠙⠂⠀⠀⠀⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠘⢿⣷⡲⣄⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠘⢮⡻⣮⣅⠳⣄⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⢥⡙⠻⣦⡈⠳⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡈⢿⣷⣦⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⣄⠀⠀⠀⠀⠀⠀⠈⠛⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠙⠦⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠈⢿⣷⡤⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⢯⡻⣮⣓⠰⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⢙⡘⠿⣧⡈⠳⡄⢀⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠙⢦⡈⢻⣶⣤⡀⠙⢦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠀⠀⠀⠀⠀⢉⠀⠻⠿⣧⡀⠀⠙⢦⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⢻⣶⣄⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠹⣿⣿⣦⡀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⠻⣦⎦

This space fragments into more sectors, which are labelled by the order of occupied spins. By defining a function that maps states to the spin order, we can split the hilbert space into those fragments. The function can be defined as

function spin_order(state, H)
    N = length(keys(H)) ÷ 2
    occupations = map(1:N) do n
        if FermionicHilbertSpaces.occupation(state, (n, :↑), H)
            :↑
        elseif FermionicHilbertSpaces.occupation(state, (n, :↓), H)
            :↓
        else
            :hole
        end
    end
    filter(x -> x != :hole, occupations)
end
spin_order (generic function with 1 method)

We can then construct the fractionalized hilbert space as

qnfrac = Base.Fix2(spin_order, H)
Hfrac = FermionicHilbertSpaces.symmetrize(H, qnfrac)
Hfrac.symmetry.qntofockstates
5-element Dictionaries.Dictionary{Vector{Symbol}, SubArray{FockNumber{UInt64}, 1, Vector{FockNumber{UInt64}}, Tuple{UnitRange{Int64}}, true}}:
 [:↑, :↑, :↑, :↑, :↓] │ FockNumber{UInt64}[FockNumber(8390528), FockNumber(8390…
 [:↑, :↑, :↑, :↓, :↑] │ FockNumber{UInt64}[FockNumber(4197248), FockNumber(4197…
 [:↑, :↑, :↓, :↑, :↑] │ FockNumber{UInt64}[FockNumber(2100608), FockNumber(2100…
 [:↑, :↓, :↑, :↑, :↑] │ FockNumber{UInt64}[FockNumber(1052288), FockNumber(1052…
 [:↓, :↑, :↑, :↑, :↑] │ FockNumber{UInt64}[FockNumber(528128), FockNumber(26598…

This splits the space of dimension 3960 into 5 fragments each of size 792.

Let's find the half-chain entanglement entropy of the ground state in each sector. We can iterate over sectors by calling sectors(H), and we can use subregion to find the hilbert space of a subsystem. 

using Arpack, LinearAlgebra
using FermionicHilbertSpaces: sectors
left_labels = Base.product(1:(N÷2), (:↑,:↓))
entropy(ρ) = sum(-λ * log(λ) for λ in eigvals(Hermitian(ρ)) if λ > 1e-12)
map(sectors(Hfrac)) do Hsec
    ham = matrix_representation(symham, Hsec; projection = true)
    vals, vecs = eigs(ham, nev = 1)
    Hleft = subregion(left_labels, Hsec)
    rholeft = partial_trace(vecs[:,1]*vecs[:,1]', Hsec => Hleft)
    entropy(rholeft)
end
5-element Vector{Float64}:
 0.7764915118824811
 0.7577846594935103
 0.7671196796227888
 0.7577846594935123
 0.7764915118824796