diff --git a/.typos.toml b/.typos.toml
index 7364886..3e4e58a 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -1,2 +1,3 @@
 [default.extend-words]
-ket = "ket"
\ No newline at end of file
+ket = "ket"
+BA = "BA"
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f2ad3a6..3a12857 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # News
 
+## v0.3.3 - 2024-07-02
+
+- Introduced `qexpand` function that manually expands expressions containing quantum objects.
+- Organized automatic scaling and flattening procedures.
+- Added `express.md` to docs.
+
 ## v0.3.2 - 2024-06-28
 
 - Added documentation for `express`.
diff --git a/Project.toml b/Project.toml
index 2ed9266..58d45f7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumSymbolics"
 uuid = "efa7fd63-0460-4890-beb7-be1bbdfbaeae"
 authors = ["QuantumSymbolics.jl contributors"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
diff --git a/docs/make.jl b/docs/make.jl
index 3afb5da..4092076 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -26,6 +26,7 @@ function main()
     pages = [
         "QuantumSymbolics.jl" => "index.md",
         "Qubit Basis Choice" => "qubit_basis.md",
+        "Express Functionality" => "express.md",
         "API" => "API.md",
     ]
     )
diff --git a/src/QSymbolicsBase/QSymbolicsBase.jl b/src/QSymbolicsBase/QSymbolicsBase.jl
index bb913bc..796cc75 100644
--- a/src/QSymbolicsBase/QSymbolicsBase.jl
+++ b/src/QSymbolicsBase/QSymbolicsBase.jl
@@ -1,20 +1,20 @@
 using Symbolics
 import Symbolics: simplify
 using SymbolicUtils
-import SymbolicUtils: Symbolic, _isone, flatten_term, isnotflat, Chain, Fixpoint, Prewalk
+import SymbolicUtils: Symbolic,_isone,flatten_term,isnotflat,Chain,Fixpoint,Prewalk
 using TermInterface
-import TermInterface: isexpr, head, iscall, children, operation, arguments, metadata, maketerm
+import TermInterface: isexpr,head,iscall,children,operation,arguments,metadata,maketerm
 
 using LinearAlgebra
-import LinearAlgebra: eigvecs, ishermitian, inv
+import LinearAlgebra: eigvecs,ishermitian,inv
 
 import QuantumInterface:
     apply!,
     tensor, ⊗,
-    basis, Basis, SpinBasis, FockBasis,
+    basis,Basis,SpinBasis,FockBasis,
     nqubits,
-    projector, dagger,
-    AbstractKet, AbstractOperator, AbstractSuperOperator, AbstractBra
+    projector,dagger,
+    AbstractBra,AbstractKet,AbstractOperator,AbstractSuperOperator
 
 export SymQObj,QObj,
        AbstractRepresentation,AbstractUse,
@@ -23,7 +23,7 @@ export SymQObj,QObj,
        apply!,
        express,
        tensor,⊗,
-       dagger,projector,commutator,anticommutator,expand,
+       dagger,projector,commutator,anticommutator,
        I,X,Y,Z,σˣ,σʸ,σᶻ,Pm,Pp,σ₋,σ₊,
        H,CNOT,CPHASE,XCX,XCY,XCZ,YCX,YCY,YCZ,ZCX,ZCY,ZCZ,
        X1,X2,Y1,Y2,Z1,Z2,X₁,X₂,Y₁,Y₂,Z₁,Z₂,L0,L1,Lp,Lm,Lpi,Lmi,L₀,L₁,L₊,L₋,L₊ᵢ,L₋ᵢ,
@@ -35,36 +35,16 @@ export SymQObj,QObj,
        SScaled,SScaledBra,SScaledOperator,SScaledKet,
        STensorBra,STensorKet,STensorOperator,
        SZeroBra,SZeroKet,SZeroOperator,
-       SProjector,MixedState,IdentityOp,SInvOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,
+       SProjector,MixedState,IdentityOp,SInvOperator,
        SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SDagger,SBraKet,SOuterKetBra,
        HGate,XGate,YGate,ZGate,CPHASEGate,CNOTGate,
        XBasisState,YBasisState,ZBasisState,
        NumberOp,CreateOp,DestroyOp,
        XCXGate,XCYGate,XCZGate,YCXGate,YCYGate,YCZGate,ZCXGate,ZCYGate,ZCZGate,
-       qsimplify,qsimplify_pauli,qsimplify_flatten,qsimplify_commutator,qsimplify_anticommutator,
+       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,
+       qexpand,
        isunitary
 
-function countmap(samples) # A simpler version of StatsBase.countmap, because StatsBase is slow to import
-    counts = Dict{Any,Any}()
-    for s in samples
-        counts[s] = get(counts, s, 0)+1
-    end
-    counts
-end
-
-function countmap_flatten(samples, flattenhead)
-    counts = Dict{Any,Any}()
-    for s in samples
-        if isexpr(s) && s isa flattenhead # TODO Could you use the TermInterface `operation` here instead of `flattenhead`?
-            coef, term = arguments(s)
-            counts[term] = get(counts, term, 0)+coef
-        else
-            counts[s] = get(counts, s, 0)+1
-        end
-    end
-    counts
-end
-
 ##
 # Metadata cache helpers
 ##
@@ -167,6 +147,13 @@ Base.isequal(::Symbolic{Complex}, ::SymQObj) = false
 # use a macro to provide specializations if that is indeed the case
 propsequal(x,y) = all(n->isequal(getproperty(x,n),getproperty(y,n)), propertynames(x))
 
+
+##
+# Utilities
+##
+
+include("utils.jl")
+
 ##
 # Most symbolic objects defined here
 ##
@@ -174,6 +161,7 @@ propsequal(x,y) = all(n->isequal(getproperty(x,n),getproperty(y,n)), propertynam
 include("literal_objects.jl")
 include("basic_ops_homogeneous.jl")
 include("basic_ops_inhomogeneous.jl")
+include("linalg.jl")
 include("predefined.jl")
 include("predefined_CPTP.jl")
 
diff --git a/src/QSymbolicsBase/basic_ops_homogeneous.jl b/src/QSymbolicsBase/basic_ops_homogeneous.jl
index 63aa58a..2714e52 100644
--- a/src/QSymbolicsBase/basic_ops_homogeneous.jl
+++ b/src/QSymbolicsBase/basic_ops_homogeneous.jl
@@ -22,7 +22,6 @@ julia> 2*A
 @withmetadata struct SScaled{T<:QObj} <: Symbolic{T}
     coeff
     obj
-    SScaled{S}(c,k) where S = _isone(c) ? k : new{S}(c,k)
 end
 isexpr(::SScaled) = true
 iscall(::SScaled) = true
@@ -31,10 +30,14 @@ operation(x::SScaled) = *
 head(x::SScaled) = :*
 children(x::SScaled) = [:*,x.coeff,x.obj]
 function Base.:(*)(c, x::Symbolic{T}) where {T<:QObj} 
-    if iszero(c) || iszero(x)
+    if (isa(c, Number) && iszero(c)) || iszero(x)
         SZero{T}()
+    elseif _isone(c)
+        x
+    elseif isa(x, SScaled)
+        SScaled{T}(c*x.coeff, x.obj)
     else 
-        x isa SScaled ? SScaled{T}(c*x.coeff, x.obj) : SScaled{T}(c, x) 
+        SScaled{T}(c, x) 
     end
 end
 Base.:(*)(x::Symbolic{T}, c) where {T<:QObj} = c*x
@@ -81,8 +84,8 @@ julia> k₁ + k₂
     _arguments_precomputed
 end
 function SAdd{S}(d) where S 
-    xs = [c*obj for (c,obj) in d]
-    length(d)==1 ? first(xs) : SAdd{S}(d,Set(xs),xs)
+    terms = [c*obj for (obj,c) in d]
+    length(d)==1 ? first(terms) : SAdd{S}(d,Set(terms),terms)
 end
 isexpr(::SAdd) = true
 iscall(::SAdd) = true
@@ -99,6 +102,11 @@ end
 Base.:(+)(xs::Vararg{Symbolic{<:QObj},0}) = 0 # to avoid undefined type parameters issue in the above method
 basis(x::SAdd) = basis(first(x.dict).first)
 
+const SAddBra = SAdd{AbstractBra}
+function Base.show(io::IO, x::SAddBra)
+    ordered_terms = sort([repr(i) for i in arguments(x)])
+    print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
+end
 const SAddKet = SAdd{AbstractKet}
 function Base.show(io::IO, x::SAddKet)
     ordered_terms = sort([repr(i) for i in arguments(x)])
@@ -109,11 +117,6 @@ function Base.show(io::IO, x::SAddOperator)
     ordered_terms = sort([repr(i) for i in arguments(x)])
     print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
 end
-const SAddBra = SAdd{AbstractBra}
-function Base.show(io::IO, x::SAddBra)
-    ordered_terms = sort([repr(i) for i in arguments(x)])
-    print(io, "("*join(ordered_terms,"+")::String*")") # type assert to help inference
-end
 
 """Symbolic application of operator on operator
 
@@ -126,10 +129,6 @@ AB
 """
 @withmetadata struct SMulOperator <: Symbolic{AbstractOperator}
     terms
-    function SMulOperator(terms)
-        coeff, cleanterms = prefactorscalings(terms)
-        coeff*new(cleanterms)
-    end
 end
 isexpr(::SMulOperator) = true
 iscall(::SMulOperator) = true
@@ -139,7 +138,13 @@ head(x::SMulOperator) = :*
 children(x::SMulOperator) = [:*;x.terms]
 function Base.:(*)(xs::Symbolic{AbstractOperator}...) 
     zero_ind = findfirst(x->iszero(x), xs)
-    isnothing(zero_ind) ? SMulOperator(collect(xs)) : SZeroOperator()
+    if isnothing(zero_ind)
+        terms = flattenop(*, collect(xs))
+        coeff, cleanterms = prefactorscalings(terms)
+        coeff * SMulOperator(cleanterms)
+    else
+        SZeroOperator()
+    end
 end
 Base.show(io::IO, x::SMulOperator) = print(io, join(map(string, arguments(x)),""))
 basis(x::SMulOperator) = basis(x.terms)
@@ -155,25 +160,27 @@ julia> k₁ ⊗ k₂
 julia> @op A; @op B;
 
 julia> A ⊗ B 
-A⊗B
+(A⊗B)
 ```
 """
 @withmetadata struct STensor{T<:QObj} <: Symbolic{T}
     terms
-    function STensor{S}(terms) where S
-        coeff, cleanterms = prefactorscalings(terms)
-        coeff * new{S}(cleanterms)
-    end
 end
 isexpr(::STensor) = true
 iscall(::STensor) = true
 arguments(x::STensor) = x.terms
 operation(x::STensor) = ⊗
 head(x::STensor) = :⊗
-children(x::STensor) = pushfirst!(x.terms,:⊗)
+children(x::STensor) = [:⊗; x.terms]
 function ⊗(xs::Symbolic{T}...) where {T<:QObj}
     zero_ind = findfirst(x->iszero(x), xs)
-    isnothing(zero_ind) ? STensor{T}(collect(xs)) : SZero{T}()
+    if isnothing(zero_ind)
+        terms = flattenop(⊗, collect(xs))
+        coeff, cleanterms = prefactorscalings(terms)
+        coeff * STensor{T}(cleanterms)
+    else
+        SZero{T}()
+    end
 end
 basis(x::STensor) = tensor(basis.(x.terms)...)
 
@@ -182,9 +189,9 @@ Base.show(io::IO, x::STensorBra) = print(io, join(map(string, arguments(x)),""))
 const STensorKet = STensor{AbstractKet}
 Base.show(io::IO, x::STensorKet) = print(io, join(map(string, arguments(x)),""))
 const STensorOperator = STensor{AbstractOperator}
-Base.show(io::IO, x::STensorOperator) = print(io, join(map(string, arguments(x)),"⊗"))
+Base.show(io::IO, x::STensorOperator) = print(io, "("*join(map(string, arguments(x)),"⊗")*")")
 const STensorSuperOperator = STensor{AbstractSuperOperator}
-Base.show(io::IO, x::STensorSuperOperator) = print(io, join(map(string, arguments(x)),"⊗"))
+Base.show(io::IO, x::STensorSuperOperator) = print(io, "("*join(map(string, arguments(x)),"⊗")*")")
 
 """Symbolic commutator of two operators
 
@@ -201,10 +208,6 @@ julia> commutator(A, A)
 @withmetadata struct SCommutator <: Symbolic{AbstractOperator}
     op1
     op2
-    function SCommutator(o1, o2) 
-        coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
-        cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff*new(cleanterms...)
-    end
 end
 isexpr(::SCommutator) = true
 iscall(::SCommutator) = true
@@ -212,13 +215,15 @@ arguments(x::SCommutator) = [x.op1, x.op2]
 operation(x::SCommutator) = commutator
 head(x::SCommutator) = :commutator
 children(x::SCommutator) = [:commutator, x.op1, x.op2]
-commutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator}) = SCommutator(o1, o2)
+function commutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
+    coeff, cleanterms = prefactorscalings([o1 o2])
+    cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff * SCommutator(cleanterms...)   
+end
 commutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 commutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
 commutator(o1::SZeroOperator, o2::SZeroOperator) = SZeroOperator()
 Base.show(io::IO, x::SCommutator) = print(io, "[$(x.op1),$(x.op2)]")
 basis(x::SCommutator) = basis(x.op1)
-expand(x::SCommutator) = x == 0 ? x : x.op1*x.op2 - x.op2*x.op1
 
 """Symbolic anticommutator of two operators
 
@@ -232,10 +237,6 @@ julia> anticommutator(A, B)
 @withmetadata struct SAnticommutator <: Symbolic{AbstractOperator}
     op1
     op2
-    function SAnticommutator(o1, o2) 
-        coeff, cleanterms = prefactorscalings([o1 o2], scalar=true)
-        coeff*new(cleanterms...)
-    end
 end
 isexpr(::SAnticommutator) = true
 iscall(::SAnticommutator) = true
@@ -243,10 +244,12 @@ arguments(x::SAnticommutator) = [x.op1, x.op2]
 operation(x::SAnticommutator) = anticommutator
 head(x::SAnticommutator) = :anticommutator
 children(x::SAnticommutator) = [:anticommutator, x.op1, x.op2]
-anticommutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator}) = SAnticommutator(o1, o2)
+function anticommutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
+    coeff, cleanterms = prefactorscalings([o1 o2])
+    coeff * SAnticommutator(cleanterms...)
+end
 anticommutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 anticommutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
 anticommutator(o1::SZeroOperator, o2::SZeroOperator) = SZeroOperator()
 Base.show(io::IO, x::SAnticommutator) = print(io, "{$(x.op1),$(x.op2)}")
 basis(x::SAnticommutator) = basis(x.op1)
-expand(x::SAnticommutator) = x == 0 ? x : x.op1*x.op2 + x.op2*x.op1
diff --git a/src/QSymbolicsBase/basic_ops_inhomogeneous.jl b/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
index 49c48ae..902d8f2 100644
--- a/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
+++ b/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
@@ -86,6 +86,8 @@ Base.:(*)(b::SZeroBra, k::Symbolic{AbstractKet}) = 0
 Base.:(*)(b::Symbolic{AbstractBra}, k::SZeroKet) = 0
 Base.:(*)(b::SZeroBra, k::SZeroKet) = 0
 Base.show(io::IO, x::SBraKet) = begin print(io,x.bra); print(io,x.ket) end
+Base.hash(x::SBraKet, h::UInt) = hash((head(x), arguments(x)), h)
+maketerm(::Type{<:SBraKet}, f, a, t, m) = f(a...)
 Base.isequal(x::SBraKet, y::SBraKet) = isequal(x.bra, y.bra) && isequal(x.ket, y.ket)
 
 """Symbolic application of a superoperator on an operator"""
diff --git a/src/QSymbolicsBase/linalg.jl b/src/QSymbolicsBase/linalg.jl
new file mode 100644
index 0000000..1cde9ee
--- /dev/null
+++ b/src/QSymbolicsBase/linalg.jl
@@ -0,0 +1,117 @@
+##
+# Linear algebra operations on quantum objects.
+##
+
+"""Projector for a given ket
+
+```jldoctest
+julia> SProjector(X1⊗X2)
+𝐏[|X₁⟩|X₂⟩]
+
+julia> express(SProjector(X2))
+Operator(dim=2x2)
+  basis: Spin(1/2)
+  0.5+0.0im  -0.5-0.0im
+ -0.5+0.0im   0.5+0.0im
+```"""
+@withmetadata struct SProjector <: Symbolic{AbstractOperator}
+    ket::Symbolic{AbstractKet} # TODO parameterize
+end
+isexpr(::SProjector) = true
+iscall(::SProjector) = true
+arguments(x::SProjector) = [x.ket]
+operation(x::SProjector) = projector
+head(x::SProjector) = :projector
+children(x::SProjector) = [:projector,x.ket]
+projector(x::Symbolic{AbstractKet}) = SProjector(x)
+projector(x::SZeroKet) = SZeroOperator()
+basis(x::SProjector) = basis(x.ket)
+function Base.show(io::IO, x::SProjector)
+    print(io,"𝐏[")
+    print(io,x.ket)
+    print(io,"]")
+end
+
+"""Dagger, i.e., adjoint of quantum objects (kets, bras, operators)
+
+```jldoctest 
+julia> @ket a; @op A;
+
+julia> dagger(2*im*A*a)
+(0 - 2im)|a⟩†A†
+
+julia> @op B;
+
+julia> dagger(A*B)
+B†A†
+
+julia> ℋ = SHermitianOperator(:ℋ); U = SUnitaryOperator(:U);
+
+julia> dagger(ℋ)
+ℋ
+
+julia> dagger(U) 
+U⁻¹
+```
+"""
+@withmetadata struct SDagger{T<:QObj} <: Symbolic{T}
+    obj
+end
+isexpr(::SDagger) = true
+iscall(::SDagger) = true
+arguments(x::SDagger) = [x.obj]
+operation(x::SDagger) = dagger
+head(x::SDagger) = :dagger
+children(x::SDagger) = [:dagger, x.obj]
+dagger(x::Symbolic{AbstractBra}) = SDagger{AbstractKet}(x)
+dagger(x::Symbolic{AbstractKet}) = SDagger{AbstractBra}(x)
+dagger(x::Symbolic{AbstractOperator}) = SDagger{AbstractOperator}(x)
+dagger(x::SScaledKet) = SScaledBra(conj(x.coeff), dagger(x.obj))
+dagger(x::SAdd) = (+)((dagger(i) for i in arguments(x))...)
+dagger(x::SScaledBra) = SScaledKet(conj(x.coeff), dagger(x.obj))
+dagger(x::SZeroOperator) = x
+dagger(x::SHermitianOperator) = x
+dagger(x::SHermitianUnitaryOperator) = x
+dagger(x::SUnitaryOperator) = inv(x)
+dagger(x::STensorBra) = STensorKet(collect(dagger(i) for i in x.terms))
+dagger(x::STensorKet) = STensorBra(collect(dagger(i) for i in x.terms))
+dagger(x::STensorOperator) = STensorOperator(collect(dagger(i) for i in x.terms))
+dagger(x::SScaledOperator) = SScaledOperator(conj(x.coeff), dagger(x.obj))
+dagger(x::SApplyKet) = dagger(x.ket)*dagger(x.op)
+dagger(x::SApplyBra) = dagger(x.op)*dagger(x.bra)
+dagger(x::SMulOperator) = SMulOperator(collect(dagger(i) for i in reverse(x.terms)))
+dagger(x::SBraKet) = SBraKet(dagger(x.ket), dagger(x.bra))
+dagger(x::SOuterKetBra) = SOuterKetBra(dagger(x.bra), dagger(x.ket))
+dagger(x::SDagger) = x.obj
+basis(x::SDagger) = basis(x.obj)
+function Base.show(io::IO, x::SDagger)
+    print(io,x.obj)
+    print(io,"†")
+end
+
+"""Inverse Operator
+
+```jldoctest
+julia> @op A;
+
+julia> inv(A)
+A⁻¹
+
+julia> inv(A)*A
+𝕀
+```
+"""
+@withmetadata struct SInvOperator <: Symbolic{AbstractOperator}
+    op::Symbolic{AbstractOperator}
+end
+isexpr(::SInvOperator) = true
+iscall(::SInvOperator) = true
+arguments(x::SInvOperator) = [x.op]
+operation(x::SInvOperator) = inv
+head(x::SInvOperator) = :inv
+children(x::SInvOperator) = [:inv, x.op]
+basis(x::SInvOperator) = basis(x.op)
+Base.show(io::IO, x::SInvOperator) = print(io, "$(x.op)⁻¹")
+Base.:(*)(invop::SInvOperator, op::SOperator) = isequal(invop.op, op) ? IdentityOp(basis(op)) : SMulOperator(invop, op)
+Base.:(*)(op::SOperator, invop::SInvOperator) = isequal(op, invop.op) ? IdentityOp(basis(op)) : SMulOperator(op, invop)
+inv(x::Symbolic{AbstractOperator}) = SInvOperator(x)
\ No newline at end of file
diff --git a/src/QSymbolicsBase/predefined.jl b/src/QSymbolicsBase/predefined.jl
index e7c28c3..7e599cc 100644
--- a/src/QSymbolicsBase/predefined.jl
+++ b/src/QSymbolicsBase/predefined.jl
@@ -200,121 +200,6 @@ const Destroy = const â = DestroyOp()
 # Other special or useful objects
 ##
 
-"""Projector for a given ket
-
-```jldoctest
-julia> SProjector(X1⊗X2)
-𝐏[|X₁⟩|X₂⟩]
-
-julia> express(SProjector(X2))
-Operator(dim=2x2)
-  basis: Spin(1/2)
-  0.5+0.0im  -0.5-0.0im
- -0.5+0.0im   0.5+0.0im
-```"""
-@withmetadata struct SProjector <: Symbolic{AbstractOperator}
-    ket::Symbolic{AbstractKet} # TODO parameterize
-end
-isexpr(::SProjector) = true
-iscall(::SProjector) = true
-arguments(x::SProjector) = [x.ket]
-operation(x::SProjector) = projector
-head(x::SProjector) = :projector
-children(x::SProjector) = [:projector,x.ket]
-projector(x::Symbolic{AbstractKet}) = SProjector(x)
-projector(x::SZeroKet) = SZeroOperator()
-basis(x::SProjector) = basis(x.ket)
-function Base.show(io::IO, x::SProjector)
-    print(io,"𝐏[")
-    print(io,x.ket)
-    print(io,"]")
-end
-
-"""Dagger, i.e., adjoint of quantum objects (kets, bras, operators)
-
-```jldoctest 
-julia> @ket a; @op A;
-
-julia> dagger(2*im*A*a)
-(0 - 2im)|a⟩†A†
-
-julia> @op B;
-
-julia> dagger(A*B)
-B†A†
-
-julia> ℋ = SHermitianOperator(:ℋ); U = SUnitaryOperator(:U);
-
-julia> dagger(ℋ)
-ℋ
-
-julia> dagger(U) 
-U⁻¹
-```
-"""
-@withmetadata struct SDagger{T<:QObj} <: Symbolic{T}
-    obj
-end
-isexpr(::SDagger) = true
-iscall(::SDagger) = true
-arguments(x::SDagger) = [x.obj]
-operation(x::SDagger) = dagger
-head(x::SDagger) = :dagger
-children(x::SDagger) = [:dagger, x.obj]
-dagger(x::Symbolic{AbstractBra}) = SDagger{AbstractKet}(x)
-dagger(x::Symbolic{AbstractKet}) = SDagger{AbstractBra}(x)
-dagger(x::Symbolic{AbstractOperator}) = SDagger{AbstractOperator}(x)
-dagger(x::SScaledKet) = SScaledBra(conj(x.coeff), dagger(x.obj))
-dagger(x::SAdd) = (+)((dagger(i) for i in arguments(x))...)
-dagger(x::SScaledBra) = SScaledKet(conj(x.coeff), dagger(x.obj))
-dagger(x::SZeroOperator) = x
-dagger(x::SHermitianOperator) = x
-dagger(x::SHermitianUnitaryOperator) = x
-dagger(x::SUnitaryOperator) = inv(x)
-dagger(x::STensorBra) = STensorKet(collect(dagger(i) for i in x.terms))
-dagger(x::STensorKet) = STensorBra(collect(dagger(i) for i in x.terms))
-dagger(x::STensorOperator) = STensorOperator(collect(dagger(i) for i in x.terms))
-dagger(x::SScaledOperator) = SScaledOperator(conj(x.coeff), dagger(x.obj))
-dagger(x::SApplyKet) = dagger(x.ket)*dagger(x.op)
-dagger(x::SApplyBra) = dagger(x.op)*dagger(x.bra)
-dagger(x::SMulOperator) = SMulOperator(collect(dagger(i) for i in reverse(x.terms)))
-dagger(x::SBraKet) = SBraKet(dagger(x.ket), dagger(x.bra))
-dagger(x::SOuterKetBra) = SOuterKetBra(dagger(x.bra), dagger(x.ket))
-dagger(x::SDagger) = x.obj
-basis(x::SDagger) = basis(x.obj)
-function Base.show(io::IO, x::SDagger)
-    print(io,x.obj)
-    print(io,"†")
-end
-symbollabel(x::SDagger) = symbollabel(x.obj)
-
-"""Inverse Operator
-
-```jldoctest
-julia> @op A;
-
-julia> inv(A)
-A⁻¹
-
-julia> inv(A)*A
-𝕀
-```
-"""
-@withmetadata struct SInvOperator <: Symbolic{AbstractOperator}
-    op::Symbolic{AbstractOperator}
-end
-isexpr(::SInvOperator) = true
-iscall(::SInvOperator) = true
-arguments(x::SInvOperator) = [x.op]
-operation(x::SInvOperator) = inv
-head(x::SInvOperator) = :inv
-children(x::SInvOperator) = [:inv, x.op]
-basis(x::SInvOperator) = basis(x.op)
-Base.show(io::IO, x::SInvOperator) = print(io, "$(x.op)⁻¹")
-Base.:(*)(invop::SInvOperator, op::SOperator) = isequal(invop.op, op) ? IdentityOp(basis(op)) : SMulOperator(invop, op)
-Base.:(*)(op::SOperator, invop::SInvOperator) = isequal(op, invop.op) ? IdentityOp(basis(op)) : SMulOperator(op, invop)
-inv(x::Symbolic{AbstractOperator}) = SInvOperator(x)
-
 """Completely depolarized state
 
 ```jldoctest
@@ -336,7 +221,8 @@ julia> express(MixedState(X1⊗X2), CliffordRepr())
 𝒵ₗ━━
 + Z_
 + _Z
-```"""
+```
+"""
 @withmetadata struct MixedState <: Symbolic{AbstractOperator}
     basis::Basis # From QuantumOpticsBase # TODO make QuantumInterface
 end
@@ -355,7 +241,8 @@ julia> IdentityOp(X1⊗X2)
 julia> express(IdentityOp(Z2))
 Operator(dim=2x2)
   basis: Spin(1/2)sparse([1, 2], [1, 2], ComplexF64[1.0 + 0.0im, 1.0 + 0.0im], 2, 2)
-```"""
+```
+"""
 @withmetadata struct IdentityOp <: Symbolic{AbstractOperator}
     basis::Basis # From QuantumOpticsBase # TODO make QuantumInterface
 end
@@ -368,4 +255,4 @@ ishermitian(::IdentityOp) = true
 isunitary(::IdentityOp) = true
 
 """Identity operator in qubit basis"""
-const I = IdentityOp(qubit_basis)
\ No newline at end of file
+const I = IdentityOp(qubit_basis)   
\ No newline at end of file
diff --git a/src/QSymbolicsBase/predefined_CPTP.jl b/src/QSymbolicsBase/predefined_CPTP.jl
index c53fde7..ea6c766 100644
--- a/src/QSymbolicsBase/predefined_CPTP.jl
+++ b/src/QSymbolicsBase/predefined_CPTP.jl
@@ -46,4 +46,4 @@ function Base.show(io::IO, x::GateCPTP)
     print(io, "[")
     print(io, x.gate)
     print(io, "]")
-end
+end
\ No newline at end of file
diff --git a/src/QSymbolicsBase/rules.jl b/src/QSymbolicsBase/rules.jl
index ef4f78c..38db184 100644
--- a/src/QSymbolicsBase/rules.jl
+++ b/src/QSymbolicsBase/rules.jl
@@ -1,5 +1,5 @@
 ##
-# This file defines automatic simplification rules for specific operations of quantum objects.
+# This file defines manual simplification and expansion rules for specific operations of quantum objects.
 ##
 
 ##
@@ -11,61 +11,12 @@ function hasscalings(xs)
     end
 end
 _isa(T) = x->isa(x,T)
-
-##
-# Determining factors for expressions containing quantum objects
-##
-
-function prefactorscalings(xs; scalar=false) # If the scalar keyword is true, then only scalar factors will be returned as coefficients
-    terms = []
-    coeff = 1::Any
-    for x in xs
-        if isexpr(x) && operation(x) == *
-            c,t = arguments(x)
-            if !scalar
-                coeff *= c
-                push!(terms,t)
-            elseif scalar && c isa Number
-                coeff *= c
-                push!(terms, t)
-            else
-                push!(terms,(*)(arguments(x)...))
-            end
-        else
-            push!(terms,x)
-        end
-    end
-    coeff, terms
-end
-
-function prefactorscalings_rule(xs)
-    coeff, terms = prefactorscalings(xs)
-    coeff * ⊗(terms...)
-end
-
-function isnotflat_precheck(*)
-    function (x)
-        operation(x) === (*) || return false
-        args = arguments(x)
-        for t in args
-            if isexpr(t) && operation(t) === (*)
-                return true
-            end
-        end
-        return false
-    end
-end
+_vecisa(T) = x->all(_isa(T), x)
 
 ##
 # Simplification rules
 ## 
 
-# Flattening expressions
-RULES_FLATTEN = [
-    @rule(~x::isnotflat_precheck(⊗) => flatten_term(⊗, ~x)),
-    @rule ⊗(~~xs::hasscalings) => prefactorscalings_rule(xs)
-]
-
 # Pauli identities
 RULES_PAULI = [
     @rule(~o1::_isa(XGate)*~o2::_isa(XGate) => I),
@@ -105,12 +56,12 @@ RULES_ANTICOMMUTATOR = [
     @rule(anticommutator(~o1::_isa(XGate), ~o2::_isa(ZGate)) => 0)
 ]
 
-RULES_ALL = [RULES_PAULI; RULES_COMMUTATOR; RULES_ANTICOMMUTATOR]
+RULES_SIMPLIFY = [RULES_PAULI; RULES_COMMUTATOR; RULES_ANTICOMMUTATOR]
 
 ##
-# Rewriters
+# Simplification rewriters
 ##
-qsimplify_flatten = Chain(RULES_FLATTEN)
+
 qsimplify_anticommutator = Chain(RULES_ANTICOMMUTATOR)
 qsimplify_pauli = Chain(RULES_PAULI)
 qsimplify_commutator = Chain(RULES_COMMUTATOR)
@@ -131,11 +82,58 @@ julia> qsimplify(anticommutator(σˣ, σˣ), rewriter=qsimplify_anticommutator)
 function qsimplify(s; rewriter=nothing)
     if QuantumSymbolics.isexpr(s)
         if isnothing(rewriter)
-            Fixpoint(Prewalk(Chain(RULES_ALL)))(s)
+            Fixpoint(Prewalk(Chain(RULES_SIMPLIFY)))(s)
         else
             Fixpoint(Prewalk(rewriter))(s)
         end
     else
         error("Object $(s) of type $(typeof(s)) is not an expression.")
     end
+end
+
+##
+# Expansion rules
+## 
+
+RULES_EXPAND = [
+    @rule(commutator(~o1, ~o2) => (~o1)*(~o2) - (~o2)*(~o1)),
+    @rule(anticommutator(~o1, ~o2) => (~o1)*(~o2) + (~o2)*(~o1)),
+    @rule(~o1 ⊗ +(~~ops) => +(map(op -> ~o1 ⊗ op, ~~ops)...)),
+    @rule(+(~~ops) ⊗ ~o1 => +(map(op -> op ⊗ ~o1, ~~ops)...)),
+    @rule(~o1 * +(~~ops) => +(map(op -> ~o1 * op, ~~ops)...)),
+    @rule(+(~~ops) * ~o1 => +(map(op -> op * ~o1, ~~ops)...)),
+    @rule(+(~~ops) * ~o1 => +(map(op -> op * ~o1, ~~ops)...)),
+    @rule(⊗(~~ops1::_vecisa(Symbolic{AbstractOperator})) * ⊗(~~ops2::_vecisa(Symbolic{AbstractOperator})) => ⊗(map(*, ~~ops1, ~~ops2)...)),
+    @rule(⊗(~~ops1::_vecisa(Symbolic{AbstractBra})) * ⊗(~~ops2::_vecisa(Symbolic{AbstractKet})) => *(map(*, ~~ops1, ~~ops2)...))
+]
+
+# 
+
+##
+# Expansion rewriter
+##
+
+"""Manually expand a symbolic expression of quantum objects. 
+
+```jldoctest
+julia> @op A; @op B; @op C;
+
+julia> qexpand(commutator(A, B))
+(-1BA+AB)
+
+julia> qexpand(A⊗(B+C))
+((A⊗B)+(A⊗C))
+
+julia> @ket k₁; @ket k₂;
+
+julia> qexpand(A*(k₁+k₂))
+(A|k₁⟩+A|k₂⟩)
+```
+"""
+function qexpand(s)
+    if QuantumSymbolics.isexpr(s)
+        Fixpoint(Prewalk(Chain(RULES_EXPAND)))(s)
+    else
+        error("Object $(s) of type $(typeof(s)) is not an expression.")
+    end
 end
\ No newline at end of file
diff --git a/src/QSymbolicsBase/utils.jl b/src/QSymbolicsBase/utils.jl
new file mode 100644
index 0000000..b14f582
--- /dev/null
+++ b/src/QSymbolicsBase/utils.jl
@@ -0,0 +1,48 @@
+function prefactorscalings(xs)
+    terms = []
+    coeff = 1::Any
+    for x in xs
+        if isa(x, SScaledOperator)
+            coeff *= x.coeff
+            push!(terms, x.obj)
+        elseif isa(x, Union{Number, Symbolic{Number}})
+            coeff *= x
+        else
+            push!(terms,x)
+        end
+    end
+    coeff, terms
+end
+
+function flattenop(f, terms)
+    newterms = []
+    for obj in terms
+        if isexpr(obj) && operation(obj) === f
+            append!(newterms, arguments(obj))
+        else
+            push!(newterms, obj)
+        end
+    end
+    newterms
+end
+
+function countmap(samples) # A simpler version of StatsBase.countmap, because StatsBase is slow to import
+    counts = Dict{Any,Any}()
+    for s in samples
+        counts[s] = get(counts, s, 0)+1
+    end
+    counts
+end
+
+function countmap_flatten(samples, flattenhead)
+    counts = Dict{Any,Any}()
+    for s in samples
+        if isexpr(s) && s isa flattenhead # TODO Could you use the TermInterface `operation` here instead of `flattenhead`?
+            coef, term = arguments(s)
+            counts[term] = get(counts, term, 0)+coef
+        else
+            counts[s] = get(counts, s, 0)+1
+        end
+    end
+    counts
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index bee2e0d..4d26e44 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,6 +35,7 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA
 @doset "anticommutator"
 @doset "dagger"
 @doset "zero_obj"
+@doset "expand"
 
 VERSION >= v"1.9" && @doset "doctests"
 get(ENV,"JET_TEST","")=="true" && @doset "jet"
diff --git a/test/test_expand.jl b/test/test_expand.jl
new file mode 100644
index 0000000..fd7fc30
--- /dev/null
+++ b/test/test_expand.jl
@@ -0,0 +1,31 @@
+using QuantumSymbolics
+using Test
+
+@bra b₁; @bra b₂; @bra b₃;
+@ket k₁; @ket k₂; @ket k₃;
+
+@op A; @op B; @op C; @op D;
+
+@testset "expand rules" begin
+    @test isequal(qexpand(commutator(A, B)), A*B - B*A)
+    @test isequal(qexpand(anticommutator(A, B)), A*B + B*A)
+
+    @test isequal(qexpand(A⊗(B+C+D)), A⊗B + A⊗C + A⊗D)
+    @test isequal(qexpand(C ⊗ commutator(A, B)), C⊗(A*B) - C⊗(B*A))
+    @test isequal(qexpand(k₁⊗(k₂+k₃)), k₁⊗k₂ + k₁⊗k₃)
+    @test isequal(qexpand(b₁⊗(b₂+b₃)), b₁⊗b₂ + b₁⊗b₃)
+
+    @test isequal(qexpand((B+C+D)⊗A), B⊗A + C⊗A + D⊗A)
+    @test isequal(qexpand(commutator(A, B) ⊗ C), (A*B)⊗C - (B*A)⊗C)
+    @test isequal(qexpand((k₂+k₃)⊗k₁), k₂⊗k₁ + k₃⊗k₁)
+    @test isequal(qexpand((b₂+b₃)⊗b₁), b₂⊗b₁ + b₃⊗b₁)
+
+    @test isequal(qexpand(A*(B+C+D)), A*B + A*C + A*D)
+    @test isequal(qexpand(C * commutator(A, B)), C*A*B - C*B*A)
+
+    @test isequal(qexpand((B+C+D)*A), B*A + C*A + D*A)
+    @test isequal(qexpand(commutator(A, B) * C), A*B*C - B*A*C)
+
+    @test isequal(qexpand((A⊗B)*(C⊗D)), (A*C)⊗(B*D))
+    @test isequal(qexpand((b₁⊗b₂)*(k₁⊗k₂)), (b₁*k₁)*(b₂*k₂))
+end
\ No newline at end of file