diff --git a/NDTensors/src/imports.jl b/NDTensors/src/imports.jl
index 04ea3d8d6b..73d8ae3771 100644
--- a/NDTensors/src/imports.jl
+++ b/NDTensors/src/imports.jl
@@ -33,9 +33,9 @@ for lib in [
   :MetalExtensions,
   :BroadcastMapConversion,
   :RankFactorization,
-  :Sectors,
   :LabelledNumbers,
   :GradedAxes,
+  :SymmetrySectors,
   :TensorAlgebra,
   :SparseArrayInterface,
   :SparseArrayDOKs,
diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml
deleted file mode 100644
index 9b1d5ccd25..0000000000
--- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml
+++ /dev/null
@@ -1,2 +0,0 @@
-[deps]
-NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl
deleted file mode 100644
index aa3056438e..0000000000
--- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-module GradedAxesSectorsExt
-using ..GradedAxes: GradedAxes
-using ...Sectors: Sectors, AbstractCategory, ⊗ # , dual
-
-GradedAxes.fuse_labels(c1::AbstractCategory, c2::AbstractCategory) = only(c1 ⊗ c2)
-
-# TODO: Decide the fate of `dual`.
-## GradedAxes.dual(c::AbstractCategory) = dual(c)
-end
diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml
deleted file mode 100644
index ef491a529c..0000000000
--- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml
+++ /dev/null
@@ -1,3 +0,0 @@
-[deps]
-NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl
deleted file mode 100644
index 371e7e57cd..0000000000
--- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-@eval module $(gensym())
-using NDTensors.GradedAxes: dual, fuse_labels
-using NDTensors.Sectors: U1, Z
-using Test: @test, @testset
-
-@testset "GradedAxesSectorsExt" begin
-  @test fuse_labels(U1(1), U1(2)) == U1(3)
-  @test dual(U1(2)) == U1(-2)
-
-  @test fuse_labels(Z{2}(1), Z{2}(1)) == Z{2}(0)
-  @test fuse_labels(Z{2}(0), Z{2}(1)) == Z{2}(1)
-  @test dual(Z{2}(1)) == Z{2}(1)
-  @test dual(Z{2}(0)) == Z{2}(0)
-end
-end
diff --git a/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl
index 7756b71575..d605289f81 100644
--- a/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl
+++ b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl
@@ -1,8 +1,7 @@
 module GradedAxes
 include("blockedunitrange.jl")
 include("gradedunitrange.jl")
-include("fusion.jl")
 include("dual.jl")
 include("unitrangedual.jl")
-include("../ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl")
+include("fusion.jl")
 end
diff --git a/NDTensors/src/lib/GradedAxes/src/dual.jl b/NDTensors/src/lib/GradedAxes/src/dual.jl
index 985ebe33cc..ff11b1103c 100644
--- a/NDTensors/src/lib/GradedAxes/src/dual.jl
+++ b/NDTensors/src/lib/GradedAxes/src/dual.jl
@@ -5,3 +5,5 @@ using NDTensors.LabelledNumbers:
 label_dual(x) = label_dual(LabelledStyle(x), x)
 label_dual(::NotLabelled, x) = x
 label_dual(::IsLabelled, x) = labelled(unlabel(x), dual(label(x)))
+
+flip(g::AbstractGradedUnitRange) = dual(gradedrange(label_dual.(blocklengths(g))))
diff --git a/NDTensors/src/lib/GradedAxes/src/fusion.jl b/NDTensors/src/lib/GradedAxes/src/fusion.jl
index 5193449f2f..3c3b58bc91 100644
--- a/NDTensors/src/lib/GradedAxes/src/fusion.jl
+++ b/NDTensors/src/lib/GradedAxes/src/fusion.jl
@@ -1,10 +1,11 @@
-using BlockArrays: AbstractBlockedUnitRange
+using BlockArrays: AbstractBlockedUnitRange, blocklengths
 
 # Represents the range `1:1` or `Base.OneTo(1)`.
 struct OneToOne{T} <: AbstractUnitRange{T} end
 OneToOne() = OneToOne{Bool}()
 Base.first(a::OneToOne) = one(eltype(a))
 Base.last(a::OneToOne) = one(eltype(a))
+BlockArrays.blockaxes(g::OneToOne) = (Block.(g),)  # BlockArrays default crashes for OneToOne{Bool}
 
 # https://github.com/ITensor/ITensors.jl/blob/v0.3.57/NDTensors/src/lib/GradedAxes/src/tensor_product.jl
 # https://en.wikipedia.org/wiki/Tensor_product
@@ -18,23 +19,25 @@ function tensor_product(
   return foldl(tensor_product, (a1, a2, a3, a_rest...))
 end
 
+flip_dual(r::AbstractUnitRange) = r
+flip_dual(r::UnitRangeDual) = flip(r)
 function tensor_product(a1::AbstractUnitRange, a2::AbstractUnitRange)
-  return error("Not implemented yet.")
+  return tensor_product(flip_dual(a1), flip_dual(a2))
 end
 
 function tensor_product(a1::Base.OneTo, a2::Base.OneTo)
   return Base.OneTo(length(a1) * length(a2))
 end
 
-function tensor_product(a1::OneToOne, a2::AbstractUnitRange)
+function tensor_product(::OneToOne, a2::AbstractUnitRange)
   return a2
 end
 
-function tensor_product(a1::AbstractUnitRange, a2::OneToOne)
+function tensor_product(a1::AbstractUnitRange, ::OneToOne)
   return a1
 end
 
-function tensor_product(a1::OneToOne, a2::OneToOne)
+function tensor_product(::OneToOne, ::OneToOne)
   return OneToOne()
 end
 
@@ -45,27 +48,28 @@ function fuse_labels(x, y)
 end
 
 function fuse_blocklengths(x::Integer, y::Integer)
-  return x * y
+  # return blocked unit range to keep non-abelian interface
+  return blockedrange([x * y])
 end
 
 using ..LabelledNumbers: LabelledInteger, label, labelled, unlabel
 function fuse_blocklengths(x::LabelledInteger, y::LabelledInteger)
-  return labelled(unlabel(x) * unlabel(y), fuse_labels(label(x), label(y)))
+  # return blocked unit range to keep non-abelian interface
+  return blockedrange([labelled(x * y, fuse_labels(label(x), label(y)))])
 end
 
 using BlockArrays: blockedrange, blocks
 function tensor_product(a1::AbstractBlockedUnitRange, a2::AbstractBlockedUnitRange)
-  blocklengths = map(vec(collect(Iterators.product(blocks(a1), blocks(a2))))) do x
-    return mapreduce(length, fuse_blocklengths, x)
+  nested = map(Iterators.flatten((Iterators.product(blocks(a1), blocks(a2)),))) do it
+    return mapreduce(length, fuse_blocklengths, it)
   end
-  return blockedrange(blocklengths)
+  new_blocklengths = mapreduce(blocklengths, vcat, nested)
+  return blockedrange(new_blocklengths)
 end
 
-function blocksortperm(a::AbstractBlockedUnitRange)
-  # TODO: Figure out how to deal with dual sectors.
-  # TODO: `rev=isdual(a)`  may not be correct for symmetries beyond `U(1)`.
-  ## return Block.(sortperm(nondual_sectors(a); rev=isdual(a)))
-  return Block.(sortperm(blocklabels(a)))
+# convention: sort UnitRangeDual according to nondual blocks
+function blocksortperm(a::AbstractUnitRange)
+  return Block.(sortperm(blocklabels(nondual(a))))
 end
 
 using BlockArrays: Block, BlockVector
@@ -82,25 +86,34 @@ end
 # Used by `TensorAlgebra.splitdims` in `BlockSparseArraysGradedAxesExt`.
 # Get the permutation for sorting, then group by common elements.
 # groupsortperm([2, 1, 2, 3]) == [[2], [1, 3], [4]]
-function blockmergesortperm(a::AbstractBlockedUnitRange)
-  # If it is dual, reverse the sorting so the sectors
-  # end up sorted in the same way whether or not the space
-  # is dual.
-  # TODO: Figure out how to deal with dual sectors.
-  # TODO: `rev=isdual(a)`  may not be correct for symmetries beyond `U(1)`.
-  ## return Block.(groupsortperm(nondual_sectors(a); rev=isdual(a)))
-  return Block.(groupsortperm(blocklabels(a)))
+function blockmergesortperm(a::AbstractUnitRange)
+  return Block.(groupsortperm(blocklabels(nondual(a))))
 end
 
 # Used by `TensorAlgebra.splitdims` in `BlockSparseArraysGradedAxesExt`.
 invblockperm(a::Vector{<:Block{1}}) = Block.(invperm(Int.(a)))
 
-# Used by `TensorAlgebra.fusedims` in `BlockSparseArraysGradedAxesExt`.
-function blockmergesortperm(a::GradedUnitRange)
-  # If it is dual, reverse the sorting so the sectors
-  # end up sorted in the same way whether or not the space
-  # is dual.
-  # TODO: Figure out how to deal with dual sectors.
-  # TODO: `rev=isdual(a)`  may not be correct for symmetries beyond `U(1)`.
-  return Block.(groupsortperm(blocklabels(a)))
+function blockmergesort(g::AbstractGradedUnitRange)
+  glabels = blocklabels(g)
+  gblocklengths = blocklengths(g)
+  new_blocklengths = map(sort(unique(glabels))) do la
+    return labelled(sum(gblocklengths[findall(==(la), glabels)]; init=0), la)
+  end
+  return gradedrange(new_blocklengths)
+end
+
+blockmergesort(g::UnitRangeDual) = flip(blockmergesort(flip(g)))
+blockmergesort(g::AbstractUnitRange) = g
+
+# fusion_product produces a sorted, non-dual GradedUnitRange
+function fusion_product(g1, g2)
+  return blockmergesort(tensor_product(g1, g2))
+end
+
+fusion_product(g::AbstractUnitRange) = blockmergesort(g)
+fusion_product(g::UnitRangeDual) = fusion_product(flip(g))
+
+# recursive fusion_product. Simpler than reduce + fix type stability issues with reduce
+function fusion_product(g1, g2, g3...)
+  return fusion_product(fusion_product(g1, g2), g3...)
 end
diff --git a/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl
index 57e9420d88..243853514b 100644
--- a/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl
+++ b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl
@@ -12,13 +12,14 @@ using BlockArrays:
   blockedrange,
   BlockIndexRange,
   blockfirsts,
-  blocklasts,
+  blockisequal,
   blocklength,
   blocklengths,
   findblock,
   findblockindex,
   mortar
 using Compat: allequal
+using FillArrays: Fill
 using ..LabelledNumbers:
   LabelledNumbers, LabelledInteger, LabelledUnitRange, label, labelled, unlabel
 
@@ -37,6 +38,18 @@ function Base.OrdinalRange{T,T}(a::GradedOneTo{<:LabelledInteger{T}}) where {T}
   return unlabel_blocks(a)
 end
 
+# == is just a range comparison that ignores labels. Need dedicated function to check equality.
+struct NoLabel end
+blocklabels(r::AbstractUnitRange) = Fill(NoLabel(), blocklength(r))
+
+function labelled_isequal(a1::AbstractUnitRange, a2::AbstractUnitRange)
+  return blockisequal(a1, a2) && (blocklabels(a1) == blocklabels(a2))
+end
+
+function space_isequal(a1::AbstractUnitRange, a2::AbstractUnitRange)
+  return (isdual(a1) == isdual(a2)) && labelled_isequal(a1, a2)
+end
+
 # This is only needed in certain Julia versions below 1.10
 # (for example Julia 1.6).
 # TODO: Delete this once we drop Julia 1.6 support.
diff --git a/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl b/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl
index aa04cc1600..8e9529e2d0 100644
--- a/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl
+++ b/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl
@@ -6,7 +6,10 @@ UnitRangeDual(a::AbstractUnitRange) = UnitRangeDual{eltype(a),typeof(a)}(a)
 dual(a::AbstractUnitRange) = UnitRangeDual(a)
 nondual(a::UnitRangeDual) = a.nondual_unitrange
 dual(a::UnitRangeDual) = nondual(a)
+flip(a::UnitRangeDual) = dual(flip(nondual(a)))
 nondual(a::AbstractUnitRange) = a
+isdual(::AbstractUnitRange) = false
+isdual(::UnitRangeDual) = true
 ## TODO: Define this to instantiate a dual unit range.
 ## materialize_dual(a::UnitRangeDual) = materialize_dual(nondual(a))
 
@@ -16,6 +19,16 @@ Base.step(a::UnitRangeDual) = label_dual(step(nondual(a)))
 
 Base.view(a::UnitRangeDual, index::Block{1}) = a[index]
 
+function Base.show(io::IO, a::UnitRangeDual)
+  return print(io, UnitRangeDual, "(", blocklasts(a), ")")
+end
+
+function Base.show(io::IO, mimetype::MIME"text/plain", a::UnitRangeDual)
+  return Base.invoke(
+    show, Tuple{typeof(io),MIME"text/plain",AbstractArray}, io, mimetype, a
+  )
+end
+
 function Base.getindex(a::UnitRangeDual, indices::AbstractUnitRange{<:Integer})
   return dual(getindex(nondual(a), indices))
 end
@@ -92,6 +105,9 @@ BlockArrays.blockaxes(a::UnitRangeDual) = blockaxes(nondual(a))
 BlockArrays.blockfirsts(a::UnitRangeDual) = label_dual.(blockfirsts(nondual(a)))
 BlockArrays.blocklasts(a::UnitRangeDual) = label_dual.(blocklasts(nondual(a)))
 BlockArrays.findblock(a::UnitRangeDual, index::Integer) = findblock(nondual(a), index)
+
+blocklabels(a::UnitRangeDual) = dual.(blocklabels(nondual(a)))
+
 function BlockArrays.combine_blockaxes(a1::UnitRangeDual, a2::UnitRangeDual)
   return dual(combine_blockaxes(dual(a1), dual(a2)))
 end
diff --git a/NDTensors/src/lib/GradedAxes/test/runtests.jl b/NDTensors/src/lib/GradedAxes/test/runtests.jl
index 09335af5e8..c0fdca21be 100644
--- a/NDTensors/src/lib/GradedAxes/test/runtests.jl
+++ b/NDTensors/src/lib/GradedAxes/test/runtests.jl
@@ -2,7 +2,7 @@
 using Test: @testset
 @testset "GradedAxes" begin
   include("test_basics.jl")
-  include("test_tensor_product.jl")
   include("test_dual.jl")
+  include("test_tensor_product.jl")
 end
 end
diff --git a/NDTensors/src/lib/GradedAxes/test/test_basics.jl b/NDTensors/src/lib/GradedAxes/test/test_basics.jl
index e1dcc67174..360c2cecff 100644
--- a/NDTensors/src/lib/GradedAxes/test/test_basics.jl
+++ b/NDTensors/src/lib/GradedAxes/test/test_basics.jl
@@ -9,7 +9,8 @@ using BlockArrays:
   blocklength,
   blocklengths,
   blocks
-using NDTensors.GradedAxes: GradedOneTo, GradedUnitRange, blocklabels, gradedrange
+using NDTensors.GradedAxes:
+  GradedOneTo, GradedUnitRange, blocklabels, labelled_isequal, gradedrange
 using NDTensors.LabelledNumbers: LabelledUnitRange, islabelled, label, labelled, unlabel
 using Test: @test, @test_broken, @testset
 @testset "GradedAxes basics" begin
@@ -40,6 +41,7 @@ using Test: @test, @test_broken, @testset
       @test label(x) == "y"
     end
     @test isnothing(iterate(a, labelled(5, "y")))
+    @test labelled_isequal(a, a)
     @test length(a) == 5
     @test step(a) == 1
     @test !islabelled(step(a))
diff --git a/NDTensors/src/lib/GradedAxes/test/test_dual.jl b/NDTensors/src/lib/GradedAxes/test/test_dual.jl
index 0fb15d31bf..dc7f740713 100644
--- a/NDTensors/src/lib/GradedAxes/test/test_dual.jl
+++ b/NDTensors/src/lib/GradedAxes/test/test_dual.jl
@@ -1,19 +1,40 @@
 @eval module $(gensym())
-using BlockArrays: Block, blockaxes, blockfirsts, blocklasts, blocks, findblock
-using NDTensors.GradedAxes: GradedAxes, UnitRangeDual, dual, gradedrange, nondual
+using BlockArrays:
+  Block, blockaxes, blockfirsts, blocklasts, blocklength, blocklengths, blocks, findblock
+using NDTensors.GradedAxes:
+  GradedAxes,
+  UnitRangeDual,
+  blocklabels,
+  blockmergesortperm,
+  blocksortperm,
+  dual,
+  flip,
+  space_isequal,
+  gradedrange,
+  isdual,
+  nondual
 using NDTensors.LabelledNumbers: LabelledInteger, label, labelled
 using Test: @test, @test_broken, @testset
 struct U1
   n::Int
 end
 GradedAxes.dual(c::U1) = U1(-c.n)
+Base.isless(c1::U1, c2::U1) = c1.n < c2.n
 @testset "dual" begin
   a = gradedrange([U1(0) => 2, U1(1) => 3])
   ad = dual(a)
   @test eltype(ad) == LabelledInteger{Int,U1}
-  @test dual(ad) == a
-  @test nondual(ad) == a
-  @test nondual(a) == a
+
+  @test space_isequal(dual(ad), a)
+  @test space_isequal(nondual(ad), a)
+  @test space_isequal(nondual(a), a)
+  @test space_isequal(ad, ad)
+  @test !space_isequal(a, ad)
+  @test !space_isequal(ad, a)
+
+  @test isdual(ad)
+  @test !isdual(a)
+
   @test blockfirsts(ad) == [labelled(1, U1(0)), labelled(3, U1(-1))]
   @test blocklasts(ad) == [labelled(2, U1(0)), labelled(5, U1(-1))]
   @test findblock(ad, 4) == Block(2)
@@ -34,5 +55,36 @@ GradedAxes.dual(c::U1) = U1(-c.n)
   @test label(ad[[Block(2), Block(1)]][Block(1)]) == U1(-1)
   @test ad[[Block(2)[1:2], Block(1)[1:2]]][Block(1)] == 3:4
   @test label(ad[[Block(2)[1:2], Block(1)[1:2]]][Block(1)]) == U1(-1)
+  @test blocksortperm(a) == [Block(1), Block(2)]
+  @test blocksortperm(ad) == [Block(1), Block(2)]
+  @test blocklength(blockmergesortperm(a)) == 2
+  @test blocklength(blockmergesortperm(ad)) == 2
+  @test blockmergesortperm(a) == [Block(1), Block(2)]
+  @test blockmergesortperm(ad) == [Block(1), Block(2)]
+end
+
+@testset "flip" begin
+  a = gradedrange([U1(0) => 2, U1(1) => 3])
+  ad = dual(a)
+  @test space_isequal(flip(a), dual(gradedrange([U1(0) => 2, U1(-1) => 3])))
+  @test space_isequal(flip(ad), gradedrange([U1(0) => 2, U1(-1) => 3]))
+
+  @test blocklabels(a) == [U1(0), U1(1)]
+  @test blocklabels(dual(a)) == [U1(0), U1(-1)]
+  @test blocklabels(flip(a)) == [U1(0), U1(1)]
+  @test blocklabels(flip(dual(a))) == [U1(0), U1(-1)]
+  @test blocklabels(dual(flip(a))) == [U1(0), U1(-1)]
+
+  @test blocklengths(a) == [2, 3]
+  @test blocklengths(dual(a)) == [2, 3]
+  @test blocklengths(flip(a)) == [2, 3]
+  @test blocklengths(flip(dual(a))) == [2, 3]
+  @test blocklengths(dual(flip(a))) == [2, 3]
+
+  @test !isdual(a)
+  @test isdual(dual(a))
+  @test isdual(flip(a))
+  @test !isdual(flip(dual(a)))
+  @test !isdual(dual(flip(a)))
 end
 end
diff --git a/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl b/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl
index c10ae4bf95..02435b5ba7 100644
--- a/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl
+++ b/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl
@@ -1,11 +1,86 @@
 @eval module $(gensym())
-using NDTensors.GradedAxes: GradedAxes, GradedOneTo, gradedrange, tensor_product
 using Test: @test, @testset
+
+using BlockArrays: blocklength, blocklengths
+
+using NDTensors.GradedAxes:
+  GradedAxes,
+  GradedOneTo,
+  OneToOne,
+  dual,
+  fusion_product,
+  flip,
+  gradedrange,
+  labelled_isequal,
+  space_isequal,
+  isdual,
+  tensor_product
+
+struct U1
+  n::Int
+end
+GradedAxes.dual(c::U1) = U1(-c.n)
+Base.isless(c1::U1, c2::U1) = c1.n < c2.n
+GradedAxes.fuse_labels(x::U1, y::U1) = U1(x.n + y.n)
+
 @testset "GradedAxes.tensor_product" begin
   GradedAxes.fuse_labels(x::String, y::String) = x * y
+
+  g0 = OneToOne()
+  @test labelled_isequal(tensor_product(g0, g0), g0)
+
   a = gradedrange(["x" => 2, "y" => 3])
   b = tensor_product(a, a)
+  @test b isa GradedOneTo
   @test length(b) == 25
+  @test blocklength(b) == 4
+  @test blocklengths(b) == [4, 6, 6, 9]
+  @test labelled_isequal(b, gradedrange(["xx" => 4, "yx" => 6, "xy" => 6, "yy" => 9]))
+
+  c = tensor_product(a, a, a)
+  @test c isa GradedOneTo
+  @test length(c) == 125
+  @test blocklength(c) == 8
+end
+
+@testset "GradedAxes.fusion_product" begin
+  g0 = OneToOne()
+  @test labelled_isequal(fusion_product(g0, g0), g0)
+
+  a = gradedrange([U1(1) => 1, U1(2) => 3, U1(1) => 1])
+
+  b = fusion_product(a)
+  @test labelled_isequal(b, gradedrange([U1(1) => 2, U1(2) => 3]))
+
+  c = fusion_product(a, a)
+  @test labelled_isequal(c, gradedrange([U1(2) => 4, U1(3) => 12, U1(4) => 9]))
+
+  d = fusion_product(a, a, a)
+  @test labelled_isequal(
+    d, gradedrange([U1(3) => 8, U1(4) => 36, U1(5) => 54, U1(6) => 27])
+  )
+end
+
+@testset "dual and tensor_product" begin
+  a = gradedrange([U1(1) => 1, U1(2) => 3, U1(1) => 1])
+  ad = dual(a)
+
+  b = fusion_product(ad)
   @test b isa GradedOneTo
+  @test !isdual(b)
+  @test space_isequal(b, gradedrange([U1(-2) => 3, U1(-1) => 2]))
+
+  c = fusion_product(ad, ad)
+  @test c isa GradedOneTo
+  @test !isdual(c)
+  @test space_isequal(c, gradedrange([U1(-4) => 9, U1(-3) => 12, U1(-2) => 4]))
+
+  d = fusion_product(ad, a)
+  @test !isdual(d)
+  @test space_isequal(d, gradedrange([U1(-1) => 6, U1(0) => 13, U1(1) => 6]))
+
+  e = fusion_product(a, ad)
+  @test !isdual(d)
+  @test space_isequal(e, d)
 end
 end
diff --git a/NDTensors/src/lib/Sectors/src/Sectors.jl b/NDTensors/src/lib/Sectors/src/Sectors.jl
deleted file mode 100644
index fc813e602e..0000000000
--- a/NDTensors/src/lib/Sectors/src/Sectors.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-module Sectors
-
-include("abstractcategory.jl")
-include("category_definitions/u1.jl")
-include("category_definitions/zn.jl")
-include("category_definitions/su.jl")
-include("category_definitions/su2.jl")
-include("category_definitions/su2k.jl")
-include("category_definitions/fib.jl")
-include("category_definitions/ising.jl")
-
-include("namedtuple_operations.jl")
-include("category_product.jl")
-
-end
diff --git a/NDTensors/src/lib/Sectors/src/abstractcategory.jl b/NDTensors/src/lib/Sectors/src/abstractcategory.jl
deleted file mode 100644
index 4aaf56c6c5..0000000000
--- a/NDTensors/src/lib/Sectors/src/abstractcategory.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-abstract type AbstractCategory end
-
-label(c::AbstractCategory) = error("method `label` not defined for type $(typeof(c))")
-
-function dimension(c::AbstractCategory)
-  return error("method `dimension` not defined for type $(typeof(c))")
-end
-
-function label_fusion_rule(category_type::Type{<:AbstractCategory}, l1, l2)
-  return error("`label_fusion_rule` not defined for type $(category_type).")
-end
-
-function fusion_rule(c1::AbstractCategory, c2::AbstractCategory)
-  category_type = typeof(c1)
-  return [category_type(l) for l in label_fusion_rule(category_type, label(c1), label(c2))]
-end
-
-⊗(c1::AbstractCategory, c2::AbstractCategory) = fusion_rule(c1, c2)
-
-⊕(c1::AbstractCategory, c2::AbstractCategory) = [c1, c2]
-⊕(cs::Vector{<:AbstractCategory}, c::AbstractCategory) = [cs; c]
-⊕(c::AbstractCategory, cs::Vector{<:AbstractCategory}) = [c; cs]
-
-function Base.show(io::IO, cs::Vector{<:AbstractCategory})
-  (length(cs) <= 1) && print(io, "[")
-  symbol = ""
-  for c in cs
-    print(io, symbol, c)
-    symbol = " ⊕ "
-  end
-  (length(cs) <= 1) && print(io, "]")
-  return nothing
-end
-
-function trivial(category_type::Type{<:AbstractCategory})
-  return error("`trivial` not defined for type $(category_type).")
-end
-
-istrivial(c::AbstractCategory) = (c == trivial(typeof(c)))
-
-function dual(category_type::Type{<:AbstractCategory})
-  return error("`dual` not defined for type $(category_type).")
-end
-
-Base.isless(c1::AbstractCategory, c2::AbstractCategory) = isless(label(c1), label(c2))
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl b/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl
deleted file mode 100644
index 1d1089046b..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-using HalfIntegers: Half, twice
-
-#
-# Ising category
-#
-# (same fusion rules as su2{2})
-#
-
-struct Ising <: AbstractCategory
-  l::Half{Int}
-end
-
-# TODO: Use `Val` dispatch here?
-function Ising(s::AbstractString)
-  for (a, v) in enumerate(("1", "σ", "ψ"))
-    (v == s) && return Ising((a - 1)//2)
-  end
-  return error("Unrecognized input \"$s\" to Ising constructor")
-end
-
-dual(i::Ising) = i
-
-label(i::Ising) = i.l
-
-trivial(::Type{Ising}) = Ising(0)
-
-dimension(i::Ising) = (label(i) == 1//2) ? √2 : 1
-
-# Fusion rules identical to su2₂
-label_fusion_rule(::Type{Ising}, l1, l2) = label_fusion_rule(su2{2}, l1, l2)
-
-# TODO: Use `Val` dispatch here?
-label_to_str(i::Ising) = ("1", "σ", "ψ")[twice(label(i)) + 1]
-
-function Base.show(io::IO, f::Ising)
-  return print(io, "Ising(", label_to_str(f), ")")
-end
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su.jl
deleted file mode 100644
index 8fc82df329..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/su.jl
+++ /dev/null
@@ -1,83 +0,0 @@
-#
-# Special unitary group SU{N}
-#
-
-struct SU{N} <: AbstractCategory
-  # l is the first row of the
-  # Gelfand-Tsetlin (GT) pattern describing
-  # an SU(N) irrep
-  #TODO: any way this could be NTuple{N-1,Int} ?
-  l::NTuple{N,Int}
-end
-
-label(s::SU) = s.l
-
-groupdim(::SU{N}) where {N} = N
-
-trivial(::Type{SU{N}}) where {N} = SU{N}(ntuple(_ -> 0, Val(N)))
-
-fundamental(::Type{SU{N}}) where {N} = SU{N}(ntuple(i -> Int(i == 1), Val(N)))
-
-adjoint(::Type{SU{N}}) where {N} = SU{N}((ntuple(i -> Int(i == 1) + Int(i < N), Val(N))))
-
-function dimension(s::SU)
-  N = groupdim(s)
-  l = label(s)
-  d = 1
-  for k1 in 1:N, k2 in (k1 + 1):N
-    d *= ((k2 - k1) + (l[k1] - l[k2]))//(k2 - k1)
-  end
-  return Int(d)
-end
-
-function dual(s::SU)
-  l = label(s)
-  nl = ((reverse(cumsum(l[begin:(end - 1)] .- l[(begin + 1):end]))..., 0))
-  return typeof(s)(nl)
-end
-
-# display SU(N) irrep as a Young tableau with utf8 box char
-function Base.show(io::IO, ::MIME"text/plain", s::SU)
-  l = label(s)
-  if l[1] == 0  # singlet = no box
-    println(io, "●")
-    return nothing
-  end
-
-  println("┌─" * "┬─"^(l[1] - 1) * "┐")
-  i = 1
-  while l[i + 1] != 0
-    println(
-      io,
-      "├─",
-      "┼─"^(l[i + 1] - 1 + (l[i] > l[i + 1])),
-      "┴─"^max(0, (l[i] - l[i + 1] - 1)),
-      "┤"^(l[i] == l[i + 1]),
-      "┘"^(l[i] > l[i + 1]),
-    )
-    i += 1
-  end
-
-  println(io, "└─", "┴─"^max(0, l[i] - 1), "┘")
-  return nothing
-end
-
-#
-# Specializations for the case SU{2}
-# Where irreps specified by dimension "d"
-#
-
-dimension(s::SU{2}) = 1 + label(s)[1]
-
-SU{2}(d::Integer) = SU{2}((d - 1, 0))
-
-dual(s::SU{2}) = s
-
-function fusion_rule(s1::SU{2}, s2::SU{2})
-  d1, d2 = dimension(s1), dimension(s2)
-  return [SU{2}(d) for d in (abs(d1 - d2) + 1):2:(d1 + d2 - 1)]
-end
-
-function Base.show(io::IO, s::SU{2})
-  return print(io, "SU{2}(", dimension(s), ")")
-end
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl
deleted file mode 100644
index 2996b63f08..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-using HalfIntegers: Half, half, twice
-
-#
-# Conventional SU2 group
-# using "J" labels
-#
-
-struct SU2 <: AbstractCategory
-  j::Half{Int}
-end
-
-dual(s::SU2) = s
-
-label(s::SU2) = s.j
-
-trivial(::Type{SU2}) = SU2(0)
-fundamental(::Type{SU2}) = SU2(half(1))
-adjoint(::Type{SU2}) = SU2(1)
-
-dimension(s::SU2) = twice(label(s)) + 1
-
-label_fusion_rule(::Type{SU2}, j1, j2) = abs(j1 - j2):(j1 + j2)
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl
deleted file mode 100644
index 86cf5f5abf..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-using HalfIntegers: Half
-
-#
-# Quantum 'group' su2ₖ
-#
-
-struct su2{k} <: AbstractCategory
-  j::Half{Int}
-end
-
-dual(s::su2) = s
-
-label(s::su2) = s.j
-
-level(s::su2{k}) where {k} = k
-
-trivial(::Type{su2{k}}) where {k} = su2{k}(0)
-
-function label_fusion_rule(::Type{su2{k}}, j1, j2) where {k}
-  return abs(j1 - j2):min(k - j1 - j2, j1 + j2)
-end
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl b/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl
deleted file mode 100644
index 7af8d22b9f..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl
+++ /dev/null
@@ -1,19 +0,0 @@
-using HalfIntegers: Half
-
-#
-# U₁ group (circle group, or particle number, total Sz etc.)
-#
-
-struct U1 <: AbstractCategory
-  n::Half{Int}
-end
-
-dual(u::U1) = U1(-u.n)
-
-label(u::U1) = u.n
-
-dimension(::U1) = 1
-
-trivial(::Type{U1}) = U1(0)
-
-label_fusion_rule(::Type{U1}, n1, n2) = (n1 + n2,)
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl b/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl
deleted file mode 100644
index 1348052d97..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-using HalfIntegers: Half
-
-#
-# Cyclic group Zₙ
-#
-
-struct Z{N} <: AbstractCategory
-  m::Half{Int}
-  Z{N}(m) where {N} = new{N}(m % N)
-end
-
-label(c::Z) = c.m
-modulus(::Type{Z{N}}) where {N} = N
-
-modulus(c::Z) = modulus(typeof(c))
-
-dimension(::Z) = 1
-
-trivial(category_type::Type{<:Z}) = category_type(0)
-
-label_fusion_rule(category_type::Type{<:Z}, n1, n2) = ((n1 + n2) % modulus(category_type),)
-
-dual(c::Z) = typeof(c)(mod(-label(c), modulus(c)))
diff --git a/NDTensors/src/lib/Sectors/src/category_product.jl b/NDTensors/src/lib/Sectors/src/category_product.jl
deleted file mode 100644
index 12baf35be1..0000000000
--- a/NDTensors/src/lib/Sectors/src/category_product.jl
+++ /dev/null
@@ -1,111 +0,0 @@
-
-struct CategoryProduct{Categories} <: AbstractCategory
-  cats::Categories
-  global _CategoryProduct(l) = new{typeof(l)}(l)
-end
-
-CategoryProduct(c::CategoryProduct) = _CategoryProduct(categories(c))
-
-categories(s::CategoryProduct) = s.cats
-
-Base.isempty(S::CategoryProduct) = isempty(categories(S))
-Base.length(S::CategoryProduct) = length(categories(S))
-Base.getindex(S::CategoryProduct, args...) = getindex(categories(S), args...)
-
-function fusion_rule(s1::CategoryProduct, s2::CategoryProduct)
-  return [
-    CategoryProduct(l) for l in categories_fusion_rule(categories(s1), categories(s2))
-  ]
-end
-
-function Base.:(==)(A::CategoryProduct, B::CategoryProduct)
-  return categories_equal(categories(A), categories(B))
-end
-
-function Base.show(io::IO, s::CategoryProduct)
-  (length(s) < 2) && print(io, "sector")
-  print(io, "(")
-  symbol = ""
-  for p in pairs(categories(s))
-    print(io, symbol)
-    category_show(io, p[1], p[2])
-    symbol = " × "
-  end
-  return print(io, ")")
-end
-
-category_show(io::IO, k, v) = print(io, v)
-
-category_show(io::IO, k::Symbol, v) = print(io, "($k=$v,)")
-
-×(c1::AbstractCategory, c2::AbstractCategory) = ×(CategoryProduct(c1), CategoryProduct(c2))
-function ×(p1::CategoryProduct, p2::CategoryProduct)
-  return CategoryProduct(categories_product(categories(p1), categories(p2)))
-end
-
-categories_product(l1::NamedTuple, l2::NamedTuple) = union_keys(l1, l2)
-
-categories_product(l1::Tuple, l2::Tuple) = (l1..., l2...)
-
-×(nt1::NamedTuple, nt2::NamedTuple) = ×(CategoryProduct(nt1), CategoryProduct(nt2))
-×(c1::NamedTuple, c2::AbstractCategory) = ×(CategoryProduct(c1), CategoryProduct(c2))
-×(c1::AbstractCategory, c2::NamedTuple) = ×(CategoryProduct(c1), CategoryProduct(c2))
-
-#
-# Dictionary-like implementation
-#
-
-function CategoryProduct(nt::NamedTuple)
-  categories = sort_keys(nt)
-  return _CategoryProduct(categories)
-end
-
-CategoryProduct(; kws...) = CategoryProduct((; kws...))
-
-function CategoryProduct(pairs::Pair...)
-  keys = ntuple(n -> Symbol(pairs[n][1]), length(pairs))
-  vals = ntuple(n -> pairs[n][2], length(pairs))
-  return CategoryProduct(NamedTuple{keys}(vals))
-end
-
-function categories_fusion_rule(A::NamedTuple, B::NamedTuple)
-  qs = [A]
-  for (la, lb) in zip(pairs(intersect_keys(A, B)), pairs(intersect_keys(B, A)))
-    @assert la[1] == lb[1]
-    fused_vals = ⊗(la[2], lb[2])
-    qs = [union_keys((; la[1] => v), q) for v in fused_vals for q in qs]
-  end
-  # Include sectors of B not in A
-  qs = [union_keys(q, B) for q in qs]
-  return qs
-end
-
-function categories_equal(A::NamedTuple, B::NamedTuple)
-  common_categories = zip(pairs(intersect_keys(A, B)), pairs(intersect_keys(B, A)))
-  common_categories_match = all(nl -> (nl[1] == nl[2]), common_categories)
-  unique_categories_zero = all(l -> istrivial(l), symdiff_keys(A, B))
-  return common_categories_match && unique_categories_zero
-end
-
-#
-# Ordered implementation
-#
-
-CategoryProduct(t::Tuple) = _CategoryProduct(t)
-CategoryProduct(cats::AbstractCategory...) = CategoryProduct((cats...,))
-
-categories_equal(o1::Tuple, o2::Tuple) = (o1 == o2)
-
-function categories_fusion_rule(o1::Tuple, o2::Tuple)
-  N = length(o1)
-  length(o2) == N ||
-    throw(DimensionMismatch("Ordered CategoryProduct must have same size in ⊗"))
-  os = [o1]
-  replace(o, n, val) = ntuple(m -> (m == n) ? val : o[m], length(o))
-  for n in 1:N
-    os = [replace(o, n, f) for f in ⊗(o1[n], o2[n]) for o in os]
-  end
-  return os
-end
-
-sector(args...; kws...) = CategoryProduct(args...; kws...)
diff --git a/NDTensors/src/lib/Sectors/test/test_category.jl b/NDTensors/src/lib/Sectors/test/test_category.jl
deleted file mode 100644
index e14582afc8..0000000000
--- a/NDTensors/src/lib/Sectors/test/test_category.jl
+++ /dev/null
@@ -1,192 +0,0 @@
-@eval module $(gensym())
-using NDTensors.Sectors:
-  Fib,
-  Ising,
-  SU,
-  SU2,
-  U1,
-  Z,
-  ⊗,
-  ⊕,
-  dimension,
-  dual,
-  istrivial,
-  trivial,
-  fundamental,
-  adjoint
-using Test: @test, @testset
-@testset "Test Category Types" begin
-  @testset "U(1)" begin
-    q1 = U1(1)
-    q2 = U1(2)
-    q3 = U1(3)
-
-    @test dimension(q1) == 1
-    @test dimension(q2) == 1
-
-    @test q1 ⊗ q1 == [q2]
-    @test q1 ⊗ q2 == [q3]
-    @test q2 ⊗ q1 == [q3]
-
-    @test trivial(U1) == U1(0)
-    @test istrivial(U1(0))
-
-    @test dual(U1(2)) == U1(-2)
-    @test isless(U1(1), U1(2))
-    @test !isless(U1(2), U1(1))
-  end
-
-  @testset "Z₂" begin
-    z0 = Z{2}(0)
-    z1 = Z{2}(1)
-
-    @test trivial(Z{2}) == Z{2}(0)
-    @test istrivial(Z{2}(0))
-
-    @test dimension(z0) == 1
-    @test dimension(z1) == 1
-
-    @test dual(z0) == z0
-    @test dual(z1) == z1
-
-    @test z0 ⊗ z0 == [z0]
-    @test z0 ⊗ z1 == [z1]
-    @test z1 ⊗ z1 == [z0]
-
-    @test dual(Z{2}(1)) == Z{2}(1)
-    @test isless(Z{2}(0), Z{2}(1))
-    @test !isless(Z{2}(1), Z{2}(0))
-  end
-
-  @testset "SU2" begin
-    j1 = SU2(0)
-    j2 = SU2(1//2)
-    j3 = SU2(1)
-    j4 = SU2(3//2)
-    j5 = SU2(2)
-
-    @test trivial(SU2) == SU2(0)
-    @test istrivial(SU2(0))
-
-    @test fundamental(SU2) == SU2(1//2)
-    @test adjoint(SU2) == SU2(1)
-
-    @test dimension(j1) == 1
-    @test dimension(j2) == 2
-    @test dimension(j3) == 3
-    @test dimension(j4) == 4
-
-    @test dual(j1) == j1
-    @test dual(j2) == j2
-    @test dual(j3) == j3
-    @test dual(j4) == j4
-
-    @test j1 ⊗ j2 == [j2]
-    @test j2 ⊗ j2 == j1 ⊕ j3
-    @test j2 ⊗ j3 == j2 ⊕ j4
-    @test j3 ⊗ j3 == j1 ⊕ j3 ⊕ j5
-  end
-
-  @testset "SU(2)" begin
-    j1 = SU{2}(1)
-    j2 = SU{2}(2)
-    j3 = SU{2}(3)
-    j4 = SU{2}(4)
-    j5 = SU{2}(5)
-
-    @test trivial(SU{2}) == SU{2}(1)
-    @test istrivial(SU{2}(1))
-
-    @test fundamental(SU{2}) == SU{2}(2)
-    @test adjoint(SU{2}) == SU{2}(3)
-
-    @test dimension(j1) == 1
-    @test dimension(j2) == 2
-    @test dimension(j3) == 3
-    @test dimension(j4) == 4
-
-    @test dual(j1) == j1
-    @test dual(j2) == j2
-    @test dual(j3) == j3
-    @test dual(j4) == j4
-
-    @test j1 ⊗ j2 == [j2]
-    @test j2 ⊗ j2 == j1 ⊕ j3
-    @test j2 ⊗ j3 == j2 ⊕ j4
-    @test j3 ⊗ j3 == j1 ⊕ j3 ⊕ j5
-  end
-
-  @testset "SU(N)" begin
-    f3 = SU{3}((1, 0, 0))
-    f4 = SU{4}((1, 0, 0, 0))
-    ad3 = SU{3}((2, 1, 0))
-    ad4 = SU{4}((2, 1, 1, 0))
-
-    @test trivial(SU{3}) == SU{3}((0, 0, 0))
-    @test istrivial(SU{3}((0, 0, 0)))
-    @test trivial(SU{4}) == SU{4}((0, 0, 0, 0))
-    @test istrivial(SU{4}((0, 0, 0, 0)))
-
-    @test fundamental(SU{3}) == f3
-    @test adjoint(SU{3}) == ad3
-    @test fundamental(SU{4}) == f4
-    @test adjoint(SU{4}) == ad4
-
-    @test dual(f3) == SU{3}((1, 1, 0))
-    @test dual(f4) == SU{4}((1, 1, 1, 0))
-    @test dual(ad3) == ad3
-    @test dual(ad4) == ad4
-
-    @test dimension(f3) == 3
-    @test dimension(f4) == 4
-    @test dimension(ad3) == 8
-    @test dimension(ad4) == 15
-    @test dimension(SU{3}((4, 2, 0))) == 27
-    @test dimension(SU{3}((3, 3, 0))) == 10
-    @test dimension(SU{3}((3, 0, 0))) == 10
-    @test dimension(SU{3}((0, 0, 0))) == 1
-  end
-
-  @testset "Fibonacci" begin
-    ı = Fib("1")
-    τ = Fib("τ")
-
-    @test trivial(Fib) == ı
-    @test istrivial(ı)
-
-    @test dual(ı) == ı
-    @test dual(τ) == τ
-
-    @test dimension(ı) == 1
-    @test dimension(τ) == ((1 + √5) / 2)
-
-    @test ı ⊗ ı == [ı]
-    @test ı ⊗ τ == [τ]
-    @test τ ⊗ ı == [τ]
-    @test τ ⊗ τ == ı ⊕ τ
-  end
-
-  @testset "Ising" begin
-    ı = Ising("1")
-    σ = Ising("σ")
-    ψ = Ising("ψ")
-
-    @test trivial(Ising) == ı
-    @test istrivial(ı)
-
-    @test dual(ı) == ı
-    @test dual(σ) == σ
-    @test dual(ψ) == ψ
-
-    @test ı ⊗ ı == [ı]
-    @test ı ⊗ σ == [σ]
-    @test σ ⊗ ı == [σ]
-    @test ı ⊗ ψ == [ψ]
-    @test ψ ⊗ ı == [ψ]
-    @test σ ⊗ σ == ı ⊕ ψ
-    @test σ ⊗ ψ == [σ]
-    @test ψ ⊗ σ == [σ]
-    @test ψ ⊗ ψ == [ı]
-  end
-end
-end
diff --git a/NDTensors/src/lib/Sectors/test/test_category_product.jl b/NDTensors/src/lib/Sectors/test/test_category_product.jl
deleted file mode 100644
index 58f095c110..0000000000
--- a/NDTensors/src/lib/Sectors/test/test_category_product.jl
+++ /dev/null
@@ -1,135 +0,0 @@
-@eval module $(gensym())
-using NDTensors.Sectors: Fib, Ising, SU, SU2, U1, Z, ⊗, ⊕, ×, sector
-using Test: @test, @testset, @test_throws
-@testset "Test Named Category Products" begin
-  @testset "Construct from × of NamedTuples" begin
-    s = (A=U1(1),) × (B=SU2(2),)
-    @test length(s) == 2
-    @test s[:A] == U1(1)
-    @test s[:B] == SU2(2)
-
-    s = s × (C=Ising("ψ"),)
-    @test length(s) == 3
-    @test s[:C] == Ising("ψ")
-  end
-
-  @testset "Construct from Pairs" begin
-    s = sector("A" => U1(2))
-    @test length(s) == 1
-    @test s[:A] == U1(2)
-    @test s == sector(; A=U1(2))
-
-    s = sector("B" => Ising("ψ"), :C => Z{2}(1))
-    @test length(s) == 2
-    @test s[:B] == Ising("ψ")
-    @test s[:C] == Z{2}(1)
-  end
-
-  @testset "Multiple U(1)'s" begin
-    q00 = sector()
-    q10 = sector(; A=U1(1))
-    q01 = sector(; B=U1(1))
-    q11 = sector(; A=U1(1), B=U1(1))
-
-    @test q00 ⊗ q00 == [q00]
-    @test q01 ⊗ q00 == [q01]
-    @test q00 ⊗ q01 == [q01]
-    @test q10 ⊗ q01 == [q11]
-  end
-
-  @testset "U(1) ⊗ SU(2) conventional" begin
-    q0 = sector()
-    q0h = sector(; J=SU2(1//2))
-    q10 = (N=U1(1),) × (J=SU2(0),)
-    # Put names in reverse order sometimes:
-    q1h = (J=SU2(1//2),) × (N=U1(1),)
-    q11 = (N=U1(1),) × (J=SU2(1),)
-    q20 = sector(; N=U1(2))
-    q2h = (N=U1(2),) × (J=SU2(1//2),)
-    q21 = (N=U1(2),) × (J=SU2(1),)
-    q22 = (N=U1(2),) × (J=SU2(2),)
-
-    @test q1h ⊗ q1h == q20 ⊕ q21
-    @test q10 ⊗ q1h == [q2h]
-    @test q0h ⊗ q1h == q10 ⊕ q11
-    @test q11 ⊗ q11 == q20 ⊕ q21 ⊕ q22
-  end
-
-  @testset "U(1) ⊗ SU(2)" begin
-    q0 = sector()
-    q0h = sector(; J=SU{2}(2))
-    q10 = (N=U1(1),) × (J=SU{2}(1),)
-    # Put names in reverse order sometimes:
-    q1h = (J=SU{2}(2),) × (N=U1(1),)
-    q11 = (N=U1(1),) × (J=SU{2}(3),)
-    q20 = sector(; N=U1(2))
-    q2h = (N=U1(2),) × (J=SU{2}(2),)
-    q21 = (N=U1(2),) × (J=SU{2}(3),)
-    q22 = (N=U1(2),) × (J=SU{2}(5),)
-
-    @test q1h ⊗ q1h == q20 ⊕ q21
-    @test q10 ⊗ q1h == [q2h]
-    @test q0h ⊗ q1h == q10 ⊕ q11
-    @test q11 ⊗ q11 == q20 ⊕ q21 ⊕ q22
-  end
-
-  @testset "Comparisons with unspecified labels" begin
-    q2 = sector(; N=U1(2))
-    q20 = (N=U1(2),) × (J=SU{2}(1),)
-    @test q20 == q2
-
-    q21 = (N=U1(2),) × (J=SU{2}(3),)
-    @test q21 != q2
-
-    a = (A=U1(0),) × (B=U1(2),)
-    b = (B=U1(2),) × (C=U1(0),)
-    @test a == b
-    c = (B=U1(2),) × (C=U1(1),)
-    @test a != c
-  end
-end
-
-@testset "Test Ordered Products" begin
-  @testset "Ordered Constructor" begin
-    s = sector(U1(1), U1(2))
-    @test length(s) == 2
-    @test s[1] == U1(1)
-    @test s[2] == U1(2)
-
-    s = U1(1) × SU2(1//2) × U1(3)
-    @test length(s) == 3
-    @test s[1] == U1(1)
-    @test s[2] == SU2(1//2)
-    @test s[3] == U1(3)
-  end
-
-  @testset "Fusion of U1 products" begin
-    p11 = U1(1) × U1(1)
-    @test p11 ⊗ p11 == [U1(2) × U1(2)]
-
-    p123 = U1(1) × U1(2) × U1(3)
-    @test p123 ⊗ p123 == [U1(2) × U1(4) × U1(6)]
-  end
-
-  @testset "Enforce same number of spaces" begin
-    p12 = U1(1) × U1(2)
-    p123 = U1(1) × U1(2) × U1(3)
-    @test_throws DimensionMismatch p12 ⊗ p123
-  end
-
-  @testset "Fusion of SU2 products" begin
-    phh = SU2(1//2) × SU2(1//2)
-    @test phh ⊗ phh ==
-      (SU2(0) × SU2(0)) ⊕ (SU2(1) × SU2(0)) ⊕ (SU2(0) × SU2(1)) ⊕ (SU2(1) × SU2(1))
-  end
-
-  @testset "Fusion of mixed U1 and SU2 products" begin
-    p2h = U1(2) × SU2(1//2)
-    p1h = U1(1) × SU2(1//2)
-    @test p2h ⊗ p1h == (U1(3) × SU2(0)) ⊕ (U1(3) × SU2(1))
-
-    p1h1 = U1(1) × SU2(1//2) × Z{2}(1)
-    @test p1h1 ⊗ p1h1 == (U1(2) × SU2(0) × Z{2}(0)) ⊕ (U1(2) × SU2(1) × Z{2}(0))
-  end
-end
-end
diff --git a/NDTensors/src/lib/Sectors/.JuliaFormatter.toml b/NDTensors/src/lib/SymmetrySectors/.JuliaFormatter.toml
similarity index 100%
rename from NDTensors/src/lib/Sectors/.JuliaFormatter.toml
rename to NDTensors/src/lib/SymmetrySectors/.JuliaFormatter.toml
diff --git a/NDTensors/src/lib/Sectors/Project.toml b/NDTensors/src/lib/SymmetrySectors/Project.toml
similarity index 100%
rename from NDTensors/src/lib/Sectors/Project.toml
rename to NDTensors/src/lib/SymmetrySectors/Project.toml
diff --git a/NDTensors/src/lib/SymmetrySectors/src/SymmetrySectors.jl b/NDTensors/src/lib/SymmetrySectors/src/SymmetrySectors.jl
new file mode 100644
index 0000000000..93ecbda986
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/SymmetrySectors.jl
@@ -0,0 +1,17 @@
+module SymmetrySectors
+
+include("symmetry_style.jl")
+include("abstractsector.jl")
+include("sector_definitions/fib.jl")
+include("sector_definitions/ising.jl")
+include("sector_definitions/o2.jl")
+include("sector_definitions/trivial.jl")
+include("sector_definitions/su.jl")
+include("sector_definitions/su2k.jl")
+include("sector_definitions/u1.jl")
+include("sector_definitions/zn.jl")
+
+include("namedtuple_operations.jl")
+include("sector_product.jl")
+
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/abstractsector.jl b/NDTensors/src/lib/SymmetrySectors/src/abstractsector.jl
new file mode 100644
index 0000000000..2257e9fb36
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/abstractsector.jl
@@ -0,0 +1,114 @@
+# This file defines the abstract type AbstractSector
+# all fusion categories (Z{2}, SU2, Ising...) are subtypes of AbstractSector
+
+using BlockArrays: blocklengths
+using ..LabelledNumbers: LabelledInteger, label, label_type, labelled, unlabel, unlabel_type
+using ..GradedAxes: GradedAxes, blocklabels, fuse_blocklengths, gradedrange, tensor_product
+
+abstract type AbstractSector end
+
+# ===================================  Base interface  =====================================
+function Base.isless(c1::C, c2::C) where {C<:AbstractSector}
+  return isless(sector_label(c1), sector_label(c2))
+end
+
+# =================================  Sectors interface  ====================================
+trivial(x) = trivial(typeof(x))
+function trivial(axis_type::Type{<:AbstractUnitRange})
+  return gradedrange([trivial(eltype(axis_type))])  # always returns nondual
+end
+function trivial(la_type::Type{<:LabelledInteger})
+  return labelled(one(unlabel_type(la_type)), trivial(label_type(la_type)))
+end
+function trivial(type::Type)
+  return error("`trivial` not defined for type $(type).")
+end
+
+istrivial(c::AbstractSector) = (c == trivial(c))
+
+function sector_label(c::AbstractSector)
+  return error("method `sector_label` not defined for type $(typeof(c))")
+end
+
+block_dimensions(g::AbstractUnitRange) = block_dimensions(SymmetryStyle(g), g)
+block_dimensions(::AbelianStyle, g) = unlabel.(blocklengths(g))
+function block_dimensions(::NotAbelianStyle, g)
+  return quantum_dimension.(blocklabels(g)) .* blocklengths(g)
+end
+
+quantum_dimension(x) = quantum_dimension(SymmetryStyle(x), x)
+
+function quantum_dimension(::NotAbelianStyle, c::AbstractSector)
+  return error("method `quantum_dimension` not defined for type $(typeof(c))")
+end
+
+quantum_dimension(::AbelianStyle, ::AbstractSector) = 1
+quantum_dimension(::AbelianStyle, g::AbstractUnitRange) = length(g)
+quantum_dimension(::NotAbelianStyle, g::AbstractUnitRange) = sum(block_dimensions(g))
+
+# ===============================  Fusion rule interface  ==================================
+⊗(c1::AbstractSector, c2::AbstractSector) = fusion_rule(c1, c2)
+
+function fusion_rule(c1::AbstractSector, c2::AbstractSector)
+  return fusion_rule(combine_styles(SymmetryStyle(c1), SymmetryStyle(c2)), c1, c2)
+end
+
+function fusion_rule(::NotAbelianStyle, c1::C, c2::C) where {C<:AbstractSector}
+  sector_degen_pairs = label_fusion_rule(C, sector_label(c1), sector_label(c2))
+  return gradedrange(sector_degen_pairs)
+end
+
+# abelian case: return Sector
+function fusion_rule(::AbelianStyle, c1::C, c2::C) where {C<:AbstractSector}
+  return label(only(fusion_rule(NotAbelianStyle(), c1, c2)))
+end
+
+function label_fusion_rule(sector_type::Type{<:AbstractSector}, l1, l2)
+  return [abelian_label_fusion_rule(sector_type, l1, l2) => 1]
+end
+
+# ================================  GradedAxes interface  ==================================
+# tensor_product interface
+function GradedAxes.fuse_blocklengths(
+  l1::LabelledInteger{<:Integer,<:AbstractSector},
+  l2::LabelledInteger{<:Integer,<:AbstractSector},
+)
+  return fuse_blocklengths(combine_styles(SymmetryStyle(l1), SymmetryStyle(l2)), l1, l2)
+end
+
+function GradedAxes.fuse_blocklengths(
+  ::NotAbelianStyle, l1::LabelledInteger, l2::LabelledInteger
+)
+  fused = label(l1) ⊗ label(l2)
+  v = labelled.(l1 * l2 .* blocklengths(fused), blocklabels(fused))
+  return gradedrange(v)
+end
+
+function GradedAxes.fuse_blocklengths(
+  ::AbelianStyle, l1::LabelledInteger, l2::LabelledInteger
+)
+  fused = label(l1) ⊗ label(l2)
+  return gradedrange([labelled(l1 * l2, fused)])
+end
+
+# cast to range
+to_gradedrange(c::AbstractSector) = to_gradedrange(labelled(1, c))
+to_gradedrange(l::LabelledInteger) = gradedrange([l])
+to_gradedrange(g::AbstractUnitRange) = g
+
+# allow to fuse a Sector with a GradedUnitRange
+function GradedAxes.tensor_product(c::AbstractSector, g::AbstractUnitRange)
+  return tensor_product(to_gradedrange(c), g)
+end
+
+function GradedAxes.tensor_product(g::AbstractUnitRange, c::AbstractSector)
+  return tensor_product(g, to_gradedrange(c))
+end
+
+function GradedAxes.tensor_product(c1::AbstractSector, c2::AbstractSector)
+  return to_gradedrange(fusion_rule(c1, c2))
+end
+
+function GradedAxes.fusion_product(c::AbstractSector)
+  return to_gradedrange(c)
+end
diff --git a/NDTensors/src/lib/Sectors/src/namedtuple_operations.jl b/NDTensors/src/lib/SymmetrySectors/src/namedtuple_operations.jl
similarity index 100%
rename from NDTensors/src/lib/Sectors/src/namedtuple_operations.jl
rename to NDTensors/src/lib/SymmetrySectors/src/namedtuple_operations.jl
diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/fib.jl
similarity index 50%
rename from NDTensors/src/lib/Sectors/src/category_definitions/fib.jl
rename to NDTensors/src/lib/SymmetrySectors/src/sector_definitions/fib.jl
index b7b8a89d49..d7fded55f2 100644
--- a/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/fib.jl
@@ -3,8 +3,9 @@
 #
 # (same fusion rules as subcategory {0,1} of su2{3})
 #
+using ..GradedAxes: GradedAxes
 
-struct Fib <: AbstractCategory
+struct Fib <: AbstractSector
   l::Int
 end
 
@@ -18,16 +19,24 @@ function Fib(s::AbstractString)
   return error("Unrecognized input \"$s\" to Fib constructor")
 end
 
-dual(f::Fib) = f
+SymmetryStyle(::Type{Fib}) = NotAbelianStyle()
 
-label(f::Fib) = f.l
+GradedAxes.dual(f::Fib) = f
+
+sector_label(f::Fib) = f.l
 
 trivial(::Type{Fib}) = Fib(0)
 
-dimension(f::Fib) = istrivial(f) ? 1 : ((1 + √5) / 2)
+quantum_dimension(::NotAbelianStyle, f::Fib) = istrivial(f) ? 1.0 : ((1 + √5) / 2)
 
 # Fusion rules identical to su2₃
-label_fusion_rule(::Type{Fib}, l1, l2) = label_fusion_rule(su2{3}, l1, l2)
+function label_fusion_rule(::Type{Fib}, l1, l2)
+  suk_sectors_degen = label_fusion_rule(su2{3}, l1, l2)
+  suk_sectors = first.(suk_sectors_degen)
+  degen = last.(suk_sectors_degen)
+  sectors = Fib.(sector_label.(suk_sectors))
+  return sectors .=> degen
+end
 
 label_to_str(f::Fib) = istrivial(f) ? "1" : "τ"
 
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/ising.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/ising.jl
new file mode 100644
index 0000000000..d4a955069e
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/ising.jl
@@ -0,0 +1,46 @@
+#
+# Ising category
+#
+# (same fusion rules as su2{2})
+#
+
+using HalfIntegers: Half, twice
+using ..GradedAxes: GradedAxes
+
+struct Ising <: AbstractSector
+  l::Half{Int}
+end
+
+# TODO: Use `Val` dispatch here?
+function Ising(s::AbstractString)
+  for (a, v) in enumerate(("1", "σ", "ψ"))
+    (v == s) && return Ising((a - 1)//2)
+  end
+  return error("Unrecognized input \"$s\" to Ising constructor")
+end
+
+SymmetryStyle(::Type{Ising}) = NotAbelianStyle()
+
+GradedAxes.dual(i::Ising) = i
+
+sector_label(i::Ising) = i.l
+
+trivial(::Type{Ising}) = Ising(0)
+
+quantum_dimension(::NotAbelianStyle, i::Ising) = (sector_label(i) == 1//2) ? √2 : 1.0
+
+# Fusion rules identical to su2₂
+function label_fusion_rule(::Type{Ising}, l1, l2)
+  suk_sectors_degen = label_fusion_rule(su2{2}, l1, l2)
+  suk_sectors = first.(suk_sectors_degen)
+  degen = last.(suk_sectors_degen)
+  sectors = Ising.(sector_label.(suk_sectors))
+  return sectors .=> degen
+end
+
+# TODO: Use `Val` dispatch here?
+label_to_str(i::Ising) = ("1", "σ", "ψ")[twice(sector_label(i)) + 1]
+
+function Base.show(io::IO, f::Ising)
+  return print(io, "Ising(", label_to_str(f), ")")
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/o2.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/o2.jl
new file mode 100644
index 0000000000..4a0dcf7646
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/o2.jl
@@ -0,0 +1,75 @@
+#
+# Orthogonal group O(2)
+# isomorphic to Z_2 ⋉ U(1)
+# isomorphic to SU(2) subgroup with Sz conservation + Sz-reversal
+#
+# O(2) has 3 kinds of irreps:
+# - trivial irrep, or "0e", corresponds to Sz=0 and even under Sz-reversal
+# - "zero odd", or "0o" irrep, corresponds to Sz=0 and odd under Sz-reversal
+# - 2-dimensional Sz=±|m| irrep, with m a half integer
+#
+
+using HalfIntegers: Half, HalfInteger
+using ..GradedAxes: GradedAxes
+
+# here we use only one half-integer as label:
+# - l=0 for trivial
+# - l=-1 for zero odd
+# - l=+|m| for Sz=±|m|
+struct O2 <: AbstractSector
+  l::Half{Int}
+end
+
+SymmetryStyle(::Type{O2}) = NotAbelianStyle()
+
+sector_label(s::O2) = s.l
+
+trivial(::Type{O2}) = O2(0)
+zero_odd(::Type{O2}) = O2(-1)
+
+is_zero_even_or_odd(s::O2) = is_zero_even_or_odd(sector_label(s))
+iszero_odd(s::O2) = iszero_odd(sector_label(s))
+
+is_zero_even_or_odd(l::HalfInteger) = iszero_even(l) || iszero_odd(l)
+iszero_even(l::HalfInteger) = l == sector_label(trivial(O2))
+iszero_odd(l::HalfInteger) = l == sector_label(zero_odd(O2))
+
+quantum_dimension(::NotAbelianStyle, s::O2) = 2 - is_zero_even_or_odd(s)
+
+GradedAxes.dual(s::O2) = s
+
+function Base.show(io::IO, s::O2)
+  if iszero_odd(s)
+    disp = "0o"
+  elseif istrivial(s)
+    disp = "0e"
+  else
+    disp = "±" * string(sector_label(s))
+  end
+  return print(io, "O(2)[", disp, "]")
+end
+
+function label_fusion_rule(::Type{O2}, l1, l2)
+  if is_zero_even_or_odd(l1)
+    degens = [1]
+    if is_zero_even_or_odd(l2)
+      labels = l1 == l2 ? [sector_label(trivial(O2))] : [sector_label(zero_odd(O2))]
+    else
+      labels = [l2]
+    end
+  else
+    if is_zero_even_or_odd(l2)
+      degens = [1]
+      labels = [l1]
+    else
+      if l1 == l2
+        degens = [1, 1, 1]
+        labels = [sector_label(zero_odd(O2)), sector_label(trivial(O2)), 2 * l1]
+      else
+        degens = [1, 1]
+        labels = [abs(l1 - l2), l1 + l2]
+      end
+    end
+  end
+  return O2.(labels) .=> degens
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su.jl
new file mode 100644
index 0000000000..b44ef4d6eb
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su.jl
@@ -0,0 +1,175 @@
+#
+# Special unitary group SU(N)
+#
+
+using HalfIntegers: HalfInteger, half, twice
+using ...GradedAxes: GradedAxes
+
+struct SU{N,M} <: AbstractSector
+  # l is the first row of the
+  # Gelfand-Tsetlin (GT) pattern describing
+  # an SU(N) irrep
+  # this first row is identical to the Young tableau of the irrep
+  l::NTuple{M,Int}
+
+  # M is there to avoid storing a N-Tuple with an extra zero.
+  # inner constructor enforces M = N - 1
+  # It does NOT check for Young Tableau validity (non-increasing positive integers)
+  function SU{N,M}(t::NTuple{M,Integer}) where {N,M}
+    return N == M + 1 && M > 0 ? new{N,M}(t) : error("Invalid tuple length")
+  end
+end
+
+SU{N}(t::Tuple) where {N} = SU{N,length(t)}(t)
+SU(t::Tuple) = SU{length(t) + 1}(t)  # infer N from tuple length
+
+SymmetryStyle(::Type{<:SU}) = NotAbelianStyle()
+
+sector_label(s::SU) = s.l
+
+groupdim(::SU{N}) where {N} = N
+
+trivial(::Type{<:SU{N}}) where {N} = SU{N}(ntuple(_ -> 0, Val(N - 1)))
+
+fundamental(::Type{<:SU{N}}) where {N} = SU{N}(ntuple(i -> i == 1, Val(N - 1)))
+
+function quantum_dimension(::NotAbelianStyle, s::SU)
+  N = groupdim(s)
+  l = (sector_label(s)..., 0)
+  d = 1
+  for k1 in 1:N, k2 in (k1 + 1):N
+    d *= ((k2 - k1) + (l[k1] - l[k2]))//(k2 - k1)
+  end
+  return Int(d)
+end
+
+function GradedAxes.dual(s::SU)
+  l = sector_label(s)
+  nl = reverse(cumsum((l[begin:(end - 1)] .- l[(begin + 1):end]..., l[end])))
+  return typeof(s)(nl)
+end
+
+function Base.show(io::IO, s::SU)
+  disp = join([string(l) for l in sector_label(s)], ", ")
+  return print(io, "SU(", groupdim(s), ")[", disp, "]")
+end
+
+# display SU(N) irrep as a Young tableau with utf8 box char
+function Base.show(io::IO, ::MIME"text/plain", s::SU)
+  if istrivial(s)   # singlet = no box
+    return print(io, "●")
+  end
+
+  N = groupdim(s)
+  l = sector_label(s)
+  println(io, "┌─" * "┬─"^(l[1] - 1) * "┐")
+  i = 1
+  while i < N - 1 && l[i + 1] != 0
+    println(
+      io,
+      "├─",
+      "┼─"^(l[i + 1] - 1 + (l[i] > l[i + 1])),
+      "┴─"^max(0, (l[i] - l[i + 1] - 1)),
+      "┤"^(l[i] == l[i + 1]),
+      "┘"^(l[i] > l[i + 1]),
+    )
+    i += 1
+  end
+
+  print(io, "└─", "┴─"^max(0, l[i] - 1), "┘")
+  return nothing
+end
+
+#
+# Specializations for the case SU{2}
+#
+
+# optimize implementation
+quantum_dimension(s::SU{2}) = sector_label(s)[1] + 1
+
+GradedAxes.dual(s::SU{2}) = s
+
+function label_fusion_rule(::Type{<:SU{2}}, s1, s2)
+  irreps = [SU{2}((i,)) for i in (abs(s1[1] - s2[1])):2:(s1[1] + s2[1])]
+  degen = ones(Int, length(irreps))
+  return irreps .=> degen
+end
+
+# define angular momentum-like interface using half-integers
+SU2(h::Number) = SU{2}((twice(HalfInteger(h)),))
+
+# display SU2 using half-integers
+function Base.show(io::IO, s::SU{2})
+  return print(io, "SU(2)[S=", half(quantum_dimension(s) - 1), "]")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", s::SU{2})
+  return print(io, "S = ", half(quantum_dimension(s) - 1))
+end
+
+#
+# Specializations for the case SU{3}
+# aimed for testing non-abelian non self-conjugate representations
+# TODO replace with generic implementation
+#
+
+function label_fusion_rule(::Type{<:SU{3}}, left, right)
+  # Compute SU(3) fusion rules using Littlewood-Richardson rule for Young tableaus.
+  # See e.g. Di Francesco, Mathieu and Sénéchal, section 13.5.3.
+  if sum(right) > sum(left)  # impose more boxes in left Young tableau
+    return label_fusion_rule(SU{3}, right, left)
+  end
+
+  if right[1] == 0  # avoid issues with singlet
+    return [SU{3}(left) => 1]
+  end
+
+  left_row1 = left[1]
+  left_row2 = left[2]
+  right_row1 = right[1]
+  right_row2 = right[2]
+
+  irreps = []
+
+  # put a23 boxes on 2nd or 3rd line
+  a23max1 = 2 * left_row1  # row2a <= row1a
+  a23max2 = right_row1  # a2 + a3 <= total number of a
+  a23max = min(a23max1, a23max2)
+  for a23 in 0:a23max
+    a3min1 = left_row2 + 2 * a23 - left_row1 - right_row1
+    a3min2 = left_row2 - left_row1 + a23  # no a below a: row2a <= row1
+    a3min = max(0, a3min1, a3min2)
+    a3max1 = left_row2  # row3a <= row2a
+    a3max2 = a23  # a3 <= a2 + a3
+    a3max3 = right_row1 - right_row2  # more a than b, right to left: b2 + b3 <= a1 + a2
+    a3max = min(a3max1, a3max2, a3max3)
+    for a3 in a3min:a3max
+      a2 = a23 - a3
+      row1a = left_row1 + right_row1 - a23
+      row2a = left_row2 + a23 - a3
+
+      # cannot put any b on 1st line: row1ab = row1a
+      b3min1 = row2a + right_row2 - row1a  # row2ab <= row1ab = row1a
+      b3min2 = right_row2 + a23 - right_row1
+      b3min = max(0, b3min1, b3min2)
+      b3max1 = right_row2  # only other.row2 b boxes
+      b3max2 = (row2a + right_row2 - a3) ÷ 2  # row3ab >= row2ab
+      b3max3 = right_row1 - a3  # more a than b, right to left: b2 <= a1
+      b3max4 = row2a - a3  # no b below b: row2a >= row3ab
+      b3max = min(b3max1, b3max2, b3max3, b3max4)
+      for b3 in b3min:b3max
+        b2 = right_row2 - b3
+        row2ab = row2a + b2
+        row3ab = a3 + b3
+        yt = (row1a - row3ab, row2ab - row3ab)
+
+        push!(irreps, yt)
+      end
+    end
+  end
+
+  unique_labels = sort(unique(irreps))
+  degen = [count(==(irr), irreps) for irr in unique_labels]
+  sectors = SU{3}.(unique_labels)
+  return sectors .=> degen
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su2k.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su2k.jl
new file mode 100644
index 0000000000..cdf9bfbb1a
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/su2k.jl
@@ -0,0 +1,27 @@
+#
+# Quantum 'group' su2ₖ
+#
+
+using HalfIntegers: Half
+using ...GradedAxes: GradedAxes
+
+struct su2{k} <: AbstractSector
+  j::Half{Int}
+end
+
+SymmetryStyle(::Type{<:su2}) = NotAbelianStyle()
+
+GradedAxes.dual(s::su2) = s
+
+sector_label(s::su2) = s.j
+
+level(::su2{k}) where {k} = k
+
+trivial(::Type{su2{k}}) where {k} = su2{k}(0)
+
+function label_fusion_rule(::Type{su2{k}}, j1, j2) where {k}
+  labels = collect(abs(j1 - j2):min(k - j1 - j2, j1 + j2))
+  degen = ones(Int, length(labels))
+  sectors = su2{k}.(labels)
+  return sectors .=> degen
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/trivial.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/trivial.jl
new file mode 100644
index 0000000000..27b3fec88e
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/trivial.jl
@@ -0,0 +1,38 @@
+#
+# Trivial sector
+# acts as a trivial sector for any AbstractSector
+#
+
+using ...GradedAxes: GradedAxes
+
+# Trivial is special as it does not have a label
+struct TrivialSector <: AbstractSector end
+
+SymmetryStyle(::Type{TrivialSector}) = AbelianStyle()
+
+trivial(::Type{TrivialSector}) = TrivialSector()
+
+GradedAxes.dual(::TrivialSector) = TrivialSector()
+
+# TrivialSector acts as trivial on any AbstractSector
+function fusion_rule(::NotAbelianStyle, ::TrivialSector, c::AbstractSector)
+  return to_gradedrange(c)
+end
+function fusion_rule(::NotAbelianStyle, c::AbstractSector, ::TrivialSector)
+  return to_gradedrange(c)
+end
+
+# abelian case: return Sector
+fusion_rule(::AbelianStyle, c::AbstractSector, ::TrivialSector) = c
+fusion_rule(::AbelianStyle, ::TrivialSector, c::AbstractSector) = c
+fusion_rule(::AbelianStyle, ::TrivialSector, ::TrivialSector) = TrivialSector()
+
+# any trivial sector equals TrivialSector
+Base.:(==)(c::AbstractSector, ::TrivialSector) = istrivial(c)
+Base.:(==)(::TrivialSector, c::AbstractSector) = istrivial(c)
+Base.:(==)(::TrivialSector, ::TrivialSector) = true
+
+# sorts as trivial for any Sector
+Base.isless(c::AbstractSector, ::TrivialSector) = c < trivial(c)
+Base.isless(::TrivialSector, c::AbstractSector) = trivial(c) < c
+Base.isless(::TrivialSector, ::TrivialSector) = false  # bypass default that calls label
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/u1.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/u1.jl
new file mode 100644
index 0000000000..2d79796c34
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/u1.jl
@@ -0,0 +1,31 @@
+#
+# U₁ group (circle group, or particle number, total Sz etc.)
+#
+
+using ...GradedAxes: GradedAxes
+
+# Parametric type to allow both integer label as well as
+# HalfInteger for easy conversion to/from SU(2)
+struct U1{T} <: AbstractSector
+  n::T
+end
+
+SymmetryStyle(::Type{<:U1}) = AbelianStyle()
+sector_label(u::U1) = u.n
+
+set_sector_label(s::U1, sector_label) = typeof(s)(sector_label)
+GradedAxes.dual(s::U1) = set_sector_label(s, -sector_label(s))
+
+trivial(::Type{U1}) = trivial(U1{Int})
+trivial(::Type{U1{T}}) where {T} = U1(zero(T))
+
+abelian_label_fusion_rule(sector_type::Type{<:U1}, n1, n2) = sector_type(n1 + n2)
+
+# hide label type in printing
+function Base.show(io::IO, u::U1)
+  return print(io, "U(1)[", sector_label(u), "]")
+end
+
+# enforce U1(Int32(1)) == U1(1)
+Base.:(==)(s1::U1, s2::U1) = sector_label(s1) == sector_label(s2)
+Base.isless(s1::U1, s2::U1) = sector_label(s1) < sector_label(s2)
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/zn.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/zn.jl
new file mode 100644
index 0000000000..8628288dc3
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_definitions/zn.jl
@@ -0,0 +1,25 @@
+#
+# Cyclic group Zₙ
+#
+
+using ...GradedAxes: GradedAxes
+
+struct Z{N} <: AbstractSector
+  m::Int
+  Z{N}(m) where {N} = new{N}(mod(m, N))
+end
+
+modulus(::Type{Z{N}}) where {N} = N
+modulus(c::Z) = modulus(typeof(c))
+
+SymmetryStyle(::Type{<:Z}) = AbelianStyle()
+sector_label(c::Z) = c.m
+
+set_sector_label(s::Z, sector_label) = typeof(s)(sector_label)
+GradedAxes.dual(s::Z) = set_sector_label(s, -sector_label(s))
+
+trivial(sector_type::Type{<:Z}) = sector_type(0)
+
+function abelian_label_fusion_rule(sector_type::Type{<:Z}, n1, n2)
+  return sector_type(n1 + n2)
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/sector_product.jl b/NDTensors/src/lib/SymmetrySectors/src/sector_product.jl
new file mode 100644
index 0000000000..12a72d022e
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/sector_product.jl
@@ -0,0 +1,238 @@
+# This files defines a structure for Cartesian product of 2 or more fusion sectors
+# e.g. U(1)×U(1), U(1)×SU2(2)×SU(3)
+
+using BlockArrays: blocklengths
+using ..LabelledNumbers: LabelledInteger, label, labelled, unlabel
+using ..GradedAxes: AbstractGradedUnitRange, GradedAxes, dual
+
+# =====================================  Definition  =======================================
+struct SectorProduct{Sectors} <: AbstractSector
+  arguments::Sectors
+  global _SectorProduct(l) = new{typeof(l)}(l)
+end
+
+SectorProduct(c::SectorProduct) = _SectorProduct(arguments(c))
+
+arguments(s::SectorProduct) = s.arguments
+
+# =================================  Sectors interface  ====================================
+function SymmetryStyle(T::Type{<:SectorProduct})
+  return arguments_symmetrystyle(arguments_type(T))
+end
+
+function quantum_dimension(::NotAbelianStyle, s::SectorProduct)
+  return mapreduce(quantum_dimension, *, arguments(s))
+end
+
+# use map instead of broadcast to support both Tuple and NamedTuple
+GradedAxes.dual(s::SectorProduct) = SectorProduct(map(dual, arguments(s)))
+
+function trivial(type::Type{<:SectorProduct})
+  return SectorProduct(arguments_trivial(arguments_type(type)))
+end
+
+# ===================================  Base interface  =====================================
+function Base.:(==)(A::SectorProduct, B::SectorProduct)
+  return arguments_isequal(arguments(A), arguments(B))
+end
+
+function Base.show(io::IO, s::SectorProduct)
+  (length(arguments(s)) < 2) && print(io, "sector")
+  print(io, "(")
+  symbol = ""
+  for p in pairs(arguments(s))
+    print(io, symbol)
+    sector_show(io, p[1], p[2])
+    symbol = " × "
+  end
+  return print(io, ")")
+end
+
+sector_show(io::IO, ::Int, v) = print(io, v)
+sector_show(io::IO, k::Symbol, v) = print(io, "($k=$v,)")
+
+function Base.isless(s1::SectorProduct, s2::SectorProduct)
+  return arguments_isless(arguments(s1), arguments(s2))
+end
+
+# =======================================  shared  =========================================
+# there are 2 implementations for SectorProduct
+# - ordered-like with a Tuple
+# - dictionary-like with a NamedTuple
+
+arguments_type(::Type{<:SectorProduct{T}}) where {T} = T
+
+arguments_maybe_insert_unspecified(s1, ::Any) = s1
+function sym_arguments_maybe_insert_unspecified(s1, s2)
+  return arguments_maybe_insert_unspecified(s1, s2),
+  arguments_maybe_insert_unspecified(s2, s1)
+end
+
+function make_empty_match(a1, b1)
+  a2 = isempty(a1) ? empty(b1) : a1
+  b2 = isempty(b1) ? empty(a2) : b1
+  return a2, b2
+end
+
+function arguments_isequal(a1, b1)
+  return ==(sym_arguments_maybe_insert_unspecified(make_empty_match(a1, b1)...)...)
+end
+
+function arguments_product(s1, s2)
+  isempty(s1) && return s2
+  isempty(s2) && return s1
+  return throw(ArgumentError("Mixing non-empty storage types is illegal"))
+end
+
+function arguments_isless(a1, b1)
+  return isless(sym_arguments_maybe_insert_unspecified(make_empty_match(a1, b1)...)...)
+end
+
+# =================================  Cartesian Product  ====================================
+×(c1::AbstractSector, c2::AbstractSector) = ×(SectorProduct(c1), SectorProduct(c2))
+function ×(p1::SectorProduct, p2::SectorProduct)
+  return SectorProduct(arguments_product(arguments(p1), arguments(p2)))
+end
+
+×(a, g::AbstractUnitRange) = ×(to_gradedrange(a), g)
+×(g::AbstractUnitRange, b) = ×(g, to_gradedrange(b))
+×(nt1::NamedTuple, nt2::NamedTuple) = ×(SectorProduct(nt1), SectorProduct(nt2))
+×(c1::NamedTuple, c2::AbstractSector) = ×(SectorProduct(c1), SectorProduct(c2))
+×(c1::AbstractSector, c2::NamedTuple) = ×(SectorProduct(c1), SectorProduct(c2))
+
+function ×(l1::LabelledInteger, l2::LabelledInteger)
+  c3 = label(l1) × label(l2)
+  m3 = unlabel(l1) * unlabel(l2)
+  return labelled(m3, c3)
+end
+
+function ×(g1::AbstractUnitRange, g2::AbstractUnitRange)
+  v = map(
+    ((l1, l2),) -> l1 × l2,
+    Iterators.flatten((Iterators.product(blocklengths(g1), blocklengths(g2)),),),
+  )
+  return gradedrange(v)
+end
+
+# ====================================  Fusion rules  ======================================
+# cast AbstractSector to SectorProduct
+function fusion_rule(style::SymmetryStyle, c1::SectorProduct, c2::AbstractSector)
+  return fusion_rule(style, c1, SectorProduct(c2))
+end
+function fusion_rule(style::SymmetryStyle, c1::AbstractSector, c2::SectorProduct)
+  return fusion_rule(style, SectorProduct(c1), c2)
+end
+
+# generic case: fusion returns a GradedAxes, even for fusion with Empty
+function fusion_rule(::NotAbelianStyle, s1::SectorProduct, s2::SectorProduct)
+  return to_gradedrange(arguments_fusion_rule(arguments(s1), arguments(s2)))
+end
+
+# Abelian case: fusion returns SectorProduct
+function fusion_rule(::AbelianStyle, s1::SectorProduct, s2::SectorProduct)
+  return label(only(fusion_rule(NotAbelianStyle(), s1, s2)))
+end
+
+# lift ambiguities for TrivialSector
+fusion_rule(::AbelianStyle, c::SectorProduct, ::TrivialSector) = c
+fusion_rule(::AbelianStyle, ::TrivialSector, c::SectorProduct) = c
+fusion_rule(::NotAbelianStyle, c::SectorProduct, ::TrivialSector) = to_gradedrange(c)
+fusion_rule(::NotAbelianStyle, ::TrivialSector, c::SectorProduct) = to_gradedrange(c)
+
+function arguments_fusion_rule(sects1, sects2)
+  isempty(sects1) && return SectorProduct(sects2)
+  isempty(sects2) && return SectorProduct(sects1)
+  shared_sect = shared_arguments_fusion_rule(arguments_common(sects1, sects2)...)
+  diff_sect = SectorProduct(arguments_diff(sects1, sects2))
+  return shared_sect × diff_sect
+end
+
+# ===============================  Ordered implementation  =================================
+SectorProduct(t::Tuple) = _SectorProduct(t)
+SectorProduct(sects::AbstractSector...) = SectorProduct(sects)
+
+function arguments_symmetrystyle(T::Type{<:Tuple})
+  return mapreduce(SymmetryStyle, combine_styles, fieldtypes(T); init=AbelianStyle())
+end
+
+arguments_product(l1::Tuple, l2::Tuple) = (l1..., l2...)
+
+arguments_trivial(T::Type{<:Tuple}) = trivial.(fieldtypes(T))
+
+function arguments_common(t1::Tuple, t2::Tuple)
+  n = min(length(t1), length(t2))
+  return t1[begin:n], t2[begin:n]
+end
+
+function arguments_diff(t1::Tuple, t2::Tuple)
+  n1 = length(t1)
+  n2 = length(t2)
+  return n1 < n2 ? t2[(n1 + 1):end] : t1[(n2 + 1):end]
+end
+
+function shared_arguments_fusion_rule(shared1::T, shared2::T) where {T<:Tuple}
+  return mapreduce(
+    to_gradedrange ∘ fusion_rule,
+    ×,
+    shared1,
+    shared2;
+    init=to_gradedrange(SectorProduct(())),
+  )
+end
+
+function arguments_maybe_insert_unspecified(t1::Tuple, t2::Tuple)
+  n1 = length(t1)
+  return (t1..., trivial.(t2[(n1 + 1):end])...)
+end
+
+# ===========================  Dictionary-like implementation  =============================
+function SectorProduct(nt::NamedTuple)
+  arguments = sort_keys(nt)
+  return _SectorProduct(arguments)
+end
+
+SectorProduct(; kws...) = SectorProduct((; kws...))
+
+function SectorProduct(pairs::Pair...)
+  keys = Symbol.(first.(pairs))
+  vals = last.(pairs)
+  return SectorProduct(NamedTuple{keys}(vals))
+end
+
+function arguments_symmetrystyle(NT::Type{<:NamedTuple})
+  return mapreduce(SymmetryStyle, combine_styles, fieldtypes(NT); init=AbelianStyle())
+end
+
+function arguments_maybe_insert_unspecified(nt1::NamedTuple, nt2::NamedTuple)
+  diff1 = arguments_trivial(typeof(setdiff_keys(nt2, nt1)))
+  return sort_keys(union_keys(nt1, diff1))
+end
+
+function arguments_product(l1::NamedTuple, l2::NamedTuple)
+  if length(intersect_keys(l1, l2)) > 0
+    throw(ArgumentError("Cannot define product of shared keys"))
+  end
+  return union_keys(l1, l2)
+end
+
+function arguments_trivial(NT::Type{<:NamedTuple{Keys}}) where {Keys}
+  return NamedTuple{Keys}(trivial.(fieldtypes(NT)))
+end
+
+function arguments_common(nt1::NamedTuple, nt2::NamedTuple)
+  # SectorProduct(nt::NamedTuple) sorts keys at init
+  @assert issorted(keys(nt1))
+  @assert issorted(keys(nt2))
+  return intersect_keys(nt1, nt2), intersect_keys(nt2, nt1)
+end
+
+arguments_diff(nt1::NamedTuple, nt2::NamedTuple) = symdiff_keys(nt1, nt2)
+
+function map_blocklabels(f, r::AbstractUnitRange)
+  return gradedrange(labelled.(unlabel.(blocklengths(r)), f.(blocklabels(r))))
+end
+
+function shared_arguments_fusion_rule(shared1::NT, shared2::NT) where {NT<:NamedTuple}
+  tuple_fused = shared_arguments_fusion_rule(values(shared1), values(shared2))
+  return map_blocklabels(SectorProduct ∘ NT ∘ arguments ∘ SectorProduct, tuple_fused)
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/src/symmetry_style.jl b/NDTensors/src/lib/SymmetrySectors/src/symmetry_style.jl
new file mode 100644
index 0000000000..f3e443dce0
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/src/symmetry_style.jl
@@ -0,0 +1,17 @@
+# This file defines SymmetryStyle, a trait to distinguish abelian groups, non-abelian groups
+# and non-group fusion categories.
+
+using ..LabelledNumbers: LabelledInteger, label_type
+
+abstract type SymmetryStyle end
+
+struct AbelianStyle <: SymmetryStyle end
+struct NotAbelianStyle <: SymmetryStyle end
+
+SymmetryStyle(x) = SymmetryStyle(typeof(x))
+SymmetryStyle(T::Type) = error("method `SymmetryStyle` not defined for type $(T)")
+SymmetryStyle(L::Type{<:LabelledInteger}) = SymmetryStyle(label_type(L))
+SymmetryStyle(G::Type{<:AbstractUnitRange}) = SymmetryStyle(eltype(G))
+
+combine_styles(::AbelianStyle, ::AbelianStyle) = AbelianStyle()
+combine_styles(::SymmetryStyle, ::SymmetryStyle) = NotAbelianStyle()
diff --git a/NDTensors/src/lib/Sectors/test/runtests.jl b/NDTensors/src/lib/SymmetrySectors/test/runtests.jl
similarity index 100%
rename from NDTensors/src/lib/Sectors/test/runtests.jl
rename to NDTensors/src/lib/SymmetrySectors/test/runtests.jl
diff --git a/NDTensors/src/lib/SymmetrySectors/test/test_fusion_rules.jl b/NDTensors/src/lib/SymmetrySectors/test/test_fusion_rules.jl
new file mode 100644
index 0000000000..bd00abef86
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/test/test_fusion_rules.jl
@@ -0,0 +1,287 @@
+@eval module $(gensym())
+using NDTensors.GradedAxes:
+  dual, fusion_product, space_isequal, gradedrange, flip, tensor_product
+using NDTensors.SymmetrySectors:
+  ⊗,
+  Fib,
+  Ising,
+  O2,
+  SU,
+  SU2,
+  TrivialSector,
+  U1,
+  Z,
+  block_dimensions,
+  quantum_dimension,
+  trivial
+using Test: @inferred, @test, @testset, @test_throws
+
+@testset "Simple SymmetrySector fusion rules" begin
+  @testset "Z{2} fusion rules" begin
+    z0 = Z{2}(0)
+    z1 = Z{2}(1)
+
+    @test z0 ⊗ z0 == z0
+    @test z0 ⊗ z1 == z1
+    @test z1 ⊗ z1 == z0
+    @test (@inferred z0 ⊗ z0) == z0  # no better way, see Julia PR 23426
+
+    q = TrivialSector()
+    @test (@inferred q ⊗ q) == q
+    @test (@inferred q ⊗ z0) == z0
+    @test (@inferred z1 ⊗ q) == z1
+
+    # using GradedAxes interface
+    @test space_isequal(fusion_product(z0, z0), gradedrange([z0 => 1]))
+    @test space_isequal(fusion_product(z0, z1), gradedrange([z1 => 1]))
+
+    # test different input number
+    @test space_isequal(fusion_product(z0), gradedrange([z0 => 1]))
+    @test space_isequal(fusion_product(z0, z0, z0), gradedrange([z0 => 1]))
+    @test space_isequal(fusion_product(z0, z0, z0, z0), gradedrange([z0 => 1]))
+    @test (@inferred block_dimensions(gradedrange([z1 => 1]))) == [1]
+  end
+  @testset "U(1) fusion rules" begin
+    q1 = U1(1)
+    q2 = U1(2)
+    q3 = U1(3)
+
+    @test q1 ⊗ q1 == U1(2)
+    @test q1 ⊗ q2 == U1(3)
+    @test q2 ⊗ q1 == U1(3)
+    @test (@inferred q1 ⊗ q2) == q3  # no better way, see Julia PR 23426
+  end
+
+  @testset "O2 fusion rules" begin
+    s0e = O2(0)
+    s0o = O2(-1)
+    s12 = O2(1//2)
+    s1 = O2(1)
+
+    q = TrivialSector()
+    @test space_isequal((@inferred s0e ⊗ q), gradedrange([s0e => 1]))
+    @test space_isequal((@inferred q ⊗ s0o), gradedrange([s0o => 1]))
+
+    @test space_isequal((@inferred s0e ⊗ s0e), gradedrange([s0e => 1]))
+    @test space_isequal((@inferred s0o ⊗ s0e), gradedrange([s0o => 1]))
+    @test space_isequal((@inferred s0o ⊗ s0e), gradedrange([s0o => 1]))
+    @test space_isequal((@inferred s0o ⊗ s0o), gradedrange([s0e => 1]))
+
+    @test space_isequal((@inferred s0e ⊗ s12), gradedrange([s12 => 1]))
+    @test space_isequal((@inferred s0o ⊗ s12), gradedrange([s12 => 1]))
+    @test space_isequal((@inferred s12 ⊗ s0e), gradedrange([s12 => 1]))
+    @test space_isequal((@inferred s12 ⊗ s0o), gradedrange([s12 => 1]))
+    @test space_isequal((@inferred s12 ⊗ s1), gradedrange([s12 => 1, O2(3//2) => 1]))
+    @test space_isequal((@inferred s12 ⊗ s12), gradedrange([s0o => 1, s0e => 1, s1 => 1]))
+
+    @test (@inferred quantum_dimension(s0o ⊗ s1)) == 2
+    @test (@inferred block_dimensions(s0o ⊗ s1)) == [2]
+  end
+
+  @testset "SU2 fusion rules" begin
+    j1 = SU2(0)
+    j2 = SU2(1//2)
+    j3 = SU2(1)
+    j4 = SU2(3//2)
+    j5 = SU2(2)
+
+    @test space_isequal(j1 ⊗ j2, gradedrange([j2 => 1]))
+    @test space_isequal(j2 ⊗ j2, gradedrange([j1 => 1, j3 => 1]))
+    @test space_isequal(j2 ⊗ j3, gradedrange([j2 => 1, j4 => 1]))
+    @test space_isequal(j3 ⊗ j3, gradedrange([j1 => 1, j3 => 1, j5 => 1]))
+    @test space_isequal((@inferred j1 ⊗ j2), gradedrange([j2 => 1]))
+    @test (@inferred quantum_dimension(j1 ⊗ j2)) == 2
+    @test (@inferred block_dimensions(j1 ⊗ j2)) == [2]
+
+    @test space_isequal(fusion_product(j2), gradedrange([j2 => 1]))
+    @test space_isequal(fusion_product(j2, j1), gradedrange([j2 => 1]))
+    @test space_isequal(fusion_product(j2, j1, j1), gradedrange([j2 => 1]))
+  end
+
+  @testset "Fibonacci fusion rules" begin
+    ı = Fib("1")
+    τ = Fib("τ")
+
+    @test space_isequal(ı ⊗ ı, gradedrange([ı => 1]))
+    @test space_isequal(ı ⊗ τ, gradedrange([τ => 1]))
+    @test space_isequal(τ ⊗ ı, gradedrange([τ => 1]))
+    @test space_isequal((@inferred τ ⊗ τ), gradedrange([ı => 1, τ => 1]))
+    @test (@inferred quantum_dimension(gradedrange([ı => 1, ı => 1]))) == 2.0
+  end
+
+  @testset "Ising fusion rules" begin
+    ı = Ising("1")
+    σ = Ising("σ")
+    ψ = Ising("ψ")
+
+    @test space_isequal(ı ⊗ ı, gradedrange([ı => 1]))
+    @test space_isequal(ı ⊗ σ, gradedrange([σ => 1]))
+    @test space_isequal(σ ⊗ ı, gradedrange([σ => 1]))
+    @test space_isequal(ı ⊗ ψ, gradedrange([ψ => 1]))
+    @test space_isequal(ψ ⊗ ı, gradedrange([ψ => 1]))
+    @test space_isequal(σ ⊗ σ, gradedrange([ı => 1, ψ => 1]))
+    @test space_isequal(σ ⊗ ψ, gradedrange([σ => 1]))
+    @test space_isequal(ψ ⊗ σ, gradedrange([σ => 1]))
+    @test space_isequal(ψ ⊗ ψ, gradedrange([ı => 1]))
+    @test space_isequal((@inferred ψ ⊗ ψ), gradedrange([ı => 1]))
+    @test (@inferred quantum_dimension(σ ⊗ σ)) == 2.0
+  end
+end
+@testset "Gradedrange fusion rules" begin
+  @testset "Trivial GradedUnitRange" begin
+    g1 = gradedrange([U1(0) => 1])
+    g2 = gradedrange([SU2(0) => 1])
+    @test space_isequal(trivial(g1), g1)
+    @test space_isequal(trivial(dual(g1)), g1)  # trivial returns nondual
+    @test space_isequal(trivial(typeof(g2)), g2)
+  end
+  @testset "GradedUnitRange abelian tensor/fusion product" begin
+    g1 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 2])
+    g2 = gradedrange([U1(-2) => 2, U1(0) => 1, U1(1) => 2])
+
+    @test space_isequal(flip(dual(g1)), gradedrange([U1(1) => 1, U1(0) => 1, U1(-1) => 2]))
+    @test (@inferred block_dimensions(g1)) == [1, 1, 2]
+
+    gt = gradedrange([
+      U1(-3) => 2,
+      U1(-2) => 2,
+      U1(-1) => 4,
+      U1(-1) => 1,
+      U1(0) => 1,
+      U1(1) => 2,
+      U1(0) => 2,
+      U1(1) => 2,
+      U1(2) => 4,
+    ])
+    gf = gradedrange([
+      U1(-3) => 2, U1(-2) => 2, U1(-1) => 5, U1(0) => 3, U1(1) => 4, U1(2) => 4
+    ])
+    @test space_isequal((@inferred tensor_product(g1, g2)), gt)
+    @test space_isequal((@inferred fusion_product(g1, g2)), gf)
+
+    gtd1 = gradedrange([
+      U1(-1) => 2,
+      U1(-2) => 2,
+      U1(-3) => 4,
+      U1(1) => 1,
+      U1(0) => 1,
+      U1(-1) => 2,
+      U1(2) => 2,
+      U1(1) => 2,
+      U1(0) => 4,
+    ])
+    gfd1 = gradedrange([
+      U1(-3) => 4, U1(-2) => 2, U1(-1) => 4, U1(0) => 5, U1(1) => 3, U1(2) => 2
+    ])
+    @test space_isequal((@inferred tensor_product(dual(g1), g2)), gtd1)
+    @test space_isequal((@inferred fusion_product(dual(g1), g2)), gfd1)
+
+    gtd2 = gradedrange([
+      U1(1) => 2,
+      U1(2) => 2,
+      U1(3) => 4,
+      U1(-1) => 1,
+      U1(0) => 1,
+      U1(1) => 2,
+      U1(-2) => 2,
+      U1(-1) => 2,
+      U1(0) => 4,
+    ])
+    gfd2 = gradedrange([
+      U1(-2) => 2, U1(-1) => 3, U1(0) => 5, U1(1) => 4, U1(2) => 2, U1(3) => 4
+    ])
+    @test space_isequal((@inferred tensor_product(g1, dual(g2))), gtd2)
+    @test space_isequal((@inferred fusion_product(g1, dual(g2))), gfd2)
+
+    gtd = gradedrange([
+      U1(3) => 2,
+      U1(2) => 2,
+      U1(1) => 4,
+      U1(1) => 1,
+      U1(0) => 1,
+      U1(-1) => 2,
+      U1(0) => 2,
+      U1(-1) => 2,
+      U1(-2) => 4,
+    ])
+    gfd = gradedrange([
+      U1(-2) => 4, U1(-1) => 4, U1(0) => 3, U1(1) => 5, U1(2) => 2, U1(3) => 2
+    ])
+    @test space_isequal((@inferred tensor_product(dual(g1), dual(g2))), gtd)
+    @test space_isequal((@inferred fusion_product(dual(g1), dual(g2))), gfd)
+
+    # test different (non-product) sectors cannot be fused
+    @test_throws MethodError fusion_product(gradedrange([Z{2}(0) => 1]), g1)
+    @test_throws MethodError tensor_product(gradedrange([Z{2}(0) => 1]), g2)
+  end
+
+  @testset "GradedUnitRange non-abelian fusion rules" begin
+    g3 = gradedrange([SU2(0) => 1, SU2(1//2) => 2, SU2(1) => 1])
+    g4 = gradedrange([SU2(1//2) => 1, SU2(1) => 2])
+    g34 = gradedrange([
+      SU2(1//2) => 1,
+      SU2(0) => 2,
+      SU2(1) => 2,
+      SU2(1//2) => 1,
+      SU2(3//2) => 1,
+      SU2(1) => 2,
+      SU2(1//2) => 4,
+      SU2(3//2) => 4,
+      SU2(0) => 2,
+      SU2(1) => 2,
+      SU2(2) => 2,
+    ])
+
+    @test space_isequal(tensor_product(g3, g4), g34)
+
+    @test space_isequal(dual(flip(g3)), g3)  # trivial for SU(2)
+    @test space_isequal(
+      (@inferred fusion_product(g3, g4)),
+      gradedrange([SU2(0) => 4, SU2(1//2) => 6, SU2(1) => 6, SU2(3//2) => 5, SU2(2) => 2]),
+    )
+    @test (@inferred block_dimensions(g3)) == [1, 4, 3]
+
+    # test dual on non self-conjugate non-abelian representations
+    s1 = SU{3}((0, 0))
+    f3 = SU{3}((1, 0))
+    c3 = SU{3}((1, 1))
+    ad8 = SU{3}((2, 1))
+
+    g5 = gradedrange([s1 => 1, f3 => 1])
+    g6 = gradedrange([s1 => 1, c3 => 1])
+    @test space_isequal(dual(flip(g5)), g6)
+    @test space_isequal(
+      fusion_product(g5, g6), gradedrange([s1 => 2, f3 => 1, c3 => 1, ad8 => 1])
+    )
+    @test space_isequal(
+      fusion_product(dual(g5), g6),
+      gradedrange([s1 => 1, f3 => 1, c3 => 2, SU{3}((2, 2)) => 1]),
+    )
+    @test space_isequal(
+      fusion_product(g5, dual(g6)),
+      gradedrange([s1 => 1, f3 => 2, c3 => 1, SU{3}((2, 0)) => 1]),
+    )
+    @test space_isequal(
+      fusion_product(dual(g5), dual(g6)), gradedrange([s1 => 2, f3 => 1, c3 => 1, ad8 => 1])
+    )
+  end
+
+  @testset "Mixed GradedUnitRange - Sector fusion rules" begin
+    g1 = gradedrange([U1(1) => 1, U1(2) => 2])
+    g2 = gradedrange([U1(2) => 1, U1(3) => 2])
+    @test space_isequal((@inferred fusion_product(g1, U1(1))), g2)
+    @test space_isequal((@inferred fusion_product(U1(1), g1)), g2)
+
+    g3 = gradedrange([SU2(0) => 1, SU2(1//2) => 2])
+    g4 = gradedrange([SU2(0) => 2, SU2(1//2) => 1, SU2(1) => 2])
+    @test space_isequal((@inferred fusion_product(g3, SU2(1//2))), g4)
+    @test space_isequal((@inferred fusion_product(SU2(1//2), g3)), g4)
+
+    # test different simple sectors cannot be fused
+    @test_throws MethodError Z{2}(0) ⊗ U1(1)
+    @test_throws MethodError SU2(1) ⊗ U1(1)
+    @test_throws MethodError fusion_product(g1, SU2(1))
+    @test_throws MethodError fusion_product(U1(1), g3)
+  end
+end
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/test/test_sector_product.jl b/NDTensors/src/lib/SymmetrySectors/test/test_sector_product.jl
new file mode 100644
index 0000000000..de046fc307
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/test/test_sector_product.jl
@@ -0,0 +1,622 @@
+@eval module $(gensym())
+using NDTensors.SymmetrySectors:
+  ×,
+  ⊗,
+  Fib,
+  Ising,
+  SectorProduct,
+  SU,
+  SU2,
+  TrivialSector,
+  U1,
+  Z,
+  block_dimensions,
+  quantum_dimension,
+  arguments,
+  trivial
+using NDTensors.GradedAxes: dual, fusion_product, space_isequal, gradedrange
+using Test: @inferred, @test, @testset, @test_throws
+
+@testset "Test Ordered Products" begin
+  @testset "Ordered Constructor" begin
+    s = SectorProduct(U1(1))
+    @test length(arguments(s)) == 1
+    @test (@inferred quantum_dimension(s)) == 1
+    @test (@inferred dual(s)) == SectorProduct(U1(-1))
+    @test arguments(s)[1] == U1(1)
+    @test (@inferred trivial(s)) == SectorProduct(U1(0))
+
+    s = SectorProduct(U1(1), U1(2))
+    @test length(arguments(s)) == 2
+    @test (@inferred quantum_dimension(s)) == 1
+    @test (@inferred dual(s)) == SectorProduct(U1(-1), U1(-2))
+    @test arguments(s)[1] == U1(1)
+    @test arguments(s)[2] == U1(2)
+    @test (@inferred trivial(s)) == SectorProduct(U1(0), U1(0))
+
+    s = U1(1) × SU2(1//2) × U1(3)
+    @test length(arguments(s)) == 3
+    @test (@inferred quantum_dimension(s)) == 2
+    @test (@inferred dual(s)) == U1(-1) × SU2(1//2) × U1(-3)
+    @test arguments(s)[1] == U1(1)
+    @test arguments(s)[2] == SU2(1//2)
+    @test arguments(s)[3] == U1(3)
+    @test (@inferred trivial(s)) == SectorProduct(U1(0), SU2(0), U1(0))
+
+    s = U1(3) × SU2(1//2) × Fib("τ")
+    @test length(arguments(s)) == 3
+    @test (@inferred quantum_dimension(s)) == 1.0 + √5
+    @test dual(s) == U1(-3) × SU2(1//2) × Fib("τ")
+    @test arguments(s)[1] == U1(3)
+    @test arguments(s)[2] == SU2(1//2)
+    @test arguments(s)[3] == Fib("τ")
+    @test (@inferred trivial(s)) == SectorProduct(U1(0), SU2(0), Fib("1"))
+
+    s = TrivialSector() × U1(3) × SU2(1 / 2)
+    @test length(arguments(s)) == 3
+    @test (@inferred quantum_dimension(s)) == 2
+    @test dual(s) == TrivialSector() × U1(-3) × SU2(1//2)
+    @test (@inferred trivial(s)) == SectorProduct(TrivialSector(), U1(0), SU2(0))
+    @test s > trivial(s)
+  end
+
+  @testset "Ordered comparisons" begin
+    # convention: missing arguments are filled with singlets
+    @test SectorProduct(U1(1), SU2(1)) == SectorProduct(U1(1), SU2(1))
+    @test SectorProduct(U1(1), SU2(0)) != SectorProduct(U1(1), SU2(1))
+    @test SectorProduct(U1(0), SU2(1)) != SectorProduct(U1(1), SU2(1))
+    @test SectorProduct(U1(1)) != U1(1)
+    @test SectorProduct(U1(1)) == SectorProduct(U1(1), U1(0))
+    @test SectorProduct(U1(1)) != SectorProduct(U1(1), U1(1))
+    @test SectorProduct(U1(0), SU2(0)) == TrivialSector()
+    @test SectorProduct(U1(0), SU2(0)) == SectorProduct(TrivialSector(), SU2(0))
+    @test SectorProduct(U1(0), SU2(0)) == SectorProduct(U1(0), TrivialSector())
+    @test SectorProduct(U1(0), SU2(0)) == SectorProduct(TrivialSector(), TrivialSector())
+
+    @test SectorProduct(U1(0)) < SectorProduct((U1(1)))
+    @test SectorProduct(U1(0), U1(2)) < SectorProduct((U1(1)), U1(0))
+    @test SectorProduct(U1(0)) < SectorProduct(U1(0), U1(1))
+    @test SectorProduct(U1(0)) > SectorProduct(U1(0), U1(-1))
+  end
+
+  @testset "Quantum dimension and GradedUnitRange" begin
+    g = gradedrange([(U1(0) × Z{2}(0)) => 1, (U1(1) × Z{2}(0)) => 2])  # abelian
+    @test (@inferred quantum_dimension(g)) == 3
+
+    g = gradedrange([  # non-abelian
+      (SU2(0) × SU2(0)) => 1,
+      (SU2(1) × SU2(0)) => 1,
+      (SU2(0) × SU2(1)) => 1,
+      (SU2(1) × SU2(1)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 16
+    @test (@inferred block_dimensions(g)) == [1, 3, 3, 9]
+
+    # mixed group
+    g = gradedrange([(U1(2) × SU2(0) × Z{2}(0)) => 1, (U1(2) × SU2(1) × Z{2}(0)) => 1])
+    @test (@inferred quantum_dimension(g)) == 4
+    @test (@inferred block_dimensions(g)) == [1, 3]
+    g = gradedrange([(SU2(0) × U1(0) × SU2(1//2)) => 1, (SU2(0) × U1(1) × SU2(1//2)) => 1])
+    @test (@inferred quantum_dimension(g)) == 4
+    @test (@inferred block_dimensions(g)) == [2, 2]
+
+    # NonGroupCategory
+    g_fib = gradedrange([(Fib("1") × Fib("1")) => 1])
+    g_ising = gradedrange([(Ising("1") × Ising("1")) => 1])
+    @test (@inferred quantum_dimension((Fib("1") × Fib("1")))) == 1.0
+    @test (@inferred quantum_dimension(g_fib)) == 1.0
+    @test (@inferred quantum_dimension(g_ising)) == 1.0
+    @test (@inferred quantum_dimension((Ising("1") × Ising("1")))) == 1.0
+    @test (@inferred block_dimensions(g_fib)) == [1.0]
+    @test (@inferred block_dimensions(g_ising)) == [1.0]
+
+    @test (@inferred quantum_dimension(U1(1) × Fib("1"))) == 1.0
+    @test (@inferred quantum_dimension(gradedrange([U1(1) × Fib("1") => 1]))) == 1.0
+
+    # mixed product Abelian / NonAbelian / NonGroup
+    g = gradedrange([
+      (U1(2) × SU2(0) × Ising(1)) => 1,
+      (U1(2) × SU2(1) × Ising(1)) => 1,
+      (U1(2) × SU2(0) × Ising("ψ")) => 1,
+      (U1(2) × SU2(1) × Ising("ψ")) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 8.0
+    @test (@inferred block_dimensions(g)) == [1.0, 3.0, 1.0, 3.0]
+
+    ϕ = (1 + √5) / 2
+    g = gradedrange([
+      (Fib("1") × SU2(0) × U1(2)) => 1,
+      (Fib("1") × SU2(1) × U1(2)) => 1,
+      (Fib("τ") × SU2(0) × U1(2)) => 1,
+      (Fib("τ") × SU2(1) × U1(2)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 4.0 + 4.0ϕ
+    @test (@inferred block_dimensions(g)) == [1.0, 3.0, 1.0ϕ, 3.0ϕ]
+  end
+
+  @testset "Fusion of Abelian products" begin
+    p1 = SectorProduct(U1(1))
+    p2 = SectorProduct(U1(2))
+    @test (@inferred p1 ⊗ TrivialSector()) == p1
+    @test (@inferred TrivialSector() ⊗ p2) == p2
+    @test (@inferred p1 ⊗ p2) == SectorProduct(U1(3))
+
+    p11 = U1(1) × U1(1)
+    @test p11 ⊗ p11 == U1(2) × U1(2)
+
+    p123 = U1(1) × U1(2) × U1(3)
+    @test p123 ⊗ p123 == U1(2) × U1(4) × U1(6)
+
+    s1 = SectorProduct(U1(1), Z{2}(1))
+    s2 = SectorProduct(U1(0), Z{2}(0))
+    @test s1 ⊗ s2 == U1(1) × Z{2}(1)
+  end
+
+  @testset "Fusion of NonAbelian products" begin
+    p0 = SectorProduct(SU2(0))
+    ph = SectorProduct(SU2(1//2))
+    @test space_isequal(
+      (@inferred p0 ⊗ TrivialSector()), gradedrange([SectorProduct(SU2(0)) => 1])
+    )
+    @test space_isequal(
+      (@inferred TrivialSector() ⊗ ph), gradedrange([SectorProduct(SU2(1//2)) => 1])
+    )
+
+    phh = SU2(1//2) × SU2(1//2)
+    @test space_isequal(
+      phh ⊗ phh,
+      gradedrange([
+        (SU2(0) × SU2(0)) => 1,
+        (SU2(1) × SU2(0)) => 1,
+        (SU2(0) × SU2(1)) => 1,
+        (SU2(1) × SU2(1)) => 1,
+      ]),
+    )
+    @test space_isequal(
+      phh ⊗ phh,
+      gradedrange([
+        (SU2(0) × SU2(0)) => 1,
+        (SU2(1) × SU2(0)) => 1,
+        (SU2(0) × SU2(1)) => 1,
+        (SU2(1) × SU2(1)) => 1,
+      ]),
+    )
+  end
+
+  @testset "Fusion of NonGroupCategory products" begin
+    ı = Fib("1")
+    τ = Fib("τ")
+    s = ı × ı
+    @test space_isequal(s ⊗ s, gradedrange([s => 1]))
+
+    s = τ × τ
+    @test space_isequal(
+      s ⊗ s, gradedrange([(ı × ı) => 1, (τ × ı) => 1, (ı × τ) => 1, (τ × τ) => 1])
+    )
+
+    σ = Ising("σ")
+    ψ = Ising("ψ")
+    s = τ × σ
+    g = gradedrange([
+      (ı × Ising("1")) => 1, (τ × Ising("1")) => 1, (ı × ψ) => 1, (τ × ψ) => 1
+    ])
+    @test space_isequal(s ⊗ s, g)
+  end
+
+  @testset "Fusion of mixed Abelian and NonAbelian products" begin
+    p2h = U1(2) × SU2(1//2)
+    p1h = U1(1) × SU2(1//2)
+    @test space_isequal(
+      p2h ⊗ p1h, gradedrange([(U1(3) × SU2(0)) => 1, (U1(3) × SU2(1)) => 1])
+    )
+
+    p1h1 = U1(1) × SU2(1//2) × Z{2}(1)
+    @test space_isequal(
+      p1h1 ⊗ p1h1,
+      gradedrange([(U1(2) × SU2(0) × Z{2}(0)) => 1, (U1(2) × SU2(1) × Z{2}(0)) => 1]),
+    )
+  end
+
+  @testset "Fusion of fully mixed products" begin
+    s = U1(1) × SU2(1//2) × Ising("σ")
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        (U1(2) × SU2(0) × Ising("1")) => 1,
+        (U1(2) × SU2(1) × Ising("1")) => 1,
+        (U1(2) × SU2(0) × Ising("ψ")) => 1,
+        (U1(2) × SU2(1) × Ising("ψ")) => 1,
+      ]),
+    )
+
+    ı = Fib("1")
+    τ = Fib("τ")
+    s = SU2(1//2) × U1(1) × τ
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        (SU2(0) × U1(2) × ı) => 1,
+        (SU2(1) × U1(2) × ı) => 1,
+        (SU2(0) × U1(2) × τ) => 1,
+        (SU2(1) × U1(2) × τ) => 1,
+      ]),
+    )
+
+    s = U1(1) × ı × τ
+    @test space_isequal(s ⊗ s, gradedrange([(U1(2) × ı × ı) => 1, (U1(2) × ı × τ) => 1]))
+  end
+
+  @testset "Fusion of different length Categories" begin
+    @test SectorProduct(U1(1) × U1(0)) ⊗ SectorProduct(U1(1)) ==
+      SectorProduct(U1(2) × U1(0))
+    @test space_isequal(
+      (@inferred SectorProduct(SU2(0) × SU2(0)) ⊗ SectorProduct(SU2(1))),
+      gradedrange([SectorProduct(SU2(1) × SU2(0)) => 1]),
+    )
+
+    @test space_isequal(
+      (@inferred SectorProduct(SU2(1) × U1(1)) ⊗ SectorProduct(SU2(0))),
+      gradedrange([SectorProduct(SU2(1) × U1(1)) => 1]),
+    )
+    @test space_isequal(
+      (@inferred SectorProduct(U1(1) × SU2(1)) ⊗ SectorProduct(U1(2))),
+      gradedrange([SectorProduct(U1(3) × SU2(1)) => 1]),
+    )
+
+    # check incompatible sectors
+    p12 = Z{2}(1) × U1(2)
+    z12 = Z{2}(1) × Z{2}(1)
+    @test_throws MethodError p12 ⊗ z12
+  end
+
+  @testset "GradedUnitRange fusion rules" begin
+    s1 = U1(1) × SU2(1//2) × Ising("σ")
+    s2 = U1(0) × SU2(1//2) × Ising("1")
+    g1 = gradedrange([s1 => 2])
+    g2 = gradedrange([s2 => 1])
+    @test space_isequal(
+      fusion_product(g1, g2),
+      gradedrange([U1(1) × SU2(0) × Ising("σ") => 2, U1(1) × SU2(1) × Ising("σ") => 2]),
+    )
+  end
+end
+
+@testset "Test Named Sector Products" begin
+  @testset "Construct from × of NamedTuples" begin
+    s = (A=U1(1),) × (B=Z{2}(0),)
+    @test length(arguments(s)) == 2
+    @test arguments(s)[:A] == U1(1)
+    @test arguments(s)[:B] == Z{2}(0)
+    @test (@inferred quantum_dimension(s)) == 1
+    @test (@inferred dual(s)) == (A=U1(-1),) × (B=Z{2}(0),)
+    @test (@inferred trivial(s)) == (A=U1(0),) × (B=Z{2}(0),)
+
+    s = (A=U1(1),) × (B=SU2(2),)
+    @test length(arguments(s)) == 2
+    @test arguments(s)[:A] == U1(1)
+    @test arguments(s)[:B] == SU2(2)
+    @test (@inferred quantum_dimension(s)) == 5
+    @test (@inferred dual(s)) == (A=U1(-1),) × (B=SU2(2),)
+    @test (@inferred trivial(s)) == (A=U1(0),) × (B=SU2(0),)
+    @test s == (B=SU2(2),) × (A=U1(1),)
+
+    s = s × (C=Ising("ψ"),)
+    @test length(arguments(s)) == 3
+    @test arguments(s)[:C] == Ising("ψ")
+    @test (@inferred quantum_dimension(s)) == 5.0
+    @test (@inferred dual(s)) == (A=U1(-1),) × (B=SU2(2),) × (C=Ising("ψ"),)
+
+    s1 = (A=U1(1),) × (B=Z{2}(0),)
+    s2 = (A=U1(1),) × (C=Z{2}(0),)
+    @test_throws ArgumentError s1 × s2
+  end
+
+  @testset "Construct from Pairs" begin
+    s = SectorProduct("A" => U1(2))
+    @test length(arguments(s)) == 1
+    @test arguments(s)[:A] == U1(2)
+    @test s == SectorProduct(; A=U1(2))
+    @test (@inferred quantum_dimension(s)) == 1
+    @test (@inferred dual(s)) == SectorProduct("A" => U1(-2))
+    @test (@inferred trivial(s)) == SectorProduct(; A=U1(0))
+
+    s = SectorProduct("B" => Ising("ψ"), :C => Z{2}(1))
+    @test length(arguments(s)) == 2
+    @test arguments(s)[:B] == Ising("ψ")
+    @test arguments(s)[:C] == Z{2}(1)
+    @test (@inferred quantum_dimension(s)) == 1.0
+  end
+
+  @testset "Comparisons with unspecified labels" begin
+    # convention: arguments evaluate as equal if unmatched labels are trivial
+    # this is different from ordered tuple convention
+    q2 = SectorProduct(; N=U1(2))
+    q20 = (N=U1(2),) × (J=SU2(0),)
+    @test q20 == q2
+    @test !(q20 < q2)
+    @test !(q2 < q20)
+
+    q21 = (N=U1(2),) × (J=SU2(1),)
+    @test q21 != q2
+    @test q20 < q21
+    @test q2 < q21
+
+    a = (A=U1(0),) × (B=U1(2),)
+    b = (B=U1(2),) × (C=U1(0),)
+    @test a == b
+    c = (B=U1(2),) × (C=U1(1),)
+    @test a != c
+  end
+
+  @testset "Quantum dimension and GradedUnitRange" begin
+    g = gradedrange([
+      SectorProduct(; A=U1(0), B=Z{2}(0)) => 1, SectorProduct(; A=U1(1), B=Z{2}(0)) => 2
+    ])  # abelian
+    @test (@inferred quantum_dimension(g)) == 3
+
+    g = gradedrange([  # non-abelian
+      SectorProduct(; A=SU2(0), B=SU2(0)) => 1,
+      SectorProduct(; A=SU2(1), B=SU2(0)) => 1,
+      SectorProduct(; A=SU2(0), B=SU2(1)) => 1,
+      SectorProduct(; A=SU2(1), B=SU2(1)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 16
+
+    # mixed group
+    g = gradedrange([
+      SectorProduct(; A=U1(2), B=SU2(0), C=Z{2}(0)) => 1,
+      SectorProduct(; A=U1(2), B=SU2(1), C=Z{2}(0)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 4
+    g = gradedrange([
+      SectorProduct(; A=SU2(0), B=Z{2}(0), C=SU2(1//2)) => 1,
+      SectorProduct(; A=SU2(0), B=Z{2}(1), C=SU2(1//2)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 4
+
+    # non group sectors
+    g_fib = gradedrange([SectorProduct(; A=Fib("1"), B=Fib("1")) => 1])
+    g_ising = gradedrange([SectorProduct(; A=Ising("1"), B=Ising("1")) => 1])
+    @test (@inferred quantum_dimension(g_fib)) == 1.0
+    @test (@inferred quantum_dimension(g_ising)) == 1.0
+
+    # mixed product Abelian / NonAbelian / NonGroup
+    g = gradedrange([
+      SectorProduct(; A=U1(2), B=SU2(0), C=Ising(1)) => 1,
+      SectorProduct(; A=U1(2), B=SU2(1), C=Ising(1)) => 1,
+      SectorProduct(; A=U1(2), B=SU2(0), C=Ising("ψ")) => 1,
+      SectorProduct(; A=U1(2), B=SU2(1), C=Ising("ψ")) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 8.0
+
+    g = gradedrange([
+      SectorProduct(; A=Fib("1"), B=SU2(0), C=U1(2)) => 1,
+      SectorProduct(; A=Fib("1"), B=SU2(1), C=U1(2)) => 1,
+      SectorProduct(; A=Fib("τ"), B=SU2(0), C=U1(2)) => 1,
+      SectorProduct(; A=Fib("τ"), B=SU2(1), C=U1(2)) => 1,
+    ])
+    @test (@inferred quantum_dimension(g)) == 4.0 + 4.0quantum_dimension(Fib("τ"))
+  end
+
+  @testset "Fusion of Abelian products" begin
+    q00 = SectorProduct(;)
+    q10 = SectorProduct(; A=U1(1))
+    q01 = SectorProduct(; B=U1(1))
+    q11 = SectorProduct(; A=U1(1), B=U1(1))
+
+    @test (@inferred q10 ⊗ q10) == SectorProduct(; A=U1(2))
+    @test (@inferred q01 ⊗ q00) == q01
+    @test (@inferred q00 ⊗ q01) == q01
+    @test (@inferred q10 ⊗ q01) == q11
+    @test q11 ⊗ q11 == SectorProduct(; A=U1(2), B=U1(2))
+
+    s11 = SectorProduct(; A=U1(1), B=Z{2}(1))
+    s10 = SectorProduct(; A=U1(1))
+    s01 = SectorProduct(; B=Z{2}(1))
+    @test (@inferred s01 ⊗ q00) == s01
+    @test (@inferred q00 ⊗ s01) == s01
+    @test (@inferred s10 ⊗ s01) == s11
+    @test s11 ⊗ s11 == SectorProduct(; A=U1(2), B=Z{2}(0))
+  end
+
+  @testset "Fusion of NonAbelian products" begin
+    p0 = SectorProduct(;)
+    pha = SectorProduct(; A=SU2(1//2))
+    phb = SectorProduct(; B=SU2(1//2))
+    phab = SectorProduct(; A=SU2(1//2), B=SU2(1//2))
+
+    @test space_isequal(
+      (@inferred pha ⊗ pha),
+      gradedrange([SectorProduct(; A=SU2(0)) => 1, SectorProduct(; A=SU2(1)) => 1]),
+    )
+    @test space_isequal((@inferred pha ⊗ p0), gradedrange([pha => 1]))
+    @test space_isequal((@inferred p0 ⊗ phb), gradedrange([phb => 1]))
+    @test space_isequal((@inferred pha ⊗ phb), gradedrange([phab => 1]))
+
+    @test space_isequal(
+      phab ⊗ phab,
+      gradedrange([
+        SectorProduct(; A=SU2(0), B=SU2(0)) => 1,
+        SectorProduct(; A=SU2(1), B=SU2(0)) => 1,
+        SectorProduct(; A=SU2(0), B=SU2(1)) => 1,
+        SectorProduct(; A=SU2(1), B=SU2(1)) => 1,
+      ]),
+    )
+  end
+
+  @testset "Fusion of NonGroupCategory products" begin
+    ı = Fib("1")
+    τ = Fib("τ")
+    s = SectorProduct(; A=ı, B=ı)
+    @test space_isequal(s ⊗ s, gradedrange([s => 1]))
+
+    s = SectorProduct(; A=τ, B=τ)
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        SectorProduct(; A=ı, B=ı) => 1,
+        SectorProduct(; A=τ, B=ı) => 1,
+        SectorProduct(; A=ı, B=τ) => 1,
+        SectorProduct(; A=τ, B=τ) => 1,
+      ]),
+    )
+
+    σ = Ising("σ")
+    ψ = Ising("ψ")
+    s = SectorProduct(; A=τ, B=σ)
+    g = gradedrange([
+      SectorProduct(; A=ı, B=Ising("1")) => 1,
+      SectorProduct(; A=τ, B=Ising("1")) => 1,
+      SectorProduct(; A=ı, B=ψ) => 1,
+      SectorProduct(; A=τ, B=ψ) => 1,
+    ])
+    @test space_isequal(s ⊗ s, g)
+  end
+
+  @testset "Fusion of mixed Abelian and NonAbelian products" begin
+    q0h = SectorProduct(; J=SU2(1//2))
+    q10 = (N=U1(1),) × (J=SU2(0),)
+    # Put names in reverse order sometimes:
+    q1h = (J=SU2(1//2),) × (N=U1(1),)
+    q11 = (N=U1(1),) × (J=SU2(1),)
+    q20 = (N=U1(2),) × (J=SU2(0),)  # julia 1.6 does not accept gradedrange without J
+    q2h = (N=U1(2),) × (J=SU2(1//2),)
+    q21 = (N=U1(2),) × (J=SU2(1),)
+    q22 = (N=U1(2),) × (J=SU2(2),)
+
+    @test space_isequal(q1h ⊗ q1h, gradedrange([q20 => 1, q21 => 1]))
+    @test space_isequal(q10 ⊗ q1h, gradedrange([q2h => 1]))
+    @test space_isequal((@inferred q0h ⊗ q1h), gradedrange([q10 => 1, q11 => 1]))
+    @test space_isequal(q11 ⊗ q11, gradedrange([q20 => 1, q21 => 1, q22 => 1]))
+  end
+
+  @testset "Fusion of fully mixed products" begin
+    s = SectorProduct(; A=U1(1), B=SU2(1//2), C=Ising("σ"))
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        SectorProduct(; A=U1(2), B=SU2(0), C=Ising("1")) => 1,
+        SectorProduct(; A=U1(2), B=SU2(1), C=Ising("1")) => 1,
+        SectorProduct(; A=U1(2), B=SU2(0), C=Ising("ψ")) => 1,
+        SectorProduct(; A=U1(2), B=SU2(1), C=Ising("ψ")) => 1,
+      ]),
+    )
+
+    ı = Fib("1")
+    τ = Fib("τ")
+    s = SectorProduct(; A=SU2(1//2), B=U1(1), C=τ)
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        SectorProduct(; A=SU2(0), B=U1(2), C=ı) => 1,
+        SectorProduct(; A=SU2(1), B=U1(2), C=ı) => 1,
+        SectorProduct(; A=SU2(0), B=U1(2), C=τ) => 1,
+        SectorProduct(; A=SU2(1), B=U1(2), C=τ) => 1,
+      ]),
+    )
+
+    s = SectorProduct(; A=τ, B=U1(1), C=ı)
+    @test space_isequal(
+      s ⊗ s,
+      gradedrange([
+        SectorProduct(; B=U1(2), A=ı, C=ı) => 1, SectorProduct(; B=U1(2), A=τ, C=ı) => 1
+      ]),
+    )
+  end
+  @testset "GradedUnitRange fusion rules" begin
+    s1 = SectorProduct(; A=U1(1), B=SU2(1//2), C=Ising("σ"))
+    s2 = SectorProduct(; A=U1(0), B=SU2(1//2), C=Ising("1"))
+    g1 = gradedrange([s1 => 2])
+    g2 = gradedrange([s2 => 1])
+    s3 = SectorProduct(; A=U1(1), B=SU2(0), C=Ising("σ"))
+    s4 = SectorProduct(; A=U1(1), B=SU2(1), C=Ising("σ"))
+    @test space_isequal(fusion_product(g1, g2), gradedrange([s3 => 2, s4 => 2]))
+
+    sA = SectorProduct(; A=U1(1))
+    sB = SectorProduct(; B=SU2(1//2))
+    sAB = SectorProduct(; A=U1(1), B=SU2(1//2))
+    gA = gradedrange([sA => 2])
+    gB = gradedrange([sB => 1])
+    @test space_isequal(fusion_product(gA, gB), gradedrange([sAB => 2]))
+  end
+end
+
+@testset "Mixing implementations" begin
+  st1 = SectorProduct(U1(1))
+  sA1 = SectorProduct(; A=U1(1))
+
+  @test sA1 != st1
+  @test_throws MethodError sA1 < st1
+  @test_throws MethodError st1 < sA1
+  @test_throws MethodError st1 ⊗ sA1
+  @test_throws MethodError sA1 ⊗ st1
+  @test_throws ArgumentError st1 × sA1
+  @test_throws ArgumentError sA1 × st1
+end
+
+@testset "Empty SymmetrySector" begin
+  st1 = SectorProduct(U1(1))
+  sA1 = SectorProduct(; A=U1(1))
+
+  for s in (SectorProduct(()), SectorProduct((;)))
+    @test s == TrivialSector()
+    @test s == SectorProduct(())
+    @test s == SectorProduct((;))
+
+    @test !(s < SectorProduct())
+    @test !(s < SectorProduct(;))
+
+    @test (@inferred s × SectorProduct(())) == s
+    @test (@inferred s × SectorProduct((;))) == s
+    @test (@inferred s ⊗ SectorProduct(())) == s
+    @test (@inferred s ⊗ SectorProduct((;))) == s
+
+    @test (@inferred dual(s)) == s
+    @test (@inferred trivial(s)) == s
+    @test (@inferred quantum_dimension(s)) == 1
+
+    g0 = gradedrange([s => 2])
+    @test space_isequal((@inferred fusion_product(g0, g0)), gradedrange([s => 4]))
+
+    @test (@inferred s × U1(1)) == st1
+    @test (@inferred U1(1) × s) == st1
+    @test (@inferred s × st1) == st1
+    @test (@inferred st1 × s) == st1
+    @test (@inferred s × sA1) == sA1
+    @test (@inferred sA1 × s) == sA1
+
+    @test (@inferred U1(1) ⊗ s) == st1
+    @test (@inferred s ⊗ U1(1)) == st1
+    @test (@inferred SU2(0) ⊗ s) == gradedrange([SectorProduct(SU2(0)) => 1])
+    @test (@inferred s ⊗ SU2(0)) == gradedrange([SectorProduct(SU2(0)) => 1])
+    @test (@inferred Fib("τ") ⊗ s) == gradedrange([SectorProduct(Fib("τ")) => 1])
+    @test (@inferred s ⊗ Fib("τ")) == gradedrange([SectorProduct(Fib("τ")) => 1])
+
+    @test (@inferred st1 ⊗ s) == st1
+    @test (@inferred SectorProduct(SU2(0)) ⊗ s) == gradedrange([SectorProduct(SU2(0)) => 1])
+    @test (@inferred SectorProduct(Fib("τ"), SU2(1), U1(2)) ⊗ s) ==
+      gradedrange([SectorProduct(Fib("τ"), SU2(1), U1(2)) => 1])
+
+    @test (@inferred sA1 ⊗ s) == sA1
+    @test (@inferred SectorProduct(; A=SU2(0)) ⊗ s) ==
+      gradedrange([SectorProduct(; A=SU2(0)) => 1])
+    @test (@inferred SectorProduct(; A=Fib("τ"), B=SU2(1), C=U1(2)) ⊗ s) ==
+      gradedrange([SectorProduct(; A=Fib("τ"), B=SU2(1), C=U1(2)) => 1])
+
+    # Empty behaves as empty NamedTuple
+    @test s != U1(0)
+    @test s == SectorProduct(U1(0))
+    @test s == SectorProduct(; A=U1(0))
+    @test SectorProduct(; A=U1(0)) == s
+    @test s != sA1
+    @test s != st1
+
+    @test s < st1
+    @test SectorProduct(U1(-1)) < s
+    @test s < sA1
+    @test s > SectorProduct(; A=U1(-1))
+    @test !(s < SectorProduct(; A=U1(0)))
+    @test !(s > SectorProduct(; A=U1(0)))
+  end
+end
+end
diff --git a/NDTensors/src/lib/SymmetrySectors/test/test_simple_sectors.jl b/NDTensors/src/lib/SymmetrySectors/test/test_simple_sectors.jl
new file mode 100644
index 0000000000..a2b243c725
--- /dev/null
+++ b/NDTensors/src/lib/SymmetrySectors/test/test_simple_sectors.jl
@@ -0,0 +1,215 @@
+@eval module $(gensym())
+using NDTensors.GradedAxes: dual
+using NDTensors.SymmetrySectors:
+  Fib,
+  Ising,
+  O2,
+  SU,
+  SU2,
+  TrivialSector,
+  U1,
+  Z,
+  quantum_dimension,
+  fundamental,
+  istrivial,
+  trivial
+using Test: @inferred, @test, @testset, @test_throws
+@testset "Test SymmetrySectors Types" begin
+  @testset "TrivialSector" begin
+    q = TrivialSector()
+
+    @test (@inferred quantum_dimension(q)) == 1
+    @test q == q
+    @test trivial(q) == q
+    @test istrivial(q)
+
+    @test dual(q) == q
+    @test !isless(q, q)
+  end
+
+  @testset "U(1)" begin
+    q1 = U1(1)
+    q2 = U1(2)
+    q3 = U1(3)
+
+    @test quantum_dimension(q1) == 1
+    @test quantum_dimension(q2) == 1
+    @test (@inferred quantum_dimension(q1)) == 1
+
+    @test trivial(q1) == U1(0)
+    @test trivial(U1) == U1(0)
+    @test istrivial(U1(0))
+
+    @test dual(U1(2)) == U1(-2)
+    @test isless(U1(1), U1(2))
+    @test !isless(U1(2), U1(1))
+    @test U1(Int8(1)) == U1(1)
+    @test U1(UInt32(1)) == U1(1)
+
+    @test U1(0) == TrivialSector()
+    @test TrivialSector() == U1(0)
+    @test U1(-1) < TrivialSector()
+    @test TrivialSector() < U1(1)
+    @test U1(Int8(1)) < U1(Int32(2))
+  end
+
+  @testset "Z₂" begin
+    z0 = Z{2}(0)
+    z1 = Z{2}(1)
+
+    @test trivial(Z{2}) == Z{2}(0)
+    @test istrivial(Z{2}(0))
+
+    @test quantum_dimension(z0) == 1
+    @test quantum_dimension(z1) == 1
+    @test (@inferred quantum_dimension(z0)) == 1
+
+    @test dual(z0) == z0
+    @test dual(z1) == z1
+
+    @test dual(Z{2}(1)) == Z{2}(1)
+    @test isless(Z{2}(0), Z{2}(1))
+    @test !isless(Z{2}(1), Z{2}(0))
+    @test Z{2}(0) == z0
+    @test Z{2}(-3) == z1
+
+    @test Z{2}(0) == TrivialSector()
+    @test TrivialSector() < Z{2}(1)
+    @test_throws MethodError U1(0) < Z{2}(1)
+    @test Z{2}(0) != Z{2}(1)
+    @test Z{2}(0) != Z{3}(0)
+    @test Z{2}(0) != U1(0)
+  end
+
+  @testset "O(2)" begin
+    s0e = O2(0)
+    s0o = O2(-1)
+    s12 = O2(1//2)
+    s1 = O2(1)
+
+    @test trivial(O2) == s0e
+    @test istrivial(s0e)
+
+    @test (@inferred quantum_dimension(s0e)) == 1
+    @test (@inferred quantum_dimension(s0o)) == 1
+    @test (@inferred quantum_dimension(s12)) == 2
+    @test (@inferred quantum_dimension(s1)) == 2
+
+    @test (@inferred dual(s0e)) == s0e
+    @test (@inferred dual(s0o)) == s0o
+    @test (@inferred dual(s12)) == s12
+    @test (@inferred dual(s1)) == s1
+
+    @test s0o < s0e < s12 < s1
+    @test s0e == TrivialSector()
+    @test s0o < TrivialSector()
+    @test TrivialSector() < s12
+  end
+
+  @testset "SU(2)" begin
+    j1 = SU2(0)
+    j2 = SU2(1//2)  # Rational will be cast to HalfInteger
+    j3 = SU2(1)
+    j4 = SU2(3//2)
+
+    # alternative constructors
+    @test j2 == SU{2}((1,))  # tuple SU(N)-like constructor
+    @test j2 == SU{2,1}((1,))  # tuple constructor with explicit {N,N-1}
+    @test j2 == SU((1,))  # infer N from tuple length
+    @test j2 == SU{2}((Int8(1),))  # any Integer type accepted
+    @test j2 == SU{2}((UInt32(1),))  # any Integer type accepted
+    @test j2 == SU2(1 / 2)  # Float will be cast to HalfInteger
+    @test_throws MethodError SU2((1,))  # avoid confusion between tuple and half-integer interfaces
+    @test_throws MethodError SU{2,1}(1)  # avoid confusion
+
+    @test trivial(SU{2}) == SU2(0)
+    @test istrivial(SU2(0))
+    @test fundamental(SU{2}) == SU2(1//2)
+
+    @test quantum_dimension(j1) == 1
+    @test quantum_dimension(j2) == 2
+    @test quantum_dimension(j3) == 3
+    @test quantum_dimension(j4) == 4
+    @test (@inferred quantum_dimension(j1)) == 1
+
+    @test dual(j1) == j1
+    @test dual(j2) == j2
+    @test dual(j3) == j3
+    @test dual(j4) == j4
+
+    @test j1 < j2 < j3 < j4
+    @test SU2(0) == TrivialSector()
+    @test !(j2 < TrivialSector())
+    @test TrivialSector() < j2
+  end
+
+  @testset "SU(N)" begin
+    f3 = SU{3}((1, 0))
+    f4 = SU{4}((1, 0, 0))
+    ad3 = SU{3}((2, 1))
+    ad4 = SU{4}((2, 1, 1))
+
+    @test trivial(SU{3}) == SU{3}((0, 0))
+    @test istrivial(SU{3}((0, 0)))
+    @test trivial(SU{4}) == SU{4}((0, 0, 0))
+    @test istrivial(SU{4}((0, 0, 0)))
+    @test SU{3}((0, 0)) == TrivialSector()
+    @test SU{4}((0, 0, 0)) == TrivialSector()
+
+    @test fundamental(SU{3}) == f3
+    @test fundamental(SU{4}) == f4
+
+    @test dual(f3) == SU{3}((1, 1))
+    @test dual(f4) == SU{4}((1, 1, 1))
+    @test dual(ad3) == ad3
+    @test dual(ad4) == ad4
+
+    @test quantum_dimension(f3) == 3
+    @test quantum_dimension(f4) == 4
+    @test quantum_dimension(ad3) == 8
+    @test quantum_dimension(ad4) == 15
+    @test quantum_dimension(SU{3}((4, 2))) == 27
+    @test quantum_dimension(SU{3}((3, 3))) == 10
+    @test quantum_dimension(SU{3}((3, 0))) == 10
+    @test quantum_dimension(SU{3}((0, 0))) == 1
+    @test (@inferred quantum_dimension(f3)) == 3
+  end
+
+  @testset "Fibonacci" begin
+    ı = Fib("1")
+    τ = Fib("τ")
+
+    @test trivial(Fib) == ı
+    @test istrivial(ı)
+    @test ı == TrivialSector()
+
+    @test dual(ı) == ı
+    @test dual(τ) == τ
+
+    @test (@inferred quantum_dimension(ı)) == 1.0
+    @test (@inferred quantum_dimension(τ)) == ((1 + √5) / 2)
+
+    @test ı < τ
+  end
+
+  @testset "Ising" begin
+    ı = Ising("1")
+    σ = Ising("σ")
+    ψ = Ising("ψ")
+
+    @test trivial(Ising) == ı
+    @test istrivial(ı)
+    @test ı == TrivialSector()
+
+    @test dual(ı) == ı
+    @test dual(σ) == σ
+    @test dual(ψ) == ψ
+
+    @test (@inferred quantum_dimension(ı)) == 1.0
+    @test (@inferred quantum_dimension(σ)) == √2
+    @test (@inferred quantum_dimension(ψ)) == 1.0
+
+    @test ı < σ < ψ
+  end
+end
+end
diff --git a/NDTensors/src/lib/TensorAlgebra/ext/TensorAlgebraGradedAxesExt/test/test_contract.jl b/NDTensors/src/lib/TensorAlgebra/ext/TensorAlgebraGradedAxesExt/test/test_contract.jl
index 68e1374531..1ada7d3393 100644
--- a/NDTensors/src/lib/TensorAlgebra/ext/TensorAlgebraGradedAxesExt/test/test_contract.jl
+++ b/NDTensors/src/lib/TensorAlgebra/ext/TensorAlgebraGradedAxesExt/test/test_contract.jl
@@ -3,8 +3,8 @@ using BlockArrays: Block, blocksize
 using Compat: Returns
 using NDTensors.BlockSparseArrays: BlockSparseArray
 using NDTensors.GradedAxes: gradedrange
-using NDTensors.Sectors: U1
 using NDTensors.SparseArrayInterface: densearray
+using NDTensors.SymmetrySectors: U1
 using NDTensors.TensorAlgebra: contract
 using Random: randn!
 using Test: @test, @testset
diff --git a/NDTensors/test/lib/runtests.jl b/NDTensors/test/lib/runtests.jl
index b1732e2620..b1b47cf866 100644
--- a/NDTensors/test/lib/runtests.jl
+++ b/NDTensors/test/lib/runtests.jl
@@ -15,10 +15,10 @@ using Test: @testset
   "LabelledNumbers",
   "MetalExtensions",
   "NamedDimsArrays",
-  "Sectors",
   "SmallVectors",
   "SortedSets",
   "SparseArrayDOKs",
+  "SymmetrySectors",
   "TagSets",
   "TensorAlgebra",
   "TypeParameterAccessors",