diff --git a/Project.toml b/Project.toml
index 85d6394..5f28d05 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,20 +1,31 @@
 name = "QuantumOperatorDefinitions"
 uuid = "826dd319-6fd5-459a-a990-3a4f214664bf"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [weakdeps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
+LabelledNumbers = "f856a3a6-4152-4ec4-b2a7-02c1a55d7993"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
+SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e"
 
 [extensions]
-QuantumOperatorDefinitionsITensorBaseExt = "ITensorBase"
+QuantumOperatorDefinitionsITensorBaseExt = ["ITensorBase", "NamedDimsArrays"]
+QuantumOperatorDefinitionsSymmetrySectorsExt = ["BlockArrays", "GradedUnitRanges", "LabelledNumbers", "SymmetrySectors"]
 
 [compat]
+BlockArrays = "1.3.0"
+GradedUnitRanges = "0.1.2"
 ITensorBase = "0.1.10"
+LabelledNumbers = "0.1.0"
 LinearAlgebra = "1.10"
+NamedDimsArrays = "0.4.0"
 Random = "1.10"
+SymmetrySectors = "0.1.3"
 julia = "1.10"
diff --git a/README.md b/README.md
index 978d575..480173d 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,6 @@ julia> Pkg.add("QuantumOperatorDefinitions")
 
 ````julia
 using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state
-using LinearAlgebra: Diagonal
 using SparseArrays: SparseMatrixCSC, SparseVector
 using Test: @test
 
@@ -64,8 +63,6 @@ using Test: @test
 @test op("Y") == [0 -im; im 0]
 @test op("Z") == [1 0; 0 -1]
 
-@test op("Z") isa Diagonal
-
 @test op(Float32, "X") == [0 1; 1 0]
 @test eltype(op(Float32, "X")) === Float32
 @test op(SparseMatrixCSC, "X") == [0 1; 1 0]
diff --git a/examples/README.jl b/examples/README.jl
index 5b8b5b1..f4f246d 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -38,7 +38,6 @@ julia> Pkg.add("QuantumOperatorDefinitions")
 # ## Examples
 
 using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state
-using LinearAlgebra: Diagonal
 using SparseArrays: SparseMatrixCSC, SparseVector
 using Test: @test
 
@@ -69,8 +68,6 @@ using Test: @test
 @test op("Y") == [0 -im; im 0]
 @test op("Z") == [1 0; 0 -1]
 
-@test op("Z") isa Diagonal
-
 @test op(Float32, "X") == [0 1; 1 0]
 @test eltype(op(Float32, "X")) === Float32
 @test op(SparseMatrixCSC, "X") == [0 1; 1 0]
diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index 280b531..84b8b4b 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -1,24 +1,58 @@
 module QuantumOperatorDefinitionsITensorBaseExt
 
-using ITensorBase: ITensor, Index, dag, gettag, prime
+using ITensorBase: ITensorBase, ITensor, Index, dag, gettag, prime, settag
+using NamedDimsArrays: dename
 using QuantumOperatorDefinitions:
-  QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
+  QuantumOperatorDefinitions,
+  @OpName_str,
+  OpName,
+  SiteType,
+  StateName,
+  has_fermion_string,
+  name
 
 function QuantumOperatorDefinitions.SiteType(r::Index)
-  return SiteType(gettag(r, "sitetype", "Qudit"); dim=Int(length(r)))
+  # We pass the axis of the (unnamed) Index because
+  # the Index may have originated from a slice, in which
+  # case the start may not be 1 (for NonContiguousIndex,
+  # which we need to add support for, it may not even
+  # be a unit range).
+  return SiteType(
+    gettag(r, "sitetype", "Qudit"); dim=Int.(length(r)), range=only(axes(dename(r)))
+  )
 end
 
+function (rangetype::Type{<:Index})(t::SiteType)
+  return settag(rangetype(AbstractUnitRange(t)), "sitetype", String(name(t)))
+end
+
+# TODO: Define in terms of `OpName` directly, and define a generic
+# forwarding method `has_fermion_string(n::String, t) = has_fermion_string(OpName(n), t)`.
 function QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index)
   return has_fermion_string(OpName(n), SiteType(r))
 end
 
-function Base.AbstractArray(n::OpName, r::Index)
-  # TODO: Define this with mapped dimnames.
-  return ITensor(AbstractArray(n, SiteType(r)), (prime(r), dag(r)))
+function Base.axes(::OpName, domain::Tuple{Vararg{Index}})
+  return (prime.(domain)..., dag.(domain)...)
 end
+## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{Index}})
+##   return (prime.(reverse(domain))..., dag.(domain)...)
+## end
 
-function Base.AbstractArray(n::StateName, r::Index)
-  return ITensor(AbstractArray(n, SiteType(r)), (r,))
+# Fix ambiguity error with generic `AbstractArray` version.
+function ITensorBase.ITensor(n::Union{OpName,StateName}, domain::Index...)
+  return ITensor(n, domain)
+end
+# Fix ambiguity error with generic `AbstractArray` version.
+function ITensorBase.ITensor(n::Union{OpName,StateName}, domain::Tuple{Vararg{Index}})
+  return ITensor(AbstractArray(n, domain), axes(n, domain))
+end
+function (arrtype::Type{<:AbstractArray})(
+  n::Union{OpName,StateName}, domain::Tuple{Vararg{Index}}
+)
+  # Convert to `SiteType` in case the Index specifies a `"sitetype"` tag.
+  # TODO: Try to build this into the generic codepath.
+  return ITensor(arrtype(n, SiteType.(domain)), axes(n, domain))
 end
 
 end
diff --git a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
new file mode 100644
index 0000000..e21f6c1
--- /dev/null
+++ b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
@@ -0,0 +1,96 @@
+module QuantumOperatorDefinitionsSymmetrySectorsExt
+
+using BlockArrays: blocklasts, blocklengths
+using GradedUnitRanges: AbstractGradedUnitRange, GradedOneTo, gradedrange
+using LabelledNumbers: label, labelled, unlabel
+using QuantumOperatorDefinitions:
+  QuantumOperatorDefinitions,
+  @SiteType_str,
+  @GradingType_str,
+  SiteType,
+  GradingType,
+  OpName,
+  name
+using SymmetrySectors: ×, dual, SectorProduct, U1, Z
+
+function Base.axes(::OpName, domain::Tuple{Vararg{AbstractGradedUnitRange}})
+  return (domain..., dual.(domain)...)
+end
+
+sortedunion(a, b) = sort(union(a, b))
+function QuantumOperatorDefinitions.combine_axes(a1::GradedOneTo, a2::GradedOneTo)
+  return gradedrange(
+    map(blocklengths(a1), blocklengths(a2)) do s1, s2
+      l1 = unlabel(s1)
+      l2 = unlabel(s2)
+      @assert l1 == l2
+      labelled(l1, label(s1) × label(s2))
+    end,
+  )
+end
+QuantumOperatorDefinitions.combine_axes(a::GradedOneTo, b::Base.OneTo) = a
+QuantumOperatorDefinitions.combine_axes(a::Base.OneTo, b::GradedOneTo) = b
+
+function Base.AbstractUnitRange(::GradingType"N", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; N=U1(i - 1))) => 1, 1:length(t)))
+end
+function Base.AbstractUnitRange(::GradingType"Sz", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; Sz=U1(i - 1))) => 1, 1:length(t)))
+end
+function Base.AbstractUnitRange(::GradingType"Sz↑", t::SiteType)
+  return AbstractUnitRange(GradingType"Sz"(), t)
+end
+function Base.AbstractUnitRange(::GradingType"Sz↓", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; Sz=U1(-(i - 1)))) => 1, 1:length(t)))
+end
+
+function sector(gradingtype::GradingType, sec)
+  sectorname = Symbol(get(gradingtype, :name, name(gradingtype)))
+  return SectorProduct(NamedTuple{(sectorname,)}((sec,)))
+end
+
+function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+# TODO: Write in terms of `GradingType"Nf"` definition.
+function Base.AbstractUnitRange(s::GradingType"NfParity", t::SiteType"Fermion")
+  return gradedrange([sector(s, Z{2}(0)) => 1, sector(s, Z{2}(1)) => 1])
+end
+function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+function Base.AbstractUnitRange(s::GradingType"Sz↑", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+function Base.AbstractUnitRange(s::GradingType"Sz↓", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(-1)) => 1])
+end
+
+# TODO: Write in terms of `SiteType"Fermion"` definitions.
+function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Electron")
+  return gradedrange([
+    sector(s, U1(0)) => 1,
+    sector(s, U1(1)) => 1,
+    sector(s, U1(1)) => 1,
+    sector(s, U1(2)) => 1,
+  ])
+end
+# TODO: Write in terms of `GradingType"Nf"` definition.
+function Base.AbstractUnitRange(s::GradingType"NfParity", t::SiteType"Electron")
+  return gradedrange([
+    sector(s, Z{2}(0)) => 1,
+    sector(s, Z{2}(1)) => 1,
+    sector(s, Z{2}(1)) => 1,
+    sector(s, Z{2}(0)) => 1,
+  ])
+end
+function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Electron")
+  return gradedrange([
+    sector(s, U1(0)) => 1,
+    sector(s, U1(1)) => 1,
+    sector(s, U1(-1)) => 1,
+    sector(s, U1(0)) => 1,
+  ])
+end
+
+end
diff --git a/src/op.jl b/src/op.jl
index c7fbf43..7aee6a0 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -8,8 +8,8 @@ struct OpName{Name,Params}
 end
 name(::OpName{Name}) where {Name} = Name
 params(n::OpName) = getfield(n, :params)
-
 Base.getproperty(n::OpName, name::Symbol) = getfield(params(n), name)
+Base.get(t::OpName, name::Symbol, default) = get(params(t), name, default)
 
 OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...))
 
@@ -54,9 +54,14 @@ end
 # Generic to `StateName` or `OpName`.
 const StateOrOpName = Union{StateName,OpName}
 alias(n::StateOrOpName) = n
-function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Integer...)
+function (arrtype::Type{<:AbstractArray})(
+  n::StateOrOpName, domain::Union{Integer,AbstractUnitRange}...
+)
   return arrtype(n, domain)
 end
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{Integer}})
+  return arrtype(n, Base.oneto.(domain))
+end
 (arrtype::Type{<:AbstractArray})(n::StateOrOpName, ts::SiteType...) = arrtype(n, ts)
 function (n::StateOrOpName)(domain...)
   # TODO: Try one alias at a time?
@@ -87,32 +92,58 @@ function nsites(n::StateOrOpName)
   return nsites(n′)
 end
 
-function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}},
-  domain::Tuple{Vararg{Integer}},
-  a::AbstractArray{<:Any,N},
-) where {N}
-  # TODO: Check the dimensions.
-  return convert(arrtype, a)
+# TODO: This does some unwanted conversions, like turning
+# `Diagonal` dense.
+function array(a::AbstractArray, ax::Tuple{Vararg{AbstractUnitRange}})
+  return a[ax...]
 end
-function op_convert(
-  arrtype::Type{<:AbstractArray}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
-)
-  # TODO: Check the dimensions.
-  return convert(arrtype, a)
+
+function Base.axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}})
+  return (domain..., domain...)
+end
+function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{Integer}})
+  return axes(n, Base.OneTo.(domain))
+end
+function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{SiteType}})
+  return axes(n, AbstractUnitRange.(domain))
+end
+
+## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{AbstractUnitRange}})
+##   return (reverse(domain)..., domain...)
+## end
+
+function reversed_sites(n::StateOrOpName, domain)
+  return reverse_sites(n, reshape(n(domain...), length.(axes(n, reverse(domain)))))
+end
+function reverse_sites(n::OpName, a::AbstractArray)
+  ndomain = Int(ndims(a)//2)
+  perm1 = reverse(ntuple(identity, ndomain))
+  perm2 = perm1 .+ ndomain
+  perm = (perm1..., perm2...)
+  return permutedims(a, perm)
 end
-function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
-) where {N}
-  size = (domain..., domain...)
-  @assert length(size) == N
-  return convert(arrtype, reshape(a, size))
+
+function state_or_op_convert(
+  n::StateOrOpName,
+  arrtype::Type{<:AbstractArray},
+  domain::Tuple{Vararg{AbstractUnitRange}},
+  a::AbstractArray,
+)
+  ax = axes(n, domain)
+  a′ = reshape(a, length.(ax))
+  a′′ = array(a′, ax)
+  return convert(arrtype, a′′)
 end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{SiteType}})
-  return op_convert(arrtype, length.(domain), n(domain...))
+
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{SiteType}})
+  domain′ = AbstractUnitRange.(domain)
+  return state_or_op_convert(n, arrtype, domain′, reversed_sites(n, domain))
 end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{Integer}})
-  return op_convert(arrtype, domain, n(Int.(domain)...))
+function (arrtype::Type{<:AbstractArray})(
+  n::StateOrOpName, domain::Tuple{Vararg{AbstractUnitRange}}
+)
+  # TODO: Make `(::OpName)(domain...)` constructor process more general inputs.
+  return state_or_op_convert(n, arrtype, domain, reversed_sites(n, Int.(length.(domain))))
 end
 
 function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
@@ -475,13 +506,13 @@ function (n::OpName"Controlled")(domain...)
   # Number of control sites.
   nc = get(params(n), :ncontrol, length(domain) - nt)
   @assert length(domain) == nc + nt
-  d_control = prod(to_dim.(domain[1:nc]))
+  d_control = prod(to_dim.(domain)) - prod(to_dim.(domain[(nc + 1):end]))
   return cat(I(d_control), n.arg(domain[(nc + 1):end]...); dims=(1, 2))
 end
-@op_alias "CNOT" "Controlled" op = OpName"X"()
-@op_alias "CX" "Controlled" op = OpName"X"()
-@op_alias "CY" "Controlled" op = OpName"Y"()
-@op_alias "CZ" "Controlled" op = OpName"Z"()
+@op_alias "CNOT" "Controlled" arg = OpName"X"()
+@op_alias "CX" "Controlled" arg = OpName"X"()
+@op_alias "CY" "Controlled" arg = OpName"Y"()
+@op_alias "CZ" "Controlled" arg = OpName"Z"()
 function alias(n::OpName"CPhase")
   return controlled(OpName"Phase"(; params(n)...))
 end
@@ -504,17 +535,17 @@ function alias(::OpName"CRn")
 end
 @op_alias "CRn̂" "CRn"
 
-@op_alias "CCNOT" "Controlled" ncontrol = 2 op = OpName"X"()
+@op_alias "CCNOT" "Controlled" ncontrol = 2 arg = OpName"X"()
 @op_alias "Toffoli" "CCNOT"
 @op_alias "CCX" "CCNOT"
 @op_alias "TOFF" "CCNOT"
 
-@op_alias "CSWAP" "Controlled" ncontrol = 2 op = OpName"SWAP"()
+@op_alias "CSWAP" "Controlled" ncontrol = 2 arg = OpName"SWAP"()
 @op_alias "Fredkin" "CSWAP"
 @op_alias "CSwap" "CSWAP"
 @op_alias "CS" "CSWAP"
 
-@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()
+@op_alias "CCCNOT" "Controlled" ncontrol = 3 arg = OpName"X"()
 
 ## # 1-qudit rotation around generic axis n̂.
 ## # exp(-im * α / 2 * n̂ ⋅ σ⃗)
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 0172f61..67936d9 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -4,29 +4,85 @@ struct SiteType{T,Params}
     return new{N,typeof(params)}(params)
   end
 end
+name(::SiteType{T}) where {T} = T
 params(t::SiteType) = getfield(t, :params)
 Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name)
-
+Base.get(t::SiteType, name::Symbol, default) = get(params(t), name, default)
+Base.haskey(t::SiteType, name::Symbol) = haskey(params(t), name)
 SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...))
-
 SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...)
 SiteType(i::Integer; kwargs...) = SiteType{Symbol(i)}(; kwargs...)
-value(::SiteType{T}) where {T} = T
 macro SiteType_str(s)
-  return SiteType{Symbol(s)}
+  return :(SiteType{$(Expr(:quote, Symbol(s)))})
 end
 
 alias(t::SiteType) = t
 alias(i::Integer) = i
 
+# Like `Base.Broadcast.axistype` (https://github.com/JuliaLang/julia/blob/v1.11.3/base/broadcast.jl#L536-L538)
+# and `BlockArrays.combine_blockaxes` (https://github.com/JuliaArrays/BlockArrays.jl/blob/v1.3.0/src/blockbroadcast.jl#L37-L38).
+combine_axes(a::T, b::T) where {T} = a
+combine_axes(a::Base.OneTo, b::Base.OneTo) = Base.OneTo{Int}(a)
+function combine_axes(a, b)
+  return UnitRange{Int}(a)
+end
+combine_axes(a) = a
+combine_axes(a, b, rest...) = combine_axes(combine_axes(a, b), rest...)
+
+struct GradingType{T,Params}
+  params::Params
+  function GradingType{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
+end
+name(::GradingType{T}) where {T} = T
+params(t::GradingType) = getfield(t, :params)
+Base.getproperty(t::GradingType, name::Symbol) = getfield(params(t), name)
+Base.get(t::GradingType, name::Symbol, default) = get(params(t), name, default)
+Base.haskey(t::GradingType, name::Symbol) = haskey(params(t), name)
+GradingType{N}(; kwargs...) where {N} = GradingType{N}((; kwargs...))
+GradingType(s::AbstractString; kwargs...) = GradingType{Symbol(s)}(; kwargs...)
+function GradingType(s::Pair{<:AbstractString,<:AbstractString}; kwargs...)
+  return GradingType(first(s); kwargs..., name=last(s))
+end
+function GradingType(s::Pair{<:AbstractString,<:NamedTuple}; kwargs...)
+  return GradingType(first(s); kwargs..., last(s)...)
+end
+macro GradingType_str(s)
+  return :(GradingType{$(Expr(:quote, Symbol(s)))})
+end
+
+function Base.AbstractUnitRange(grading::GradingType, t::SiteType)
+  return error("Not implemented.")
+end
+function Base.AbstractUnitRange(grading::GradingType"Trivial", t::SiteType)
+  return Base.OneTo(length(t))
+end
+
 function Base.length(t::SiteType)
   t′ = alias(t)
   if t == t′
-    return t.length
+    return t.dim
   end
   return length(t′)
 end
-Base.AbstractUnitRange(t::SiteType) = Base.OneTo(length(t))
+# TODO: Use a shorthand `(t::SiteType)() = AbstractUnitRange(t)`,
+# i.e. make `SiteType` callable like `OpName` and `StateName`
+# are right now.
+function Base.AbstractUnitRange(t::SiteType)
+  # This logic allows specifying a range with extra properties,
+  # like ones with symmetry sectors.
+  haskey(t, :range) && return t.range
+  if haskey(t, :gradings)
+    rs = map(grading -> AbstractUnitRange(GradingType(grading), t), t.gradings)
+    return combine_axes(Base.OneTo(length(t)), rs...)
+  end
+  return Base.OneTo(length(t))
+end
+# kwargs are passed for fancier constructors, like `ITensors.Index`.
+function (rangetype::Type{<:AbstractUnitRange})(t::SiteType)
+  return rangetype(AbstractUnitRange(t))
+end
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]
 Base.axes(t::SiteType) = (AbstractUnitRange(t),)
diff --git a/src/state.jl b/src/state.jl
index d9904a5..56262c3 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -6,14 +6,15 @@ struct StateName{Name,Params}
     return new{N,typeof(params)}(params)
   end
 end
+name(::StateName{N}) where {N} = N
 params(n::StateName) = getfield(n, :params)
 Base.getproperty(n::StateName, name::Symbol) = getfield(params(n), name)
+Base.get(t::StateName, name::Symbol, default) = get(params(t), name, default)
 
 StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...))
 
 StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...)
 StateName(s::Symbol; kwargs...) = StateName{s}(; kwargs...)
-name(::StateName{N}) where {N} = N
 macro StateName_str(s)
   return StateName{Symbol(s)}
 end
@@ -31,17 +32,6 @@ macro state_alias(name1, name2, params...)
   return state_alias_expr(name1, name2)
 end
 
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain::Tuple{Vararg{SiteType}})
-  # TODO: Define `state_convert` to handle reshaping multisite states
-  # to higher order arrays.
-  return convert(arrtype, n(domain...))
-end
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain::Tuple{Vararg{Integer}})
-  # TODO: Define `state_convert` to handle reshaping multisite states
-  # to higher order arrays.
-  return convert(arrtype, n(Int.(domain)...))
-end
-
 # This compiles operator expressions, such as:
 # ```julia
 # stateexpr("0 + 1") == StateName("0") + StateName("1")
@@ -126,6 +116,15 @@ function state_or_op_expr(ntype::Type, ex::Expr, depth::Int)
   return error("Can't parse expression $ex.")
 end
 
+function Base.axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
+  return domain
+end
+
+function reverse_sites(n::StateName, a::AbstractArray)
+  perm = reverse(ntuple(identity, ndims(a)))
+  return permutedims(a, perm)
+end
+
 function state(arrtype::Type{<:AbstractArray}, n::Union{Int,String}, domain...; kwargs...)
   return arrtype(stateexpr(n; kwargs...), domain...)
 end
diff --git a/test/Project.toml b/test/Project.toml
index 32b3272..69e5ba0 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,11 +1,16 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 53f2595..dfafddc 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -33,22 +33,30 @@ const elts = (real_elts..., complex_elts...)
         [0 -im; im 0],
         [1 0; 0 -1],
         [0 0; 0 1],
-        [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1],
-        [1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1],
-        (_, θ) -> [
-          cos(θ / 2) 0 0 -im*sin(θ / 2)
-          0 cos(θ / 2) -im*sin(θ / 2) 0
-          0 -im*sin(θ / 2) cos(θ / 2) 0
-          -im*sin(θ / 2) 0 0 cos(θ / 2)
-        ],
-        (_, θ) -> [
-          cos(θ / 2) 0 0 im*sin(θ / 2)
-          0 cos(θ / 2) -im*sin(θ / 2) 0
-          0 -im*sin(θ / 2) cos(θ / 2) 0
-          im*sin(θ / 2) 0 0 cos(θ / 2)
-        ],
-        (_, θ) ->
+        reshape([1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1], (2, 2, 2, 2)),
+        reshape([1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1], (2, 2, 2, 2)),
+        (_, θ) -> reshape(
+          [
+            cos(θ / 2) 0 0 -im*sin(θ / 2)
+            0 cos(θ / 2) -im*sin(θ / 2) 0
+            0 -im*sin(θ / 2) cos(θ / 2) 0
+            -im*sin(θ / 2) 0 0 cos(θ / 2)
+          ],
+          (2, 2, 2, 2),
+        ),
+        (_, θ) -> reshape(
+          [
+            cos(θ / 2) 0 0 im*sin(θ / 2)
+            0 cos(θ / 2) -im*sin(θ / 2) 0
+            0 -im*sin(θ / 2) cos(θ / 2) 0
+            im*sin(θ / 2) 0 0 cos(θ / 2)
+          ],
+          (2, 2, 2, 2),
+        ),
+        (_, θ) -> reshape(
           Diagonal([exp(-im * θ / 2), exp(im * θ / 2), exp(im * θ / 2), exp(-im * θ / 2)]),
+          (2, 2, 2, 2),
+        ),
         [1 0; 0 0],
         [0 0; 0 1],
         [0 1; 0 0],
@@ -60,31 +68,37 @@ const elts = (real_elts..., complex_elts...)
         √2 * [0 -im 0; im 0 -im; 0 im 0],
         2 * [1 0 0; 0 0 0; 0 0 -1],
         [0 0 0; 0 1 0; 0 0 2],
-        [
-          1 0 0 0 0 0 0 0 0
-          0 0 0 1 0 0 0 0 0
-          0 0 0 0 0 0 1 0 0
-          0 1 0 0 0 0 0 0 0
-          0 0 0 0 1 0 0 0 0
-          0 0 0 0 0 0 0 1 0
-          0 0 1 0 0 0 0 0 0
-          0 0 0 0 0 1 0 0 0
-          0 0 0 0 0 0 0 0 1
-        ],
-        [
-          1 0 0 0 0 0 0 0 0
-          0 0 0 im 0 0 0 0 0
-          0 0 0 0 0 0 im 0 0
-          0 im 0 0 0 0 0 0 0
-          0 0 0 0 1 0 0 0 0
-          0 0 0 0 0 0 0 im 0
-          0 0 im 0 0 0 0 0 0
-          0 0 0 0 0 im 0 0 0
-          0 0 0 0 0 0 0 0 1
-        ],
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
+        reshape(
+          [
+            1 0 0 0 0 0 0 0 0
+            0 0 0 1 0 0 0 0 0
+            0 0 0 0 0 0 1 0 0
+            0 1 0 0 0 0 0 0 0
+            0 0 0 0 1 0 0 0 0
+            0 0 0 0 0 0 0 1 0
+            0 0 1 0 0 0 0 0 0
+            0 0 0 0 0 1 0 0 0
+            0 0 0 0 0 0 0 0 1
+          ],
+          (3, 3, 3, 3),
+        ),
+        reshape(
+          [
+            1 0 0 0 0 0 0 0 0
+            0 0 0 im 0 0 0 0 0
+            0 0 0 0 0 0 im 0 0
+            0 im 0 0 0 0 0 0 0
+            0 0 0 0 1 0 0 0 0
+            0 0 0 0 0 0 0 im 0
+            0 0 im 0 0 0 0 0 0
+            0 0 0 0 0 im 0 0 0
+            0 0 0 0 0 0 0 0 1
+          ],
+          (3, 3, 3, 3),
+        ),
+        (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)),
+        (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)),
+        (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)),
         [1 0 0; 0 0 0; 0 0 0],
         [0 0 0; 0 1 0; 0 0 0],
         [0 1 0; 0 0 0; 0 0 0],
@@ -114,9 +128,9 @@ const elts = (real_elts..., complex_elts...)
         (OpName("Ry"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Ymat)),
         (OpName("Rz"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Zmat)),
         (OpName("SWAP"), 2, elts, SWAPmat),
-        (OpName("√SWAP"), 2, complex_elts, √SWAPmat),
+        # (OpName("√SWAP"), 2, complex_elts, √SWAPmat),
         (OpName("iSWAP"), 2, complex_elts, iSWAPmat),
-        (OpName("√iSWAP"), 2, complex_elts, √iSWAPmat),
+        # (OpName("√iSWAP"), 2, complex_elts, √iSWAPmat),
         (OpName("Rxx"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)),
         (OpName("RXX"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)),
         (OpName("Ryy"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)),
@@ -128,7 +142,7 @@ const elts = (real_elts..., complex_elts...)
         (OpName("StandardBasis"; index=(1, 2)), 1, elts, StandardBasis12mat),
       )
         @test nsites(o) == nbits
-        for arraytype in (AbstractArray, AbstractMatrix, Array, Matrix)
+        for arraytype in (AbstractArray, Array)
           for elt in elts
             ts = ntuple(Returns(t), nbits)
             lens = ntuple(Returns(len), nbits)
@@ -148,8 +162,10 @@ const elts = (real_elts..., complex_elts...)
     @test op("X * Y + Z") == op("X") * op("Y") + op("Z")
     @test op("X * Y + 2 * Z") == op("X") * op("Y") + 2 * op("Z")
     @test op("exp(im * (X * Y + 2 * Z))") == exp(im * (op("X") * op("Y") + 2 * op("Z")))
-    @test op("exp(im * (X ⊗ Y + Z ⊗ Z))") ==
-      exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z"))))
+    @test op("exp(im * (X ⊗ Y + Z ⊗ Z))") == permutedims(
+      reshape(exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))), (2, 2, 2, 2)),
+      (2, 1, 4, 3),
+    )
     @test op("Ry{θ=π/2}") == op("Ry"; θ=π / 2)
     # Awkward parsing corner cases.
     @test op("S+") == Matrix(OpName("S+"))
@@ -187,7 +203,19 @@ const elts = (real_elts..., complex_elts...)
     @test state("2", 3) == [0, 0, 1]
 
     @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+")
-    @test state("|0⟩ ⊗ |+⟩") == kron(state("0"), state("+"))
+    @test state("|0⟩ ⊗ |+⟩") ==
+      permutedims(reshape(kron(state("0"), state("+")), (2, 2)), (2, 1))
+  end
+  @testset "Ranges" begin
+    @test AbstractUnitRange(SiteType("S=1/2")) == Base.OneTo(2)
+    @test UnitRange{Int32}(SiteType("S=1/2")) == Base.OneTo(2)
+    @test UnitRange{Int32}(SiteType("S=1/2")) isa UnitRange{Int32}
+    @test AbstractUnitRange(SiteType("Qudit"; dim=3)) == Base.OneTo(3)
+    @test UnitRange{Int32}(SiteType("Qudit"; dim=3)) == Base.OneTo(3)
+    @test UnitRange{Int32}(SiteType("Qudit"; dim=3)) isa UnitRange{Int32}
+
+    @test op("X", 2) == op("X", Base.OneTo(2))
+    @test op("X", 3) == op("X", Base.OneTo(3))
   end
   @testset "Electron/tJ" begin
     for (ns, x) in (
diff --git a/test/test_itensorbaseext.jl b/test/test_itensorbaseext.jl
index e71fe46..cb274b5 100644
--- a/test/test_itensorbaseext.jl
+++ b/test/test_itensorbaseext.jl
@@ -1,13 +1,64 @@
-using ITensorBase: ITensor, Index, prime
-using QuantumOperatorDefinitions: op, state
+using ITensorBase: ITensor, Index, gettag, prime, settag
+using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state
 using Test: @test, @testset
 
-@testset "QuantumOperatorDefinitionsITensorBaseExt" begin
+@testset "ITensorBaseExt" begin
+  i = Index(SiteType("S=1/2"))
+  @test gettag(i, "sitetype") == "S=1/2"
+
   i = Index(2)
+  a = state("0", i)
+  @test a == ITensor(state("0", 2), i)
+  @test a == state("0", (i,))
+  @test a == ITensor(StateName("0"), i)
+  @test a == ITensor(StateName("0"), (i,))
+
+  i = settag(Index(2), "sitetype", "S=1/2")
+  a = state("X+", i)
+  @test a == ITensor(state("X+", SiteType("S=1/2")), i)
 
+  i = Index(2)
   a = op("X", i)
-  @test a == ITensor([0 1; 1 0], (prime(i), i))
+  @test a == ITensor(op("X", 2), (prime(i), i))
+  @test a == op("X", (i,))
+  @test a == ITensor(OpName("X"), i)
+  @test a == ITensor(OpName("X"), (i,))
+
+  i1 = Index(2)
+  i2 = Index(2)
+  i1′ = prime(i1)
+  i2′ = prime(i2)
+  a = ITensor(OpName("CX"), (i1, i2))
+  @test a == ITensor(op("CX", (2, 2)), (i1′, i2′, i1, i2))
+  @test a[i1′[1], i2′, i1[1], i2] == op("I", i2)
+  @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2)
+  @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2)
+  @test a[i1′[2], i2′, i1[2], i2] == op("X", i2)
+
+  i1 = Index(3)
+  i2 = Index(2)
+  i1′ = prime(i1)
+  i2′ = prime(i2)
+  a = ITensor(OpName("CX"), (i1, i2))
+  @test a == ITensor(op("CX", (3, 2)), (i1′, i2′, i1, i2))
+  @test a[i1′[1], i2′, i1[1], i2] == op("I", i2)
+  @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2)
+  @test a[i1′[3], i2′, i1[1], i2] == zeros(i2′, i2)
+  @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2)
+  @test a[i1′[2], i2′, i1[2], i2] == op("I", i2)
+  @test a[i1′[3], i2′, i1[2], i2] == zeros(i2′, i2)
+  @test a[i1′[1], i2′, i1[3], i2] == zeros(i2′, i2)
+  @test a[i1′[2], i2′, i1[3], i2] == zeros(i2′, i2)
+  @test a[i1′[3], i2′, i1[3], i2] == op("X", i2)
 
-  a = state(1, i)
-  @test a == ITensor([1, 0], (i,))
+  i1 = Index(2)
+  i2 = Index(3)
+  i1′ = prime(i1)
+  i2′ = prime(i2)
+  a = ITensor(OpName("CX"), (i1, i2))
+  @test a == ITensor(op("CX", (2, 3)), (i1′, i2′, i1, i2))
+  @test a[i1′[1], i2′, i1[1], i2] == op("I", i2)
+  @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2)
+  @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2)
+  @test a[i1′[2], i2′, i1[2], i2] == op("X", i2)
 end
diff --git a/test/test_symmetrysectorsext.jl b/test/test_symmetrysectorsext.jl
new file mode 100644
index 0000000..e7da36d
--- /dev/null
+++ b/test/test_symmetrysectorsext.jl
@@ -0,0 +1,113 @@
+using BlockArrays: AbstractBlockArray, blocklengths
+using BlockSparseArrays: BlockSparseArray
+using GradedUnitRanges: blocklabels
+using ITensorBase: ITensor, Index, gettag, prime, settag
+using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state
+using SymmetrySectors: SectorProduct, U1, Z
+using NamedDimsArrays: dename
+using Test: @test, @test_broken, @testset
+
+@testset "SymmetrySectorsExt" begin
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  r = AbstractUnitRange(t)
+  @test r == 1:2
+  @test blocklabels(r) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r) == [1, 1]
+
+  t = SiteType("Electron"; gradings=("Nf", "Sz"))
+  r = AbstractUnitRange(t)
+  @test r == 1:4
+  @test blocklabels(r) == [
+    SectorProduct((; Nf=U1(0), Sz=U1(0))),
+    SectorProduct((; Nf=U1(1), Sz=U1(1))),
+    SectorProduct((; Nf=U1(1), Sz=U1(-1))),
+    SectorProduct((; Nf=U1(2), Sz=U1(0))),
+  ]
+  @test blocklengths(r) == [1, 1, 1, 1]
+
+  t = SiteType("Electron"; gradings=("Nf" => "NfA", "Sz" => "SzA"))
+  r = AbstractUnitRange(t)
+  @test r == 1:4
+  @test blocklabels(r) == [
+    SectorProduct((; NfA=U1(0), SzA=U1(0))),
+    SectorProduct((; NfA=U1(1), SzA=U1(1))),
+    SectorProduct((; NfA=U1(1), SzA=U1(-1))),
+    SectorProduct((; NfA=U1(2), SzA=U1(0))),
+  ]
+  @test blocklengths(r) == [1, 1, 1, 1]
+
+  t = SiteType("Electron"; gradings=("NfParity", "Sz"))
+  r = AbstractUnitRange(t)
+  @test r == 1:4
+  @test blocklabels(r) == [
+    SectorProduct((; NfParity=Z{2}(0), Sz=U1(0))),
+    SectorProduct((; NfParity=Z{2}(1), Sz=U1(1))),
+    SectorProduct((; NfParity=Z{2}(1), Sz=U1(-1))),
+    SectorProduct((; NfParity=Z{2}(0), Sz=U1(0))),
+  ]
+  @test blocklengths(r) == [1, 1, 1, 1]
+
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  (r1, r2) = axes(OpName("σ⁺"), (t,))
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  @test blocklabels(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+
+  # TODO: There is a bug slicing `BitVector` by `GradedOneTo` in Julia 1.11,
+  # investigate. See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  @test state("0", t) == [1, 0] broken = VERSION ≥ v"1.11"
+
+  # Force conversion to `Vector{Float64}` before conversion,
+  # since there is a bug slicing `BitVector` by `GradedOneTo`.
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  a = AbstractArray(2.0 * StateName("0"), t)
+  @test a == [2, 0]
+  @test a isa AbstractBlockArray
+  # TODO: Currently slicing a dense array by graded ranges outputs a `BlockedArray`
+  # rather than a `BlockSparseArray`.
+  # See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
+  @test_broken a isa BlockSparseArray
+  (r1,) = axes(a)
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  a = op("σ⁺", t)
+  @test a == [0 2; 0 0]
+  @test a isa AbstractBlockArray
+  @test_broken a isa BlockSparseArray
+  (r1, r2) = axes(a)
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  # TODO: This is a bug in indexing with GradedUnitRangeDual, fix this.
+  @test_broken blocklabels(r2) ==
+    [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+end
+
+@testset "SymmetrySectorsExt + ITensorBaseExt" begin
+  i = Index(SiteType("S=1/2"; gradings=("Sz",)))
+  @test gettag(i, "sitetype") == "S=1/2"
+  # TODO: Test without denaming.
+  @test dename(i) == 1:2
+  @test blocklabels(dename(i)) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(dename(i)) == [1, 1]
+
+  i′ = prime(i)
+  a = op("σ⁺", i)
+  # TODO: The indices should be `(i′, dual(i))`.
+  @test a == ITensor([0 2; 0 0], (i′, i))
+  a′ = dename(a)
+  @test a′ isa AbstractBlockArray
+  @test_broken a′ isa BlockSparseArray
+  # TODO: Test these without denaming `a`.
+  (r1, r2) = axes(a′)
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  # TODO: This is a bug in indexing with GradedUnitRangeDual, fix this.
+  @test_broken blocklabels(r2) ==
+    [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+end