Constrained Hilbert spaces and conserved quantum numbers

Conservation laws let you work in a small part of an exponentially large Hilbert space. This page shows how a constraint selects basis states and groups them into sectors of fixed quantum numbers, lists the available constraints, and ends with worked examples where the full space would be far too large to build.

Constraints and sectors

A constraint selects basis states of a space. Constraints for conserved quantities also group the selected states into sectors, one for each value of the quantum number. A constraint is given when building a space, either to hilbert_space or as the constraint keyword of tensor_product, and constrain_space applies one to an existing space:

using FermionicHilbertSpaces, LinearAlgebra@fermions fHleft = hilbert_space(f, 1:2)Hright = hilbert_space(f, 3:4)Hone = tensor_product(Hleft, Hright; constraint = NumberConservation(1))Heven = constrain_space(hilbert_space(f, 1:4), ParityConservation(1))H = hilbert_space(f, 1:4, NumberConservation())

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

The sector labels of a space are its quantum numbers. A sector is a space of its own, and an operator that conserves the quantity can be represented directly on it. embed places a state on a sector back into the full space, where indices locates the sector's states:

quantumnumbers(H)H2 = sector(2, H)hopping = sum(f[k]'f[k+1] + hc for k in 1:3)ψ = eigvecs(Hermitian(representation(hopping, H2, :dense)))[:, 1]@assert embed(ψ, H2 => H)[indices(2, H)] ≈ ψindices(2, H)

6:11

Passing a constraint while building a space prunes states as they are generated, so the full space is never built; only FilterConstraint and SectorConstraint test the full space afterwards. constrain_space always filters a space that already exists.

Available constraints

NumberConservation

Conserves the particle number, or a weighted sum of the particle numbers of selected spaces. The allowed total is one value or a collection; leaving it out keeps all sectors.

@boson bHf = hilbert_space(f, 1:3)Hb = hilbert_space(b, 10)# fermions in Hf minus bosons in Hb is 0 or 1tensor_product(Hf, Hb; constraint = NumberConservation(0:1, [Hf, Hb], [1, -1]))

15-dimensional SectorHilbertSpace
Parent: ProductSpace(80-dim, 2 factors)
Sectors: [0 (8-dim), 1 (7-dim)]

ParityConservation

Conserves fermion parity, 1 for even and -1 for odd states; without arguments it keeps both sectors.

quantumnumbers(hilbert_space(f, 1:4, ParityConservation()))

2-element Vector{Int64}:
 -1
  1

Products of constraints

Multiplying constraints imposes all of them, and the quantum number of a sector has one entry per factor. Here each pair of modes holds at most one fermion, and there are two fermions in total:

pairs = [hilbert_space(f, 2k-1:2k) for k in 1:3]pair_constraint = prod(NumberConservation(0:1, Hk) for Hk in pairs) * NumberConservation(2)tensor_product(pairs; constraint = pair_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)]

AdditiveConstraint

Conserves a sum of contributions from the factors of a tensor product, given by maps applied to the state of each factor. This covers additive quantum numbers that are not particle numbers, such as the magnetization of spins:

@spins S 1//2spins = [hilbert_space(S[k]) for k in 1:4]Hmag = tensor_product(spins; constraint = AdditiveConstraint(-2:2; maps = s -> s.m))quantumnumbers(Hmag)

5-element Vector{HalfIntegers.Half{Int64}}:
 -2
 -1
  0
  1
  2

FilterConstraint

Keeps the states for which a function returns true. With maps and spaces, the function receives the value of each map on the state of its space. Here the particle numbers of two bosonic modes must be sorted:

import FermionicHilbertSpaces: particle_number@bosons dHbb = hilbert_space(d, 1:2, 3)sorted = FilterConstraint(issorted; maps = particle_number, spaces = factors(Hbb))basisstates(constrain_space(Hbb, sorted))

6-element basis of ConstrainedSpace(ProductSpace(9-dim, 2 factors), 6-dim):
 (0, 0)
 (0, 1)
 (1, 1)
 (0, 2)
 (1, 2)
 (2, 2)

SectorConstraint

Like FilterConstraint, but the function returns a sector label, or missing to discard the state. Here the two modes must hold different numbers of bosons, and the sectors are labelled by the total:

distinct = SectorConstraint(n -> allunique(n) ? sum(n) : missing;    maps = particle_number, spaces = factors(Hbb))Hdistinct = constrain_space(Hbb, distinct)map((qn, Hsec) -> qn => basisstates(Hsec), quantumnumbers(Hdistinct), sectors(Hdistinct))

3-element Vector{Pair{Int64, FermionicHilbertSpaces.BasisStates{ProductState{Tuple{BosonicState, BosonicState}}, FermionicHilbertSpaces.ConstrainedSpace{ProductState{Tuple{BosonicState, BosonicState}}, FermionicHilbertSpaces.ProductSpace{ProductState{Tuple{BosonicState, BosonicState}}, Tuple{FermionicHilbertSpaces.TruncatedBosonicHilbertSpace{Int64, BosonField{FermionicHilbertSpaces.Tags{Nothing}}}, FermionicHilbertSpaces.TruncatedBosonicHilbertSpace{Int64, BosonField{FermionicHilbertSpaces.Tags{Nothing}}}}, FermionicHilbertSpaces.TruncatedBosonicHilbertSpace{Int64, BosonField{FermionicHilbertSpaces.Tags{Nothing}}}, Int64, LinearIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}, Vector{ProductState{Tuple{BosonicState, BosonicState}}}, FermionicHilbertSpaces.IndexTable}, Vector{ProductState{Tuple{BosonicState, BosonicState}}}}}}:
 1 => [(1, 0), (0, 1)]
 2 => [(2, 0), (0, 2)]
 3 => [(2, 1), (1, 2)]

BranchConstraint

Builds the states one factor at a time and calls f(partial_state, depth, spaces) after each factor is assigned. Returning false discards every state that starts this way, so a sector much smaller than the full space is found without enumerating it. Here no two neighbouring modes are both occupied:

no_neighbours = BranchConstraint((p, depth, spaces) -> depth == 1 || count_ones(p[depth]) + count_ones(p[depth-1]) < 2)dim(hilbert_space(f, 1:20, no_neighbours))

Examples

Spinful fermions

The package does not special-case spinful fermions: spin is an extra label. If the numbers of spin up and spin down fermions are separately conserved, the product of their number conservations divides the space into sectors labelled by both numbers:

using FermionicHilbertSpaces@fermions fHup = hilbert_space(f, [(i, :↑) for i in 1:4])Hdn = hilbert_space(f, [(i, :↓) for i in 1:4])H = tensor_product(Hup, Hdn)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)]

The sector with two spin up and no spin down fermions can be picked out of this space, or built directly, which is cheaper. Building it in the tensor product is cheaper still, since it never enumerates the full space:

constraint = NumberConservation(2, Hup) * NumberConservation(0, Hdn)@assert basisstates(sector([2, 0], Hsectors)) == basisstates(constrain_space(H, constraint))tensor_product(Hup, Hdn; constraint)

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

Hubbard model

The Hubbard model

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

conserves the number of spin up and spin down fermions separately. On 20 sites the full space has dimension $4^{20} \approx 10^{12}$, but the sector with two spin up fermions and one spin down fermion is small, and is generated directly:

using FermionicHilbertSpacesfunction 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@fermions fN = 20Hup = hilbert_space(f, [(i, :↑) for i in 1:N], NumberConservation(2))Hdn = hilbert_space(f, [(i, :↓) for i in 1:N], NumberConservation(1))H = tensor_product(Hup, Hdn)ham = representation(hubbard_hamiltonian(f, N, 1.0, 4.0), H);size(ham)

(3800, 3800)

No double occupation

For strong interactions, double occupation can be forbidden by allowing at most one fermion on each site. The hamiltonian maps some states of this subspace outside it, and projection = true drops those terms:

local_modes(k) = [f[k, :↑], f[k, :↓]]no_double_occ = prod(NumberConservation(0:1, local_modes(k)) for k in 1:N)H_ndo = constrain_space(H, no_double_occ)ham_ndo = representation(hubbard_hamiltonian(f, N, 1, 0), H_ndo; projection = true);size(ham_ndo)

(3420, 3420)

Tavis–Cummings model

$N$ spin-1/2 emitters coupled to a cavity mode,

\[H = \omega a^\dagger a + \omega_0 \sum_i S^z_i + g \sum_i \left(a^\dagger S^-_i + a S^+_i\right),\]

conserve the number of excitations, the boson number plus $\sum_i (S^z_i + 1/2)$. Spins have no particle number, so an AdditiveConstraint counts their excitations with one map per factor:

using FermionicHilbertSpacesimport FermionicHilbertSpaces: particle_numberN = 3nmax = 4@spins S 1//2@boson aspaces = [hilbert_space(a, nmax + 1); [hilbert_space(S[i]) for i in 1:N]]excitations = AdditiveConstraint(0:nmax+N; spaces, maps = [particle_number; fill(s -> Int(s.m + 1//2), N)])H = tensor_product(spaces; constraint = excitations)

40-dimensional SectorHilbertSpace
Parent: ProductSpace(40-dim, 4 factors)
Sectors: 8 total [0 (1-dim), ..., 7 (1-dim)]

The Hamiltonian is block diagonal in these sectors, one per excitation number.

Fractionalized Hilbert space of the t-Jz model

The $t-J_z$ model

\[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},\]

with $S^z_i = n_{i,↑} - n_{i,↓}$ and no double occupation, has a Hilbert space that splits into exponentially many dynamically disconnected sectors^[1]. Each sector is labelled by a spin pattern such as ↑↑↓↑↓: its states have the spins of the occupied sites in exactly that order, with holes in between. A BranchConstraint builds one such sector without the full space:

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# partials holds the fock numbers of the 2N modes, of which the first `depth` are assignedfunction valid_state(partials, depth, spaces, expected_order)    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            site == prev_site && return false # the spaces are ordered by site            prev_site = site            expected_order[exc_count] == spin || return false        end    end    depth == length(spaces) && return exc_count == length(expected_order)    return trueendspin_order_constraint(order) = BranchConstraint((p, d, s) -> valid_state(p, d, s, order))N = 16labels = [(i, s) for i in 1:N for s in (:↑,:↓)]Hfrac = hilbert_space(f, labels, spin_order_constraint([:↑, :↑, :↓, :↑, :↓]))ham = representation(tjz(f, N, 1, 1/4), Hfrac; projection = true);dim(Hfrac)

The full space of these modes has dimension $2^{32}$. subregion gives the space of the left half compatible with the constraint, which is what partial_trace(ρ, Hfrac => Hsub_frac) needs:

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

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

1T. Rakovszky, P. Sala, R. Verresen, M. Knap and F. Pollmann, "Statistical localization: from strong fragmentation to strong edge modes", Phys. Rev. B 101, 125126 (2020), arXiv:1910.06341.