diff --git a/CHANGELOG.md b/CHANGELOG.md
index f336dbf..b2b1ba5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # News
 
+## v0.3.3 - dev
+
+- Add single qubit simplification rules.
+
 ## v0.3.2 - 2024-07-02
 
 - Added documentation for `express`.
diff --git a/Project.toml b/Project.toml
index 2ed9266..758d68d 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-dev"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
diff --git a/docs/src/index.md b/docs/src/index.md
index 2ce5bd4..70dcbd2 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -202,4 +202,4 @@ express(MixedState(X1)/2+SProjector(Z1)/2, CliffordRepr())
 
 !!! warning "Stabilizer state expressions"
 
-    The state written as $\frac{|Z₁⟩⊗|Z₁⟩+|Z₂⟩⊗|Z₂⟩}{√2}$ is a well known stabilizer state, namely a Bell state. However, automatically expressing it as a stabilizer is a prohibitively expensive computational operation in general. We do not perform that computation automatically. If you want to ensure that states you define can be automatically converted to tableaux for Clifford simulations, avoid using sumation of kets. On the other hand, in all of our Clifford Monte-Carlo simulations, `⊗` is fully supported, as well as [`SProjector`](@ref), [`MixedState`](@ref), [`StabilizerState`](@ref), and sumation of density matrices.
+    The state written as $\frac{|Z₁⟩⊗|Z₁⟩+|Z₂⟩⊗|Z₂⟩}{√2}$ is a well known stabilizer state, namely a Bell state. However, automatically expressing it as a stabilizer is a prohibitively expensive computational operation in general. We do not perform that computation automatically. If you want to ensure that states you define can be automatically converted to tableaux for Clifford simulations, avoid using summation of kets. On the other hand, in all of our Clifford Monte-Carlo simulations, `⊗` is fully supported, as well as [`SProjector`](@ref), [`MixedState`](@ref), [`StabilizerState`](@ref), and summation of density matrices.
diff --git a/src/QSymbolicsBase/rules.jl b/src/QSymbolicsBase/rules.jl
index 38db184..890fea0 100644
--- a/src/QSymbolicsBase/rules.jl
+++ b/src/QSymbolicsBase/rules.jl
@@ -11,6 +11,7 @@ function hasscalings(xs)
     end
 end
 _isa(T) = x->isa(x,T)
+_isequal(obj) = x->(x==obj)
 _vecisa(T) = x->all(_isa(T), x)
 
 ##
@@ -30,7 +31,38 @@ RULES_PAULI = [
     @rule(~o1::_isa(XGate)*~o2::_isa(ZGate) => -im*Y),
     @rule(~o1::_isa(HGate)*~o2::_isa(XGate)*~o3::_isa(HGate) => Z),
     @rule(~o1::_isa(HGate)*~o2::_isa(YGate)*~o3::_isa(HGate) => -Y),
-    @rule(~o1::_isa(HGate)*~o2::_isa(ZGate)*~o3::_isa(HGate) => X)
+    @rule(~o1::_isa(HGate)*~o2::_isa(ZGate)*~o3::_isa(HGate) => X),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(X1) => X1),
+    @rule(~o::_isa(YGate)*~k::_isequal(X1) => -im*X2),
+    @rule(~o::_isa(ZGate)*~k::_isequal(X1) => X2),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(X2) => -X2),
+    @rule(~o::_isa(YGate)*~k::_isequal(X2) => im*X1),
+    @rule(~o::_isa(ZGate)*~k::_isequal(X2) => X1),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(Y1) => im*Y2),
+    @rule(~o::_isa(YGate)*~k::_isequal(Y1) => Y1),
+    @rule(~o::_isa(ZGate)*~k::_isequal(Y1) => Y2),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(Y2) => -im*Y1),
+    @rule(~o::_isa(YGate)*~k::_isequal(Y2) => -Y2),
+    @rule(~o::_isa(ZGate)*~k::_isequal(Y2) => Y1),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(Z1) => Z2),
+    @rule(~o::_isa(YGate)*~k::_isequal(Z1) => im*Z2),
+    @rule(~o::_isa(ZGate)*~k::_isequal(Z1) => Z1),
+
+    @rule(~o::_isa(XGate)*~k::_isequal(Z2) => Z1),
+    @rule(~o::_isa(YGate)*~k::_isequal(Z2) => -im*Z1),
+    @rule(~o::_isa(ZGate)*~k::_isequal(Z2) => -Z2),
+
+    @rule(~o::_isa(HGate)*~k::_isequal(X1) => Z1),
+    @rule(~o::_isa(HGate)*~k::_isequal(X2) => Z2),
+    @rule(~o::_isa(HGate)*~k::_isequal(Y1) => (X1+im*X2)/sqrt(2)),
+    @rule(~o::_isa(HGate)*~k::_isequal(Y2) => (X1-im*X2)/sqrt(2)),
+    @rule(~o::_isa(HGate)*~k::_isequal(Z1) => X1),
+    @rule(~o::_isa(HGate)*~k::_isequal(Z2) => X2)
 ]
 
 # Commutator identities
diff --git a/test/runtests.jl b/test/runtests.jl
index 4d26e44..ad0dce6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -36,6 +36,7 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA
 @doset "dagger"
 @doset "zero_obj"
 @doset "expand"
+@doset "pauli"
 
 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
index fd7fc30..5514709 100644
--- a/test/test_expand.jl
+++ b/test/test_expand.jl
@@ -6,6 +6,10 @@ using Test
 
 @op A; @op B; @op C; @op D;
 
+@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)
diff --git a/test/test_pauli.jl b/test/test_pauli.jl
new file mode 100644
index 0000000..65a912f
--- /dev/null
+++ b/test/test_pauli.jl
@@ -0,0 +1,50 @@
+using QuantumSymbolics
+using Test
+
+@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 "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*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*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