diff --git a/ext/QuantumOpticsExt/QuantumOpticsExt.jl b/ext/QuantumOpticsExt/QuantumOpticsExt.jl
index 0046ff6..c1a592e 100644
--- a/ext/QuantumOpticsExt/QuantumOpticsExt.jl
+++ b/ext/QuantumOpticsExt/QuantumOpticsExt.jl
@@ -4,7 +4,7 @@ using QuantumInterface, QuantumOpticsBase
 using QuantumInterface: samebases
 using QuantumSymbolics
 using QuantumSymbolics:
-    HGate, XGate, YGate, ZGate, CPHASEGate, CNOTGate, PauliP, PauliM,
+    HGate, XGate, YGate, ZGate, RotXGate, RotYGate, RotZGate, CPHASEGate, CNOTGate, PauliP, PauliM,
     XCXGate, XCYGate, XCZGate, YCXGate, YCYGate, YCZGate, ZCXGate, ZCYGate, ZCZGate,
     XBasisState, YBasisState, ZBasisState,
     NumberOp, CreateOp, DestroyOp,
@@ -31,6 +31,9 @@ const _z = sigmaz(_b2)
 const _x = sigmax(_b2)
 const _y = sigmay(_b2)
 const _hadamard = (sigmaz(_b2)+sigmax(_b2))/√2
+_rotx(θ) = cos(θ/2)*_id - im*sin(θ/2)*_x
+_roty(θ) = cos(θ/2)*_id - im*sin(θ/2)*_y
+_rotz(θ) = cos(θ/2)*_id - im*sin(θ/2)*_z
 const _cnot = _l00⊗_id + _l11⊗_x
 const _cphase = _l00⊗_id + _l11⊗_z
 const _phase = _l00 + im*_l11
@@ -47,6 +50,9 @@ express_nolookup(::HGate, ::QuantumOpticsRepr) = _hadamard
 express_nolookup(::XGate, ::QuantumOpticsRepr) = _x
 express_nolookup(::YGate, ::QuantumOpticsRepr) = _y
 express_nolookup(::ZGate, ::QuantumOpticsRepr) = _z
+express_nolookup(g::RotXGate, ::QuantumOpticsRepr) = _rotx(g.θ)
+express_nolookup(g::RotYGate, ::QuantumOpticsRepr) = _roty(g.θ)
+express_nolookup(g::RotZGate, ::QuantumOpticsRepr) = _rotz(g.θ)
 express_nolookup(::CPHASEGate, ::QuantumOpticsRepr) = _cphase
 express_nolookup(::CNOTGate, ::QuantumOpticsRepr) = _cnot
 
diff --git a/src/QSymbolicsBase/QSymbolicsBase.jl b/src/QSymbolicsBase/QSymbolicsBase.jl
index 6f73ffc..8b7a85b 100644
--- a/src/QSymbolicsBase/QSymbolicsBase.jl
+++ b/src/QSymbolicsBase/QSymbolicsBase.jl
@@ -29,7 +29,7 @@ export SymQObj,QObj,
        tensor,⊗,
        dagger,projector,commutator,anticommutator,conj,transpose,inv,exp,vec,tr,ptrace,
        I,X,Y,Z,σˣ,σʸ,σᶻ,Pm,Pp,σ₋,σ₊,
-       H,CNOT,CPHASE,XCX,XCY,XCZ,YCX,YCY,YCZ,ZCX,ZCY,ZCZ,
+       H,Rx,Ry,Rz,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,inf_fock_basis,
        N,n̂,Create,âꜛ,Destroy,â,basis,SpinBasis,FockBasis,
@@ -42,12 +42,12 @@ export SymQObj,QObj,
        SConjugate,STranspose,SProjector,SDagger,SInvOperator,SExpOperator,SVec,STrace,SPartialTrace,
        MixedState,IdentityOp,
        SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SBraKet,SOuterKetBra,
-       HGate,XGate,YGate,ZGate,CPHASEGate,CNOTGate,
+       HGate,XGate,YGate,ZGate,RotXGate,RotYGate,RotZGate,CPHASEGate,CNOTGate,
        XBasisState,YBasisState,ZBasisState,FockState,CoherentState,SqueezedState,TwoSqueezedState,BosonicThermalState,
        NumberOp,CreateOp,DestroyOp,PhaseShiftOp,DisplaceOp,SqueezeOp,
        TwoSqueezeOp,BeamSplitterOp,
        XCXGate,XCYGate,XCZGate,YCXGate,YCYGate,YCZGate,ZCXGate,ZCYGate,ZCZGate,
-       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,qsimplify_fock,
+       qsimplify,qsimplify_pauli,qsimplify_commutator,qsimplify_anticommutator,qsimplify_fock,qsimplify_rot,
        qexpand,
        isunitary,
        KrausRepr,kraus
diff --git a/src/QSymbolicsBase/predefined.jl b/src/QSymbolicsBase/predefined.jl
index 8b77b80..5f8182f 100644
--- a/src/QSymbolicsBase/predefined.jl
+++ b/src/QSymbolicsBase/predefined.jl
@@ -106,6 +106,43 @@ symbollabel(::HGate) = "H"
 ishermitian(::HGate) = true
 isunitary(::HGate) = true
 
+abstract type AbstractRotGate <: AbstractSingleQubitGate end
+ishermitian(::AbstractRotGate) = false
+isunitary(::AbstractRotGate) = true
+isexpr(::AbstractRotGate) = true
+iscall(::AbstractRotGate) = true
+arguments(x::AbstractRotGate) = [x.θ]
+
+@withmetadata struct RotXGate <: AbstractRotGate
+    θ
+end
+operation(::RotXGate) = Rx
+head(::RotXGate) = :Rx
+children(x::RotXGate) = [:Rx, x.θ]
+symbollabel(x::RotXGate) = "Rx($(x.θ))"
+"""Rotation around the X axis"""
+Rx(θ) = (θ == 0) ? I : RotXGate(θ)
+
+@withmetadata struct RotYGate <: AbstractRotGate
+    θ
+end
+operation(::RotYGate) = Ry
+head(::RotYGate) = :Ry
+children(x::RotYGate) = [:Ry, x.θ]
+symbollabel(x::RotYGate) = "Ry($(x.θ))"
+"""Rotation around the Y axis"""
+Ry(θ) = (θ == 0) ? I : RotYGate(θ)
+
+@withmetadata struct RotZGate <: AbstractRotGate
+    θ
+end
+operation(::RotZGate) = Rz
+head(::RotZGate) = :Rz
+children(x::RotZGate) = [:Rz, x.θ]
+symbollabel(x::RotZGate) = "Rz($(x.θ))"
+"""Rotation around the Z axis"""
+Rz(θ) = (θ == 0) ? I : RotZGate(θ)
+
 @withmetadata struct CNOTGate <: AbstractTwoQubitGate end
 symbollabel(::CNOTGate) = "CNOT"
 ishermitian(::CNOTGate) = true
diff --git a/src/QSymbolicsBase/rules.jl b/src/QSymbolicsBase/rules.jl
index 1664ffe..6fd6be4 100644
--- a/src/QSymbolicsBase/rules.jl
+++ b/src/QSymbolicsBase/rules.jl
@@ -109,7 +109,25 @@ RULES_FOCK = [
     @rule(~o::_isa(TwoSqueezeOp) * ~k::isequal(vac ⊗ vac) => TwoSqueezedState((~o).z))
 ]
 
-RULES_SIMPLIFY = [RULES_PAULI; RULES_COMMUTATOR; RULES_ANTICOMMUTATOR; RULES_FOCK]
+RULES_ROT = [
+    @rule(~r::_isa(RotXGate) => (-im*X) where (~r).θ in [π, 1π]),
+    @rule(~r::_isa(RotYGate) => (-im*Y) where (~r).θ in [π, 1π]),
+    @rule(~r::_isa(RotZGate) => (-im*Z) where (~r).θ in [π, 1π]),
+    @rule(~r::_isa(RotXGate) => try Rx(mod((~r).θ, 4π)) catch end),
+    @rule(~r::_isa(RotYGate) => try Ry(mod((~r).θ, 4π)) catch end),
+    @rule(~r::_isa(RotZGate) => try Rz(mod((~r).θ, 4π)) catch end),
+    @rule(~r::_isa(RotXGate) => try (~r).θ ≥ 2π ? -Rx((~r).θ - 2π) : nothing catch end),
+    @rule(~r::_isa(RotYGate) => try (~r).θ ≥ 2π ? -Ry((~r).θ - 2π) : nothing catch end),
+    @rule(~r::_isa(RotZGate) => try (~r).θ ≥ 2π ? -Rz((~r).θ - 2π) : nothing catch end),
+    @rule(~r1::_isa(RotXGate) * ~r2::_isa(RotXGate) => try Rx((~r1).θ + (~r2).θ) catch end),
+    @rule(~r1::_isa(RotYGate) * ~r2::_isa(RotYGate) => try Ry((~r1).θ + (~r2).θ) catch end),
+    @rule(~r1::_isa(RotZGate) * ~r2::_isa(RotZGate) => try Rz((~r1).θ + (~r2).θ) catch end),
+    @rule(exp(~α * ~x::_isa(XGate)) => try if real(~α) == 0 Rx(-2imag(~α)) end catch end),
+    @rule(exp(~α * ~x::_isa(YGate)) => try if real(~α) == 0 Ry(-2imag(~α)) end catch end),
+    @rule(exp(~α * ~x::_isa(ZGate)) => try if real(~α) == 0 Rz(-2imag(~α)) end catch end)
+]
+
+RULES_SIMPLIFY = [RULES_PAULI; RULES_COMMUTATOR; RULES_ANTICOMMUTATOR; RULES_FOCK; RULES_ROT]
 
 ##
 # Simplification rewriters
@@ -119,6 +137,7 @@ qsimplify_pauli = Chain(RULES_PAULI)
 qsimplify_commutator = Chain(RULES_COMMUTATOR)
 qsimplify_anticommutator = Chain(RULES_ANTICOMMUTATOR)
 qsimplify_fock = Chain(RULES_FOCK)
+qsimplify_rot = Chain(RULES_ROT)
 
 """
     qsimplify(s; rewriter=nothing)
diff --git a/test/test_rotation.jl b/test/test_rotation.jl
new file mode 100644
index 0000000..f8df318
--- /dev/null
+++ b/test/test_rotation.jl
@@ -0,0 +1,67 @@
+@testitem "Test Rotation" begin
+    @testset "Identity tests" begin
+        @test isequal(Rx(0), I)
+        @test isequal(Ry(0), I)
+        @test isequal(Rz(0), I)
+    end
+
+    @testset "Pauli tests" begin
+        @test isequal(qsimplify(Rx(π), rewriter=qsimplify_rot), -im*X)
+        @test isequal(qsimplify(Ry(π), rewriter=qsimplify_rot), -im*Y)
+        @test isequal(qsimplify(Rz(π), rewriter=qsimplify_rot), -im*Z)
+
+        @test isequal(qsimplify(Rx(1π), rewriter=qsimplify_rot), -im*X)
+        @test isequal(qsimplify(Ry(1π), rewriter=qsimplify_rot), -im*Y)
+        @test isequal(qsimplify(Rz(1π), rewriter=qsimplify_rot), -im*Z)
+    end
+
+    @testset "Fusion tests" begin 
+        @test isequal(qsimplify(Rx(π/3) * Rx(π/3), rewriter=qsimplify_rot), Rx(2π/3))
+        @test isequal(qsimplify(Ry(π/3) * Ry(π/3), rewriter=qsimplify_rot), Ry(2π/3))
+        @test isequal(qsimplify(Rz(π/3) * Rz(π/3), rewriter=qsimplify_rot), Rz(2π/3))
+
+        @test isequal(qsimplify(Rx(π/2) * Rx(-π/2), rewriter=qsimplify_rot), I)
+        @test isequal(qsimplify(Ry(π/2) * Ry(-π/2), rewriter=qsimplify_rot), I)
+        @test isequal(qsimplify(Rz(π/2) * Rz(-π/2), rewriter=qsimplify_rot), I)
+
+        @test isequal(qsimplify(2 * Rx(π/2) * Rx(π/2), rewriter=qsimplify_rot), -2im*X)
+        @test isequal(qsimplify(2 * Ry(π/2) * Ry(π/2), rewriter=qsimplify_rot), -2im*Y)
+        @test isequal(qsimplify(2 * Rz(π/2) * Rz(π/2), rewriter=qsimplify_rot), -2im*Z)
+    end
+
+    @testset "Modulo tests" begin
+        @test isequal(qsimplify(Rx(2π), rewriter=qsimplify_rot), -I)
+        @test isequal(qsimplify(Ry(2π), rewriter=qsimplify_rot), -I)
+        @test isequal(qsimplify(Rz(2π), rewriter=qsimplify_rot), -I)
+
+        @test isequal(qsimplify(Rx(3π), rewriter=qsimplify_rot), im*X)
+        @test isequal(qsimplify(Ry(3π), rewriter=qsimplify_rot), im*Y)
+        @test isequal(qsimplify(Rz(3π), rewriter=qsimplify_rot), im*Z)
+
+        @test isequal(qsimplify(Rx(4π), rewriter=qsimplify_rot), I)
+        @test isequal(qsimplify(Ry(4π), rewriter=qsimplify_rot), I)
+        @test isequal(qsimplify(Rz(4π), rewriter=qsimplify_rot), I)
+
+        @test isequal(qsimplify(Rx(5π), rewriter=qsimplify_rot), -im*X)
+        @test isequal(qsimplify(Ry(5π), rewriter=qsimplify_rot), -im*Y)
+        @test isequal(qsimplify(Rz(5π), rewriter=qsimplify_rot), -im*Z)
+
+        @test isequal(qsimplify(Rx(π) * Rx(π), rewriter=qsimplify_rot), -I)
+        @test isequal(qsimplify(Ry(π) * Ry(π), rewriter=qsimplify_rot), -I)
+        @test isequal(qsimplify(Rz(π) * Rz(π), rewriter=qsimplify_rot), -I)
+
+        @test isequal(qsimplify(2 * Rx(3π/2) * Rx(2π/2), rewriter=qsimplify_rot), -2Rx(π/2))
+        @test isequal(qsimplify(2 * Ry(3π/2) * Ry(2π/2), rewriter=qsimplify_rot), -2Ry(π/2))
+        @test isequal(qsimplify(2 * Rz(3π/2) * Rz(2π/2), rewriter=qsimplify_rot), -2Rz(π/2))
+    end
+
+    @testset "Exponential tests" begin
+        @test isequal(qsimplify(exp(-im * π/2 * X), rewriter=qsimplify_rot), -im*X)
+        @test isequal(qsimplify(exp(-im * π/2 * Y), rewriter=qsimplify_rot), -im*Y)
+        @test isequal(qsimplify(exp(-im * π/2 * Z), rewriter=qsimplify_rot), -im*Z)
+
+        @test isequal(qsimplify(2 * exp(-im * 2π/2 * X) * exp(-im * 5π/2 * X), rewriter=qsimplify_rot), 2im*X)
+        @test isequal(qsimplify(2 * exp(-im * 2π/2 * Y) * exp(-im * 5π/2 * Y), rewriter=qsimplify_rot), 2im*Y)
+        @test isequal(qsimplify(2 * exp(-im * 2π/2 * Z) * exp(-im * 5π/2 * Z), rewriter=qsimplify_rot), 2im*Z)
+    end
+end
\ No newline at end of file