diff --git a/Project.toml b/Project.toml
index 430819a..c882343 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,13 +1,20 @@
 name = "QuantumOperatorDefinitions"
 uuid = "826dd319-6fd5-459a-a990-3a4f214664bf"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[weakdeps]
+ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
+
+[extensions]
+QuantumOperatorDefinitionsITensorBaseExt = "ITensorBase"
+
 [compat]
+ITensorBase = "0.1.10"
 LinearAlgebra = "1.10"
 Random = "1.10"
 julia = "1.10"
diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
new file mode 100644
index 0000000..280b531
--- /dev/null
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -0,0 +1,24 @@
+module QuantumOperatorDefinitionsITensorBaseExt
+
+using ITensorBase: ITensor, Index, dag, gettag, prime
+using QuantumOperatorDefinitions:
+  QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
+
+function QuantumOperatorDefinitions.SiteType(r::Index)
+  return SiteType(gettag(r, "sitetype", "Qudit"); dim=Int(length(r)))
+end
+
+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)))
+end
+
+function Base.AbstractArray(n::StateName, r::Index)
+  return ITensor(AbstractArray(n, SiteType(r)), (r,))
+end
+
+end
diff --git a/src/has_fermion_string.jl b/src/has_fermion_string.jl
index 9620299..614f940 100644
--- a/src/has_fermion_string.jl
+++ b/src/has_fermion_string.jl
@@ -1 +1,7 @@
-has_fermion_string(::OpName, ::SiteType) = nothing
+function has_fermion_string(n::OpName, t::SiteType)
+  n′ = alias(n)
+  if n == n′
+    return false
+  end
+  return has_fermion_string(n′, t)
+end
diff --git a/src/op.jl b/src/op.jl
index 57c9396..13b23b7 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -1,60 +1,25 @@
-using Random: randstring
+using LinearAlgebra: diag, qr
 
 struct OpName{Name,Params}
-  function OpName{Name,Params}(params::NamedTuple) where {Name,Params}
-    return new{Name,(; Params..., params...)}()
+  params::Params
+  function OpName{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
   end
 end
 name(::OpName{Name}) where {Name} = Name
-params(::OpName{<:Any,Params}) where {Params} = Params
+params(n::OpName) = getfield(n, :params)
 
 Base.getproperty(n::OpName, name::Symbol) = getfield(params(n), name)
 
-OpName{Name,Params}(; kwargs...) where {Name,Params} = OpName{Name,Params}((; kwargs...))
-
-OpName{N}(params::NamedTuple) where {N} = OpName{N,params}()
 OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...))
 
-const DAGGER_STRING = randstring()
-const UPARROW_STRING = randstring()
-const DOWNARROW_STRING = randstring()
-const OPEXPR_REPLACEMENTS = (
-  "†" => DAGGER_STRING, "↑" => UPARROW_STRING, "↓" => DOWNARROW_STRING
-)
-
 # This compiles operator expressions, such as:
 # ```julia
 # opexpr("X + Y") == OpName("X") + OpName("Y")
 # opexpr("Ry{θ=π/2}") == OpName("Ry"; θ=π/2)
 # ```
 function opexpr(n::String; kwargs...)
-  n = replace(n, OPEXPR_REPLACEMENTS...)
-  return opexpr(Meta.parse(n); kwargs...)
-end
-opexpr(n::Number) = n
-function opexpr(n::Symbol; kwargs...)
-  n === :im && return im
-  n === :π && return π
-  n = Symbol(replace(String(n), reverse.(OPEXPR_REPLACEMENTS)...))
-  return OpName{n}(; kwargs...)
-end
-function opexpr(ex::Expr)
-  if Meta.isexpr(ex, :call)
-    return eval(ex.args[1])(opexpr.(ex.args[2:end])...)
-  end
-  if Meta.isexpr(ex, :curly)
-    # Syntax for parametrized gates, i.e.
-    # `opexpr("Ry{θ=π/2}")`.
-    params = ex.args[2:end]
-    kwargs = Dict(
-      map(params) do param
-        @assert Meta.isexpr(param, :(=))
-        return param.args[1] => eval(param.args[2])
-      end,
-    )
-    return OpName{ex.args[1]}(; kwargs...)
-  end
-  return error("Can't parse expression $ex.")
+  return state_or_op_expr(OpName, n; kwargs...)
 end
 
 # TODO: Should this parse the string?
@@ -73,91 +38,9 @@ end
 # macro opexpr_str(s)
 #   return :(typeof(opexpr($s)))
 # end
-
-for f in (
-  :(Base.sqrt),
-  :(Base.real),
-  :(Base.imag),
-  :(Base.complex),
-  :(Base.exp),
-  :(Base.cis),
-  :(Base.cos),
-  :(Base.sin),
-  :(Base.adjoint),
-  :(Base.:+),
-  :(Base.:-),
-)
-  @eval begin
-    $f(n::OpName) = OpName"f"(; f=$f, op=n)
-  end
-end
-
-# Unary operations
-nsites(n::OpName"f") = nsites(n.op)
-function Base.AbstractArray(n::OpName"f", domain_size::Tuple{Vararg{Int}})
-  return n.f(AbstractArray(n.op, domain_size))
-end
-
-nsites(n::OpName"^") = nsites(n.op)
-function Base.AbstractArray(n::OpName"^", domain_size::Tuple{Vararg{Int}})
-  return AbstractArray(n.op, domain_size)^n.exponent
-end
-Base.:^(n::OpName, exponent) = OpName"^"(; op=n, exponent)
-
-nsites(n::OpName"kron") = nsites(n.op1) + nsites(n.op2)
-function Base.AbstractArray(n::OpName"kron", domain_size::Tuple{Vararg{Int}})
-  domain_size1 = domain_size[1:nsites(n.op1)]
-  domain_size2 = domain_size[(nsites(n.op1) + 1):end]
-  @assert length(domain_size2) == nsites(n.op2)
-  return kron(AbstractArray(n.op1, domain_size1), AbstractArray(n.op2, domain_size2))
-end
-Base.kron(n1::OpName, n2::OpName) = OpName"kron"(; op1=n1, op2=n2)
-⊗(n1::OpName, n2::OpName) = kron(n1, n2)
-
-function nsites(n::OpName"+")
-  @assert nsites(n.op1) == nsites(n.op2)
-  return nsites(n.op1)
-end
-function Base.AbstractArray(n::OpName"+", domain_size::Tuple{Vararg{Int}})
-  return AbstractArray(n.op1, domain_size) + AbstractArray(n.op2, domain_size)
-end
-Base.:+(n1::OpName, n2::OpName) = OpName"+"(; op1=n1, op2=n2)
-Base.:-(n1::OpName, n2::OpName) = n1 + (-n2)
-
-function nsites(n::OpName"*")
-  @assert nsites(n.op1) == nsites(n.op2)
-  return nsites(n.op1)
-end
-function Base.AbstractArray(n::OpName"*", domain_size::Tuple{Vararg{Int}})
-  return AbstractArray(n.op1, domain_size) * AbstractArray(n.op2, domain_size)
-end
-Base.:*(n1::OpName, n2::OpName) = OpName"*"(; op1=n1, op2=n2)
-
-nsites(n::OpName"scaled") = nsites(n.op)
-function Base.AbstractArray(n::OpName"scaled", domain_size::Tuple{Vararg{Int}})
-  return AbstractArray(n.op, domain_size) * n.c
-end
-function Base.:*(c::Number, n::OpName)
-  return OpName"scaled"(; op=n, c)
-end
-function Base.:*(n::OpName, c::Number)
-  return OpName"scaled"(; op=n, c)
-end
-function Base.:/(n::OpName, c::Number)
-  return OpName"scaled"(; op=n, c=inv(c))
-end
-
-function Base.:*(c::Number, n::OpName"scaled")
-  return OpName"scaled"(; op=n.op, c=(c * n.c))
-end
-function Base.:*(n::OpName"scaled", c::Number)
-  return OpName"scaled"(; op=n.op, c=(n.c * c))
-end
-function Base.:/(n::OpName"scaled", c::Number)
-  return OpName"scaled"(; op=n.op, c=(n.c / c))
-end
-
-controlled(n::OpName; ncontrol=1) = OpName"Controlled"(; ncontrol, op=n)
+# macro statexpr_str(s)
+#   return :(typeof(stateexpr($s)))
+# end
 
 function op_alias_expr(name1, name2, pars...)
   return :(function alias(n::OpName{Symbol($name1)})
@@ -168,63 +51,68 @@ macro op_alias(name1, name2, pars...)
   return op_alias_expr(name1, name2, pars...)
 end
 
-alias(n::OpName) = n
+# Generic to `StateName` or `OpName`.
+const StateOrOpName = Union{StateName,OpName}
+alias(n::StateOrOpName) = n
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Integer...)
+  return arrtype(n, domain)
+end
+(arrtype::Type{<:AbstractArray})(n::StateOrOpName, ts::SiteType...) = arrtype(n, ts)
+function (n::StateOrOpName)(domain...)
+  # TODO: Try one alias at a time?
+  # TODO: First call `alias(n, domain...)`
+  # to allow for aliases specific to certain
+  # SiteTypes?
+  n′ = alias(n)
+  domain′ = alias.(domain)
+  if n′ == n && domain′ == domain
+    error("Not implemented.")
+  end
+  return n′(domain′...)
+end
+# TODO: Decide on this.
+function (n::StateOrOpName)()
+  return n(ntuple(Returns(default_sitetype()), nsites(n))...)
+end
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName)
+  return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n)))
+end
+# `Int(ndims // 2)`, i.e. `ndims_domain`/`ndims_codomain`.
+function nsites(n::StateOrOpName)
+  n′ = alias(n)
+  if n′ == n
+    # Default value, assume 1-site operator/state.
+    return 1
+  end
+  return nsites(n′)
+end
+
 function op_convert(
   arrtype::Type{<:AbstractArray{<:Any,N}},
-  domain_size::Tuple{Vararg{Integer}},
+  domain::Tuple{Vararg{Integer}},
   a::AbstractArray{<:Any,N},
 ) where {N}
   # TODO: Check the dimensions.
   return convert(arrtype, a)
 end
 function op_convert(
-  arrtype::Type{<:AbstractArray}, domain_size::Tuple{Vararg{Integer}}, a::AbstractArray
+  arrtype::Type{<:AbstractArray}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
 )
   # TODO: Check the dimensions.
   return convert(arrtype, a)
 end
 function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}},
-  domain_size::Tuple{Vararg{Integer}},
-  a::AbstractArray,
+  arrtype::Type{<:AbstractArray{<:Any,N}}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
 ) where {N}
-  size = (domain_size..., domain_size...)
+  size = (domain..., domain...)
   @assert length(size) == N
   return convert(arrtype, reshape(a, size))
 end
-function (arrtype::Type{<:AbstractArray})(n::OpName, ts::Tuple{Vararg{SiteType}})
-  return op_convert(arrtype, length.(ts), AbstractArray(n, ts))
+function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{SiteType}})
+  return op_convert(arrtype, length.(domain), n(domain...))
 end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain_size::Tuple{Vararg{Integer}})
-  return op_convert(arrtype, domain_size, AbstractArray(n, Int.(domain_size)))
-end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain_size::Integer...)
-  return arrtype(n, domain_size)
-end
-(arrtype::Type{<:AbstractArray})(n::OpName, ts::SiteType...) = arrtype(n, ts)
-Base.AbstractArray(n::OpName, ts::SiteType...) = AbstractArray(n, ts)
-function Base.AbstractArray(n::OpName, ts::Tuple{Vararg{SiteType}})
-  n′ = alias(n)
-  ts′ = alias.(ts)
-  if n′ == n && ts′ == ts
-    return AbstractArray(n′, length.(ts′))
-  end
-  return AbstractArray(n′, ts′)
-end
-function Base.AbstractArray(n::OpName, domain_size::Tuple{Vararg{Int}})
-  n′ = alias(n)
-  if n′ == n
-    error("Not implemented.")
-  end
-  return AbstractArray(n′, domain_size)
-end
-
-# TODO: Decide on this.
-function Base.AbstractArray(n::OpName)
-  return AbstractArray(n, ntuple(Returns(default_sitetype()), nsites(n)))
-end
-function (arrtype::Type{<:AbstractArray})(n::OpName)
-  return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n)))
+function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{Integer}})
+  return op_convert(arrtype, domain, n(Int.(domain)...))
 end
 
 function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
@@ -237,22 +125,136 @@ function op(n::String, domain...; kwargs...)
   return op(AbstractArray, n, domain...; kwargs...)
 end
 
-# `Int(ndims // 2)`, i.e. `ndims_domain`/`ndims_codomain`.
-function nsites(n::Union{StateName,OpName})
-  n′ = alias(n)
-  if n′ == n
-    # Default value, assume 1-site operator/state.
-    return 1
+# Unary operations
+for nametype in (:StateName, :OpName)
+  applied = :($(nametype){:applied})
+  @eval begin
+    nsites(n::$(applied)) = nsites(n.arg)
+    function (n::$(applied))(domain...)
+      return n.f(n.arg(domain...))
+    end
+  end
+  for f in (
+    :(Base.real),
+    :(Base.imag),
+    :(Base.complex),
+    # :(Base.adjoint), # Decide on this, since this becomes a matrix.
+    :(Base.:+),
+    :(Base.:-),
+  )
+    @eval begin
+      $f(n::$(nametype)) = $(applied)(; f=$f, arg=n)
+    end
+  end
+end
+
+for nametype in (:StateName, :OpName)
+  kronned = :($(nametype){:kronned})
+  @eval begin
+    nsites(n::$(kronned)) = sum(nsites, n.args)
+    function (n::$(kronned))(domain...)
+      # TODO: Generalize beyond two arguments.
+      # Use `cumsum(nsites.(n.args))`.
+      stops = cumsum(nsites.(n.args))
+      starts = [1, stops[1:(end - 1)] .+ 1...]
+      domains = map((start, stop) -> domain[start:stop], starts, stops)
+      return kron(map((arg, domain) -> arg(domain...), n.args, domains)...)
+    end
+    Base.kron(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) =
+      $(kronned)(; args=(n1, n2, n_rest...))
+    ⊗(n1::$(nametype), n2::$(nametype)) = kron(n1, n2)
+    ⊗(n1::$(kronned), n2::$(kronned)) = kron(n1.args..., n2.args...)
+    ⊗(n1::$(nametype), n2::$(kronned)) = kron(n1, n2.args...)
+    ⊗(n1::$(kronned), n2::$(nametype)) = kron(n1.args..., n2)
+  end
+end
+
+for nametype in (:StateName, :OpName)
+  summed = :($(nametype){:summed})
+  @eval begin
+    function nsites(n::$(summed))
+      # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support.
+      @assert allequal(nsites.(n.args))
+      return nsites(first(n.args))
+    end
+    function (n::$(summed))(domain...)
+      return mapreduce(a -> a(domain...), +, n.args)
+    end
+    Base.:+(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) =
+      $(summed)(; args=(n1, n2, n_rest...))
+    Base.:+(n1::$(summed), n2::$(summed)) = +(n1.args..., n2.args...)
+    Base.:+(n1::$(summed), n2::$(nametype)) = +(n1.args..., n2)
+    Base.:+(n1::$(nametype), n2::$(summed)) = +(n1, n2.args...)
+    Base.:-(n1::$(nametype), n2::$(nametype)) = n1 + (-n2)
+  end
+end
+
+for nametype in (:StateName, :OpName)
+  scaled = :($(nametype){:scaled})
+  @eval begin
+    nsites(n::$(scaled)) = nsites(n.arg)
+    function (n::$(scaled))(domain...)
+      return n.arg(domain...) * n.c
+    end
+    function Base.:*(c::Number, n::$(nametype))
+      return $(scaled)(; arg=n, c)
+    end
+    function Base.:*(n::$(nametype), c::Number)
+      return $(scaled)(; arg=n, c)
+    end
+    function Base.:/(n::$(nametype), c::Number)
+      return $(scaled)(; arg=n, c=inv(c))
+    end
+
+    function Base.:*(c::Number, n::$(scaled))
+      return $(scaled)(; arg=n.arg, c=(c * n.c))
+    end
+    function Base.:*(n::$(scaled), c::Number)
+      return $(scaled)(; arg=n.arg, c=(n.c * c))
+    end
+    function Base.:/(n::$(scaled), c::Number)
+      return $(scaled)(; arg=n.arg, c=(n.c / c))
+    end
+  end
+end
+
+# Unary operations unique to operators.
+for f in (:(Base.sqrt), :(Base.exp), :(Base.cis), :(Base.cos), :(Base.sin), :(Base.adjoint))
+  @eval begin
+    $f(n::OpName) = OpName"applied"(; f=$f, arg=n)
   end
-  return nsites(n′)
 end
 
+nsites(n::OpName"exponentiated") = nsites(n.arg)
+function (n::OpName"exponentiated")(domain...)
+  return n.arg(domain...)^n.exponent
+end
+Base.:^(n::OpName, exponent) = OpName"exponentiated"(; arg=n, exponent)
+
+function nsites(n::OpName"producted")
+  # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support.
+  @assert allequal(nsites.(n.args))
+  return nsites(first(n.args))
+end
+function (n::OpName"producted")(domain...)
+  return mapreduce(a -> a(domain...), *, n.args)
+end
+function Base.:*(n1::OpName, n2::OpName, n_rest::OpName...)
+  return OpName"producted"(; args=(n1, n2, n_rest...))
+end
+Base.:*(n1::StateName"producted", n2::StateName"producted") = *(n1.args..., n2.args...)
+Base.:*(n1::StateName, n2::StateName"producted") = *(n1, n2.args...)
+Base.:*(n1::StateName"producted", n2::StateName) = *(n1.args..., n2)
+
+controlled(n::OpName; ncontrol=1) = OpName"Controlled"(; ncontrol, arg=n)
+
 using LinearAlgebra: Diagonal
-function Base.AbstractArray(::OpName"Id", domain_size::Tuple{Int})
-  return Diagonal(trues(only(domain_size)))
+function (::OpName"Id")(domain)
+  return Diagonal(trues(to_dim(domain)))
 end
-function Base.AbstractArray(n::OpName"Id", domain_size::Tuple{Int,Vararg{Int}})
-  return Base.AbstractArray(kron(ntuple(Returns(n), length(domain_size))...), domain_size)
+function (n::OpName"Id")(domain1, domain_rest...)
+  domain = (domain1, domain_rest)
+  return kron(ntuple(Returns(n), length(domain))...)(domain...)
 end
 @op_alias "I" "Id"
 @op_alias "σ0" "Id"
@@ -261,9 +263,20 @@ end
 # TODO: Is this a good definition?
 @op_alias "F" "Id"
 
-# TODO: Define as `::Tuple{Int}`.
-function Base.AbstractArray(n::OpName"a†", domain_size::Tuple{Int})
-  d = only(domain_size)
+function (n::OpName"StandardBasis")(domain)
+  d = to_dim(domain)
+  a = falses(d, d)
+  a[n.index...] = one(Bool)
+  return a
+end
+
+function alias(n::OpName"Proj")
+  return OpName"StandardBasis"(; index=(n.index, n.index))
+end
+@op_alias "proj" "Proj"
+
+function (n::OpName"a†")(domain)
+  d = to_dim(domain)
   a = zeros(d, d)
   for k in 1:(d - 1)
     a[k + 1, k] = √k
@@ -299,19 +312,21 @@ alias(::OpName"a†a†") = OpName("a†") ⊗ OpName("a†")
 # https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices
 # https://en.wikipedia.org/wiki/Generalized_Clifford_algebra
 # https://github.com/QuantumKitHub/MPSKitModels.jl/blob/v0.4.0/src/operators/spinoperators.jl
-function Base.AbstractArray(n::OpName"σ⁺", domain_size::Tuple{Int})
-  d = only(domain_size)
+function (n::OpName"σ⁺")(domain)
+  d = to_dim(domain)
   s = (d - 1) / 2
   return [2 * δ(i + 1, j) * √((s + 1) * (i + j - 1) - i * j) for i in 1:d, j in 1:d]
 end
-alias(::OpName"S⁺") = OpName("σ⁺") / 2
-@op_alias "Splus" "S⁺"
-@op_alias "Sp" "S⁺"
+alias(::OpName"S+") = OpName("σ⁺") / 2
+@op_alias "S⁺" "S+"
+@op_alias "Splus" "S+"
+@op_alias "Sp" "S+"
 
 alias(::OpName"σ⁻") = OpName"σ⁺"()'
-alias(::OpName"S⁻") = OpName("σ⁻") / 2
-@op_alias "Sminus" "S⁻"
-@op_alias "Sm" "S⁻"
+alias(::OpName"S-") = OpName("σ⁻") / 2
+@op_alias "S⁻" "S-"
+@op_alias "Sminus" "S-"
+@op_alias "Sm" "S-"
 
 alias(::OpName"X") = (OpName"σ⁺"() + OpName"σ⁻"()) / 2
 @op_alias "σx" "X"
@@ -329,9 +344,13 @@ alias(n::OpName"Sx") = OpName"X"() / 2
 @op_alias "Sₓ" "Sx"
 alias(::OpName"Sx2") = OpName"Sx"()^2
 
+# Generic rotation.
+# exp(-im * θ / 2 * O)
+alias(n::OpName"R") = cis(-(n.θ / 2) * n.arg)
+
 # Rotation around X-axis
-# exp(-im * n.θ / 2 * X)
-alias(n::OpName"Rx") = cis(-(n.θ / 2) * OpName"X"())
+# exp(-im * θ / 2 * X)
+alias(n::OpName"Rx") = OpName"R"(; params(n)..., arg=OpName"X"())
 
 alias(::OpName"Y") = -im * (OpName"σ⁺"() - OpName"σ⁻"()) / 2
 # TODO: No subsript `\_y` available
@@ -355,22 +374,22 @@ alias(n::OpName"iSy") = OpName"iY"() / 2
 alias(::OpName"Sy2") = -OpName"iSy"()^2
 
 # Rotation around Y-axis
-# exp(-im * n.θ / 2 * Y)
-alias(n::OpName"Ry") = exp(-(n.θ / 2) * OpName"iY"())
+# exp(-im * θ / 2 * Y) = exp(-θ / 2 * iY)
+alias(n::OpName"Ry") = OpName"R"(; params(n)..., arg=OpName"Y"())
 
 # Ising (XX) coupling gate
 # exp(-im * θ/2 * X ⊗ X)
-alias(n::OpName"Rxx") = exp(-im * (n.θ / 2) * OpName"X"() ⊗ OpName"X"())
+alias(n::OpName"Rxx") = OpName"R"(; params(n)..., arg=OpName"X"() ⊗ OpName"X"())
 @op_alias "RXX" "Rxx"
 
 # Ising (YY) coupling gate
 # exp(-im * θ/2 * Y ⊗ Y)
-alias(n::OpName"Ryy") = exp(-im * (n.θ / 2) * OpName"Y"() ⊗ OpName"Y"())
+alias(n::OpName"Ryy") = OpName"R"(; params(n)..., arg=OpName"Y"() ⊗ OpName"Y"())
 @op_alias "RYY" "Ryy"
 
 # Ising (ZZ) coupling gate
 # exp(-im * θ/2 * Z ⊗ Z)
-alias(n::OpName"Rzz") = exp(-im * (n.θ / 2) * OpName"Z"() ⊗ OpName"Z"())
+alias(n::OpName"Rzz") = OpName"R"(; params(n)..., arg=OpName"Z"() ⊗ OpName"Z"())
 @op_alias "RZZ" "Rzz"
 
 ## TODO: Check this definition and see if it is worth defining this.
@@ -379,8 +398,8 @@ alias(n::OpName"Rzz") = exp(-im * (n.θ / 2) * OpName"Z"() ⊗ OpName"Z"())
 ## alias(n::OpName"Rxy") = exp(-im * (n.θ / 2) * OpName"X"() ⊗ OpName"Y"())
 ## @op_alias "RXY" "Rxy"
 
-function Base.AbstractArray(n::OpName"σᶻ", domain_size::Tuple{Int})
-  d = only(domain_size)
+function (n::OpName"σᶻ")(domain)
+  d = to_dim(domain)
   s = (d - 1) / 2
   return Diagonal([2 * (s + 1 - i) for i in 1:d])
 end
@@ -407,23 +426,24 @@ alias(n::OpName"S2") = OpName"Sx2"() + OpName"Sy2"() + OpName"Sz2"()
 alias(::OpName"Sz2") = OpName"Sz"()^2
 
 # Rotation around Z-axis
-# exp(-im * n.θ / 2 * Z)
-alias(n::OpName"Rz") = exp(-im * (n.θ / 2) * OpName"Z"())
+# exp(-im * θ / 2 * Z)
+alias(n::OpName"Rz") = OpName"R"(; params(n)..., arg=OpName"Z"())
 
 using LinearAlgebra: eigen
-function Base.AbstractArray(n::OpName"H", domain_size::Tuple{Int})
-  Λ, H = eigen(AbstractArray(OpName("X"), domain_size))
+function (n::OpName"H")(domain)
+  Λ, H = eigen(OpName("X")(domain))
   p = sortperm(Λ; rev=true)
   return H[:, p]
 end
 
 using LinearAlgebra: Diagonal
 nsites(::OpName"SWAP") = 2
-function Base.AbstractArray(::OpName"SWAP", domain_size::Tuple{Int,Int})
-  I_matrix = Diagonal(trues(prod(domain_size)))
-  I_array = reshape(I_matrix, (domain_size..., domain_size...))
+function (::OpName"SWAP")(domain1, domain2)
+  domain = to_dim.((domain1, domain2))
+  I_matrix = Diagonal(trues(prod(domain)))
+  I_array = reshape(I_matrix, (domain..., domain...))
   SWAP_array = permutedims(I_array, (2, 1, 3, 4))
-  SWAP_matrix = reshape(SWAP_array, (prod(domain_size), prod(domain_size)))
+  SWAP_matrix = reshape(SWAP_array, (prod(domain), prod(domain)))
   return SWAP_matrix
 end
 @op_alias "Swap" "SWAP"
@@ -432,8 +452,9 @@ alias(::OpName"√SWAP") = √(OpName"SWAP"())
 
 using LinearAlgebra: diagind
 nsites(::OpName"iSWAP") = 2
-function Base.AbstractArray(::OpName"iSWAP", domain_size::Tuple{Int,Int})
-  swap = AbstractArray(OpName"SWAP"(), domain_size)
+function (::OpName"iSWAP")(domain1, domain2)
+  domain = (domain1, domain2)
+  swap = OpName"SWAP"()(domain...)
   iswap = im * swap
   iswap[diagind(iswap)] .*= -im
   return iswap
@@ -442,35 +463,84 @@ end
 alias(::OpName"√iSWAP") = √(OpName"iSWAP"())
 @op_alias "√iSwap" "√iSWAP"
 
-## TODO: Bring back these definitions.
-## function default_random_matrix(eltype::Type, s::Index...)
-##   n = prod(dim.(s))
-##   return randn(eltype, n, n)
-## end
-##
-## # Haar-random unitary
-## #
-## # Reference:
-## # Section 4.6
-## # http://math.mit.edu/~edelman/publications/random_matrix_theory.pdf
-## function op!(
-##   o::ITensor,
-##   ::OpName"RandomUnitary",
-##   ::SiteType"Generic",
-##   s1::Index,
-##   sn::Index...;
-##   eltype=ComplexF64,
-##   random_matrix=default_random_matrix(eltype, s1, sn...),
-## )
-##   s = (s1, sn...)
-##   Q, _ = NDTensors.qr_positive(random_matrix)
-##   t = itensor(Q, prime.(s)..., dag.(s)...)
-##   return settensor!(o, tensor(t))
-## end
-##
-## function op!(o::ITensor, ::OpName"randU", st::SiteType"Generic", s::Index...; kwargs...)
-##   return op!(o, OpName("RandomUnitary"), st, s...; kwargs...)
-## end
+# TODO: Generalize to more controlled operators, right now
+# there is only one controlled operator that is enabled when
+# all of the last values of each sites/Qudit is set. See:
+# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UCGate
+# https://arxiv.org/abs/2410.05122
+nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.arg)
+function (n::OpName"Controlled")(domain...)
+  # Number of target sites.
+  nt = nsites(n.arg)
+  # 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]))
+  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"()
+function alias(n::OpName"CPhase")
+  return controlled(OpName"Phase"(; params(n)...))
+end
+@op_alias "CPHASE" "CPhase"
+@op_alias "Cphase" "CPhase"
+function alias(n::OpName"CRx")
+  return controlled(OpName"Rx"(; params(n)...))
+end
+@op_alias "CRX" "CRx"
+function alias(::OpName"CRy")
+  return controlled(OpName"Ry"(; params(n)...))
+end
+@op_alias "CRY" "CRy"
+function alias(::OpName"CRz")
+  return controlled(OpName"Rz"(; params(n)...))
+end
+@op_alias "CRZ" "CRz"
+function alias(::OpName"CRn")
+  return controlled(; arg=OpName"Rn"(; params(n)...))
+end
+@op_alias "CRn̂" "CRn"
+
+@op_alias "CCNOT" "Controlled" ncontrol = 2 op = 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 "Fredkin" "CSWAP"
+@op_alias "CSwap" "CSWAP"
+@op_alias "CS" "CSWAP"
+
+@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()
+
+# Version of `sign` that returns one
+# if `x == 0`.
+function nonzero_sign(x)
+  iszero(x) && return one(x)
+  return sign(x)
+end
+
+function qr_positive(M::AbstractMatrix)
+  Q, R = qr(M)
+  Q′ = typeof(R)(Q)
+  signs = nonzero_sign.(diag(R))
+  Q′ = Q′ * Diagonal(signs)
+  R = Diagonal(conj.(signs)) * R
+  return Q′, R
+end
+
+# TODO: Store a random matrix or seed as a parameter
+# of the `OpName`?
+function (n::OpName"RandomUnitary")(domain...)
+  d = prod(to_dim.(domain))
+  elt = get(params(n), :eltype, Complex{Float64})
+  Q, _ = qr_positive(randn(elt, (d, d)))
+  return Q
+end
+@op_alias "randU" "RandomUnitary"
 
 # Expand the operator in a new basis.
 using LinearAlgebra: ⋅
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 2d06f9b..0172f61 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -1,11 +1,12 @@
 struct SiteType{T,Params}
   params::Params
+  function SiteType{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
 end
 params(t::SiteType) = getfield(t, :params)
-
 Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name)
 
-SiteType{N}(params) where {N} = SiteType{N,typeof(params)}(params)
 SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...))
 
 SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...)
@@ -16,8 +17,25 @@ macro SiteType_str(s)
 end
 
 alias(t::SiteType) = t
+alias(i::Integer) = i
+
+function Base.length(t::SiteType)
+  t′ = alias(t)
+  if t == t′
+    return t.length
+  end
+  return length(t′)
+end
 Base.AbstractUnitRange(t::SiteType) = Base.OneTo(length(t))
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]
 Base.axes(t::SiteType) = (AbstractUnitRange(t),)
 Base.axes(t::SiteType, dim::Integer) = axes(t)[dim]
+
+to_dim(d::Base.OneTo) = length(d)
+to_dim(d::SiteType) = length(d)
+to_dim(d::Integer) = d
+
+# TODO: Decide on this.
+# TODO: Move to `sitetype.jl`.
+default_sitetype() = SiteType"Qubit"()
diff --git a/src/sitetypes/electron.jl b/src/sitetypes/electron.jl
index b33519b..0de8438 100644
--- a/src/sitetypes/electron.jl
+++ b/src/sitetypes/electron.jl
@@ -1,303 +1,98 @@
 Base.length(::SiteType"Electron") = 4
 
-Base.AbstractArray(::StateName"Emp", ::Tuple{SiteType"Electron"}) = [1.0, 0, 0, 0]
-Base.AbstractArray(::StateName"Up", ::Tuple{SiteType"Electron"}) = [0.0, 1, 0, 0]
-Base.AbstractArray(::StateName"Dn", ::Tuple{SiteType"Electron"}) = [0.0, 0, 1, 0]
-Base.AbstractArray(::StateName"UpDn", ::Tuple{SiteType"Electron"}) = [0.0, 0, 0, 1]
-# TODO: Use aliasing.
-function Base.AbstractArray(::StateName"0", st::Tuple{SiteType"Electron"})
-  return AbstractArray(StateName("Emp"), st)
-end
-function Base.AbstractArray(::StateName"↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(StateName("Up"), st)
-end
-function Base.AbstractArray(::StateName"↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(StateName("Dn"), st)
-end
-function Base.AbstractArray(::StateName"↑↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(StateName("UpDn"), st)
-end
+(::StateName"Emp")(domain::SiteType"Electron") = StateName"0"()(domain)
 
-# I ⊗ n
-# TODO: Define as `AbstractArray(OpName"I"() ⊗ OpName"n"(), (SiteType("Fermion"), SiteType("Fermion")))`?
-function Base.AbstractArray(::OpName"Nup", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 1.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 1.0
-  ]
-end
-function Base.AbstractArray(on::OpName"n↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+(::StateName"↑")(::SiteType"Electron") = (StateName"0"() ⊗ StateName"1"())(2, 2)
+(::StateName"Up")(domain::SiteType"Electron") = StateName"↑"()(domain)
 
-# n ⊗ I
-function Base.AbstractArray(::OpName"Ndn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 1.0
-  ]
-end
-function Base.AbstractArray(on::OpName"n↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+(::StateName"↓")(domain::SiteType"Electron") = (StateName"1"() ⊗ StateName"0"())(2, 2)
+(::StateName"Dn")(domain::SiteType"Electron") = StateName"↓"()(domain)
 
-# n ⊗ n
-function Base.AbstractArray(::OpName"Nupdn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 1.0
-  ]
-end
-function Base.AbstractArray(on::OpName"n↑↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+(::StateName"↑↓")(domain::SiteType"Electron") = (StateName"1"() ⊗ StateName"1"())(2, 2)
+(::StateName"UpDn")(domain::SiteType"Electron") = StateName"↑↓"()(domain)
 
-# I ⊗ n + n ⊗ I
-function Base.AbstractArray(::OpName"Ntot", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 1.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 2.0
-  ]
-end
-function Base.AbstractArray(on::OpName"ntot", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+# I ⊗ n
+(::OpName"n↑")(::SiteType"Electron") = (OpName"I"() ⊗ OpName"n"())(2, 2)
+@op_alias "Nup" "n↑"
 
-# I ⊗ c
-function Base.AbstractArray(::OpName"Cup", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 1.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 1.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(on::OpName"c↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+# n ⊗ I
+(::OpName"n↓")(::SiteType"Electron") = (OpName"n"() ⊗ OpName"I"())(2, 2)
+@op_alias "Ndn" "n↓"
 
-# I ⊗ c†
-function Base.AbstractArray(::OpName"Cdagup", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    1.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-  ]
-end
-function Base.AbstractArray(on::OpName"c†↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
+# n ⊗ n
+(::OpName"n↑↓")(::SiteType"Electron") = (OpName"n"() ⊗ OpName"n"())(2, 2)
+@op_alias "Nupdn" "n↑↓"
 
-# c ⊗ F
-function Base.AbstractArray(::OpName"Cdn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 -1.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(on::OpName"c↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(alias(on), st)
-end
-
-# c† ⊗ F
-function Base.AbstractArray(::OpName"Cdagdn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    1.0 0.0 0.0 0.0
-    0.0 -1.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"c†↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Cdagdn"), st)
-end
+# I ⊗ n + n ⊗ I = n↑ + n↓
+alias(::OpName"ntot") = OpName"n↑"() + OpName"n↓"()
+@op_alias "Ntot" "ntot"
 
 # I ⊗ a
-function Base.AbstractArray(::OpName"Aup", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 1.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 1.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"a↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Aup"), st)
-end
+(::OpName"a↑")(::SiteType"Electron") = (OpName"I"() ⊗ OpName"a"())(2, 2)
+@op_alias "Aup" "a↑"
+# I ⊗ c
+@op_alias "c↑" "a↑"
+@op_alias "Cup" "c↑"
 
 # I ⊗ a†
-function Base.AbstractArray(::OpName"Adagup", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    1.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"a†↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Adagup"), st)
-end
+(::OpName"a†↑")(::SiteType"Electron") = (OpName"I"() ⊗ OpName"a†"())(2, 2)
+@op_alias "Adagup" "a†↑"
+# I ⊗ c†
+@op_alias "c†↑" "a†↑"
+@op_alias "Cdagup" "c†↑"
 
 # a ⊗ I
-function Base.AbstractArray(::OpName"Adn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 1.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"a↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Adn"), st)
+(::OpName"a↓")(::SiteType"Electron") = (OpName"a"() ⊗ OpName"I"())(2, 2)
+@op_alias "Adn" "a↓"
+# c ⊗ F
+function (::OpName"c↓")(::SiteType"Electron")
+  return (OpName"C"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"())
 end
+@op_alias "Cdn" "c↓"
 
 # a† ⊗ I
-function Base.AbstractArray(::OpName"Adagdn", ::Tuple{SiteType"Electron"})
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    1.0 0.0 0.0 0.0
-    0.0 1.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"a†↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Adagdn"), st)
+(::OpName"a†↓")(::SiteType"Electron") = (OpName"a†"() ⊗ OpName"I"())(2, 2)
+@op_alias "Adagdn" "a†↓"
+# c† ⊗ F
+function (::OpName"c†↓")(::SiteType"Electron")
+  return (OpName"Cdag"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"())
 end
+@op_alias "Cdagdn" "c†↓"
 
 # F ⊗ F
-function Base.AbstractArray(::OpName"F", ::Tuple{SiteType"Electron"})
-  return [
-    1.0 0.0 0.0 0.0
-    0.0 -1.0 0.0 0.0
-    0.0 0.0 -1.0 0.0
-    0.0 0.0 0.0 1.0
-  ]
+function (::OpName"F")(::SiteType"Electron")
+  return (OpName"F"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"())
 end
 
 # I ⊗ F
-function Base.AbstractArray(::OpName"Fup", ::Tuple{SiteType"Electron"})
-  return [
-    1.0 0.0 0.0 0.0
-    0.0 -1.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 -1.0
-  ]
-end
-function Base.AbstractArray(::OpName"F↑", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Fup"), st)
+function (::OpName"F↑")(::SiteType"Electron")
+  return (OpName"I"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"())
 end
+@op_alias "Fup" "F↑"
 
 # F ⊗ I
-function Base.AbstractArray(::OpName"Fdn", ::Tuple{SiteType"Electron"})
-  return [
-    1.0 0.0 0.0 0.0
-    0.0 1.0 0.0 0.0
-    0.0 0.0 -1.0 0.0
-    0.0 0.0 0.0 -1.0
-  ]
-end
-function Base.AbstractArray(::OpName"F↓", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Fdn"), st)
+function (::OpName"F↓")(::SiteType"Electron")
+  return (OpName"F"() ⊗ OpName"I"())(SiteType"Fermion"(), SiteType"Fermion"())
 end
+@op_alias "Fdn" "F↓"
 
-function Base.AbstractArray(::OpName"Sz", ::Tuple{SiteType"Electron"})
-  # cat(falses(1, 1), Matrix(OpName("Sz")), falses(1, 1); dims=(1, 2))
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.5 0.0 0.0
-    0.0 0.0 -0.5 0.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"Sᶻ", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Sz"), st)
-end
-
-function Base.AbstractArray(::OpName"Sx", ::Tuple{SiteType"Electron"})
-  # cat(falses(1, 1), Matrix(OpName("Sx")), falses(1, 1); dims=(1, 2))
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.5 0.0
-    0.0 0.5 0.0 0.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-
-function Base.AbstractArray(::OpName"Sˣ", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("Sx"), st)
-end
-
-function Base.AbstractArray(::OpName"S⁺", ::Tuple{SiteType"Electron"})
-  # cat(falses(1, 1), Matrix(OpName("S⁺")), falses(1, 1); dims=(1, 2))
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 1.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-  ]
-end
-function Base.AbstractArray(::OpName"Sp", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("S⁺"), st)
-end
-function Base.AbstractArray(::OpName"Splus", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("S⁺"), st)
+function (n::OpName"SpinOp")(::SiteType"Electron")
+  return cat(falses(1, 1), n.arg(2), falses(1, 1); dims=(1, 2))
 end
 
-function Base.AbstractArray(::OpName"S⁻", ::Tuple{SiteType"Electron"})
-  # cat(falses(1, 1), Matrix(OpName("S⁻")), falses(1, 1); dims=(1, 2))
-  return [
-    0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-    0.0 1.0 0.0 0.0
-    0.0 0.0 0.0 0.0
-  ]
+# These implicitly define other spin operators.
+# TODO: Maybe require calling it as `OpName("SpinOp"; arg=OpName("Sz"))`?
+function (n::OpName"σ⁺")(domain::SiteType"Electron")
+  return OpName"SpinOp"(; arg=n)(domain)
 end
-function Base.AbstractArray(::OpName"Sm", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("S⁻"), st)
+function (n::OpName"σᶻ")(domain::SiteType"Electron")
+  return OpName"SpinOp"(; arg=n)(domain)
 end
-function Base.AbstractArray(::OpName"Sminus", st::Tuple{SiteType"Electron"})
-  return AbstractArray(OpName("S⁻"), st)
+function (n::OpName"R")(domain::SiteType"Electron")
+  return OpName"SpinOp"(; arg=n)(domain)
 end
 
-@op_alias "a↑" "Aup"
-@op_alias "a↓" "Adn"
-@op_alias "a†↓" "Adagdn"
-@op_alias "a†↑" "Adagup"
-@op_alias "c↑" "Cup"
-@op_alias "c↓" "Cdn"
-@op_alias "c†↑" "Cdagup"
-@op_alias "c†↓" "Cdagdn"
-@op_alias "n↑" "Nup"
-@op_alias "n↓" "Ndn"
-@op_alias "n↑↓" "Nupdn"
-@op_alias "ntot" "Ntot"
-@op_alias "F↑" "Fup"
-@op_alias "F↓" "Fdn"
-
-has_fermion_string(::OpName"Cup", ::Tuple{SiteType"Electron"}) = true
-function has_fermion_string(on::OpName"c↑", st::Tuple{SiteType"Electron"})
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdagup", ::Tuple{SiteType"Electron"}) = true
-function has_fermion_string(on::OpName"c†↑", st::Tuple{SiteType"Electron"})
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdn", ::Tuple{SiteType"Electron"}) = true
-function has_fermion_string(on::OpName"c↓", st::Tuple{SiteType"Electron"})
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdagdn", ::Tuple{SiteType"Electron"}) = true
-function has_fermion_string(on::OpName"c†↓", st::Tuple{SiteType"Electron"})
-  return has_fermion_string(alias(on), st)
-end
+has_fermion_string(::OpName"c↑", ::SiteType"Electron") = true
+has_fermion_string(::OpName"c†↑", ::SiteType"Electron") = true
+has_fermion_string(::OpName"c↓", ::SiteType"Electron") = true
+has_fermion_string(::OpName"c†↓", ::SiteType"Electron") = true
diff --git a/src/sitetypes/fermion.jl b/src/sitetypes/fermion.jl
index 797b584..d5a942b 100644
--- a/src/sitetypes/fermion.jl
+++ b/src/sitetypes/fermion.jl
@@ -1,20 +1,9 @@
-Base.length(::SiteType"Fermion") = 2
+# TODO: Make an alias of `"Qubit"` to inherit
+# more operator and state definitions?
 
-Base.AbstractArray(::StateName"Emp", ::SiteType"Fermion") = [1.0 0.0]
-Base.AbstractArray(::StateName"Occ", ::SiteType"Fermion") = [0.0 1.0]
-function Base.AbstractArray(::StateName"0", st::SiteType"Fermion")
-  return AbstractArray(StateName("Emp"), st)
-end
-function Base.AbstractArray(::StateName"1", st::SiteType"Fermion")
-  return AbstractArray(StateName("Occ"), st)
-end
+Base.length(::SiteType"Fermion") = 2
 
-function Base.AbstractArray(::OpName"F", ::Tuple{SiteType"Fermion"})
-  return [
-    1 0
-    0 -1
-  ]
-end
+(::OpName"F")(::SiteType"Fermion") = Diagonal([1, -1])
 
 @op_alias "c" "a"
 @op_alias "C" "c"
@@ -23,14 +12,5 @@ end
 @op_alias "Cdag" "c†"
 @op_alias "cdag" "c†"
 
-has_fermion_string(::OpName"C", ::SiteType"Fermion") = true
-function has_fermion_string(on::OpName"c", st::SiteType"Fermion")
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdag", ::SiteType"Fermion") = true
-function has_fermion_string(on::OpName"c†", st::SiteType"Fermion")
-  return has_fermion_string(alias(on), st)
-end
-function has_fermion_string(on::OpName"cdag", st::SiteType"Fermion")
-  return has_fermion_string(alias(on), st)
-end
+has_fermion_string(::OpName"c", ::SiteType"Fermion") = true
+has_fermion_string(::OpName"c†", ::SiteType"Fermion") = true
diff --git a/src/sitetypes/qubit.jl b/src/sitetypes/qubit.jl
index fef6d8f..ceccae4 100644
--- a/src/sitetypes/qubit.jl
+++ b/src/sitetypes/qubit.jl
@@ -6,86 +6,71 @@ alias(::SiteType"SpinHalf=1/2") = SiteType"Qubit"()
 
 Base.length(::SiteType"Qubit") = 2
 
-# `eigvecs(Z)`
-Base.AbstractArray(::StateName"0", ::Tuple{SiteType"Qubit"}) = [1, 0]
-Base.AbstractArray(::StateName"1", ::Tuple{SiteType"Qubit"}) = [0, 1]
-@state_alias "Up" "0"
-@state_alias "↑" "0"
-@state_alias "Z+" "0"
-@state_alias "Zp" "0"
-@state_alias "↓" "1"
-@state_alias "Dn" "1"
-@state_alias "Z-" "1"
-@state_alias "Zm" "1"
+(::StateName"↑")(::SiteType"Qubit") = StateName"0"()(2)
+(::StateName"Up")(::SiteType"Qubit") = StateName"0"()(2)
+(::StateName"Z+")(::SiteType"Qubit") = StateName"0"()(2)
+(::StateName"Zp")(::SiteType"Qubit") = StateName"0"()(2)
+(::StateName"Emp")(::SiteType"Qubit") = StateName"0"()(2)
+
+(::StateName"↓")(::SiteType"Qubit") = StateName"1"()(2)
+(::StateName"Dn")(::SiteType"Qubit") = StateName"1"()(2)
+(::StateName"Z-")(::SiteType"Qubit") = StateName"1"()(2)
+(::StateName"Zm")(::SiteType"Qubit") = StateName"1"()(2)
+(::StateName"Occ")(::SiteType"Qubit") = StateName"1"()(2)
 
 # `eigvecs(X)`
-Base.AbstractArray(::StateName"+", ::Tuple{SiteType"Qubit"}) = [1, 1] / √2
-Base.AbstractArray(::StateName"-", ::Tuple{SiteType"Qubit"}) = [1, -1] / √2
+alias(::StateName"+") = (StateName"0"() + StateName"1"()) / √2
 @state_alias "X+" "+"
 @state_alias "Xp" "+"
+
+alias(::StateName"-") = (StateName"0"() - StateName"1"()) / √2
 @state_alias "X-" "-"
 @state_alias "Xm" "-"
 
 # `eigvecs(Y)`
-Base.AbstractArray(::StateName"i", ::Tuple{SiteType"Qubit"}) = [1, im] / √2
-Base.AbstractArray(::StateName"-i", ::Tuple{SiteType"Qubit"}) = [1, -im] / √2
+alias(::StateName"i") = (StateName"0"() + im * StateName"1"()) / √2
 @state_alias "Y+" "i"
 @state_alias "Yp" "i"
+
+alias(::StateName"-i") = (StateName"0"() - im * StateName"1"()) / √2
 @state_alias "Y-" "-i"
 @state_alias "Ym" "-i"
 
 # SIC-POVMs
-Base.AbstractArray(::StateName"Tetra0", ::Tuple{SiteType"Qubit"}) = [
+(::StateName"Tetra0")(::SiteType"Qubit") = [
   1
   0
 ]
-Base.AbstractArray(::StateName"Tetra2", ::Tuple{SiteType"Qubit"}) = [
+(::StateName"Tetra2")(::SiteType"Qubit") = [
   1 / √3
   √2 / √3
 ]
-function Base.AbstractArray(::StateName"Tetra3", ::Tuple{SiteType"Qubit"})
-  return [
-    1 / √3
-    √2 / √3 * exp(im * 2π / 3)
-  ]
-end
-function Base.AbstractArray(::StateName"Tetra4", ::Tuple{SiteType"Qubit"})
-  return [
-    1 / √3
-    √2 / √3 * exp(im * 4π / 3)
-  ]
-end
-
-# TODO: Write as `(I + σᶻ) / 2`?
-Base.AbstractArray(::OpName"ProjUp", ::Tuple{SiteType"Qubit"}) = [
-  1 0
-  0 0
+(::StateName"Tetra3")(::SiteType"Qubit") = [
+  1 / √3
+  √2 / √3 * exp(im * 2π / 3)
 ]
+(::StateName"Tetra4")(::SiteType"Qubit") = [
+  1 / √3
+  √2 / √3 * exp(im * 4π / 3)
+]
+
+# TODO: Define as `(I + σᶻ) / 2`?
+alias(::OpName"ProjUp") = OpName"Proj"(; index=1)
 @op_alias "projUp" "ProjUp"
 @op_alias "Proj↑" "ProjUp"
 @op_alias "proj↑" "ProjUp"
 @op_alias "Proj0" "ProjUp"
 @op_alias "proj0" "ProjUp"
 
-# TODO: Define as `σ⁺ * σ−`?
-# TODO: Write as `(I - σᶻ) / 2`?
-Base.AbstractArray(::OpName"ProjDn", ::Tuple{SiteType"Qubit"}) = [
-  0 0
-  0 1
-]
+# TODO: Define as `σ⁺ * σ⁻`?
+# TODO: Define as `(I - σᶻ) / 2`?
+alias(::OpName"ProjDn") = OpName"Proj"(; index=2)
 @op_alias "projDn" "ProjDn"
 @op_alias "Proj↓" "ProjDn"
 @op_alias "proj↓" "ProjDn"
 @op_alias "Proj1" "ProjDn"
 @op_alias "proj1" "ProjDn"
 
-# TODO: Determine a general spin definition, such as
-# `eigvecs(X)`.
-Base.AbstractArray(::OpName"H", ::Tuple{SiteType"Qubit"}) = [
-  1/√2 1/√2
-  1/√2 -1/√2
-]
-
 # Rotation around generic axis n̂
 # exp(-im * n.θ / 2 * n̂ ⋅ σ⃗)
 #=
@@ -103,62 +88,15 @@ or:
 alias(n::OpName"R") = OpName"Rn"(; θ=n.θ, ϕ=n.ϕ, λ=-n.ϕ)
 =#
 # https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UGate
-function Base.AbstractArray(n::OpName"Rn", ::Tuple{SiteType"Qubit"})
+# TODO: Generalize to `"Qudit"`, see:
+# https://quantumcomputing.stackexchange.com/questions/16251/how-does-a-general-rotation-r-hatn-theta-related-to-u-3-gate
+# https://quantumcomputing.stackexchange.com/questions/4249/decomposition-of-an-arbitrary-1-qubit-gate-into-a-specific-gateset
+# https://en.wikipedia.org/wiki/List_of_quantum_logic_gates#Other_named_gates
+# https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins
+function (n::OpName"Rn")(::SiteType"Qubit")
   return [
     cos(n.θ / 2) -exp(im * n.λ)*sin(n.θ / 2)
     exp(im * n.ϕ)*sin(n.θ / 2) exp(im * (n.ϕ + n.λ))*cos(n.θ / 2)
   ]
 end
 @op_alias "Rn̂" "Rn"
-
-# TODO: Generalize to `"Qudit"` and other SiteTypes.
-# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UCGate
-nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.op)
-function Base.AbstractArray(n::OpName"Controlled", ts::Tuple{Vararg{SiteType"Qubit"}})
-  # Number of target qubits.
-  nt = nsites(n.op)
-  # Number of control qubits.
-  nc = get(params(n), :ncontrol, length(ts) - nt)
-  @assert length(ts) == nc + nt
-  return [
-    I(2^nc) falses(2^nc, 2^nt)
-    falses(2^nt, 2^nc) AbstractArray(n.op, ts[(nc + 1):end])
-  ]
-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"()
-function alias(n::OpName"CPhase")
-  return controlled(OpName"Phase"(; params(n)...))
-end
-@op_alias "CPHASE" "CPhase"
-@op_alias "Cphase" "CPhase"
-function alias(n::OpName"CRx")
-  return controlled(OpName"Rx"(; params(n)...))
-end
-@op_alias "CRX" "CRx"
-function Base.AbstractArray(::OpName"CRy")
-  return controlled(OpName"Ry"(; params(n)...))
-end
-@op_alias "CRY" "CRy"
-function Base.AbstractArray(::OpName"CRz")
-  return controlled(OpName"Rz"(; params(n)...))
-end
-@op_alias "CRZ" "CRz"
-function Base.AbstractArray(::OpName"CRn")
-  return controlled(; op=OpName"Rn"(; params(n)...))
-end
-@op_alias "CRn̂" "CRn"
-
-@op_alias "CCNOT" "Controlled" ncontrol = 2 op = 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 "Fredkin" "CSWAP"
-@op_alias "CSwap" "CSWAP"
-@op_alias "CS" "CSWAP"
-
-@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()
diff --git a/src/sitetypes/qudit.jl b/src/sitetypes/qudit.jl
index ad8463d..7205c98 100644
--- a/src/sitetypes/qudit.jl
+++ b/src/sitetypes/qudit.jl
@@ -1,2 +1,3 @@
-Base.length(t::SiteType"Qudit") = t.length
-alias(::SiteType"Boson") = SiteType"Qudit"()
+Base.length(t::SiteType"Qudit") = t.dim
+alias(t::SiteType"Boson") = SiteType"Qudit"(; dim=t.dim)
+alias(t::SiteType"S") = SiteType"Qudit"(; dim=Int(2t.spin + 1))
diff --git a/src/sitetypes/spinone.jl b/src/sitetypes/spinone.jl
index 87d612c..84f12e8 100644
--- a/src/sitetypes/spinone.jl
+++ b/src/sitetypes/spinone.jl
@@ -2,27 +2,27 @@ Base.length(::SiteType"S=1") = 3
 
 alias(::SiteType"SpinOne") = SiteType"S=1"()
 
-# TODO: Decide on these names, use `alias`.
-# TODO: Make a more general `SiteType"Spin`/`SiteType"S"`
-# with a `spin` field that can be set to a rational number
-# `1//2, `2//2`, `3//2`, etc. that maps to `Qudit`
-# of length `2 * spin + 1`.
-Base.AbstractArray(::StateName"Up", ::SiteType"S=1") = [1, 0, 0]
-Base.AbstractArray(::StateName"Z0", ::SiteType"S=1") = [0, 1, 0]
-Base.AbstractArray(::StateName"Dn", ::SiteType"S=1") = [0, 0, 1]
+# TODO: Make these more general, define as something like:
+# `(n::StateName"Z")(::SiteType"S=1") = eigvecs(OpName"Z"())[n.eigval]`
+(::StateName"Z+")(::SiteType"S=1") = [1, 0, 0]
+(::StateName"Z0")(::SiteType"S=1") = [0, 1, 0]
+(::StateName"Z-")(::SiteType"S=1") = [0, 0, 1]
 
-Base.AbstractArray(::StateName"↑", st::SiteType"S=1") = [1, 0, 0]
-Base.AbstractArray(::StateName"0", st::SiteType"S=1") = [0, 1, 0]
-Base.AbstractArray(::StateName"↓", st::SiteType"S=1") = [0, 0, 1]
+## TODO: Decide on these aliases.
+(::StateName"↑")(::SiteType"S=1") = StateName"Z+"()(domain)
+(::StateName"Up")(::SiteType"S=1") = StateName"Z+"()(domain)
+(::StateName"0")(::SiteType"S=1") = StateName"Z0"()(domain)
+(::StateName"↓")(::SiteType"S=1") = StateName"Z-"()(domain)
+(::StateName"Dn")(::SiteType"S=1") = StateName"Z-"()(domain)
 
-Base.AbstractArray(::StateName"Z+", st::SiteType"S=1") = [1.0, 0.0, 0.0]
-# -- Z0 is already defined above --
-Base.AbstractArray(::StateName"Z-", st::SiteType"S=1") = [0.0, 0.0, 1.0]
+# TODO: Make these more general, define as something like:
+# `(n::StateName"X")(::SiteType"S=1") = eigvecs(OpName"X"())[n.eigval]`
+(::StateName"X+")(::SiteType"S=1") = [1 / 2, 1 / sqrt(2), 1 / 2]
+(::StateName"X0")(::SiteType"S=1") = [-1 / sqrt(2), 0, 1 / sqrt(2)]
+(::StateName"X-")(::SiteType"S=1") = [1 / 2, -1 / sqrt(2), 1 / 2]
 
-Base.AbstractArray(::StateName"X+", ::SiteType"S=1") = [1 / 2, 1 / sqrt(2), 1 / 2]
-Base.AbstractArray(::StateName"X0", ::SiteType"S=1") = [-1 / sqrt(2), 0, 1 / sqrt(2)]
-Base.AbstractArray(::StateName"X-", ::SiteType"S=1") = [1 / 2, -1 / sqrt(2), 1 / 2]
-
-Base.AbstractArray(::StateName"Y+", ::SiteType"S=1") = [-1 / 2, -im / sqrt(2), 1 / 2]
-Base.AbstractArray(::StateName"Y0", ::SiteType"S=1") = [1 / sqrt(2), 0, 1 / sqrt(2)]
-Base.AbstractArray(::StateName"Y-", ::SiteType"S=1") = [-1 / 2, im / sqrt(2), 1 / 2]
+# TODO: Make these more general, define as something like:
+# `(n::StateName"Y")(::SiteType"S=1") = eigvecs(OpName"Y"())[n.eigval]`
+(::StateName"Y+")(::SiteType"S=1") = [-1 / 2, -im / sqrt(2), 1 / 2]
+(::StateName"Y0")(::SiteType"S=1") = [1 / sqrt(2), 0, 1 / sqrt(2)]
+(::StateName"Y-")(::SiteType"S=1") = [-1 / 2, im / sqrt(2), 1 / 2]
diff --git a/src/sitetypes/tj.jl b/src/sitetypes/tj.jl
index 92bebce..4efa7ef 100644
--- a/src/sitetypes/tj.jl
+++ b/src/sitetypes/tj.jl
@@ -1,174 +1,14 @@
 Base.length(::SiteType"tJ") = 3
 
-Base.AbstractArray(::StateName"Emp", ::SiteType"tJ") = [1.0, 0, 0]
-Base.AbstractArray(::StateName"Up", ::SiteType"tJ") = [0.0, 1, 0]
-Base.AbstractArray(::StateName"Dn", ::SiteType"tJ") = [0.0, 0, 1]
-Base.AbstractArray(::StateName"0", st::SiteType"tJ") = AbstractArray(StateName("Emp"), st)
-Base.AbstractArray(::StateName"↑", st::SiteType"tJ") = AbstractArray(StateName("Up"), st)
-Base.AbstractArray(::StateName"↓", st::SiteType"tJ") = AbstractArray(StateName("Dn"), st)
-
-# TODO: Update these to the latest syntax.
-## function op!(Op::ITensor, ::OpName"Nup", ::SiteType"tJ", s::Index)
-##   return Op[s' => 2, s => 2] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"n↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Ndn", ::SiteType"tJ", s::Index)
-##   return Op[s' => 3, s => 3] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"n↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Ntot", ::SiteType"tJ", s::Index)
-##   Op[s' => 2, s => 2] = 1.0
-##   return Op[s' => 3, s => 3] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"ntot", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Cup", ::SiteType"tJ", s::Index)
-##   return Op[s' => 1, s => 2] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"c↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Cdagup", ::SiteType"tJ", s::Index)
-##   return Op[s' => 2, s => 1] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"c†↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Cdn", ::SiteType"tJ", s::Index)
-##   return Op[s' => 1, s => 3] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"c↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Cdagdn", ::SiteType"tJ", s::Index)
-##   return Op[s' => 3, s => 1] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"c†↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Aup", ::SiteType"tJ", s::Index)
-##   return Op[s' => 1, s => 2] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"a↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Adagup", ::SiteType"tJ", s::Index)
-##   return Op[s' => 2, s => 1] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"a†↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Adn", ::SiteType"tJ", s::Index)
-##   return Op[s' => 1, s => 3] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"a↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Adagdn", ::SiteType"tJ", s::Index)
-##   return Op[s' => 3, s => 1] = 1.0
-## end
-## function op!(Op::ITensor, on::OpName"a†↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"F", ::SiteType"tJ", s::Index)
-##   Op[s' => 1, s => 1] = +1.0
-##   Op[s' => 2, s => 2] = -1.0
-##   return Op[s' => 3, s => 3] = -1.0
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Fup", ::SiteType"tJ", s::Index)
-##   Op[s' => 1, s => 1] = +1.0
-##   Op[s' => 2, s => 2] = -1.0
-##   return Op[s' => 3, s => 3] = +1.0
-## end
-## function op!(Op::ITensor, on::OpName"F↑", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Fdn", ::SiteType"tJ", s::Index)
-##   Op[s' => 1, s => 1] = +1.0
-##   Op[s' => 2, s => 2] = +1.0
-##   return Op[s' => 3, s => 3] = -1.0
-## end
-## function op!(Op::ITensor, on::OpName"F↓", st::SiteType"tJ", s::Index)
-##   return op!(Op, alias(on), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Sz", ::SiteType"tJ", s::Index)
-##   Op[s' => 2, s => 2] = +0.5
-##   return Op[s' => 3, s => 3] = -0.5
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Sᶻ", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("Sz"), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Sx", ::SiteType"tJ", s::Index)
-##   Op[s' => 2, s => 3] = 0.5
-##   return Op[s' => 3, s => 2] = 0.5
-## end
-## 
-## function op!(Op::ITensor, ::OpName"Sˣ", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("Sx"), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"S+", ::SiteType"tJ", s::Index)
-##   return Op[s' => 2, s => 3] = 1.0
-## end
-## 
-## function op!(Op::ITensor, ::OpName"S⁺", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S+"), st, s)
-## end
-## function op!(Op::ITensor, ::OpName"Sp", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S+"), st, s)
-## end
-## function op!(Op::ITensor, ::OpName"Splus", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S+"), st, s)
-## end
-## 
-## function op!(Op::ITensor, ::OpName"S-", ::SiteType"tJ", s::Index)
-##   return Op[s' => 3, s => 2] = 1.0
-## end
-## 
-## function op!(Op::ITensor, ::OpName"S⁻", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S-"), st, s)
-## end
-## function op!(Op::ITensor, ::OpName"Sm", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S-"), st, s)
-## end
-## function op!(Op::ITensor, ::OpName"Sminus", st::SiteType"tJ", s::Index)
-##   return op!(Op, OpName("S-"), st, s)
-## end
-
-has_fermion_string(::OpName"Cup", ::SiteType"tJ") = true
-function has_fermion_string(on::OpName"c↑", st::SiteType"tJ")
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdagup", ::SiteType"tJ") = true
-function has_fermion_string(on::OpName"c†↑", st::SiteType"tJ")
-  return has_fermion_string(alias(on), st)
+function (n::StateName)(domain::SiteType"tJ")
+  return n(SiteType"Electron"())[1:length(domain)]
 end
-has_fermion_string(::OpName"Cdn", ::SiteType"tJ") = true
-function has_fermion_string(on::OpName"c↓", st::SiteType"tJ")
-  return has_fermion_string(alias(on), st)
-end
-has_fermion_string(::OpName"Cdagdn", ::SiteType"tJ") = true
-function has_fermion_string(on::OpName"c†↓", st::SiteType"tJ")
-  return has_fermion_string(alias(on), st)
+
+function (n::OpName)(domain::SiteType"tJ")
+  return n(SiteType"Electron"())[1:length(domain), 1:length(domain)]
 end
+
+has_fermion_string(::OpName"c↑", ::SiteType"tJ") = true
+has_fermion_string(::OpName"c†↑", ::SiteType"tJ") = true
+has_fermion_string(::OpName"c↓", ::SiteType"tJ") = true
+has_fermion_string(::OpName"c†↓", ::SiteType"tJ") = true
diff --git a/src/state.jl b/src/state.jl
index d0cb08b..969ead4 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -1,11 +1,14 @@
+using Random: randstring
+
 struct StateName{Name,Params}
   params::Params
+  function StateName{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
 end
 params(n::StateName) = getfield(n, :params)
-
 Base.getproperty(n::StateName, name::Symbol) = getfield(params(n), name)
 
-StateName{N}(params) where {N} = StateName{N,typeof(params)}(params)
 StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...))
 
 StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...)
@@ -15,6 +18,10 @@ macro StateName_str(s)
   return StateName{Symbol(s)}
 end
 
+# Handles special case `state(1) == [1, 0]`. Note the
+# one-based indexing, as opposed to `state("0") == [1, 0]`.
+StateName(i::Integer; kwargs...) = StateName"StandardBasis"(; index=i, kwargs...)
+
 function state_alias_expr(name1, name2, pars...)
   return :(function alias(n::StateName{Symbol($name1)})
     return StateName{Symbol($name2)}(; params(n)..., $(esc.(pars)...))
@@ -24,68 +31,121 @@ macro state_alias(name1, name2, params...)
   return state_alias_expr(name1, name2)
 end
 
-# TODO: Decide on this.
-default_sitetype() = SiteType"Qubit"()
-
-alias(n::StateName) = n
-function (arrtype::Type{<:AbstractArray})(n::StateName, ts::Tuple{Vararg{SiteType}})
+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, AbstractArray(n, ts))
+  return convert(arrtype, n(domain...))
 end
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain_size::Tuple{Vararg{Integer}})
+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, AbstractArray(n, Int.(domain_size)))
+  return convert(arrtype, n(Int.(domain)...))
 end
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain_size::Integer...)
-  return arrtype(n, domain_size)
+
+# This compiles operator expressions, such as:
+# ```julia
+# stateexpr("0 + 1") == StateName("0") + StateName("1")
+# ```
+function stateexpr(n::String; kwargs...)
+  return state_or_op_expr(StateName, n; kwargs...)
 end
-(arrtype::Type{<:AbstractArray})(n::StateName, ts::SiteType...) = arrtype(n, ts)
-function Base.AbstractArray(n::StateName, ts::Tuple{Vararg{SiteType}})
-  n′ = alias(n)
-  ts′ = alias.(ts)
-  if n′ == n && ts′ == ts
-    return AbstractArray(n′, length.(ts′))
-  end
-  return AbstractArray(n′, ts′)
+
+# Handles special case `state(1) == [1, 0]`. Note the
+# one-based indexing, as opposed to `state("0") == [1, 0]`.
+stateexpr(n::Integer; kwargs...) = StateName(n; kwargs...)
+
+randcharstring() = randstring(['A':'Z'; 'a':'z'], 12)
+const DAGGER_STRING = randcharstring()
+const UPARROW_STRING = randcharstring()
+const DOWNARROW_STRING = randcharstring()
+const PLUS_STRING = randcharstring()
+const MINUS_STRING = randcharstring()
+const VERTICALBAR_STRING = randcharstring()
+const RANGLE_STRING = randcharstring()
+const EXPR_REPLACEMENTS_1 = (
+  "†" => DAGGER_STRING,
+  "↑" => UPARROW_STRING,
+  "↓" => DOWNARROW_STRING,
+  # Replace trailing plus and minus characters
+  # in operators, which don't parse properly.
+  r"(\S)\+" => SubstitutionString("\\1$(PLUS_STRING)"),
+  r"(\S)\-" => SubstitutionString("\\1$(MINUS_STRING)"),
+)
+const EXPR_REPLACEMENTS_2 = (
+  r"\|(\S*)⟩" => SubstitutionString("$(VERTICALBAR_STRING)\\1$(RANGLE_STRING)"),
+)
+const INVERSE_EXPR_REPLACEMENTS = (
+  DAGGER_STRING => "†",
+  UPARROW_STRING => "↑",
+  DOWNARROW_STRING => "↓",
+  PLUS_STRING => "+",
+  MINUS_STRING => "-",
+  # We remove the bra-ket notation and search for states
+  # with names stored inside.
+  VERTICALBAR_STRING => "",
+  RANGLE_STRING => "",
+)
+
+function state_or_op_expr(ntype::Type, n::String; kwargs...)
+  # Do this in two rounds since for some reason
+  # one round doesn't work for expressions
+  # like `"|+⟩"`.
+  n = replace(n, EXPR_REPLACEMENTS_1...)
+  n = replace(n, EXPR_REPLACEMENTS_2...)
+  depth = 1
+  return state_or_op_expr(ntype, Meta.parse(n), depth; kwargs...)
 end
-function Base.AbstractArray(n::StateName, domain_size::Tuple{Vararg{Int}})
-  n′ = alias(n)
-  if n′ == n
-    error("Not implemented.")
+function state_or_op_expr(ntype::Type, n::Number, depth::Int; kwargs...)
+  if depth == 1
+    return ntype{Symbol(n)}(; kwargs...)
   end
-  return AbstractArray(n′, domain_size)
+  return n
 end
-
-# TODO: Decide on this.
-function Base.AbstractArray(n::StateName)
-  return AbstractArray(n, ntuple(Returns(default_sitetype()), nsites(n)))
+function state_or_op_expr(ntype::Type, n::Symbol, depth::Int; kwargs...)
+  n === :im && return im
+  n === :π && return π
+  n = Symbol(replace(String(n), INVERSE_EXPR_REPLACEMENTS...))
+  return ntype{n}(; kwargs...)
 end
-function (arrtype::Type{<:AbstractArray})(n::StateName)
-  return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n)))
+function state_or_op_expr(ntype::Type, ex::Expr, depth::Int)
+  if Meta.isexpr(ex, :call)
+    return eval(ex.args[1])(state_or_op_expr.(ntype, ex.args[2:end], depth + 1)...)
+  end
+  if Meta.isexpr(ex, :curly)
+    # Syntax for parametrized gates, i.e.
+    # `state_or_op_expr("Ry{θ=π/2}")`.
+    params = ex.args[2:end]
+    kwargs = Dict(
+      map(params) do param
+        @assert Meta.isexpr(param, :(=))
+        return param.args[1] => eval(param.args[2])
+      end,
+    )
+    return ntype{ex.args[1]}(; kwargs...)
+  end
+  return error("Can't parse expression $ex.")
 end
 
-function state(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
-  return arrtype(StateName(n; kwargs...), domain...)
+function state(arrtype::Type{<:AbstractArray}, n::Union{Int,String}, domain...; kwargs...)
+  return arrtype(stateexpr(n; kwargs...), domain...)
 end
-function state(elt::Type{<:Number}, n::String, domain...; kwargs...)
+function state(elt::Type{<:Number}, n::Union{Int,String}, domain...; kwargs...)
   return state(AbstractArray{elt}, n, domain...; kwargs...)
 end
-function state(n::String, domain...; kwargs...)
+function state(n::Union{Int,String}, domain...; kwargs...)
   return state(AbstractArray, n, domain...; kwargs...)
 end
 
-# TODO: Add this.
-## function Base.Integer(::StateName{N}) where {N}
-##   return parse(Int, String(N))
-## end
-
-# TODO: Define as `::Tuple{Int}`.
-function Base.AbstractArray(n::StateName{N}, domain_size::Tuple{Int}) where {N}
-  # TODO: Use `Integer(n)`.
-  n = parse(Int, String(N))
-  a = falses(only(domain_size))
-  a[n + 1] = one(Bool)
+function (n::StateName"StandardBasis")(domain)
+  a = falses(to_dim(domain))
+  a[n.index] = one(Bool)
   return a
 end
+function (n::StateName{N})(domain) where {N}
+  n′ = alias(n)
+  if n == n′
+    index = parse(Int, String(N)) + 1
+    return StateName"StandardBasis"(; index)(domain)
+  end
+  return n′(domain)
+end
diff --git a/test/test_basics.jl b/test/test_basics.jl
index b2160bf..b1d5eeb 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -1,5 +1,5 @@
 using QuantumOperatorDefinitions: OpName, SiteType, ⊗, expand, op, opexpr, state, nsites
-using LinearAlgebra: Diagonal
+using LinearAlgebra: Diagonal, I
 using Test: @test, @testset
 
 const real_elts = (Float32, Float64)
@@ -10,7 +10,22 @@ const elts = (real_elts..., complex_elts...)
   @testset "Qubit/Qudit" begin
     # https://en.wikipedia.org/wiki/Pauli_matrices
     # https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins
-    for (t, len, Xmat, Ymat, Zmat, Nmat, SWAPmat, iSWAPmat, RXXmat, RYYmat, RZZmat) in (
+    for (
+      t,
+      len,
+      Xmat,
+      Ymat,
+      Zmat,
+      Nmat,
+      SWAPmat,
+      iSWAPmat,
+      RXXmat,
+      RYYmat,
+      RZZmat,
+      Proj1mat,
+      Proj2mat,
+      StandardBasis12mat,
+    ) in (
       (
         SiteType("Qubit"),
         2,
@@ -34,9 +49,12 @@ const elts = (real_elts..., complex_elts...)
         ],
         (_, θ) ->
           Diagonal([exp(-im * θ / 2), exp(im * θ / 2), exp(im * θ / 2), exp(-im * θ / 2)]),
+        [1 0; 0 0],
+        [0 0; 0 1],
+        [0 1; 0 0],
       ),
       (
-        SiteType("Qudit"; length=3),
+        SiteType("Qudit"; dim=3),
         3,
         √2 * [0 1 0; 1 0 1; 0 1 0],
         √2 * [0 -im 0; im 0 -im; 0 im 0],
@@ -67,6 +85,9 @@ const elts = (real_elts..., complex_elts...)
         (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
         (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
         (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
+        [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],
       ),
     )
       @test length(t) == len
@@ -102,6 +123,9 @@ const elts = (real_elts..., complex_elts...)
         (OpName("RYY"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)),
         (OpName("Rzz"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)),
         (OpName("RZZ"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)),
+        (OpName("Proj"; index=1), 1, elts, Proj1mat),
+        (OpName("Proj"; index=2), 1, elts, Proj2mat),
+        (OpName("StandardBasis"; index=(1, 2)), 1, elts, StandardBasis12mat),
       )
         @test nsites(o) == nbits
         for arraytype in (AbstractArray, AbstractMatrix, Array, Matrix)
@@ -118,7 +142,7 @@ const elts = (real_elts..., complex_elts...)
       end
     end
   end
-  @testset "Parsing" begin
+  @testset "op parsing" begin
     @test Matrix(opexpr("X * Y")) == op("X") * op("Y")
     @test op("X * Y") == op("X") * op("Y")
     @test op("X * Y + Z") == op("X") * op("Y") + op("Z")
@@ -127,5 +151,81 @@ const elts = (real_elts..., complex_elts...)
     @test op("exp(im * (X ⊗ Y + Z ⊗ Z))") ==
       exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z"))))
     @test op("Ry{θ=π/2}") == op("Ry"; θ=π / 2)
+    # Awkward parsing corner cases.
+    @test op("S+") == Matrix(OpName("S+"))
+    @test op("S-") == Matrix(OpName("S-"))
+    @test op("S+ + S-") == Matrix(OpName("S+") + OpName("S-"))
+    @test op("S+ - S-") == Matrix(OpName("S+") - OpName("S-"))
+    @test op("a†") == Matrix(OpName("a†"))
+    for name in ("c↑", "c†↑", "c↓", "c†↓")
+      @test op(name, SiteType("Electron")) == Matrix(OpName(name), SiteType("Electron"))
+    end
+  end
+  @testset "Random unitary" begin
+    U = op("RandomUnitary")
+    @test eltype(U) === Complex{Float64}
+    @test U * U' ≈ I(2)
+
+    U = op("RandomUnitary"; eltype=Float32)
+    @test eltype(U) === Float32
+    @test U * U' ≈ I(2)
+
+    U = op("RandomUnitary", 3)
+    @test eltype(U) === Complex{Float64}
+    @test U * U' ≈ I(3)
+  end
+  @testset "state" begin
+    @test state(1) == [1, 0]
+    @test state("0") == [1, 0]
+    @test state(2) == [0, 1]
+    @test state("1") == [0, 1]
+    @test state(1, 3) == [1, 0, 0]
+    @test state("0", 3) == [1, 0, 0]
+    @test state(2, 3) == [0, 1, 0]
+    @test state("1", 3) == [0, 1, 0]
+    @test state(3, 3) == [0, 0, 1]
+    @test state("2", 3) == [0, 0, 1]
+
+    @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+")
+    @test state("|0⟩ ⊗ |+⟩") == kron(state("0"), state("+"))
+  end
+  @testset "Electron/tJ" begin
+    for (ns, x) in (
+      (("0", "Emp"), [1, 0, 0, 0]),
+      (("↑", "Up"), [0, 1, 0, 0]),
+      (("↓", "Dn"), [0, 0, 1, 0]),
+      (("↑↓", "UpDn"), [0, 0, 0, 1]),
+    )
+      for n in ns
+        @test state(n, SiteType("Electron")) == x
+        @test state(n, SiteType("tJ")) == x[1:3]
+      end
+    end
+    for (ns, x) in (
+      (("n↑", "Nup"), [0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]),
+      (("n↓", "Ndn"), [0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]),
+      (("n↑↓", "Nupdn"), [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1]),
+      (("ntot", "Ntot"), [0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2]),
+      (("c↑", "Cup"), [0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]),
+      (("Sx",), [0 0 0 0; 0 0 1/2 0; 0 1/2 0 0; 0 0 0 0]),
+      # TODO: Add more tests.
+    )
+      for n in ns
+        @test op(n, SiteType("Electron")) == x
+        @test op(n, SiteType("tJ")) == x[1:3, 1:3]
+      end
+    end
+  end
+  @testset "Boson and spin" begin
+    for t in (SiteType("Qudit"; dim=3), SiteType("Boson"; dim=3), SiteType("S"; spin=1))
+      @test length(t) == 3
+      @test op("X", t) == op("X", 3)
+    end
+
+    for t in (SiteType("S=1"), SiteType("SpinOne"))
+      @test state("Z+", SiteType("S=1")) == [1, 0, 0]
+      @test state("Z0", SiteType("S=1")) == [0, 1, 0]
+      @test state("Z-", SiteType("S=1")) == [0, 0, 1]
+    end
   end
 end
diff --git a/test/test_itensorbaseext.jl b/test/test_itensorbaseext.jl
new file mode 100644
index 0000000..e71fe46
--- /dev/null
+++ b/test/test_itensorbaseext.jl
@@ -0,0 +1,13 @@
+using ITensorBase: ITensor, Index, prime
+using QuantumOperatorDefinitions: op, state
+using Test: @test, @testset
+
+@testset "QuantumOperatorDefinitionsITensorBaseExt" begin
+  i = Index(2)
+
+  a = op("X", i)
+  @test a == ITensor([0 1; 1 0], (prime(i), i))
+
+  a = state(1, i)
+  @test a == ITensor([1, 0], (i,))
+end