diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4978657..ff068e4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,10 @@
 # News
 
+## v0.3.4 - dev
+
+- Added `tr` and `ptrace` functionalities.
+- New symbolic superoperator representations.
+  
 ## v0.3.3 - 2024-07-12
 
 - Added single qubit simplification rules.
diff --git a/Project.toml b/Project.toml
index 58d45f7..fcd3a93 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.3"
+version = "0.3.4-dev"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
diff --git a/src/QSymbolicsBase/QSymbolicsBase.jl b/src/QSymbolicsBase/QSymbolicsBase.jl
index b53808f..10e2f93 100644
--- a/src/QSymbolicsBase/QSymbolicsBase.jl
+++ b/src/QSymbolicsBase/QSymbolicsBase.jl
@@ -1,7 +1,7 @@
 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,sorted_arguments
 using TermInterface
 import TermInterface: isexpr,head,iscall,children,operation,arguments,metadata,maketerm
 
@@ -11,9 +11,9 @@ import LinearAlgebra: eigvecs,ishermitian,inv
 import QuantumInterface:
     apply!,
     tensor, ⊗,
-    basis,Basis,SpinBasis,FockBasis,
+    basis,Basis,samebases,IncompatibleBases,SpinBasis,FockBasis,CompositeBasis,
     nqubits,
-    projector,dagger,
+    projector,dagger,tr,ptrace,
     AbstractBra,AbstractKet,AbstractOperator,AbstractSuperOperator
 
 export SymQObj,QObj,
@@ -23,28 +23,30 @@ export SymQObj,QObj,
        apply!,
        express,
        tensor,⊗,
-       dagger,projector,commutator,anticommutator,
+       dagger,projector,commutator,anticommutator,tr,ptrace,
        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₋ᵢ,
        vac,F₀,F0,F₁,F1,
-       N,n̂,Create,âꜛ,Destroy,â,SpinBasis,FockBasis,
-       SBra,SKet,SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,
-       @ket,@bra,@op,
+       N,n̂,Create,âꜛ,Destroy,â,basis,SpinBasis,FockBasis,
+       SBra,SKet,SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,SSuperOperator,
+       @ket,@bra,@op,@superop,
        SAdd,SAddBra,SAddKet,SAddOperator,
        SScaled,SScaledBra,SScaledOperator,SScaledKet,
        STensorBra,STensorKet,STensorOperator,
        SZeroBra,SZeroKet,SZeroOperator,
-       SProjector,MixedState,IdentityOp,SInvOperator,
-       SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SDagger,SBraKet,SOuterKetBra,
+       MixedState,IdentityOp,
+       SProjector,SDagger,STrace,SPartialTrace,SInvOperator,
+       SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,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_commutator,qsimplify_anticommutator,
        qexpand,
-       isunitary
-
+       isunitary,
+       KrausRepr,kraus
+       
 ##
 # Metadata cache helpers
 ##
@@ -161,6 +163,7 @@ include("utils.jl")
 include("literal_objects.jl")
 include("basic_ops_homogeneous.jl")
 include("basic_ops_inhomogeneous.jl")
+include("basic_superops.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 2714e52..97322a5 100644
--- a/src/QSymbolicsBase/basic_ops_homogeneous.jl
+++ b/src/QSymbolicsBase/basic_ops_homogeneous.jl
@@ -139,15 +139,19 @@ children(x::SMulOperator) = [:*;x.terms]
 function Base.:(*)(xs::Symbolic{AbstractOperator}...) 
     zero_ind = findfirst(x->iszero(x), xs)
     if isnothing(zero_ind)
-        terms = flattenop(*, collect(xs))
-        coeff, cleanterms = prefactorscalings(terms)
-        coeff * SMulOperator(cleanterms)
+        if any(x->!(samebases(basis(x),basis(first(xs)))),xs)
+            throw(IncompatibleBases())
+        else
+            terms = flattenop(*, collect(xs))
+            coeff, cleanterms = prefactorscalings(terms)
+            coeff * SMulOperator(cleanterms)
+        end
     else
         SZeroOperator()
     end
 end
 Base.show(io::IO, x::SMulOperator) = print(io, join(map(string, arguments(x)),""))
-basis(x::SMulOperator) = basis(x.terms)
+basis(x::SMulOperator) = basis(first(x.terms))
 
 """Tensor product of quantum objects (kets, operators, or bras)
 
@@ -216,8 +220,12 @@ operation(x::SCommutator) = commutator
 head(x::SCommutator) = :commutator
 children(x::SCommutator) = [:commutator, x.op1, x.op2]
 function commutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
-    coeff, cleanterms = prefactorscalings([o1 o2])
-    cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff * SCommutator(cleanterms...)   
+    if !(samebases(basis(o1),basis(o2)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([o1 o2])
+        cleanterms[1] === cleanterms[2] ? SZeroOperator() : coeff * SCommutator(cleanterms...)  
+    end 
 end
 commutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 commutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
@@ -245,8 +253,12 @@ operation(x::SAnticommutator) = anticommutator
 head(x::SAnticommutator) = :anticommutator
 children(x::SAnticommutator) = [:anticommutator, x.op1, x.op2]
 function anticommutator(o1::Symbolic{AbstractOperator}, o2::Symbolic{AbstractOperator})
-    coeff, cleanterms = prefactorscalings([o1 o2])
-    coeff * SAnticommutator(cleanterms...)
+    if !(samebases(basis(o1),basis(o2)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([o1 o2])
+        coeff * SAnticommutator(cleanterms...)
+    end
 end
 anticommutator(o1::SZeroOperator, o2::Symbolic{AbstractOperator}) = SZeroOperator()
 anticommutator(o1::Symbolic{AbstractOperator}, o2::SZeroOperator) = SZeroOperator()
diff --git a/src/QSymbolicsBase/basic_ops_inhomogeneous.jl b/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
index 902d8f2..616639d 100644
--- a/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
+++ b/src/QSymbolicsBase/basic_ops_inhomogeneous.jl
@@ -14,10 +14,6 @@ A|k⟩
 @withmetadata struct SApplyKet <: Symbolic{AbstractKet}
     op
     ket
-    function SApplyKet(o, k)
-        coeff, cleanterms = prefactorscalings([o k])
-        coeff*new(cleanterms...)
-    end
 end
 isexpr(::SApplyKet) = true
 iscall(::SApplyKet) = true
@@ -25,7 +21,14 @@ arguments(x::SApplyKet) = [x.op,x.ket]
 operation(x::SApplyKet) = *
 head(x::SApplyKet) = :*
 children(x::SApplyKet) = [:*,x.op,x.ket]
-Base.:(*)(op::Symbolic{AbstractOperator}, k::Symbolic{AbstractKet}) = SApplyKet(op,k)
+function Base.:(*)(op::Symbolic{AbstractOperator}, k::Symbolic{AbstractKet})
+    if !(samebases(basis(op),basis(k)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([op k])
+        coeff*SApplyKet(cleanterms...)
+    end
+end
 Base.:(*)(op::SZeroOperator, k::Symbolic{AbstractKet}) = SZeroKet()
 Base.:(*)(op::Symbolic{AbstractOperator}, k::SZeroKet) = SZeroKet()
 Base.:(*)(op::SZeroOperator, k::SZeroKet) = SZeroKet()
@@ -44,10 +47,6 @@ julia> b*A
 @withmetadata struct SApplyBra <: Symbolic{AbstractBra}
     bra
     op
-    function SApplyBra(b, o)
-        coeff, cleanterms = prefactorscalings([b o])
-        coeff*new(cleanterms...)
-    end
 end
 isexpr(::SApplyBra) = true
 iscall(::SApplyBra) = true
@@ -55,7 +54,14 @@ arguments(x::SApplyBra) = [x.bra,x.op]
 operation(x::SApplyBra) = *
 head(x::SApplyBra) = :*
 children(x::SApplyBra) = [:*,x.bra,x.op]
-Base.:(*)(b::Symbolic{AbstractBra}, op::Symbolic{AbstractOperator}) = SApplyBra(b,op)
+function Base.:(*)(b::Symbolic{AbstractBra}, op::Symbolic{AbstractOperator}) 
+    if !(samebases(basis(b),basis(op)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([b op])
+        coeff*SApplyBra(cleanterms...)
+    end
+end
 Base.:(*)(b::SZeroBra, op::Symbolic{AbstractOperator}) = SZeroBra()
 Base.:(*)(b::Symbolic{AbstractBra}, op::SZeroOperator) = SZeroBra()
 Base.:(*)(b::SZeroBra, op::SZeroOperator) = SZeroBra()
@@ -81,33 +87,22 @@ arguments(x::SBraKet) = [x.bra,x.ket]
 operation(x::SBraKet) = *
 head(x::SBraKet) = :*
 children(x::SBraKet) = [:*,x.bra,x.ket]
-Base.:(*)(b::Symbolic{AbstractBra}, k::Symbolic{AbstractKet}) = SBraKet(b,k)
+function Base.:(*)(b::Symbolic{AbstractBra}, k::Symbolic{AbstractKet}) 
+    if !(samebases(basis(b),basis(k)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([b k])
+        coeff == 1 ? SBraKet(cleanterms...) : coeff*SBraKet(cleanterms...)
+    end
+end
 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...)
+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"""
-@withmetadata struct SSuperOpApply <: Symbolic{AbstractOperator}
-    sop
-    op
-end
-isexpr(::SSuperOpApply) = true
-iscall(::SSuperOpApply) = true
-arguments(x::SSuperOpApply) = [x.sop,x.op]
-operation(x::SSuperOpApply) = *
-head(x::SSuperOpApply) = :*
-children(x::SSuperOpApply) = [:*,x.sop,x.op]
-Base.:(*)(sop::Symbolic{AbstractSuperOperator}, op::Symbolic{AbstractOperator}) = SSuperOpApply(sop,op)
-Base.:(*)(sop::Symbolic{AbstractSuperOperator}, op::SZeroOperator) = SZeroOperator()
-Base.:(*)(sop::Symbolic{AbstractSuperOperator}, k::Symbolic{AbstractKet}) = SSuperOpApply(sop,SProjector(k))
-Base.:(*)(sop::Symbolic{AbstractSuperOperator}, k::SZeroKet) = SZeroKet()
-Base.show(io::IO, x::SSuperOpApply) = begin print(io, x.sop); print(io, x.op) end
-basis(x::SSuperOpApply) = basis(x.op)
-
 """Symbolic outer product of a ket and a bra
 ```jldoctest 
 julia> @bra b; @ket k;
@@ -119,10 +114,6 @@ julia> k*b
 @withmetadata struct SOuterKetBra <: Symbolic{AbstractOperator}
     ket
     bra
-    function SOuterKetBra(k, b)
-        coeff, cleanterms = prefactorscalings([k b])
-        coeff*new(cleanterms...)
-    end
 end
 isexpr(::SOuterKetBra) = true
 iscall(::SOuterKetBra) = true
@@ -130,7 +121,14 @@ arguments(x::SOuterKetBra) = [x.ket,x.bra]
 operation(x::SOuterKetBra) = *
 head(x::SOuterKetBra) = :*
 children(x::SOuterKetBra) = [:*,x.ket,x.bra]
-Base.:(*)(k::Symbolic{AbstractKet}, b::Symbolic{AbstractBra}) = SOuterKetBra(k,b)
+function Base.:(*)(k::Symbolic{AbstractKet}, b::Symbolic{AbstractBra})
+    if !(samebases(basis(k),basis(b)))
+        throw(IncompatibleBases())
+    else
+        coeff, cleanterms = prefactorscalings([k b])
+        coeff*SOuterKetBra(cleanterms...)
+    end
+end
 Base.:(*)(k::SZeroKet, b::Symbolic{AbstractBra}) = SZeroOperator()
 Base.:(*)(k::Symbolic{AbstractKet}, b::SZeroBra) = SZeroOperator()
 Base.:(*)(k::SZeroKet, b::SZeroBra) = SZeroOperator()
diff --git a/src/QSymbolicsBase/basic_superops.jl b/src/QSymbolicsBase/basic_superops.jl
new file mode 100644
index 0000000..d6f0363
--- /dev/null
+++ b/src/QSymbolicsBase/basic_superops.jl
@@ -0,0 +1,62 @@
+##
+# Superoperator representations
+##
+
+"""Kraus representation of a quantum channel
+
+```jldoctest
+julia> @op A₁; @op A₂; @op A₃;
+
+julia> K = kraus(A₁, A₂, A₃)
+𝒦(A₁,A₂,A₃)
+
+julia> @op ρ;
+
+julia> K*ρ
+(A₁ρA₁†+A₂ρA₂†+A₃ρA₃†)
+```
+"""
+@withmetadata struct KrausRepr <: Symbolic{AbstractSuperOperator}
+    krausops
+end
+isexpr(::KrausRepr) = true
+iscall(::KrausRepr) = true
+arguments(x::KrausRepr) = x.krausops
+operation(x::KrausRepr) = kraus
+head(x::KrausRepr) = :kraus
+children(x::KrausRepr) = [:kraus; x.krausops]
+kraus(xs::Symbolic{AbstractOperator}...) = KrausRepr(collect(xs))
+basis(x::KrausRepr) = basis(first(x.krausops))
+Base.show(io::IO, x::KrausRepr) = print(io, "𝒦("*join([symbollabel(i) for i in x.krausops], ",")*")")
+
+##
+# Superoperator operations
+##
+
+"""Symbolic application of a superoperator on an operator
+
+```jldoctest
+julia> @op A; @superop S;
+
+julia> S*A
+S[A]
+```
+"""
+@withmetadata struct SSuperOpApply <: Symbolic{AbstractOperator}
+    sop
+    op
+end
+isexpr(::SSuperOpApply) = true
+iscall(::SSuperOpApply) = true
+arguments(x::SSuperOpApply) = [x.sop,x.op]
+operation(x::SSuperOpApply) = *
+head(x::SSuperOpApply) = :*
+children(x::SSuperOpApply) = [:*,x.sop,x.op]
+Base.:(*)(sop::Symbolic{AbstractSuperOperator}, op::Symbolic{AbstractOperator}) = SSuperOpApply(sop,op)
+Base.:(*)(sop::Symbolic{AbstractSuperOperator}, op::SZeroOperator) = SZeroOperator()
+Base.:(*)(sop::Symbolic{AbstractSuperOperator}, k::Symbolic{AbstractKet}) = SSuperOpApply(sop,SProjector(k))
+Base.:(*)(sop::Symbolic{AbstractSuperOperator}, k::SZeroKet) = SZeroOperator()
+Base.:(*)(sop::KrausRepr, op::Symbolic{AbstractOperator}) = (+)((i*op*dagger(i) for i in sop.krausops)...)
+Base.:(*)(sop::KrausRepr, k::Symbolic{AbstractKet}) = (+)((i*SProjector(k)*dagger(i) for i in sop.krausops)...)
+Base.show(io::IO, x::SSuperOpApply) = print(io, "$(x.sop)[$(x.op)]")
+basis(x::SSuperOpApply) = basis(x.op)
\ No newline at end of file
diff --git a/src/QSymbolicsBase/linalg.jl b/src/QSymbolicsBase/linalg.jl
index 1cde9ee..c4cd36e 100644
--- a/src/QSymbolicsBase/linalg.jl
+++ b/src/QSymbolicsBase/linalg.jl
@@ -89,6 +89,160 @@ function Base.show(io::IO, x::SDagger)
     print(io,"†")
 end
 
+"""Trace of an operator
+
+```jldoctest
+julia> @op A; @op B;
+
+julia> tr(A)
+tr(A)
+
+julia> tr(commutator(A, B))
+0
+
+julia> @bra b; @ket k;
+
+julia> tr(k*b)
+⟨b||k⟩
+```
+"""
+@withmetadata struct STrace <: Symbolic{Complex}
+    op::Symbolic{AbstractOperator}
+end
+isexpr(::STrace) = true
+iscall(::STrace) = true
+arguments(x::STrace) = [x.op]
+sorted_arguments(x::STrace) = arguments(x)
+operation(x::STrace) = tr
+head(x::STrace) = :tr
+children(x::STrace) = [:tr, x.op]
+Base.show(io::IO, x::STrace) = print(io, "tr($(x.op))")
+tr(x::Symbolic{AbstractOperator}) = STrace(x)
+tr(x::SScaled{AbstractOperator}) = x.coeff*tr(x.obj)
+tr(x::SAdd{AbstractOperator}) = (+)((tr(i) for i in arguments(x))...)
+tr(x::SOuterKetBra) = x.bra*x.ket
+tr(x::SCommutator) = 0
+tr(x::STensorOperator) = (*)((tr(i) for i in arguments(x))...)
+Base.hash(x::STrace, h::UInt) = hash((head(x), arguments(x)), h)
+Base.isequal(x::STrace, y::STrace) = isequal(x.op, y.op)
+
+"""Partial trace over system i of a composite quantum system
+
+```jldoctest
+julia> @op 𝒪 SpinBasis(1//2)⊗SpinBasis(1//2);
+
+julia> op = ptrace(𝒪, 1)
+tr1(𝒪)
+
+julia> QuantumSymbolics.basis(op)
+Spin(1/2)
+
+julia> @op A; @op B;
+
+julia> ptrace(A⊗B, 1)
+(tr(A))B
+
+julia> @ket k; @bra b;
+
+julia> factorizable = A ⊗ (k*b)
+(A⊗|k⟩⟨b|)
+
+julia> ptrace(factorizable, 1)
+(tr(A))|k⟩⟨b|
+
+julia> ptrace(factorizable, 2)
+(⟨b||k⟩)A
+
+julia> mixed_state = (A⊗(k*b)) + ((k*b)⊗B)
+((A⊗|k⟩⟨b|)+(|k⟩⟨b|⊗B))
+
+julia> ptrace(mixed_state, 1)
+((0 + ⟨b||k⟩)B+(tr(A))|k⟩⟨b|)
+
+julia> ptrace(mixed_state, 2)
+((0 + ⟨b||k⟩)A+(tr(B))|k⟩⟨b|)
+```
+"""
+@withmetadata struct SPartialTrace <: Symbolic{AbstractOperator}
+    obj
+    sys::Int
+end
+isexpr(::SPartialTrace) = true
+iscall(::SPartialTrace) = true
+arguments(x::SPartialTrace) = [x.obj, x.sys]
+operation(x::SPartialTrace) = ptrace
+head(x::SPartialTrace) = :ptrace
+children(x::SPartialTrace) = [:ptrace, x.obj, x.sys]
+function basis(x::SPartialTrace)
+    obj_bases = collect(basis(x.obj).bases)
+    new_bases = deleteat!(copy(obj_bases), x.sys)
+    tensor(new_bases...)
+end
+Base.show(io::IO, x::SPartialTrace) = print(io, "tr$(x.sys)($(x.obj))")
+function ptrace(x::Symbolic{AbstractOperator}, s) 
+    ex = isexpr(x) ? qexpand(x) : x
+    if isa(ex, typeof(x))
+        if isa(basis(x), CompositeBasis)
+            SPartialTrace(x, s)
+        elseif s==1
+            tr(x)
+        else
+            throw(ArgumentError("cannot take partial trace of a single quantum system"))
+        end
+    else
+        ptrace(ex, s)
+    end
+end
+function ptrace(x::SAddOperator, s)
+    add_terms = []
+    if isa(basis(x), CompositeBasis)
+        for i in arguments(x)
+            if isexpr(i)
+                if isa(i, SScaledOperator) && operation(i.obj) === ⊗  # scaled tensor product
+                    prod_terms = arguments(i.obj)
+                    coeff = i.coeff
+                elseif operation(i) === ⊗  # tensor product
+                    prod_terms = arguments(i)
+                    coeff = 1
+                else  # multiplication of operators with composite basis
+                    return SPartialTrace(x,s)
+                end
+            else  # operator with composite basis
+                return SPartialTrace(x,s)  
+            end
+            if any(j -> isa(basis(j), CompositeBasis), prod_terms)  # tensor product of operators with composite bases
+                return SPartialTrace(x,s)
+            else  # tensor product without composite basis
+                sys_op = coeff*prod_terms[s]
+                new_terms = deleteat!(copy(prod_terms), s)
+                trace = tr(sys_op)
+                isone(length(new_terms)) ? push!(add_terms, trace*first(new_terms)) : push!(add_terms, trace*(⊗)(new_terms...))
+            end  
+        end
+        (+)(add_terms...)
+    elseif s==1 # partial trace must be over the first system if sum does not have a composite basis
+        tr(x)
+    else
+        throw(ArgumentError("cannot take partial trace of a single quantum system"))
+    end
+end
+function ptrace(x::STensorOperator, s)
+    ex = qexpand(x)
+    if isa(ex, SAddOperator)
+        ptrace(ex, s)
+    else
+        terms = arguments(ex)
+        newterms = []
+        if any(i -> isa(basis(i), CompositeBasis), terms)
+            SPartialTrace(ex, s)
+        else 
+            sys_op = terms[s]
+            new_terms = deleteat!(copy(terms), s)
+            isone(length(new_terms)) ? tr(sys_op)*first(new_terms) : tr(sys_op)*STensorOperator(new_terms)
+        end
+    end
+end
+
 """Inverse Operator
 
 ```jldoctest
diff --git a/src/QSymbolicsBase/literal_objects.jl b/src/QSymbolicsBase/literal_objects.jl
index d6ebfbc..a65c850 100644
--- a/src/QSymbolicsBase/literal_objects.jl
+++ b/src/QSymbolicsBase/literal_objects.jl
@@ -67,7 +67,21 @@ SHermitianUnitaryOperator(name) = SHermitianUnitaryOperator(name, qubit_basis)
 ishermitian(::SHermitianUnitaryOperator) = true
 isunitary(::SHermitianUnitaryOperator) = true
 
-const SymQ = Union{SKet,SBra,SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator}
+struct SSuperOperator <: Symbolic{AbstractSuperOperator}
+    name::Symbol
+    basis::Basis
+end
+SSuperOperator(name) = SSuperOperator(name, qubit_basis)
+macro superop(name, basis)
+    :($(esc(name)) = SSuperOperator($(Expr(:quote, name)), $(basis)))
+end
+macro superop(name)
+    :($(esc(name)) = SSuperOperator($(Expr(:quote, name))))
+end
+ishermitian(x::SSuperOperator) = false
+isunitary(x::SSuperOperator) = false
+
+const SymQ = Union{SKet,SBra,SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,SSuperOperator}
 isexpr(::SymQ) = false
 metadata(::SymQ) = nothing
 symbollabel(x::SymQ) = x.name
@@ -75,7 +89,7 @@ basis(x::SymQ) = x.basis
 
 Base.show(io::IO, x::SKet) = print(io, "|$(symbollabel(x))⟩")
 Base.show(io::IO, x::SBra) = print(io, "⟨$(symbollabel(x))|")
-Base.show(io::IO, x::Union{SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator}) = print(io, "$(symbollabel(x))")
+Base.show(io::IO, x::Union{SOperator,SHermitianOperator,SUnitaryOperator,SHermitianUnitaryOperator,SSuperOperator}) = print(io, "$(symbollabel(x))")
 Base.show(io::IO, x::SymQObj) = print(io, symbollabel(x)) # fallback that probably is not great
 
 struct SZero{T<:QObj} <: Symbolic{T} end
diff --git a/src/QSymbolicsBase/predefined.jl b/src/QSymbolicsBase/predefined.jl
index 7e599cc..2eb194d 100644
--- a/src/QSymbolicsBase/predefined.jl
+++ b/src/QSymbolicsBase/predefined.jl
@@ -180,7 +180,7 @@ const CPHASE = CPHASEGate()
 abstract type AbstractSingleBosonOp <: Symbolic{AbstractOperator} end
 abstract type AbstractSingleBosonGate <: AbstractSingleBosonOp end # TODO maybe an IsUnitaryTrait is a better choice
 isexpr(::AbstractSingleBosonGate) = false
-basis(::AbstractSingleBosonGate) = inf_fock_basis
+basis(x::AbstractSingleBosonOp) = inf_fock_basis
 
 @withmetadata struct NumberOp <: AbstractSingleBosonOp end
 symbollabel(::NumberOp) = "n"
diff --git a/src/QSymbolicsBase/rules.jl b/src/QSymbolicsBase/rules.jl
index 890fea0..b7d1a79 100644
--- a/src/QSymbolicsBase/rules.jl
+++ b/src/QSymbolicsBase/rules.jl
@@ -134,9 +134,8 @@ RULES_EXPAND = [
     @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{AbstractBra})) * ⊗(~~ops2::_vecisa(Symbolic{AbstractKet})) => *(map(*, ~~ops1, ~~ops2)...)),
     @rule(⊗(~~ops1::_vecisa(Symbolic{AbstractOperator})) * ⊗(~~ops2::_vecisa(Symbolic{AbstractOperator})) => ⊗(map(*, ~~ops1, ~~ops2)...)),
-    @rule(⊗(~~ops1::_vecisa(Symbolic{AbstractBra})) * ⊗(~~ops2::_vecisa(Symbolic{AbstractKet})) => *(map(*, ~~ops1, ~~ops2)...))
 ]
 
 # 
diff --git a/src/QSymbolicsBase/utils.jl b/src/QSymbolicsBase/utils.jl
index b14f582..5afe9f1 100644
--- a/src/QSymbolicsBase/utils.jl
+++ b/src/QSymbolicsBase/utils.jl
@@ -2,7 +2,7 @@ function prefactorscalings(xs)
     terms = []
     coeff = 1::Any
     for x in xs
-        if isa(x, SScaledOperator)
+        if isa(x, SScaled)
             coeff *= x.coeff
             push!(terms, x.obj)
         elseif isa(x, Union{Number, Symbolic{Number}})
diff --git a/test/runtests.jl b/test/runtests.jl
index ad0dce6..5a13c9e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,7 +35,9 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA
 @doset "anticommutator"
 @doset "dagger"
 @doset "zero_obj"
+@doset "trace"
 @doset "expand"
+@doset "throws"
 @doset "pauli"
 
 VERSION >= v"1.9" && @doset "doctests"
diff --git a/test/test_basis_consistency.jl b/test/test_basis_consistency.jl
index df4da34..9862da1 100644
--- a/test/test_basis_consistency.jl
+++ b/test/test_basis_consistency.jl
@@ -1,5 +1,6 @@
 using Test
 using QuantumSymbolics
+using QuantumOptics
 
 @test express(Z*Z1) == express(Z1)
 @test express(Z*Z2) == -express(Z2)
@@ -11,3 +12,8 @@ using QuantumSymbolics
 @test express(Pp*Z2) == express(Z1)
 @test express(Pm*L0) == express(L1)
 @test express(Pp*L1) == express(L0)
+
+@op A; @op B; @op C; @op O; @ket k;
+@superop S; K = kraus(A, B, C);
+
+@test basis(K) == basis(A)
\ No newline at end of file
diff --git a/test/test_superop.jl b/test/test_superop.jl
index 07e4b20..f9afe8a 100644
--- a/test/test_superop.jl
+++ b/test/test_superop.jl
@@ -44,6 +44,19 @@ noisy_pair_fromdm = (noiseop ⊗ noiseop) * pure_pair_dm
 @test tr(express(noisy_pair)) ≈ tr(express(noisy_pair_fromdm)) ≈ 1
 @test express(noisy_pair) ≈ express(noisy_pair_fromdm)
 
+@op A; @op B; @op C; @op O; @ket k;
+@superop S; K = kraus(A, B, C);
+
+@testset "symbolic superoperator tests" begin
+    @test ishermitian(S) == false
+    @test isunitary(S) == false
+    @test isequal(S*SZeroOperator(), SZeroOperator())
+    @test isequal(S*SZeroKet(), SZeroOperator())
+    @test isequal(S*k, S*projector(k))
+    @test isequal(K*O, A*O*dagger(A) + B*O*dagger(B) + C*O*dagger(C))
+    @test isequal(K*k, A*projector(k)*dagger(A) + B*projector(k)*dagger(B) + C*projector(k)*dagger(C))
+end
+
 # TODO
 # test against depolarization
 # Depolarization over two qubits is different from depolarizing each separately (see related tutorial)
diff --git a/test/test_throws.jl b/test/test_throws.jl
new file mode 100644
index 0000000..fec1ca3
--- /dev/null
+++ b/test/test_throws.jl
@@ -0,0 +1,17 @@
+using Test
+using QuantumSymbolics
+using QuantumOptics
+using QuantumInterface: IncompatibleBases
+
+@op A SpinBasis(1//2) ⊗ SpinBasis(1//2); @op B; @op C;
+@ket k; @bra b; @ket l SpinBasis(1//2) ⊗ SpinBasis(1//2);
+
+@test_throws IncompatibleBases A*B
+@test_throws IncompatibleBases commutator(A, B)
+@test_throws IncompatibleBases anticommutator(A, B)
+@test_throws IncompatibleBases A*k
+@test_throws IncompatibleBases b*A
+@test_throws IncompatibleBases l*b
+@test_throws IncompatibleBases b*l
+@test_throws ArgumentError ptrace(B, 2)
+@test_throws ArgumentError ptrace(B+C, 2)
\ No newline at end of file
diff --git a/test/test_trace.jl b/test/test_trace.jl
new file mode 100644
index 0000000..746da93
--- /dev/null
+++ b/test/test_trace.jl
@@ -0,0 +1,57 @@
+using QuantumSymbolics
+using Test
+
+@bra b₁; @bra b₂;
+@ket k₁; @ket k₂;
+@op A; @op B; @op C; @op D; @op E; @op F;
+@op 𝒪 SpinBasis(1//2)⊗SpinBasis(1//2); @op 𝒫 SpinBasis(1//2)⊗SpinBasis(1//2);
+@op ℒ SpinBasis(1//2)⊗SpinBasis(1//2);
+
+@testset "trace tests" begin
+    @test isequal(tr(2*A), 2*tr(A))
+    @test isequal(tr(A+B), tr(A)+tr(B))
+    @test isequal(tr(k₁*b₁), b₁*k₁)
+    @test isequal(tr(commutator(A, B)), 0)
+    @test isequal(tr(A⊗B⊗C), tr(A)*tr(B)*tr(C))
+end
+
+exp1 = (k₁*b₁)⊗A + (k₂*b₂)⊗B
+exp2 = A⊗(B⊗C + D⊗E)
+@testset "partial trace tests" begin
+
+    # tests for ptrace(x::Symbolic{AbstractOperator}, s)
+    @test isequal(ptrace(𝒪, 1), SPartialTrace(𝒪, 1))
+    @test isequal(ptrace(𝒪, 2), SPartialTrace(𝒪, 2))
+    @test isequal(ptrace(A, 1), tr(A))
+    @test isequal(ptrace(A*(B+C), 1), tr(A*B)+tr(A*C))
+    @test isequal(QuantumSymbolics.basis(ptrace(𝒪, 1)), SpinBasis(1//2))
+    @test isequal(QuantumSymbolics.basis(ptrace(𝒪, 2)), SpinBasis(1//2))
+
+    # tests for ptrace(x::SAddOperator, s)
+    @test isequal(ptrace(A+B, 1), tr(A+B))
+    @test isequal(ptrace(2*(A⊗B)+(C⊗D), 1), 2*tr(A)*B + tr(C)*D)
+    @test isequal(ptrace((A⊗B)+(C⊗D), 1), tr(A)*B + tr(C)*D)
+    @test isequal(ptrace((A⊗B⊗C)+(D⊗E⊗F), 1), tr(A)*(B⊗C) + tr(D)*(E⊗F))
+    @test isequal(ptrace(𝒪 + 𝒫, 1), SPartialTrace(𝒪 + 𝒫, 1))
+    @test isequal(ptrace(𝒪*ℒ + 𝒫*ℒ, 1), SPartialTrace(𝒪*ℒ + 𝒫*ℒ, 1))
+    @test isequal(ptrace(𝒪⊗ℒ + 𝒫⊗ℒ, 1), SPartialTrace(𝒪⊗ℒ + 𝒫⊗ℒ, 1))
+
+    # tests for ptrace(x::STensorOperator, s)
+    @test isequal(ptrace(A⊗(B⊗C + D⊗E), 1),  tr(A)*(B⊗C) + tr(A)*(D⊗E)) 
+    @test isequal(ptrace(𝒪⊗A, 1), SPartialTrace(𝒪⊗A, 1))
+    @test isequal(ptrace(A⊗B, 1), tr(A)*B)
+    @test isequal(ptrace(A⊗B⊗C, 1), tr(A)*(B⊗C))
+
+    # additional tests 
+    @test isequal(ptrace(exp1, 1), (b₁*k₁)*A + (b₂*k₂)*B)
+    @test isequal(basis(ptrace(exp1, 1)), SpinBasis(1//2))
+    @test isequal(ptrace(exp1, 2), tr(A)*(k₁*b₁) + tr(B)*(k₂*b₂))
+    @test isequal(basis(ptrace(exp1, 2)), SpinBasis(1//2))
+
+    @test isequal(ptrace(exp2, 1), tr(A)*(B⊗C) + tr(A)*(D⊗E))
+    @test isequal(basis(ptrace(exp2, 1)), SpinBasis(1//2)⊗SpinBasis(1//2))
+    @test isequal(ptrace(exp2, 2), tr(B)*(A⊗C) + tr(D)*(A⊗E))
+    @test isequal(basis(ptrace(exp2, 2)), SpinBasis(1//2)⊗SpinBasis(1//2))
+    @test isequal(ptrace(exp2, 3), tr(C)*(A⊗B) + tr(E)*(A⊗D))
+    @test isequal(basis(ptrace(exp2, 3)), SpinBasis(1//2)⊗SpinBasis(1//2))
+end