diff --git a/Project.toml b/Project.toml
index cd10327..ce7c92e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FusionTensors"
 uuid = "e16ca583-1f51-4df0-8e12-57d32947d33e"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.5.13"
+version = "0.5.14"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
 TensorProducts = "decf83d6-1968-43f4-96dc-fdb3fe15fc6d"
 TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
@@ -22,13 +23,14 @@ WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 Accessors = "0.1.42"
 BlockArrays = "1.7"
 BlockSparseArrays = "0.10"
-GradedArrays = "0.4"
+GradedArrays = "0.5.2"
 HalfIntegers = "1.6"
 LRUCache = "1.6"
 LinearAlgebra = "1.10"
 Random = "1.10"
 Strided = "2.3"
 TensorAlgebra = "0.4"
+TensorKitSectors = "0.1, 0.2"
 TensorProducts = "0.1.7"
 TypeParameterAccessors = "0.4"
 WignerSymbols = "2.0.0"
diff --git a/src/fusion_trees/clebsch_gordan_tensors.jl b/src/fusion_trees/clebsch_gordan_tensors.jl
index 3815564..09245a0 100644
--- a/src/fusion_trees/clebsch_gordan_tensors.jl
+++ b/src/fusion_trees/clebsch_gordan_tensors.jl
@@ -7,11 +7,12 @@ using WignerSymbols: clebschgordan
 
 using GradedArrays:
     AbelianStyle,
-    AbstractSector,
     NotAbelianStyle,
     O2,
-    SU,
+    SU2,
+    SectorRange,
     SymmetryStyle,
+    TrivialSector,
     dual,
     istrivial,
     quantum_dimension,
@@ -21,16 +22,17 @@ using GradedArrays:
     zero_odd
 using TensorAlgebra: contract
 using TensorProducts: ⊗
+import TensorKitSectors as TKS
 
-function symbol_1j(s::AbstractSector)
+function symbol_1j(s::SectorRange)
     cgt = clebsch_gordan_tensor(s, dual(s), trivial(s), 1)
     return sqrt(quantum_dimension(s)) * cgt[:, :, 1]
 end
 
 function clebsch_gordan_tensor(
-        s1::AbstractSector,
-        s2::AbstractSector,
-        s3::AbstractSector,
+        s1::SectorRange,
+        s2::SectorRange,
+        s3::SectorRange,
         arrow1::Bool,
         arrow2::Bool,
         inner_mult_index::Int,
@@ -47,91 +49,19 @@ function clebsch_gordan_tensor(
     return cgt
 end
 
-function clebsch_gordan_tensor(s1::S, s2::S, s3::S, outer_mult_index::Int) where {S}
-    return clebsch_gordan_tensor(SymmetryStyle(S), s1, s2, s3, outer_mult_index)
-end
-
-function clebsch_gordan_tensor(
-        ::AbelianStyle, s1::S, s2::S, s3::S, outer_mult_index::Int
-    ) where {S}
-    @assert outer_mult_index == 1
-    return s1 ⊗ s2 == s3 ? ones((1, 1, 1)) : zeros((1, 1, 1))
-end
-
-function clebsch_gordan_tensor(::NotAbelianStyle, s1::O2, s2::O2, s3::O2, ::Int)
-    return clebsch_gordan_tensor(s1, s2, s3)  # no outer multiplicity
-end
-
-function clebsch_gordan_tensor(s1::O2, s2::O2, s3::O2)
-    d1 = quantum_dimension(s1)
-    d2 = quantum_dimension(s2)
-    d3 = quantum_dimension(s3)
-    cgt = zeros((d1, d2, d3))
-    s3 ∉ sectors(s1 ⊗ s2) && return cgt
-
-    # adapted from TensorKit
-    l1 = sector_label(s1)
-    l2 = sector_label(s2)
-    l3 = sector_label(s3)
-    if l3 <= 0  # 0even or 0odd
-        if l1 <= 0 && l2 <= 0
-            cgt[1, 1, 1, 1] = 1.0
-        else
-            if istrivial(s3)
-                cgt[1, 2, 1, 1] = 1.0 / sqrt(2)
-                cgt[2, 1, 1, 1] = 1.0 / sqrt(2)
-            else
-                cgt[1, 2, 1, 1] = 1.0 / sqrt(2)
-                cgt[2, 1, 1, 1] = -1.0 / sqrt(2)
-            end
-        end
-    elseif l1 <= 0  # 0even or 0odd
-        cgt[1, 1, 1, 1] = 1.0
-        cgt[1, 2, 2, 1] = s1 == zero_odd(O2) ? -1.0 : 1.0
-    elseif l2 == 0
-        cgt[1, 1, 1, 1] = 1.0
-        cgt[2, 1, 2, 1] = s2 == zero_odd(O2) ? -1.0 : 1.0
-    elseif l3 == l1 + l2
-        cgt[1, 1, 1, 1] = 1.0
-        cgt[2, 2, 2, 1] = 1.0
-    elseif l3 == l1 - l2
-        cgt[1, 2, 1, 1] = 1.0
-        cgt[2, 1, 2, 1] = 1.0
-    elseif l3 == l2 - l1
-        cgt[2, 1, 1, 1] = 1.0
-        cgt[1, 2, 2, 1] = 1.0
+function clebsch_gordan_tensor(s1::S, s2::S, s3::S, outer_mult_index::Int = 1) where {S}
+    CGC = TKS.fusiontensor(GradedArrays.label.((s1, s2, s3))...)
+    outer_mult_index ∈ axes(CGC, 4) || throw(ArgumentError("invalid outer multiplicity index"))
+    if TKS.FusionStyle(S) === TKS.GenericFusion()
+        # TODO: do we want a view here?
+        return CGC[:, :, :, outer_mult_index]
+    else
+        return dropdims(CGC; dims = 4)
     end
-    return cgt
-end
-
-function clebsch_gordan_tensor(::NotAbelianStyle, s1::SU{2}, s2::SU{2}, s3::SU{2}, ::Int)
-    return clebsch_gordan_tensor(s1, s2, s3)  # no outer multiplicity
-end
-
-function clebsch_gordan_tensor(s1::SU{2}, s2::SU{2}, s3::SU{2})
-    d1 = quantum_dimension(s1)
-    d2 = quantum_dimension(s2)
-    d3 = quantum_dimension(s3)
-    j1 = half(d1 - 1)
-    j2 = half(d2 - 1)
-    j3 = half(d3 - 1)
-    cgtensor = Array{Float64, 3}(undef, (d1, d2, d3))
-    for (i, j, k) in Iterators.product(1:d1, 1:d2, 1:d3)
-        m1 = j1 - i + 1
-        m2 = j2 - j + 1
-        m3 = j3 - k + 1
-        cgtensor[i, j, k] = clebschgordan(j1, m1, j2, m2, j3, m3)
-    end
-    return cgtensor
 end
 
-function clebsch_gordan_tensor(
-        ::NotAbelianStyle, s1::SU{3}, s2::SU{3}, s3::SU{3}, outer_mult_index::Int
-    )
-    d1 = quantum_dimension(s1)
-    d2 = quantum_dimension(s2)
-    d3 = quantum_dimension(s3)
-    cgtensor = zeros(d1, d2, d3)
-    # dummy
-    return cgtensor
+# TODO: remove once TensorKitSectors fixes this
+function clebsch_gordan_tensor(s1::TrivialSector, s2::TrivialSector, s3::TrivialSector, outer_mult_index::Int = 1)
+    outer_mult_index == 1 || throw(ArgumentError("invalid outer multiplicity index"))
+    return fill(1, (1, 1, 1))
 end
diff --git a/src/fusion_trees/fusiontree.jl b/src/fusion_trees/fusiontree.jl
index b448d6f..c8b8830 100644
--- a/src/fusion_trees/fusiontree.jl
+++ b/src/fusion_trees/fusiontree.jl
@@ -5,10 +5,10 @@
 
 using GradedArrays:
     GradedArrays,
-    ×,
     AbstractGradedUnitRange,
-    AbstractSector,
-    SectorProduct,
+    SectorProductRange,
+    SectorRange,
+    ×,
     arguments,
     flip,
     flip_dual,
@@ -16,6 +16,7 @@ using GradedArrays:
     nsymbol,
     sector_multiplicities,
     sector_type,
+    sectorproduct,
     sectors,
     to_gradedrange,
     trivial
@@ -125,7 +126,7 @@ function GradedArrays.:×(f1::SectorFusionTree, f2::SectorFusionTree)
     )
 end
 
-function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProductRange})
     transposed_indices = outer_multiplicity_split.(
         Base.tail(leaves(f)),
         branch_sectors(f),
@@ -136,7 +137,7 @@ function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct})
     arguments_leaves = arguments.(leaves(f))
     arguments_branch_sectors = arguments.(branch_sectors(f))
     # TODO way to avoid explicit ntuple?
-    # works fine for Tuple and NamedTuple SectorProduct
+    # works fine for Tuple and NamedTuple SectorProductRange
     return ntuple(
         i -> SectorFusionTree(
             getindex.(arguments_leaves, i),
@@ -149,21 +150,21 @@ function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct})
     )
 end
 
-function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct, 0})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProductRange, 0})
     return map(arg -> SectorFusionTree((), (), arg, (), ()), arguments(root_sector(f)))
 end
 
-function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct, 1})
+function GradedArrays.arguments(f::SectorFusionTree{<:SectorProductRange, 1})
     arguments_root = arguments(root_sector(f))
     arguments_leave = arguments(only(leaves(f)))
-    # use map(keys) to stay agnostic with respect to SectorProduct implementation
+    # use map(keys) to stay agnostic with respect to SectorProductRange implementation
     return map(keys(arguments_root)) do k
         return SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ())
     end
 end
 
 # TBD change type depending on AbelianStyle?
-fusiontree_eltype(::Type{<:AbstractSector}) = Float64
+fusiontree_eltype(::Type{<:SectorRange}) = Float64
 
 # constructors
 function build_trees(legs::Vararg{AbstractGradedUnitRange})
@@ -175,7 +176,7 @@ function build_trees(legs::Vararg{AbstractGradedUnitRange})
 end
 
 function build_trees(
-        sectors_to_fuse::NTuple{N, <:AbstractSector}, arrows_to_fuse::NTuple{N, Bool}
+        sectors_to_fuse::NTuple{N, <:SectorRange}, arrows_to_fuse::NTuple{N, Bool}
     ) where {N}
     # construct all authorized trees with fixed outer sectors
     trees = [SectorFusionTree(first(sectors_to_fuse), first(arrows_to_fuse))]
@@ -186,18 +187,18 @@ end
 # =====================================  Internals  ========================================
 #
 
-# --------------- SectorProduct helper functions  ---------------
+# --------------- SectorProductRange helper functions  ---------------
 function outer_multiplicity_kron(
         sec1, sec2, fused, outer_multiplicity1, outer_multiplicity2
     )
-    n = nsymbol(sec1, sec2, fused)
+    n = Int(nsymbol(sec1, sec2, fused))
     linear_inds = LinearIndices((n, outer_multiplicity2))
     return linear_inds[outer_multiplicity1, outer_multiplicity2]
 end
 
 function outer_multiplicity_split(
         sec1::S, sec2::S, fused::S, outer_mult_index::Integer
-    ) where {S <: SectorProduct}
+    ) where {S <: SectorProductRange}
     args1 = arguments(sec1)
     args2 = arguments(sec2)
     args12 = arguments(fused)
@@ -207,15 +208,15 @@ end
 
 # --------------- Build trees  ---------------
 # zero leg: need S to get sector type information
-function SectorFusionTree{S}() where {S <: AbstractSector}
+function SectorFusionTree{S}() where {S <: SectorRange}
     return SectorFusionTree((), (), trivial(S), (), ())
 end
-function SectorFusionTree{S}(::Tuple{}, ::Tuple{}) where {S <: AbstractSector}
+function SectorFusionTree{S}(::Tuple{}, ::Tuple{}) where {S <: SectorRange}
     return SectorFusionTree((), (), trivial(S), (), ())
 end
 
 # one leg
-function SectorFusionTree(sect::AbstractSector, arrow::Bool)
+function SectorFusionTree(sect::SectorRange, arrow::Bool)
     return SectorFusionTree((sect,), (arrow,), sect, (), ())
 end
 
@@ -227,7 +228,7 @@ end
 
 function append_tree_leave(
         parent_tree::SectorFusionTree,
-        branch_sector::AbstractSector,
+        branch_sector::SectorRange,
         level_arrow::Bool,
         child_root_sector,
         outer_mult,
@@ -242,7 +243,7 @@ function append_tree_leave(
 end
 
 function fuse_next_sector(
-        parent_tree::SectorFusionTree, branch_sector::AbstractSector, level_arrow::Bool
+        parent_tree::SectorFusionTree, branch_sector::SectorRange, level_arrow::Bool
     )
     new_space = to_gradedrange(root_sector(parent_tree) ⊗ branch_sector)
     return mapreduce(
@@ -282,7 +283,7 @@ function to_array(f::SectorFusionTree)
     return grow_tensor_tree(tree_tensor, f)
 end
 
-function to_array(f::SectorFusionTree{<:SectorProduct})
+function to_array(f::SectorFusionTree{<:SectorProductRange})
     args = convert.(Array, arguments(f))
     return reduce(_tensor_kron, args)
 end
diff --git a/src/fusiontensor/fusiontensor.jl b/src/fusiontensor/fusiontensor.jl
index fcc232e..51c8c9d 100644
--- a/src/fusiontensor/fusiontensor.jl
+++ b/src/fusiontensor/fusiontensor.jl
@@ -6,7 +6,7 @@ using BlockSparseArrays:
     AbstractBlockSparseMatrix, BlockSparseArray, eachblockstoredindex, to_block_indices
 using GradedArrays:
     AbstractGradedUnitRange,
-    SectorProduct,
+    SymmetryStyle,
     TrivialSector,
     dual,
     findfirstblock,
@@ -35,7 +35,7 @@ function flip_domain(nonflipped_col_axis, nonflipped_trees_to_ranges)
     return col_axis, domain_trees_to_ranges_mapping
 end
 
-function fuse_axes(::Type{S}, ::Tuple{}) where {S <: AbstractSector}
+function fuse_axes(::Type{S}, ::Tuple{}) where {S <: SectorRange}
     fused_axis = trivial_axis(S)
     trees_to_ranges_mapping = Dict([SectorFusionTree{S}() => Block(1)[1:1]])
     return fused_axis, trees_to_ranges_mapping
@@ -253,7 +253,9 @@ function GradedArrays.sector_type(::Type{FT}) where {FT <: FusionTensor}
     return sector_type(type_parameters(FT, 3))
 end
 
-SymmetryStyle(::Type{FT}) where {FT <: FusionTensor} = SymmetryStyle(sector_type(FT))
+function GradedArrays.SymmetryStyle(::Type{FT}) where {FT <: FusionTensor}
+    return SymmetryStyle(sector_type(FT))
+end
 
 # ==============================  FusionTensor interface  ==================================
 
diff --git a/src/fusiontensor/fusiontensoraxes.jl b/src/fusiontensor/fusiontensoraxes.jl
index 0e999f5..efc9799 100644
--- a/src/fusiontensor/fusiontensoraxes.jl
+++ b/src/fusiontensor/fusiontensoraxes.jl
@@ -2,7 +2,6 @@ using BlockArrays: BlockArrays
 using GradedArrays:
     GradedArrays,
     AbstractGradedUnitRange,
-    AbstractSector,
     SymmetryStyle,
     TrivialSector,
     dual,
diff --git a/test/Project.toml b/test/Project.toml
index 268e23a..3fb032e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ FusionTensors = "e16ca583-1f51-4df0-8e12-57d32947d33e"
 GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SUNRepresentations = "1a50b95c-7aac-476d-a9ce-2bfc675fc617"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
@@ -17,9 +18,10 @@ Aqua = "0.8.11"
 BlockArrays = "1.6"
 BlockSparseArrays = "0.10"
 FusionTensors = "0.5"
-GradedArrays = "0.4"
+GradedArrays = "0.5"
 LinearAlgebra = "1.10.0"
 Random = "1.10"
+SUNRepresentations = "0.3.1"
 SafeTestsets = "0.1.0"
 Suppressor = "0.2.8"
 TensorAlgebra = "0.4"
diff --git a/test/test_array_cast.jl b/test/test_array_cast.jl
index b3fddbc..ed3f633 100644
--- a/test/test_array_cast.jl
+++ b/test/test_array_cast.jl
@@ -4,7 +4,7 @@ using Test: @test, @test_throws, @testset
 using BlockArrays: Block, BlockedArray, blocksize
 
 using FusionTensors: FusionTensor, data_matrix, to_fusiontensor
-using GradedArrays: O2, SectorProduct, SU2, TrivialSector, U1, dual, gradedrange
+using GradedArrays: O2, SU2, TrivialSector, U1, dual, gradedrange, sectorproduct
 using TensorProducts: tensor_product
 
 include("setup.jl")
@@ -389,20 +389,21 @@ end
             tRVB[i, 1, 1, 1, i + 1] = 1.0
         end
 
-        gd = gradedrange([SectorProduct(s, U1(3)) => 1])
+        gd = gradedrange([sectorproduct(s, U1(3)) => 1])
         codomain_legs = (dual(gd),)
-        gD = gradedrange([SectorProduct(SU2(0), U1(1)) => 1, SectorProduct(s, U1(0)) => 1])
+        gD = gradedrange([sectorproduct(SU2(0), U1(1)) => 1, sectorproduct(s, U1(0)) => 1])
         domain_legs = (gD, gD, gD, gD)
         ft = to_fusiontensor(tRVB, codomain_legs, domain_legs)
         @test isnothing(check_sanity(ft))
         @test Array(ft) ≈ tRVB
 
         # same with NamedTuples
-        gd_nt = gradedrange([SectorProduct(; S = s, N = U1(3)) => 1])
+        gd_nt = gradedrange([sectorproduct((; S = s, N = U1(3))) => 1])
         codomain_legs_nt = (dual(gd_nt),)
         gD_nt = gradedrange(
             [
-                SectorProduct(; S = SU2(0), N = U1(1)) => 1, SectorProduct(; S = s, N = U1(0)) => 1,
+                sectorproduct((; S = SU2(0), N = U1(1))) => 1,
+                sectorproduct((; S = s, N = U1(0))) => 1,
             ]
         )
         domain_legs_nt = (gD_nt, gD_nt, gD_nt, gD_nt)
diff --git a/test/test_basics.jl b/test/test_basics.jl
index be017ce..e9327b1 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -17,7 +17,6 @@ using FusionTensors:
 using GradedArrays:
     U1,
     SU2,
-    SectorProduct,
     TrivialSector,
     Z,
     checkspaces,
@@ -26,6 +25,7 @@ using GradedArrays:
     flip,
     gradedrange,
     sector_type,
+    sectorproduct,
     space_isequal
 using TensorAlgebra: tuplemortar
 using TensorProducts: tensor_product
@@ -68,15 +68,15 @@ include("setup.jl")
     @test space_isequal(domain_axis(ft1), g2)
     @test isnothing(check_sanity(ft0))
     @test isnothing(check_sanity(ft1))
-    @test sector_type(ft1) === U1{Int}
-    @test sector_type(typeof(ft1)) === U1{Int}
+    @test sector_type(ft1) ≡ U1
+    @test sector_type(typeof(ft1)) ≡ U1
 
     m1 = BlockSparseArray{Float64}(undef, g1, g2)
     m1[Block(2, 1)] = ones(2, 2)  # forbidden
     @test_throws ArgumentError FusionTensor(m1, (g1,), (g2,))
 
     # Base methods
-    @test eltype(ft1) === Float64
+    @test eltype(ft1) ≡ Float64
     @test length(ft1) == 30
     @test ndims(ft1) == 2
     @test size(ft1) == tuplemortar(((6,), (5,)))
@@ -167,7 +167,7 @@ end
     @test size(ft1) == tuplemortar(((6,), ()))
     @test size(data_matrix(ft1)) == (6, 1)
     @test isnothing(check_sanity(ft1))
-    @test sector_type(ft1) === sector_type(g1)
+    @test sector_type(ft1) ≡ sector_type(g1)
 
     # one column axis
     ft2 = FusionTensor{Float64}(undef, (), (g1,))
@@ -179,7 +179,7 @@ end
     @test size(ft2) == tuplemortar(((), (6,)))
     @test size(data_matrix(ft2)) == (1, 6)
     @test isnothing(check_sanity(ft2))
-    @test sector_type(ft2) === sector_type(g1)
+    @test sector_type(ft2) ≡ sector_type(g1)
 
     # zero axis
     ft3 = FusionTensor{Float64}(undef, (), ())
@@ -191,7 +191,7 @@ end
     @test size(ft3) == tuplemortar(((), ()))
     @test size(data_matrix(ft3)) == (1, 1)
     @test isnothing(check_sanity(ft3))
-    @test sector_type(ft3) === TrivialSector
+    @test sector_type(ft3) ≡ TrivialSector
 
     ft4 = FusionTensor{Float64}(undef, (g1, g1), ())
     @test ft4 isa FusionTensor
@@ -251,34 +251,34 @@ end
     @test isnothing(check_sanity(ft3))
 
     ft4 = +ft3
-    @test ft4 === ft3  # same object
+    @test ft4 ≡ ft3  # same object
 
     ft4 = -ft3
     @test isnothing(check_sanity(ft4))
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test norm(ft4) ≈ norm(ft3)
     @test norm(ft4 + ft3) ≈ 0.0
 
     ft4 = ft3 + ft3
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft4), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
     @test norm(ft4) ≈ 2norm(ft3)
 
     ft4 = ft3 - ft3
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft4), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
     @test norm(ft4) ≈ 0.0
 
     ft4 = 2 * ft3
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft4), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
@@ -286,16 +286,16 @@ end
     @test norm(ft4) ≈ 2norm(ft3)
 
     ft4 = 2.0 * ft3
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft4), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
     @test eltype(ft4) == Float64
 
     ft4 = ft3 / 2.0
-    @test codomain_axes(ft4) === codomain_axes(ft3)
-    @test domain_axes(ft4) === domain_axes(ft3)
+    @test codomain_axes(ft4) ≡ codomain_axes(ft3)
+    @test domain_axes(ft4) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft4), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft4), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
@@ -303,16 +303,16 @@ end
     @test norm(ft4) ≈ norm(ft3) / 2.0
 
     ft5 = (1.0 + 2.0im) * ft3
-    @test codomain_axes(ft5) === codomain_axes(ft3)
-    @test domain_axes(ft5) === domain_axes(ft3)
+    @test codomain_axes(ft5) ≡ codomain_axes(ft3)
+    @test domain_axes(ft5) ≡ domain_axes(ft3)
     @test space_isequal(codomain_axis(ft5), codomain_axis(ft3))
     @test space_isequal(domain_axis(ft5), domain_axis(ft3))
     @test isnothing(check_sanity(ft4))
     @test eltype(ft5) == ComplexF64
     @test norm(ft5) ≈ √5 * norm(ft3)
 
-    @test conj(ft3) === ft3  # same object
-    @test real(ft3) === ft3
+    @test conj(ft3) ≡ ft3  # same object
+    @test real(ft3) ≡ ft3
     @test norm(imag(ft3)) == 0
 
     @test conj(ft5) isa FusionTensor{ComplexF64, 4}
@@ -325,8 +325,8 @@ end
     ft6 = conj(ft5)
     @test ft6 !== ft5  # different object
     @test isnothing(check_sanity(ft6))
-    @test codomain_axes(ft6) === codomain_axes(ft5)
-    @test domain_axes(ft6) === domain_axes(ft5)
+    @test codomain_axes(ft6) ≡ codomain_axes(ft5)
+    @test domain_axes(ft6) ≡ domain_axes(ft5)
     @test space_isequal(codomain_axis(ft6), codomain_axis(ft5))
     @test space_isequal(domain_axis(ft6), domain_axis(ft5))
     @test eltype(ft6) == ComplexF64
@@ -348,31 +348,31 @@ end
     @test_throws ArgumentError ft7 * ft3
 end
 
-@testset "missing SectorProduct" begin
-    g1 = gradedrange([SectorProduct(U1(1)) => 1])
-    g2 = gradedrange([SectorProduct(U1(1), SU2(1 // 2)) => 1])
-    g3 = gradedrange([SectorProduct(U1(1), SU2(1 // 2), Z{2}(1)) => 1])
+@testset "missing sectorproduct" begin
+    g1 = gradedrange([sectorproduct(U1(1)) => 1])
+    g2 = gradedrange([sectorproduct(U1(1), SU2(1 // 2)) => 1])
+    g3 = gradedrange([sectorproduct(U1(1), SU2(1 // 2), Z{2}(1)) => 1])
     S = sector_type(g3)
 
     ft = FusionTensor{Float64}(undef, (g1,), (dual(g2), dual(g3)))
-    @test sector_type(ft) === S
-    gr = gradedrange([SectorProduct(U1(1), SU2(0), Z{2}(0)) => 1])
+    @test sector_type(ft) ≡ S
+    gr = gradedrange([sectorproduct(U1(1), SU2(0), Z{2}(0)) => 1])
     @test space_isequal(codomain_axis(ft), gr)
     gc = gradedrange(
         [
-            SectorProduct(U1(2), SU2(0), Z{2}(1)) => 1, SectorProduct(U1(2), SU2(1), Z{2}(1)) => 1,
+            sectorproduct(U1(2), SU2(0), Z{2}(1)) => 1, sectorproduct(U1(2), SU2(1), Z{2}(1)) => 1,
         ]
     )
     @test space_isequal(domain_axis(ft), dual(gc))
 
-    gA = gradedrange([SectorProduct(; A = U1(1)) => 1])
-    gB = gradedrange([SectorProduct(; B = SU2(1 // 2)) => 1])
-    gC = gradedrange([SectorProduct(; C = Z{2}(0)) => 1])
+    gA = gradedrange([sectorproduct((; A = U1(1))) => 1])
+    gB = gradedrange([sectorproduct((; B = SU2(1 // 2))) => 1])
+    gC = gradedrange([sectorproduct((; C = Z{2}(0))) => 1])
     gABC = tensor_product(gA, gB, gC)
     S = sector_type(gABC)
 
     ft = FusionTensor{Float64}(undef, (gA, gB), (dual(gA), dual(gB), gC))
-    @test sector_type(ft) === S
+    @test sector_type(ft) ≡ S
     @test space_isequal(codomain_axis(ft), gABC)
     @test space_isequal(domain_axis(ft), dual(gABC))
 end
diff --git a/test/test_fusion_trees.jl b/test/test_fusion_trees.jl
index 7e2aa90..d63d0f7 100644
--- a/test/test_fusion_trees.jl
+++ b/test/test_fusion_trees.jl
@@ -3,17 +3,13 @@ using LinearAlgebra: I
 
 using BlockArrays: BlockArrays
 
-using FusionTensors:
-    SectorFusionTree,
-    arrows,
-    branch_sectors,
-    build_trees,
-    fusiontree_eltype,
-    leaves,
-    outer_multiplicity_indices,
-    root_sector
-using GradedArrays:
-    ×, SectorProduct, SU, SU2, TrivialSector, arguments, dual, flip, gradedrange, sector_type
+using FusionTensors: SectorFusionTree, arrows, branch_sectors, build_trees,
+    fusiontree_eltype, leaves, outer_multiplicity_indices, root_sector
+using GradedArrays: ×, SectorProduct, SectorRange, SU2, TrivialSector, arguments, dual,
+    flip, gradedrange, sector_type
+using SUNRepresentations: SUNRepresentations
+
+const SU{N} = SectorRange{SUNRepresentations.SUNIrrep{N}}
 
 @testset "Trivial fusion trees" begin
     q = TrivialSector()
@@ -116,7 +112,7 @@ end
     @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))]
 
     trees = build_trees((c3, c3), (false, false))
-    @test root_sector.(trees) == [SU{3}((1, 0)), SU{3}((2, 2))]
+    @test root_sector.(trees) == [SU{3}((2, 2)), SU{3}((1, 0))]
 
     # test GradedUnitRange interface
     g = gradedrange([SU{3}((1, 0)) => 1])
@@ -137,8 +133,8 @@ end
     @test root_sector(f) == SU{3}((0, 0))
     @test sector_type(f) == typeof(a8)
 
-    f8a = trees[2]
-    f8b = trees[3]
+    f8a = trees[4]
+    f8b = trees[5]
     @test root_sector(f8a) == a8
     @test root_sector(f8b) == a8
     @test branch_sectors(f8a) == (a8,)
diff --git a/test/test_fusiontensoraxes.jl b/test/test_fusiontensoraxes.jl
index d0c66b7..eb302a9 100644
--- a/test/test_fusiontensoraxes.jl
+++ b/test/test_fusiontensoraxes.jl
@@ -17,12 +17,12 @@ using FusionTensors:
 using GradedArrays:
     ×,
     U1,
-    SectorProduct,
     TrivialSector,
     SU2,
     checkspaces,
     dual,
     gradedrange,
+    sectorproduct,
     sector_type,
     space_isequal
 
@@ -33,7 +33,7 @@ using GradedArrays:
     @test promote_sector_type(g1, g1) === typeof(U1(1))
     @test promote_sector_type((g1, g1)) === typeof(U1(1))
 
-    sNS = SectorProduct(; N = U1(1), S = SU2(1 / 2))
+    sNS = sectorproduct((; N = U1(1), S = SU2(1 / 2)))
     gN = gradedrange([(; N = U1(1)) => 1])
     gS = gradedrange([(; S = SU2(1 / 2)) => 1])
     @test promote_sector_type(gN, gS) == typeof(sNS)