diff --git a/test/Project.toml b/test/Project.toml
index 449d3c3..1527d69 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,9 +13,9 @@ QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
diff --git a/test/runtests.jl b/test/runtests.jl
index 4644d7b..b699694 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,46 +1,21 @@
-using SafeTestsets
+using TestItemRunner
 using QuantumSymbolics
 
-function doset(descr)
-    if length(ARGS) == 0
-        return true
-    end
-    for a in ARGS
-        if occursin(lowercase(a), lowercase(descr))
-            return true
-        end
-    end
-    return false
-end
+# filter for the test
+testfilter = ti -> begin
+  exclude = Symbol[:jet]
+  if !(VERSION >= v"1.10")
+    push!(exclude, :doctests)
+    push!(exclude, :aqua)
+  end
 
-macro doset(descr)
-    quote
-        if doset($descr)
-            @safetestset $descr begin include("test_"*$descr*".jl") end
-        end
-    end
+  return all(!in(exclude), ti.tags)
 end
 
 println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
 
-@doset "sym_expressions"
-@doset "express_opt"
-@doset "express_cliff"
-@doset "qo"
-@doset "qo_qc_interop"
-@doset "basis_consistency"
-@doset "superop"
-@doset "conditional_cliffords"
-@doset "commutator"
-@doset "anticommutator"
-@doset "dagger"
-@doset "zero_obj"
-@doset "trace"
-@doset "expand"
-@doset "misc_linalg"
-@doset "throws"
-@doset "pauli"
+@run_package_tests filter=testfilter
 
-VERSION >= v"1.9" && @doset "doctests"
-get(ENV,"JET_TEST","")=="true" && @doset "jet"
-VERSION >= v"1.9" && @doset "aqua"
+if get(ENV,"JET_TEST","")=="true"
+    @run_package_tests filter=(ti -> :jet in ti.tags)
+end
diff --git a/test/test_anticommutator.jl b/test/test_anticommutator.jl
index 98100b8..ca3b365 100644
--- a/test/test_anticommutator.jl
+++ b/test/test_anticommutator.jl
@@ -1,20 +1,20 @@
-using QuantumSymbolics
-using Test
+@testitem "Anticommutator" begin
 
-@op A; @op B;
+    @op A; @op B;
 
-@testset "symbolic anticommutator tests" begin
-    @test isequal(anticommutator(2*A, B), anticommutator(A, 2*B)) && isequal(2*anticommutator(A, B), anticommutator(2*A, B)) && isequal(2*anticommutator(A, B), anticommutator(2*A, B))
-end
+    @testset "symbolic anticommutator tests" begin
+        @test isequal(anticommutator(2*A, B), anticommutator(A, 2*B)) && isequal(2*anticommutator(A, B), anticommutator(2*A, B)) && isequal(2*anticommutator(A, B), anticommutator(2*A, B))
+    end
 
-@testset "anticommutator Pauli tests" begin
-    @test isequal(qsimplify(anticommutator(X, X), rewriter=qsimplify_anticommutator), 2*I)
-    @test isequal(qsimplify(anticommutator(Y, Y), rewriter=qsimplify_anticommutator), 2*I)
-    @test isequal(qsimplify(anticommutator(Z, Z), rewriter=qsimplify_anticommutator), 2*I)
-    @test isequal(qsimplify(anticommutator(X, Y), rewriter=qsimplify_anticommutator), 0)
-    @test isequal(qsimplify(anticommutator(Y, X), rewriter=qsimplify_anticommutator), 0)
-    @test isequal(qsimplify(anticommutator(Y, Z), rewriter=qsimplify_anticommutator), 0)
-    @test isequal(qsimplify(anticommutator(Z, Y), rewriter=qsimplify_anticommutator), 0)
-    @test isequal(qsimplify(anticommutator(Z, X), rewriter=qsimplify_anticommutator), 0)
-    @test isequal(qsimplify(anticommutator(X, Z), rewriter=qsimplify_anticommutator), 0)
+    @testset "anticommutator Pauli tests" begin
+        @test isequal(qsimplify(anticommutator(X, X), rewriter=qsimplify_anticommutator), 2*I)
+        @test isequal(qsimplify(anticommutator(Y, Y), rewriter=qsimplify_anticommutator), 2*I)
+        @test isequal(qsimplify(anticommutator(Z, Z), rewriter=qsimplify_anticommutator), 2*I)
+        @test isequal(qsimplify(anticommutator(X, Y), rewriter=qsimplify_anticommutator), 0)
+        @test isequal(qsimplify(anticommutator(Y, X), rewriter=qsimplify_anticommutator), 0)
+        @test isequal(qsimplify(anticommutator(Y, Z), rewriter=qsimplify_anticommutator), 0)
+        @test isequal(qsimplify(anticommutator(Z, Y), rewriter=qsimplify_anticommutator), 0)
+        @test isequal(qsimplify(anticommutator(Z, X), rewriter=qsimplify_anticommutator), 0)
+        @test isequal(qsimplify(anticommutator(X, Z), rewriter=qsimplify_anticommutator), 0)
+    end
 end
\ No newline at end of file
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 8d60ebf..24f6440 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -1,6 +1,7 @@
-using Aqua, QuantumSymbolics, Test
-
-Aqua.test_all(QuantumSymbolics,
-ambiguities=(;broken=true),
-piracies=(;broken=true),
-)
+@testitem "Aqua" tags=[:aqua] begin
+    using Aqua
+    Aqua.test_all(QuantumSymbolics,
+            ambiguities=(;broken=true),
+            piracies=(;broken=true),
+    )
+end
diff --git a/test/test_basis_consistency.jl b/test/test_basis_consistency.jl
index 9862da1..d9d7c92 100644
--- a/test/test_basis_consistency.jl
+++ b/test/test_basis_consistency.jl
@@ -1,19 +1,19 @@
-using Test
-using QuantumSymbolics
-using QuantumOptics
+@testitem "Basis consistency" begin
+    using QuantumOptics
 
-@test express(Z*Z1) == express(Z1)
-@test express(Z*Z2) == -express(Z2)
-@test express(X*X1) == express(X1)
-@test express(X*X2) == -express(X2)
-@test express(Y*Y1) == express(Y1)
-@test express(Y*Y2) == -express(Y2)
-@test express(Pm*Z1) == express(Z2)
-@test express(Pp*Z2) == express(Z1)
-@test express(Pm*L0) == express(L1)
-@test express(Pp*L1) == express(L0)
+    @test express(Z*Z1) == express(Z1)
+    @test express(Z*Z2) == -express(Z2)
+    @test express(X*X1) == express(X1)
+    @test express(X*X2) == -express(X2)
+    @test express(Y*Y1) == express(Y1)
+    @test express(Y*Y2) == -express(Y2)
+    @test express(Pm*Z1) == express(Z2)
+    @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);
+    @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
+    @test basis(K) == basis(A)
+end
diff --git a/test/test_commutator.jl b/test/test_commutator.jl
index 23d7520..73b9546 100644
--- a/test/test_commutator.jl
+++ b/test/test_commutator.jl
@@ -1,18 +1,17 @@
-using QuantumSymbolics
-using Test
+@testitem "Commutator" begin
+    @op A; @op B;
 
-@op A; @op B;
+    @testset "symbolic commutator tests" begin
+        @test isequal(commutator(2*A, B), commutator(A, 2*B)) && isequal(2*commutator(A, B), commutator(2*A, B)) && isequal(commutator(A, 2*B), 2*commutator(A, B))
+        @test commutator(A, A) == SZeroOperator()
+    end
 
-@testset "symbolic commutator tests" begin
-    @test isequal(commutator(2*A, B), commutator(A, 2*B)) && isequal(2*commutator(A, B), commutator(2*A, B)) && isequal(commutator(A, 2*B), 2*commutator(A, B))
-    @test commutator(A, A) == SZeroOperator()
+    @testset "commutator Pauli tests" begin
+        @test isequal(qsimplify(commutator(X, Y), rewriter=qsimplify_commutator), 2*im*Z)
+        @test isequal(qsimplify(commutator(Y, X), rewriter=qsimplify_commutator), -2*im*Z)
+        @test isequal(qsimplify(commutator(Y, Z), rewriter=qsimplify_commutator), 2*im*X)
+        @test isequal(qsimplify(commutator(Z, Y), rewriter=qsimplify_commutator), -2*im*X)
+        @test isequal(qsimplify(commutator(Z, X), rewriter=qsimplify_commutator), 2*im*Y)
+        @test isequal(qsimplify(commutator(X, Z), rewriter=qsimplify_commutator), -2*im*Y)
+    end
 end
-
-@testset "commutator Pauli tests" begin
-    @test isequal(qsimplify(commutator(X, Y), rewriter=qsimplify_commutator), 2*im*Z)
-    @test isequal(qsimplify(commutator(Y, X), rewriter=qsimplify_commutator), -2*im*Z)
-    @test isequal(qsimplify(commutator(Y, Z), rewriter=qsimplify_commutator), 2*im*X)
-    @test isequal(qsimplify(commutator(Z, Y), rewriter=qsimplify_commutator), -2*im*X)
-    @test isequal(qsimplify(commutator(Z, X), rewriter=qsimplify_commutator), 2*im*Y)
-    @test isequal(qsimplify(commutator(X, Z), rewriter=qsimplify_commutator), -2*im*Y)
-end
\ No newline at end of file
diff --git a/test/test_conditional_cliffords.jl b/test/test_conditional_cliffords.jl
index bfa9950..4b19245 100644
--- a/test/test_conditional_cliffords.jl
+++ b/test/test_conditional_cliffords.jl
@@ -1,15 +1,15 @@
-using Test
-using QuantumSymbolics
-using QuantumClifford
-using QuantumOpticsBase
-using LinearAlgebra
+@testitem "Conditional Cliffords" begin
+    using QuantumClifford
+    using QuantumOpticsBase
+    using LinearAlgebra
 
-for control in (:X, :Y, :Z)
-    for target in (:X, :Y, :Z)
-        structname = Symbol(control,"C",target,"Gate")
-        gate = eval(structname)()
-        gate_qo = express(gate, QuantumOpticsRepr())
-        gate_qc = Operator(CliffordOperator(express(gate, CliffordRepr(), UseAsOperation())(1,2),2))
-        @test gate_qo ≈ gate_qc
+    for control in (:X, :Y, :Z)
+        for target in (:X, :Y, :Z)
+            structname = Symbol(control,"C",target,"Gate")
+            gate = eval(structname)()
+            gate_qo = express(gate, QuantumOpticsRepr())
+            gate_qc = Operator(CliffordOperator(express(gate, CliffordRepr(), UseAsOperation())(1,2),2))
+            @test gate_qo ≈ gate_qc
+        end
     end
 end
diff --git a/test/test_dagger.jl b/test/test_dagger.jl
index dbabafc..88f9c9c 100644
--- a/test/test_dagger.jl
+++ b/test/test_dagger.jl
@@ -1,30 +1,30 @@
-using QuantumSymbolics
-using QuantumInterface: AbstractOperator
-using Test
+@testitem "Dagger" begin
+    using QuantumInterface: AbstractOperator
 
-@bra b₁; @bra b₂;
-@ket k₁; @ket k₂;
+    @bra b₁; @bra b₂;
+    @ket k₁; @ket k₂;
 
-@op A; @op B; @op C;
-U = SUnitaryOperator(:U); ℋ = SHermitianOperator(:ℋ);
+    @op A; @op B; @op C;
+    U = SUnitaryOperator(:U); ℋ = SHermitianOperator(:ℋ);
 
-@testset "symbolic dagger tests" begin
-    @test isequal(dagger(im*k₁), -im*dagger(k₁))
-    @test isequal(dagger(k₁+k₂), dagger(k₁)+dagger(k₂))
-    @test isequal(dagger(im*b₁), -im*dagger(b₁))
-    @test isequal(dagger(b₁+b₂), dagger(b₁)+dagger(b₂))
-    @test isequal(dagger(A+B), dagger(A) + dagger(B))
-    @test isequal(dagger(ℋ), ℋ)
-    @test isequal(dagger(U), inv(U))
-    @test isequal(dagger(b₁⊗b₂), dagger(b₁)⊗dagger(b₂))
-    @test isequal(dagger(k₁⊗k₂), dagger(k₁)⊗dagger(k₂))
-    @test isequal(dagger(A⊗B), dagger(A)⊗dagger(B))
-    @test isequal(dagger(im*A), -im*dagger(A))
-    @test isequal(dagger(A*k₁), dagger(k₁)*dagger(A))
-    @test isequal(dagger(b₁*A), dagger(A)*dagger(b₁))
-    @test isequal(dagger(A*B*C), dagger(C)*dagger(B)*dagger(A))
-    @test isequal(dagger(b₁*k₁), dagger(k₁)*dagger(b₁))
-    @test isequal(dagger(k₁*b₁), dagger(b₁)*  dagger(k₁))
-    @test isequal(dagger(dagger(A)), A)
-    @test isequal(dagger(dagger(A)), A)
-end
\ No newline at end of file
+    @testset "symbolic dagger tests" begin
+        @test isequal(dagger(im*k₁), -im*dagger(k₁))
+        @test isequal(dagger(k₁+k₂), dagger(k₁)+dagger(k₂))
+        @test isequal(dagger(im*b₁), -im*dagger(b₁))
+        @test isequal(dagger(b₁+b₂), dagger(b₁)+dagger(b₂))
+        @test isequal(dagger(A+B), dagger(A) + dagger(B))
+        @test isequal(dagger(ℋ), ℋ)
+        @test isequal(dagger(U), inv(U))
+        @test isequal(dagger(b₁⊗b₂), dagger(b₁)⊗dagger(b₂))
+        @test isequal(dagger(k₁⊗k₂), dagger(k₁)⊗dagger(k₂))
+        @test isequal(dagger(A⊗B), dagger(A)⊗dagger(B))
+        @test isequal(dagger(im*A), -im*dagger(A))
+        @test isequal(dagger(A*k₁), dagger(k₁)*dagger(A))
+        @test isequal(dagger(b₁*A), dagger(A)*dagger(b₁))
+        @test isequal(dagger(A*B*C), dagger(C)*dagger(B)*dagger(A))
+        @test isequal(dagger(b₁*k₁), dagger(k₁)*dagger(b₁))
+        @test isequal(dagger(k₁*b₁), dagger(b₁)*  dagger(k₁))
+        @test isequal(dagger(dagger(A)), A)
+        @test isequal(dagger(dagger(A)), A)
+    end
+end
diff --git a/test/test_doctests.jl b/test/test_doctests.jl
index 8ac845f..4850ad3 100644
--- a/test/test_doctests.jl
+++ b/test/test_doctests.jl
@@ -1,14 +1,8 @@
-using Documenter
-using QuantumSymbolics
-using QuantumOptics
-using QuantumClifford
-using Test
+@testitem "Doctests" tags=[:doctests] begin
+    using Documenter
+    using QuantumOptics
+    using QuantumClifford
 
-function doctests()
-    @testset "Doctests" begin
-        DocMeta.setdocmeta!(QuantumSymbolics, :DocTestSetup, :(using QuantumSymbolics, QuantumOptics, QuantumClifford); recursive=true)
-        doctest(QuantumSymbolics)
-    end
-end
-
-doctests()
+    DocMeta.setdocmeta!(QuantumSymbolics, :DocTestSetup, :(using QuantumSymbolics, QuantumOptics, QuantumClifford); recursive=true)
+    doctest(QuantumSymbolics)
+end
\ No newline at end of file
diff --git a/test/test_expand.jl b/test/test_expand.jl
index 5514709..443011a 100644
--- a/test/test_expand.jl
+++ b/test/test_expand.jl
@@ -1,35 +1,34 @@
-using QuantumSymbolics
-using Test
+@testitem "Expand" begin
+    @bra b₁; @bra b₂; @bra b₃;
+    @ket k₁; @ket k₂; @ket k₃;
 
-@bra b₁; @bra b₂; @bra b₃;
-@ket k₁; @ket k₂; @ket k₃;
+    @op A; @op B; @op C; @op D;
 
-@op A; @op B; @op C; @op D;
+    @testset "expand errors" begin
+        @test_throws ErrorException qexpand(X)
+    end
 
-@testset "expand errors" begin
-    @test_throws ErrorException qexpand(X)
-end
-
-@testset "expand rules" begin
-    @test isequal(qexpand(commutator(A, B)), A*B - B*A)
-    @test isequal(qexpand(anticommutator(A, B)), A*B + B*A)
+    @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(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((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(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((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
+        @test isequal(qexpand((A⊗B)*(C⊗D)), (A*C)⊗(B*D))
+        @test isequal(qexpand((b₁⊗b₂)*(k₁⊗k₂)), (b₁*k₁)*(b₂*k₂))
+    end
+end
diff --git a/test/test_express_cliff.jl b/test/test_express_cliff.jl
index 3375283..4193020 100644
--- a/test/test_express_cliff.jl
+++ b/test/test_express_cliff.jl
@@ -1,45 +1,46 @@
-using Test
-using QuantumSymbolics
-using QuantumClifford
+@testitem "Clifford" begin
+    using QuantumClifford
+    using QuantumOptics
 
-state = StabilizerState(S"ZZ XX")
-state = SProjector(state)*0.5 + 0.5*MixedState(state)
-state2 = deepcopy(state)
-express(state, CliffordRepr())
-express(state, CliffordRepr())
-express(state2)
-express(state2)
-nocache = @timed express(state2, CliffordRepr())
-withcache = @timed express(state2, CliffordRepr())
-@test nocache.time > 2*withcache.time
-@test withcache.bytes <= 200
-results = Set([express(state2, CliffordRepr()) for i in 1:20])
-@test length(results)==2
+    state = StabilizerState(S"ZZ XX")
+    state = SProjector(state)*0.5 + 0.5*MixedState(state)
+    state2 = deepcopy(state)
+    express(state, CliffordRepr())
+    express(state, CliffordRepr())
+    express(state2)
+    express(state2)
+    nocache = @timed express(state2, CliffordRepr())
+    withcache = @timed express(state2, CliffordRepr())
+    @test nocache.time > 2*withcache.time
+    @test withcache.bytes <= 200
+    results = Set([express(state2, CliffordRepr()) for i in 1:20])
+    @test length(results)==2
 
-CR = CliffordRepr()
-UseOp = UseAsOperation()
-UseObs = UseAsObservable()
+    CR = CliffordRepr()
+    UseOp = UseAsOperation()
+    UseObs = UseAsObservable()
 
-@testset "Clifford representations for basis states" begin
-    isequal(express(X1, CR), MixedDestabilizer(S"X"))
-    isequal(express(X2, CR), MixedDestabilizer(S"-X"))
-    isequal(express(Y1, CR), MixedDestabilizer(S"Y"))
-    isequal(express(Y2, CR), MixedDestabilizer(S"-Y"))
-    isequal(express(Z1, CR), MixedDestabilizer(S"Z"))
-    isequal(express(Z2, CR), MixedDestabilizer(S"-Z"))
-end
+    @testset "Clifford representations for basis states" begin
+        isequal(express(X1, CR), MixedDestabilizer(S"X"))
+        isequal(express(X2, CR), MixedDestabilizer(S"-X"))
+        isequal(express(Y1, CR), MixedDestabilizer(S"Y"))
+        isequal(express(Y2, CR), MixedDestabilizer(S"-Y"))
+        isequal(express(Z1, CR), MixedDestabilizer(S"Z"))
+        isequal(express(Z2, CR), MixedDestabilizer(S"-Z"))
+    end
 
-@testset "Clifford representations as observables" begin
-    isequal(express(σˣ, CR, UseObs), P"X")
-    isequal(express(σʸ, CR, UseObs), P"Y")
-    isequal(express(σᶻ, CR, UseObs), P"Z")
-    isequal(express(im*σˣ, CR, UseObs), im*P"X")
-    isequal(express(σˣ⊗σʸ⊗σᶻ), P"X"⊗P"Y"⊗P"Z")
-    isequal(express(σˣ*σʸ*σᶻ), P"X"*P"Y"*P"Z")
-end
+    @testset "Clifford representations as observables" begin
+        isequal(express(σˣ, CR, UseObs), P"X")
+        isequal(express(σʸ, CR, UseObs), P"Y")
+        isequal(express(σᶻ, CR, UseObs), P"Z")
+        isequal(express(im*σˣ, CR, UseObs), im*P"X")
+        isequal(express(σˣ⊗σʸ⊗σᶻ), P"X"⊗P"Y"⊗P"Z")
+        isequal(express(σˣ*σʸ*σᶻ), P"X"*P"Y"*P"Z")
+    end
 
-@testset "Clifford representations as operations" begin
-    isequal(express(σˣ, CR, UseOp), sX)
-    isequal(express(σʸ, CR, UseOp), sY)
-    isequal(express(σᶻ, CR, UseOp), sZ)
-end
\ No newline at end of file
+    @testset "Clifford representations as operations" begin
+        isequal(express(σˣ, CR, UseOp), sX)
+        isequal(express(σʸ, CR, UseOp), sY)
+        isequal(express(σᶻ, CR, UseOp), sZ)
+    end
+end
diff --git a/test/test_express_opt.jl b/test/test_express_opt.jl
index bb06d0a..9d86060 100644
--- a/test/test_express_opt.jl
+++ b/test/test_express_opt.jl
@@ -1,38 +1,38 @@
-using Test
-using QuantumSymbolics
-import QuantumClifford
-import QuantumOptics
+@testitem "Express opt" begin
+    import QuantumClifford
+    import QuantumOptics
 
-state = 1im*X2⊗Z1+2*Y1⊗(Z2+X2)+StabilizerState("XZ YY")
-express(state)
-express(state)
-state = 1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ")
-nocache = @timed express(state)
-withcache = @timed express(state)
-@test nocache.time > 10*withcache.time
-@test withcache.bytes == 0
-@test nocache.value ≈ withcache.value ≈ express(1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ"))
+    state = 1im*X2⊗Z1+2*Y1⊗(Z2+X2)+StabilizerState("XZ YY")
+    express(state)
+    express(state)
+    state = 1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ")
+    nocache = @timed express(state)
+    withcache = @timed express(state)
+    @test nocache.time > 10*withcache.time
+    @test withcache.bytes == 0
+    @test nocache.value ≈ withcache.value ≈ express(1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ"))
 
-state = 1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ")
-state = SProjector(state)+2*X⊗(Z+Y)/3im
-state = state+MixedState(state)
-state2 = deepcopy(state)
-express(state)
-express(state)
-nocache = @timed express(state2)
-withcache = @timed express(state2)
-@test nocache.time > 50*withcache.time
-@test withcache.bytes == 0
-@test nocache.value ≈ withcache.value ≈ express(state2)
+    state = 1im*X1⊗Z2+2*Y2⊗(Z1+X1)+StabilizerState("YX ZZ")
+    state = SProjector(state)+2*X⊗(Z+Y)/3im
+    state = state+MixedState(state)
+    state2 = deepcopy(state)
+    express(state)
+    express(state)
+    nocache = @timed express(state2)
+    withcache = @timed express(state2)
+    @test nocache.time > 50*withcache.time
+    @test withcache.bytes == 0
+    @test nocache.value ≈ withcache.value ≈ express(state2)
 
-state = 1im*F1⊗F0
-state1 = N⊗Create * state
-@test express(state1) ≈ 1im*express(F1)⊗express(F1)
-@test express(IdentityOp(F1)⊗Destroy)*express(state1) ≈ express((IdentityOp(F1)⊗Destroy)*state1) ≈ express(state)
+    state = 1im*F1⊗F0
+    state1 = N⊗Create * state
+    @test express(state1) ≈ 1im*express(F1)⊗express(F1)
+    @test express(IdentityOp(F1)⊗Destroy)*express(state1) ≈ express((IdentityOp(F1)⊗Destroy)*state1) ≈ express(state)
 
-state = F0⊗X1 + F1⊗Z1
-op = N⊗X
-@test express(op*state) ≈ express(op)*express(state)
-@test express(op*state) ≈ express(F1⊗Z2)
+    state = F0⊗X1 + F1⊗Z1
+    op = N⊗X
+    @test express(op*state) ≈ express(op)*express(state)
+    @test express(op*state) ≈ express(F1⊗Z2)
 
-state = (3im*(2*dagger(Z1)+dagger(Y1))) * (3im*(2*X1+X2))
+    state = (3im*(2*dagger(Z1)+dagger(Y1))) * (3im*(2*X1+X2))
+end
diff --git a/test/test_jet.jl b/test/test_jet.jl
index 60d19c5..e50688f 100644
--- a/test/test_jet.jl
+++ b/test/test_jet.jl
@@ -1,4 +1,5 @@
-using QuantumSymbolics, JET
+@testitem "JET checks" tags=[:jet] begin
+using JET
 using QuantumOptics, QuantumClifford # to load the extensions
 
 using JET: ReportPass, BasicPass, InferenceErrorReport, UncaughtExceptionReport
@@ -30,3 +31,4 @@ rep = report_package("QuantumSymbolics";
 @show rep
 @test_broken length(JET.get_reports(rep)) == 0
 @test length(JET.get_reports(rep)) <= 7
+end
diff --git a/test/test_misc_linalg.jl b/test/test_misc_linalg.jl
index 4517664..5473af5 100644
--- a/test/test_misc_linalg.jl
+++ b/test/test_misc_linalg.jl
@@ -1,29 +1,28 @@
-using QuantumSymbolics
-using Test
+@testitem "Misc linear algebra" begin
+    @op A; @op B;
+    O = SZeroOperator()
 
-@op A; @op B;
-O = SZeroOperator()
+    @testset "Complex Conjugate" begin
+        @test isequal(conj(O), O)
+        @test isequal(conj(conj(A)), A)
+    end
 
-@testset "Complex Conjugate" begin
-    @test isequal(conj(O), O)
-    @test isequal(conj(conj(A)), A)
-end
+    @testset "Transpose" begin
+        @test isequal(transpose(2*A), 2*transpose(A))
+        @test isequal(transpose(A+B), transpose(A)+transpose(B))
+        @test isequal(transpose(A*B), transpose(B)*transpose(A))
+        @test isequal(transpose(A⊗B), transpose(A)⊗transpose(B))
+        @test isequal(transpose(O), O)
+        @test isequal(transpose(transpose(A)), A)
+    end
 
-@testset "Transpose" begin
-    @test isequal(transpose(2*A), 2*transpose(A))
-    @test isequal(transpose(A+B), transpose(A)+transpose(B))
-    @test isequal(transpose(A*B), transpose(B)*transpose(A))
-    @test isequal(transpose(A⊗B), transpose(A)⊗transpose(B))
-    @test isequal(transpose(O), O)
-    @test isequal(transpose(transpose(A)), A)
-end
+    @testset "Vectorization" begin
+        @test isequal(vec(2*A), 2*vec(A))
+        @test isequal(vec(A+B), vec(A)+vec(B))
+        @test isequal(basis(vec(A)), SpinBasis(1//2)⊗SpinBasis(1//2))
+    end
 
-@testset "Vectorization" begin
-    @test isequal(vec(2*A), 2*vec(A))
-    @test isequal(vec(A+B), vec(A)+vec(B))
-    @test isequal(basis(vec(A)), SpinBasis(1//2)⊗SpinBasis(1//2))
+    @testset "Exponential" begin
+        @test isequal(exp(A), SExpOperator(A))
+    end
 end
-
-@testset "Exponential" begin
-    @test isequal(exp(A), SExpOperator(A))
-end
\ No newline at end of file
diff --git a/test/test_pauli.jl b/test/test_pauli.jl
index 65a912f..4644a93 100644
--- a/test/test_pauli.jl
+++ b/test/test_pauli.jl
@@ -1,50 +1,49 @@
-using QuantumSymbolics
-using Test
+@testitem "Test Pauli" begin
+    @testset "simplify errors" begin
+        @test_throws ErrorException qsimplify(X)
+    end
 
-@testset "simplify errors" begin
-    @test_throws ErrorException qsimplify(X)
-end
+    @testset "MulOperator tests" begin
+        @test isequal(qsimplify(X*X,rewriter=qsimplify_pauli), I)
+        @test isequal(qsimplify(Y*Y,rewriter=qsimplify_pauli), I)
+        @test isequal(qsimplify(Z*Z,rewriter=qsimplify_pauli), I)
+        @test isequal(qsimplify(X*Y,rewriter=qsimplify_pauli), im*Z)
+        @test isequal(qsimplify(Y*Z,rewriter=qsimplify_pauli), im*X)
+        @test isequal(qsimplify(Z*X,rewriter=qsimplify_pauli), im*Y)
+        @test isequal(qsimplify(Y*X,rewriter=qsimplify_pauli), -im*Z)
+        @test isequal(qsimplify(Z*Y,rewriter=qsimplify_pauli), -im*X)
+        @test isequal(qsimplify(X*Z,rewriter=qsimplify_pauli), -im*Y)
+        @test isequal(qsimplify(H*X*H,rewriter=qsimplify_pauli), Z)
+        @test isequal(qsimplify(H*Y*H,rewriter=qsimplify_pauli), -Y)
+        @test isequal(qsimplify(H*Z*H,rewriter=qsimplify_pauli), X)
+    end
 
-@testset "MulOperator tests" begin
-    @test isequal(qsimplify(X*X,rewriter=qsimplify_pauli), I)
-    @test isequal(qsimplify(Y*Y,rewriter=qsimplify_pauli), I)
-    @test isequal(qsimplify(Z*Z,rewriter=qsimplify_pauli), I)
-    @test isequal(qsimplify(X*Y,rewriter=qsimplify_pauli), im*Z)
-    @test isequal(qsimplify(Y*Z,rewriter=qsimplify_pauli), im*X)
-    @test isequal(qsimplify(Z*X,rewriter=qsimplify_pauli), im*Y)
-    @test isequal(qsimplify(Y*X,rewriter=qsimplify_pauli), -im*Z)
-    @test isequal(qsimplify(Z*Y,rewriter=qsimplify_pauli), -im*X)
-    @test isequal(qsimplify(X*Z,rewriter=qsimplify_pauli), -im*Y)
-    @test isequal(qsimplify(H*X*H,rewriter=qsimplify_pauli), Z)
-    @test isequal(qsimplify(H*Y*H,rewriter=qsimplify_pauli), -Y)
-    @test isequal(qsimplify(H*Z*H,rewriter=qsimplify_pauli), X)
-end
+    @testset "ApplyKet tests" begin
+        @test isequal(qsimplify(X*X1,rewriter=qsimplify_pauli), X1)
+        @test isequal(qsimplify(Y*X1,rewriter=qsimplify_pauli), -im*X2)
+        @test isequal(qsimplify(Z*X1,rewriter=qsimplify_pauli), X2)
 
-@testset "ApplyKet tests" begin
-   @test isequal(qsimplify(X*X1,rewriter=qsimplify_pauli), X1)
-   @test isequal(qsimplify(Y*X1,rewriter=qsimplify_pauli), -im*X2)
-   @test isequal(qsimplify(Z*X1,rewriter=qsimplify_pauli), X2)
+        @test isequal(qsimplify(X*X2,rewriter=qsimplify_pauli), -X2)
+        @test isequal(qsimplify(Y*X2,rewriter=qsimplify_pauli), im*X1)
+        @test isequal(qsimplify(Z*X2,rewriter=qsimplify_pauli), X1)
 
-   @test isequal(qsimplify(X*X2,rewriter=qsimplify_pauli), -X2)
-   @test isequal(qsimplify(Y*X2,rewriter=qsimplify_pauli), im*X1)
-   @test isequal(qsimplify(Z*X2,rewriter=qsimplify_pauli), X1)
-    
-   @test isequal(qsimplify(X*Y1,rewriter=qsimplify_pauli), im*Y2)
-   @test isequal(qsimplify(Y*Y1,rewriter=qsimplify_pauli), Y1)
-   @test isequal(qsimplify(Z*Y1,rewriter=qsimplify_pauli), Y2)
+        @test isequal(qsimplify(X*Y1,rewriter=qsimplify_pauli), im*Y2)
+        @test isequal(qsimplify(Y*Y1,rewriter=qsimplify_pauli), Y1)
+        @test isequal(qsimplify(Z*Y1,rewriter=qsimplify_pauli), Y2)
 
-   @test isequal(qsimplify(X*Z1,rewriter=qsimplify_pauli), Z2)
-   @test isequal(qsimplify(Y*Z1,rewriter=qsimplify_pauli), im*Z2)
-   @test isequal(qsimplify(Z*Z1,rewriter=qsimplify_pauli), Z1)
+        @test isequal(qsimplify(X*Z1,rewriter=qsimplify_pauli), Z2)
+        @test isequal(qsimplify(Y*Z1,rewriter=qsimplify_pauli), im*Z2)
+        @test isequal(qsimplify(Z*Z1,rewriter=qsimplify_pauli), Z1)
 
-   @test isequal(qsimplify(X*Z2,rewriter=qsimplify_pauli), Z1)
-   @test isequal(qsimplify(Y*Z2,rewriter=qsimplify_pauli), -im*Z1)
-   @test isequal(qsimplify(Z*Z2,rewriter=qsimplify_pauli), -Z2)
+        @test isequal(qsimplify(X*Z2,rewriter=qsimplify_pauli), Z1)
+        @test isequal(qsimplify(Y*Z2,rewriter=qsimplify_pauli), -im*Z1)
+        @test isequal(qsimplify(Z*Z2,rewriter=qsimplify_pauli), -Z2)
 
-   @test isequal(qsimplify(H*X1,rewriter=qsimplify_pauli), Z1)
-   @test isequal(qsimplify(H*X2,rewriter=qsimplify_pauli), Z2)
-   @test isequal(qsimplify(H*Y1,rewriter=qsimplify_pauli), (X1+im*X2)/sqrt(2))
-   @test isequal(qsimplify(H*Y2,rewriter=qsimplify_pauli), (X1-im*X2)/sqrt(2))
-   @test isequal(qsimplify(H*Z1,rewriter=qsimplify_pauli), X1)
-   @test isequal(qsimplify(H*Z2,rewriter=qsimplify_pauli), X2)
-end
\ No newline at end of file
+        @test isequal(qsimplify(H*X1,rewriter=qsimplify_pauli), Z1)
+        @test isequal(qsimplify(H*X2,rewriter=qsimplify_pauli), Z2)
+        @test isequal(qsimplify(H*Y1,rewriter=qsimplify_pauli), (X1+im*X2)/sqrt(2))
+        @test isequal(qsimplify(H*Y2,rewriter=qsimplify_pauli), (X1-im*X2)/sqrt(2))
+        @test isequal(qsimplify(H*Z1,rewriter=qsimplify_pauli), X1)
+        @test isequal(qsimplify(H*Z2,rewriter=qsimplify_pauli), X2)
+    end
+end
diff --git a/test/test_qo.jl b/test/test_qo.jl
index 727f518..2727ff3 100644
--- a/test/test_qo.jl
+++ b/test/test_qo.jl
@@ -1,25 +1,26 @@
-using Test
-using QuantumOptics
-using QuantumSymbolics
-#using QuantumOpticsExt: LazyPrePost
-const LazyPrePost = Base.get_extension(QuantumSymbolics, :QuantumOpticsExt).LazyPrePost
+@testitem "Test qo" begin
+    using QuantumOptics
+    using QuantumSymbolics
+    #using QuantumOpticsExt: LazyPrePost
+    const LazyPrePost = Base.get_extension(QuantumSymbolics, :QuantumOpticsExt).LazyPrePost
 
-bs = GenericBasis(2),GenericBasis(2)
-op0 = Operator(bs...,rand(2,2))
-op21 = Operator(bs...,rand(2,2))
-op22 = Operator(bs...,rand(2,2))
-op31 = Operator(bs...,rand(2,2))
-op32 = Operator(bs...,rand(2,2))
-l2 = LazyPrePost(op21,op22)
-l3 = LazyPrePost(op31,op32)
-@test spre(op21)*spost(op22) ≈ spost(op22)*spre(op21)
-@test spre(op21)*spost(op22)*op0 ≈ l2*op0
-@test spre(op31)*spost(op32)*spre(op21)*spost(op22)*op0 ≈ (l3*l2)*op0 ≈ l3*(l2*op0)
-@test (l2+l3) * op0 ≈ spre(op21)*spost(op22)*op0 + spre(op31)*spost(op32)*op0
+    bs = GenericBasis(2),GenericBasis(2)
+    op0 = Operator(bs...,rand(2,2))
+    op21 = Operator(bs...,rand(2,2))
+    op22 = Operator(bs...,rand(2,2))
+    op31 = Operator(bs...,rand(2,2))
+    op32 = Operator(bs...,rand(2,2))
+    l2 = LazyPrePost(op21,op22)
+    l3 = LazyPrePost(op31,op32)
+    @test spre(op21)*spost(op22) ≈ spost(op22)*spre(op21)
+    @test spre(op21)*spost(op22)*op0 ≈ l2*op0
+    @test spre(op31)*spost(op32)*spre(op21)*spost(op22)*op0 ≈ (l3*l2)*op0 ≈ l3*(l2*op0)
+    @test (l2+l3) * op0 ≈ spre(op21)*spost(op22)*op0 + spre(op31)*spost(op32)*op0
 
-op0a = Operator(bs...,rand(2,2))
-op0b = Operator(bs...,rand(2,2))
-opt0 = op0⊗op0a⊗op0b
-b = basis(opt0)
-@test embed(b,b,[1],l2)*opt0 ≈ (spre(op21)*spost(op22)*op0)⊗op0a⊗op0b
-@test embed(b,b,[1],l2+l3)*opt0 ≈ (spre(op21)*spost(op22)*op0 + spre(op31)*spost(op32)*op0)⊗op0a⊗op0b
+    op0a = Operator(bs...,rand(2,2))
+    op0b = Operator(bs...,rand(2,2))
+    opt0 = op0⊗op0a⊗op0b
+    b = basis(opt0)
+    @test embed(b,b,[1],l2)*opt0 ≈ (spre(op21)*spost(op22)*op0)⊗op0a⊗op0b
+    @test embed(b,b,[1],l2+l3)*opt0 ≈ (spre(op21)*spost(op22)*op0 + spre(op31)*spost(op32)*op0)⊗op0a⊗op0b
+end
diff --git a/test/test_qo_qc_interop.jl b/test/test_qo_qc_interop.jl
index a00109b..b2a9d56 100644
--- a/test/test_qo_qc_interop.jl
+++ b/test/test_qo_qc_interop.jl
@@ -1,30 +1,30 @@
-using Test
-using QuantumSymbolics
-using QuantumClifford: @S_str, random_stabilizer
-using QuantumOpticsBase
-using LinearAlgebra
-#using QuantumOpticsExt: _l0, _l1, _s₊, _s₋, _i₊, _i₋
-const qo = Base.get_extension(QuantumSymbolics, :QuantumOpticsExt)
-const _l0 = qo._l0
-const _l1 = qo._l1
-const _s₊ = qo._s₊
-const _s₋ = qo._s₋
-const _i₊ = qo._i₊
-const _i₋ = qo._i₋
+@testitem "QuantumOptics QuantumClifford Interop" begin
+    using QuantumClifford: @S_str, random_stabilizer
+    using QuantumOpticsBase
+    using LinearAlgebra
+    #using QuantumOpticsExt: _l0, _l1, _s₊, _s₋, _i₊, _i₋
+    const qo = Base.get_extension(QuantumSymbolics, :QuantumOpticsExt)
+    const _l0 = qo._l0
+    const _l1 = qo._l1
+    const _s₊ = qo._s₊
+    const _s₋ = qo._s₋
+    const _i₊ = qo._i₊
+    const _i₋ = qo._i₋
 
-for n in 1:5
-    stabs = [random_stabilizer(1) for _ in 1:n]
-    stab = tensor(stabs...)
-    translate = Dict(S"X"=>_s₊,S"-X"=>_s₋,S"Z"=>_l0,S"-Z"=>_l1,S"Y"=>_i₊,S"-Y"=>_i₋)
-    kets = [translate[s] for s in stabs]
-    ket = tensor(kets...)
-    @test ket.data ≈ Ket(stab).data
+    for n in 1:5
+        stabs = [random_stabilizer(1) for _ in 1:n]
+        stab = tensor(stabs...)
+        translate = Dict(S"X"=>_s₊,S"-X"=>_s₋,S"Z"=>_l0,S"-Z"=>_l1,S"Y"=>_i₊,S"-Y"=>_i₋)
+        kets = [translate[s] for s in stabs]
+        ket = tensor(kets...)
+        @test ket.data ≈ Ket(stab).data
 
-    rstab = random_stabilizer(n)
-    lstab = random_stabilizer(n)
-    lket = Ket(rstab)
-    rket = Ket(lstab)
-    dotket = abs(lket'*rket)
-    dotstab = abs(dot(lstab,rstab))
-    @test (dotket<=1e-10 && dotstab==0) || dotket≈dotstab
+        rstab = random_stabilizer(n)
+        lstab = random_stabilizer(n)
+        lket = Ket(rstab)
+        rket = Ket(lstab)
+        dotket = abs(lket'*rket)
+        dotstab = abs(dot(lstab,rstab))
+        @test (dotket<=1e-10 && dotstab==0) || dotket≈dotstab
+    end
 end
diff --git a/test/test_superop.jl b/test/test_superop.jl
index f9afe8a..98a09a2 100644
--- a/test/test_superop.jl
+++ b/test/test_superop.jl
@@ -1,62 +1,63 @@
-using Test
-using QuantumSymbolics
-using LinearAlgebra
+@testitem "Superop" begin
+    using LinearAlgebra
+    using QuantumOptics
 
-#
-# single qubit
-#
+    #
+    # single qubit
+    #
 
-state = Y₁
-state_dm = SProjector(state)
+    state = Y₁
+    state_dm = SProjector(state)
 
-p = 0.1
+    p = 0.1
 
-noiseop = PauliNoiseCPTP(p,p,p)
-noisy_state = noiseop * state
-noisy_state_fromdm = noiseop * state_dm
+    noiseop = PauliNoiseCPTP(p,p,p)
+    noisy_state = noiseop * state
+    noisy_state_fromdm = noiseop * state_dm
 
-# for one qubit
-# Fρ + (1-F)I/4 = Fρ + (1-F)(ρ + XρX + YρY + ZρZ)/4 = (1-3p)ρ + pXρX + pYρY + pZρZ
-# therefore 3/4F+1/4 = 1-3p
-# therefore 3/4F = 3/4-3p
-F = 1-4p
+    # for one qubit
+    # Fρ + (1-F)I/4 = Fρ + (1-F)(ρ + XρX + YρY + ZρZ)/4 = (1-3p)ρ + pXρX + pYρY + pZρZ
+    # therefore 3/4F+1/4 = 1-3p
+    # therefore 3/4F = 3/4-3p
+    F = 1-4p
 
-mixed_dm = MixedState(state_dm)
-noisy_state_depol = F*state_dm + (1-F)*mixed_dm # TODO make a depolarization helper
+    mixed_dm = MixedState(state_dm)
+    noisy_state_depol = F*state_dm + (1-F)*mixed_dm # TODO make a depolarization helper
 
-@test tr(express(noisy_state)) ≈ tr(express(noisy_state_fromdm)) ≈ tr(express(noisy_state_depol)) ≈ 1
-@test express(noisy_state) ≈ express(noisy_state_fromdm)
-@test express(noisy_state) ≈ express(noisy_state_depol)
+    @test tr(express(noisy_state)) ≈ tr(express(noisy_state_fromdm)) ≈ tr(express(noisy_state_depol)) ≈ 1
+    @test express(noisy_state) ≈ express(noisy_state_fromdm)
+    @test express(noisy_state) ≈ express(noisy_state_depol)
 
-#
-# two qubits
-#
+    #
+    # two qubits
+    #
 
-pure_pair = (Z₁⊗Z₁ + Z₂⊗Z₂) / √2
-pure_pair_dm = SProjector(pure_pair)
+    pure_pair = (Z₁⊗Z₁ + Z₂⊗Z₂) / √2
+    pure_pair_dm = SProjector(pure_pair)
 
-p = 0.1
+    p = 0.1
 
-noiseop = PauliNoiseCPTP(p,p,p)
-noisy_pair = (noiseop ⊗ noiseop) * pure_pair
-noisy_pair_fromdm = (noiseop ⊗ noiseop) * pure_pair_dm
+    noiseop = PauliNoiseCPTP(p,p,p)
+    noisy_pair = (noiseop ⊗ noiseop) * pure_pair
+    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)
+    @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);
+    @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
+    @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)
+    # TODO
+    # test against depolarization
+    # Depolarization over two qubits is different from depolarizing each separately (see related tutorial)
+end
diff --git a/test/test_sym_expressions.jl b/test/test_sym_expressions.jl
index 6ded747..1fb9c7c 100644
--- a/test/test_sym_expressions.jl
+++ b/test/test_sym_expressions.jl
@@ -1,15 +1,14 @@
-using Test
-using QuantumSymbolics
+@testitem "Sym expressions" begin
+    @test +(Z1) == Z1
+    @test +(Z) == Z
+    @test isequal(Z1 - Z2, Z1 + (-Z2))
+    @test_broken isequal(Z1 - 2*Z2 + 2*X1, -2*Z2 + Z1 + 2*X1)
+    @test_broken isequal(Z1 - 2*Z2 + 2*X1, Z1 + 2*(-Z2+X1))
 
-@test +(Z1) == Z1
-@test +(Z) == Z
-@test isequal(Z1 - Z2, Z1 + (-Z2))
-@test_broken isequal(Z1 - 2*Z2 + 2*X1, -2*Z2 + Z1 + 2*X1)
-@test_broken isequal(Z1 - 2*Z2 + 2*X1, Z1 + 2*(-Z2+X1))
+    state1 = XBasisState(1, SpinBasis(1//2))
+    state2 = XBasisState(1, SpinBasis(1//2))
+    state3 = XBasisState(2, SpinBasis(1//2))
 
-state1 = XBasisState(1, SpinBasis(1//2))
-state2 = XBasisState(1, SpinBasis(1//2))
-state3 = XBasisState(2, SpinBasis(1//2))
-
-@test isequal(state1, state2)
-@test !isequal(state1, state3)
+    @test isequal(state1, state2)
+    @test !isequal(state1, state3)
+end
diff --git a/test/test_throws.jl b/test/test_throws.jl
index fec1ca3..c2f9cc6 100644
--- a/test/test_throws.jl
+++ b/test/test_throws.jl
@@ -1,17 +1,17 @@
-using Test
-using QuantumSymbolics
-using QuantumOptics
-using QuantumInterface: IncompatibleBases
+@testitem "Test throwing" begin
+    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);
+    @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
+    @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)
+end
\ No newline at end of file
diff --git a/test/test_trace.jl b/test/test_trace.jl
index 746da93..f525dcd 100644
--- a/test/test_trace.jl
+++ b/test/test_trace.jl
@@ -1,57 +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))
+@testitem "Trace" begin
+
+    @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
 end
diff --git a/test/test_zero_obj.jl b/test/test_zero_obj.jl
index 2225297..188acdc 100644
--- a/test/test_zero_obj.jl
+++ b/test/test_zero_obj.jl
@@ -1,39 +1,36 @@
-using QuantumSymbolics
-using Test
+@testitem "Zero operator" begin
+    A = SOperator(:A, SpinBasis(1//2))
+    k = SKet(:k, SpinBasis(1//2))
+    b = SBra(:b, SpinBasis(1//2))
 
-A = SOperator(:A, SpinBasis(1//2))
-k = SKet(:k, SpinBasis(1//2))
-b = SBra(:b, SpinBasis(1//2))
+    Oop = SZeroOperator()
+    Ok = SZeroKet()
+    Ob = SZeroBra()
 
-Oop = SZeroOperator()
-Ok = SZeroKet()
-Ob = SZeroBra()
+    @testset "zero operator tests" begin
+        @test isequal(0*A, Oop) && isequal(A*0, Oop)
+        @test isequal(2*Oop, Oop) && isequal(Oop*2, Oop) && isequal(Oop/2, Oop)
+        @test isequal(Oop + A, A) && isequal(A + Oop, A) && isequal(Oop + Oop, Oop)
+        @test isequal(Oop*A, Oop) && isequal(A*Oop, Oop)
+        @test isequal(Oop ⊗ A, Oop) && isequal(A ⊗ Oop, Oop) && isequal(Oop*Oop, Oop)
+        @test isequal(commutator(A, Oop), Oop) && isequal(commutator(Oop, A), Oop) && isequal(commutator(Oop, Oop), Oop)
+        @test isequal(anticommutator(A, Oop), Oop) && isequal(anticommutator(Oop, A), Oop) && isequal(anticommutator(Oop, Oop), Oop)
+        @test isequal(projector(Ok), Oop)
+        @test isequal(dagger(Oop), Oop)
+    end
 
-@testset "zero operator tests" begin
-    @test isequal(0*A, Oop) && isequal(A*0, Oop)
-    @test isequal(2*Oop, Oop) && isequal(Oop*2, Oop) && isequal(Oop/2, Oop)
-    @test isequal(Oop + A, A) && isequal(A + Oop, A) && isequal(Oop + Oop, Oop)
-    @test isequal(Oop*A, Oop) && isequal(A*Oop, Oop)
-    @test isequal(Oop ⊗ A, Oop) && isequal(A ⊗ Oop, Oop) && isequal(Oop*Oop, Oop)
-    @test isequal(commutator(A, Oop), Oop) && isequal(commutator(Oop, A), Oop) && isequal(commutator(Oop, Oop), Oop)
-    @test isequal(anticommutator(A, Oop), Oop) && isequal(anticommutator(Oop, A), Oop) && isequal(anticommutator(Oop, Oop), Oop)
-    @test isequal(projector(Ok), Oop)
-    @test isequal(dagger(Oop), Oop)
+    @testset "zero bra and ket tests" begin
+        @test isequal(0*k, Ok) && isequal(k*0, Ok)
+        @test isequal(2*Ok, Ok) && isequal(Ok*2, Ok) && isequal(Ok/2, Ok)
+        @test isequal(Ok + k, k) && isequal(k + Ok, k) && isequal(Ok + Ok, Ok)
+        @test isequal(Ok ⊗ k, Ok) && isequal(k ⊗ Ok, Ok)
+        @test isequal(0*b, Ob) && isequal(b*0, Ob)
+        @test isequal(2*Ob, Ob) && isequal(Ob*2, Ob) && isequal(Ob/2, Ob)
+        @test isequal(Ob + b, b) && isequal(b + Ob, b) && isequal(Ob + Ob, Ob)
+        @test isequal(Ob ⊗ b, Ob) && isequal(b ⊗ Ob, Ob)
+        @test isequal(Oop*k, Ok) && isequal(A*Ok, Ok) && isequal(Oop*Ok, Ok)
+        @test isequal(b*Oop, Ob) && isequal(Ob*A, Ob) && isequal(Ob*Oop, Ob)
+        @test isequal(Ok*b, Oop) && isequal(k*Ob, Oop) && isequal(Ok*Ob, Oop)
+        @test isequal(Ob*k, 0) && isequal(b*Ok, 0) && isequal(Ob*Ok, 0)
+    end
 end
-
-@testset "zero bra and ket tests" begin
-    @test isequal(0*k, Ok) && isequal(k*0, Ok)
-    @test isequal(2*Ok, Ok) && isequal(Ok*2, Ok) && isequal(Ok/2, Ok)
-    @test isequal(Ok + k, k) && isequal(k + Ok, k) && isequal(Ok + Ok, Ok)
-    @test isequal(Ok ⊗ k, Ok) && isequal(k ⊗ Ok, Ok)
-    @test isequal(0*b, Ob) && isequal(b*0, Ob)
-    @test isequal(2*Ob, Ob) && isequal(Ob*2, Ob) && isequal(Ob/2, Ob)
-    @test isequal(Ob + b, b) && isequal(b + Ob, b) && isequal(Ob + Ob, Ob)
-    @test isequal(Ob ⊗ b, Ob) && isequal(b ⊗ Ob, Ob)
-    @test isequal(Oop*k, Ok) && isequal(A*Ok, Ok) && isequal(Oop*Ok, Ok)
-    @test isequal(b*Oop, Ob) && isequal(Ob*A, Ob) && isequal(Ob*Oop, Ob)
-    @test isequal(Ok*b, Oop) && isequal(k*Ob, Oop) && isequal(Ok*Ob, Oop)
-    @test isequal(Ob*k, 0) && isequal(b*Ok, 0) && isequal(Ob*Ok, 0)
-end
-
-