cvsvensson · cvsvensson · Apr 8, 2026 · Mar 9, 2026 · Mar 10, 2026 · Mar 11, 2026
diff --git a/Project.toml b/Project.toml
@@ -1,33 +1,27 @@
 name = "FermionicHilbertSpaces"
 uuid = "82e5dd95-83c3-46b2-a303-acc930c81756"
-version = "0.9.3"
+version = "0.10.0"
 authors = ["Viktor Svensson", "William Samuelson"]
 
 [deps]
 BitPermutations = "81206824-6155-40a8-95f7-18696067be4d"
-Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
-FlexiGroups = "1e56b746-2900-429a-8028-5ec1f00612ec"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonCommutativeProducts = "a9a0a982-c603-4485-bf80-6b61b825abb4"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
-UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [compat]
 BitPermutations = "0.3"
-Dictionaries = "0.4.5"
 FillArrays = "1"
-FlexiGroups = "0.1"
 LinearAlgebra = "1.11"
-NonCommutativeProducts = "0.3"
+NonCommutativeProducts = "0.4.1"
 OrderedCollections = "1"
 PrecompileTools = "1.2"
 SparseArrays = "1"
 TestItems = "1.0"
 TupleTools = "1.6"
-UUIDs = "1.11"
 julia = "1.10"
diff --git a/README.md b/README.md
@@ -40,7 +40,7 @@ entanglement_entropy = sum(λ -> -λ * log(λ), eigvals(ρsub))
 The hamiltonian above conserves the number of fermions, which we can exploit as
 
 ````julia
-Hcons = hilbert_space(1:4, number_conservation(2))
+Hcons = hilbert_space(1:4, NumberConservation(2))
 ````
 
 ````

diff --git a/docs/make.jl b/docs/make.jl
@@ -5,7 +5,7 @@ using Literate
 
 DocMeta.setdocmeta!(FermionicHilbertSpaces, :DocTestSetup, :(using FermionicHilbertSpaces); recursive=true)
 
-literate_files = ["examples/kitaev_chain.jl", "examples/free_fermions.jl"]
+literate_files = ["examples/kitaev_chain.jl", "examples/free_fermions.jl", "examples/floquet_tutorial.jl", "examples/open_system_lindblad.jl"]
 output_directory = "docs/src/literate_output"
 for file in literate_files
     Literate.markdown(file, output_directory; documenter=true, execute=false)
@@ -24,9 +24,13 @@ makedocs(;
         "Home" => "index.md",
         "Conserved quantities" => "conservation.md",
         "Non interacting systems" => "non_interacting.md",
+        "Tutorials" => [       
+                 "Defining your own algebra: Floquet" => "literate_output/floquet_tutorial.md",
+        ],
         "Examples" => [
             "Interacting Kitaev chain" => "literate_output/kitaev_chain.md",
-            "Free fermions" => "literate_output/free_fermions.md"
+            "Free fermions" => "literate_output/free_fermions.md",
+            "Open systems" => "literate_output/open_system_lindblad.md"
         ],
         "Misc" => "misc.md",
         "Functions" => "docstrings.md",

diff --git a/docs/src/conservation.md b/docs/src/conservation.md
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -13,17 +13,14 @@ CurrentModule = FermionicHilbertSpaces
 The following example demonstrates how to define a fermionic Hilbert space and construct a simple Hamiltonian.
 ```@example intro
 using FermionicHilbertSpaces, LinearAlgebra
-N = 2 # number of fermions
-spatial_labels = 1:N 
-internal_labels = (:↑,:↓)
-labels = Base.product(spatial_labels, internal_labels) 
-H = hilbert_space(labels) 
+@fermions c
+labels = [(i, σ) for i in 1:2 for σ in (:↑,:↓)]
+H = hilbert_space(c, labels) 
 ```
 We now have a Hilbert space representing 2N fermions. To construct operators on this space we first call `@fermions c` which makes `c` represent symbolic fermions. Indexing into `c` returns an operator representing an annihilation operator e.g. `c[1,:↑]`. The creation operator is given by the adjoint `c[1,:↑]'`. These can be multiplied and added together. Here is a simple Hamiltonian with hopping and Coulomb interaction.
 ```@example intro
-@fermions c
 hopping = c[1,:↑]'c[2,:↑] + c[1,:↓]'c[2,:↓] + hc 
-coulomb = sum(c[n,:↑]'c[n,:↑]c[n,:↓]'c[n,:↓] for n in spatial_labels)
+coulomb = sum(c[n,:↑]'c[n,:↑]c[n,:↓]'c[n,:↓] for n in 1:2)
 ham = hopping + coulomb
 ```
 To get the matrix representation of this operator on the Hilbert space, do
@@ -34,8 +31,8 @@ mat = matrix_representation(ham, H)
 ## Tensor product and partial trace
 This package also includes functionality for combining Hilbert spaces and operators on them, and taking partial traces, in a way that is consistent with fermionic anticommutation relations. 
 ```@example intro
-H1 = hilbert_space(1:2)
-H2 = hilbert_space(3:4)
+H1 = hilbert_space(c, 1:2)
+H2 = hilbert_space(c, 3:4)
 H = tensor_product(H1, H2)
 c1 = matrix_representation(c[1], H1)
 c3 = matrix_representation(c[3], H2)
@@ -47,25 +44,25 @@ tensor_product((c1, c3), (H1, H2) => H) ≈ c1c3
 ```
 Let's partial trace to sites 1 and 3. Let's get a Hilbert space for those sites by using the function `subregion`, and then we can partial trace to that space with `partial_trace`.
 ```@example intro
-Hsub = subregion([1,3], H)
+Hsub = subregion([c[1], c[3]], H)
 dim(Hsub) / dim(H) * partial_trace(c1c3, H => Hsub) ≈ matrix_representation(c[1] * c[3], Hsub)
 ```
 
 ## Conserved quantum numbers
 This package also includes some functionality for working with conserved quantum numbers. If we have for example number conservation, we might want to get a block structure of the hamiltonian. Here's how one can do that:
 ```@example intro
-H = hilbert_space(labels, number_conservation())
+H = hilbert_space(c, labels, NumberConservation())
 matrix_representation(ham, H)
 ```
 This has a block structure corresponding to the different sectors. To only look at some sectors, for example the sectors with 0, 2 and 4 particles, use
 ```@example intro
-H = hilbert_space(labels, number_conservation([0, 2, 4]))
+H = hilbert_space(c, labels, NumberConservation([0, 2, 4]))
 matrix_representation(ham, H)
 ```
 
 Those sectors have even fermionic parity, which can alternatively be specified with `ParityConservation`.
 ```@example intro
-H = hilbert_space(labels, ParityConservation(1))
+H = hilbert_space(c, labels, ParityConservation(1))
 matrix_representation(ham, H)
 ```
 

diff --git a/docs/src/misc.md b/docs/src/misc.md
@@ -2,16 +2,22 @@
 # Misc
 
 ## Subregion
-The function subregion can be used to extract the Hilbert space of a subregion.
+When the Hilbert space has a restricted set of fock states, the Hilbert space of the subregion will only include fock states compatible with this restriction. In the example below, the total number of particles 1, and the subregion will have three possible states: (1,0), (0,1), and (0,0).
 ```@example subregion
 using FermionicHilbertSpaces
-H = hilbert_space(1:4)
-Hsub = subregion(1:2, H)
+@fermions f
+H = hilbert_space(f, 1:4, NumberConservation(1))
+Hsub = subregion(hilbert_space(f, 1:2), H)
+basisstates(Hsub)
 ``` 
 
-When the Hilbert space has a restricted set of fock states, the Hilbert space of the subregion will only include fock states compatible with this restriction. In the example below, the total number of particles 1, and the subregion will have three possible states: (1,0), (0,1), and (0,0).
-```@example subregion
-H = hilbert_space(1:4, number_conservation(1))
-Hsub = subregion(1:2, H)
-basisstates(Hsub)
-``` 
+## State mapper interface
+
+Internal tensor/reshape/partial-trace routines use a common mapper protocol:
+
+- `state_mapper(H, Hs)` returns a mapper object.
+- `split_state(state, mapper)` returns a tuple with one entry per target subsystem.
+- Each tuple entry is a weighted collection `((substate, weight), ...)`.
+- `combine_states(substates, mapper)` returns a weighted collection `((state, weight), ...)`.
+
+This package does not require a single concrete container type for weighted collections; callers should treat them as iterable collections of `(state, weight)` outcomes.
diff --git a/docs/src/non_interacting.md b/docs/src/non_interacting.md
@@ -14,9 +14,9 @@ If ``H`` is hermitian so is ``h`` and ``E`` is real, but we don't have to restri
 One can manually define the hilbert space using only the single particle states as
 ```@example single_particle_hilbert_space
 using FermionicHilbertSpaces, LinearAlgebra
-N = 2
-H = hilbert_space(1:N, FermionicHilbertSpaces.SingleParticleState.(1:N))
 @fermions c
+N = 2
+H = hilbert_space(c, 1:N, FermionicHilbertSpaces.SingleParticleState.(1:N))
 h = rand(ComplexF64, N, N)
 E = rand(ComplexF64)
 op = E * I + sum(c[i]' * h[i, j] * c[j] for i in 1:N, j in 1:N)
@@ -27,7 +27,7 @@ Often, $h_{nm}$ is of interest because diagonalizing it gives information on the
 !!! tip "Use `single_particle_hilbert_space` instead"
     For convenience, `single_particle_hilbert_space` can be used define the hilbert space which will give only the single particle states, and will remove the contribution of the identity operator when calling `matrix_representation`:
     ```julia
-    H = single_particle_hilbert_space(1:N)
+    H = single_particle_hilbert_space(c, 1:N)
     matrix_representation(op,H) == h # true
     ```
 
@@ -53,9 +53,9 @@ We can get this representation manually as follows:
 ```@example bdg_particle_hilbert_space
 using FermionicHilbertSpaces, LinearAlgebra
 N = 2
-states = [FermionicHilbertSpaces.NambuState(n, hole) for (n, hole) in Base.product(1:N, (true, false))]
-H = hilbert_space(1:N, states)
+states = vec([FermionicHilbertSpaces.NambuState(n, hole) for (n, hole) in Base.product(1:N, (true, false))])
 @fermions c
+H = hilbert_space(c, 1:N, states)
 h = Hermitian(rand(N, N))
 Δ = rand(N, N) |> m -> m - transpose(m)
 E = rand()
@@ -67,7 +67,7 @@ The matrix returned by `matrix_representation` will depend on the ordering of op
 !!! tip "Use `bdg_hilbert_space` instead"
     By defining the hilbert space using `bdg_hilbert_space`, one automatically gets the Nambu states and `matrix_representation` will return a matrix of the form above without the need to manually convert it:
     ```julia
-    Hbdg = bdg_hilbert_space(1:N)
+    Hbdg = bdg_hilbert_space(c, 1:N)
     matrix_representation(op, Hbdg)
     #= example output
     4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries:

diff --git a/examples/bose_hubbard.jl b/examples/bose_hubbard.jl
@@ -0,0 +1,38 @@
+using FermionicHilbertSpaces
+using Arpack, LinearAlgebra
+
+N = 6
+local_dimension = 4
+total_particles = N
+@bosons b
+
+H = hilbert_space(b, 1:6, local_dimension, NumberConservation(total_particles))
+
+number_ops = [matrix_representation(b[i]'b[i], H) for i in 1:N]
+hopping_ops = [matrix_representation(b[i]'b[i+1], H) for i in 1:(N-1)]
+
+function bose_hubbard_observables(U, t)
+    ham = -t * sum(b[i]'b[i+1] + hc for i in 1:(N-1)) +
+          U * sum(b[i]'b[i] * (b[i]'b[i] - 1) for i in 1:N)
+    M = matrix_representation(ham, H)
+    _, vecs = eigs(M; nev=1, which=:SR)
+    ψ = vecs[:, 1]
+    occupation = sum(dot(ψ, ni, ψ) for ni in number_ops) / N
+    fluctuation = sum(dot(ψ, ni^2, ψ) - dot(ψ, ni, ψ)^2 for ni in number_ops) / N
+    coherence = sum(dot(ψ, hij, ψ) for hij in hopping_ops) / (N - 1)
+    return (; occupation, fluctuation, coherence)
+end
+
+##
+t = 1.0
+Us = range(0, 20, length=10)
+data = [bose_hubbard_observables(U, t) for U in Us]
+##
+occs = getproperty.(data, :occupation)
+flucts = getproperty.(data, :fluctuation)
+cohs = getproperty.(data, :coherence)
+using Plots
+plot(Us, flucts, xlabel="U", label="Average Fluctuations")
+plot!(Us, occs, label="Average Occupation")
+plot!(Us, cohs, label="Average Coherence")
+
diff --git a/examples/braiding.jl b/examples/braiding.jl
@@ -1,8 +1,8 @@
-using FermionicHilbertSpaces, LinearAlgebra, Plots, OrdinaryDiffEqTsit5
-using Symbolics
+using FermionicHilbertSpaces
+using Symbolics, LinearAlgebra, Plots, OrdinaryDiffEqTsit5
 
-H = majorana_hilbert_space([0, 1, 2, 3, 22, 33], ParityConservation())
 @majoranas γ
+H = hilbert_space(γ, [0, 1, 2, 3, 22, 33], ParityConservation())
 @variables Δ[1:3]::Real
 symbolic_ham = sum(1im * Δ[i] * γ[0] * γ[i] for i in 1:3)
 ham = collect(matrix_representation(symbolic_ham, H))
@@ -58,7 +58,7 @@ isapprox(abs2(sol[end]' * exchange_gate^2 * u0), 1, atol=1e-2)
 isapprox(abs2(sol(T)' * exchange_gate * u0), 1, atol=1e-2)
 ## lets measure the parities
 measurements = [matrix_representation(1.0im * γ[i] * γ[j], H) for (i, j) in [(2, 22), (3, 33)]]
-plot(ts, [real(sol(t)'m * sol(t)) for m in measurements2, t in ts]', xlabel="t", label=["P2" "P3"], frame=:box, size=(400, 250), lw=2)
+plot(ts, [real(sol(t)'m * sol(t)) for m in measurements, t in ts]', xlabel="t", label=["P2" "P3"], frame=:box, size=(400, 250), lw=2)
 
 ## Let's calculate the Non-abelian berry pase with the Kato method
 const totalparity = parityoperator(H)
@@ -83,7 +83,7 @@ U0 = Matrix{ComplexF64}(I, 8, 8)
 kato_prob = ODEProblem(kato_ode!, U0, tspan, p)
 ts = range(0, tspan[2], 100)
 ## Solve the ODE and check the norm
-@time kato_sol = solve(kato_prob, Tsit5(); saveat=[0, T, 2T], reltol=1e-10, abstol=1e-10); #saveat=ts, reltol=1e-4)
+@time kato_sol = solve(kato_prob, Tsit5(); saveat=[0, T, 2T], reltol=1e-10, abstol=1e-10);
 kato_sol[end]' * kato_sol[end] ≈ I
 1 ≈ dot(kato_sol[end], exchange_gate^2) / (dot(exchange_gate, exchange_gate)) |> abs
 1 ≈ dot(kato_sol(T), exchange_gate) / (dot(exchange_gate, exchange_gate)) |> abs

diff --git a/examples/cavity_DQD.jl b/examples/cavity_DQD.jl
@@ -0,0 +1,60 @@
+using FermionicHilbertSpaces
+using LinearAlgebra, Plots
+
+@fermions f
+@boson a
+
+Ha = hilbert_space(a, 8)
+Hf = hilbert_space(f, [(j, σ) for j in 1:2 for σ in [:↑, :↓]], NumberConservation(1))
+H = tensor_product(Hf, Ha)
+n_photon = matrix_representation(a'a, H)
+n_f(j, σ) = f[j, σ]'f[j, σ]
+
+h_dqd(; δϵ, t, α, U, B) = sum(δϵ * (n_f(1, σ) - n_f(2, σ)) for σ in [:↑, :↓]) +
+                        B * sum(n_f(j, :↑) - n_f(j, :↓) for j in 1:2) +
+                        t * sum(f[1, σ]'f[2, σ] + hc for σ in [:↑, :↓]) +
+                        t * α * (f[1, :↓]'f[2, :↑] - f[1, :↑]'f[2, :↓] + hc) +
+                        U * sum(n_f(j, :↑) * n_f(j, :↓) for j in 1:2)
+
+h_cavity_dqd(; ωr, δϵ, t, α, U, B, g) = ωr * a' * a + h_dqd(; δϵ, t, α, U, B) +
+                                            g * (a' + a) * sum(n_f(1, σ) - n_f(2, σ) for σ in [:↑, :↓])
+
+ωr = 1.0
+t  = 0.2
+α = 0.2
+U  = 0.0
+g  = 0.1
+B = 0.2
+fix_params = (; ωr, t, α, U, B, g)
+
+function cavity_dqd_eigensystem(; ωr, δϵ, t, α, U, B, g)
+    Hsym = h_cavity_dqd(; ωr, δϵ, t, α, U, B, g)
+    M = matrix_representation(Hsym, H)
+    vals, vecs = eigen(Hermitian(Matrix(M)))
+    return vals, vecs
+end
+
+function branch_data(; ωr, δϵ, t, α, U, B, g, branches=3:6)
+    vals, vecs = cavity_dqd_eigensystem(; ωr, δϵ, t, α, U, B, g)
+    splittings = vals[branches] .- vals[1]
+    photon_numbers = [dot(ψ, n_photon, ψ) for ψ in eachcol(vecs[:, branches])]
+    return (; splittings, photon_numbers)
+end
+
+δϵ_res = sqrt(ωr^2/4 - t^2)
+δϵs = range(-3g, 3g, length=50)
+data = [branch_data(;ωr, δϵ=δϵ_res + δϵ, t, α, U, B, g) for δϵ in δϵs]
+splittings = reduce(hcat, getproperty.(data, :splittings))'
+photon_numbers = reduce(hcat, getproperty.(data, :photon_numbers))'
+
+plot(
+    δϵ_res .+ δϵs,
+    splittings,
+    label = "",
+    line_z = photon_numbers,
+    color = :viridis,
+    xlabel = "δϵ",
+    ylabel = "Energy splitting",
+    linewidth = 3,
+    colorbar_title = "Photon number",
+)
diff --git a/examples/dicke_model.jl b/examples/dicke_model.jl
@@ -0,0 +1,22 @@
+using FermionicHilbertSpaces
+using FermionicHilbertSpaces: FilterConstraint, parity
+
+N = 3
+@spins S
+@boson a
+
+ωc = 1.0
+ωz = 0.5
+λ = 0.1
+symham = ωc * a' * a + ωz * sum(S[k][:z] for k in 1:N) + 2λ / sqrt(N) * (a' + a) * sum(S[k][:x] for k in 1:N)
+
+Hs = hilbert_space(S, 1:N, 1 // 2)
+Ha = hilbert_space(a, 10)
+H = tensor_product(Hs, Ha)
+spin_parity(p) = (-1)^(Int(N // 2 + sum(s -> s.m, p.states)))
+constraint = FilterConstraint([Ha, Hs], [parity, spin_parity], ==(1) ∘ prod)
+H = tensor_product(Hs, Ha; constraint=constraint)
+H2 = constrain_space(tensor_product(Hs, Ha), constraint)
+H == H2
+ham = matrix_representation(symham, H)
+# One should use sectors with permutationally invariant states, but we need to add methods to constrain product spaces to do that