diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1851f3f..4503079 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # News
 
+## v0.4.11 - dev
+
+- Add new symbolic bosonic states and operators (`TwoSqueezeOp`, `TwoSqueezedState`, `BosonicThermalState`, `BeamSplitterOp`).
+
 ## v0.4.10 - 2025-05-11
 
 - Polish `Base.show` methods for application products and scaled quantum objects.
diff --git a/Project.toml b/Project.toml
index ba5d880..377e580 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.4.10"
+version = "0.4.11-dev"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
diff --git a/src/QSymbolicsBase/QSymbolicsBase.jl b/src/QSymbolicsBase/QSymbolicsBase.jl
index 84a886f..6f73ffc 100644
--- a/src/QSymbolicsBase/QSymbolicsBase.jl
+++ b/src/QSymbolicsBase/QSymbolicsBase.jl
@@ -43,8 +43,9 @@ export SymQObj,QObj,
        MixedState,IdentityOp,
        SApplyKet,SApplyBra,SMulOperator,SSuperOpApply,SCommutator,SAnticommutator,SBraKet,SOuterKetBra,
        HGate,XGate,YGate,ZGate,CPHASEGate,CNOTGate,
-       XBasisState,YBasisState,ZBasisState,FockState,CoherentState,SqueezedState,
+       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,
        qexpand,
diff --git a/src/QSymbolicsBase/predefined_fock.jl b/src/QSymbolicsBase/predefined_fock.jl
index b3294d5..c1b256d 100644
--- a/src/QSymbolicsBase/predefined_fock.jl
+++ b/src/QSymbolicsBase/predefined_fock.jl
@@ -2,34 +2,40 @@
 # Predefined objects in the Fock space.
 ##
 
+const inf_fock_basis = FockBasis(Inf,0.0) # Inf is Float, so the second parameter has to be Float too
+
+abstract type AbstractSingleBosonState <: SpecialKet end
+abstract type AbstractTwoBosonState <: SpecialKet end
+basis(::AbstractSingleBosonState) = inf_fock_basis
+basis(::AbstractTwoBosonState) = inf_fock_basis^2
+
 """Fock state in defined Fock basis."""
-@withmetadata struct FockState <: SpecialKet
+@withmetadata struct FockState <: AbstractSingleBosonState
     idx::Int
-    basis::FockBasis
 end
-FockState(idx::Int) = FockState(idx, inf_fock_basis)
 symbollabel(x::FockState) = "$(x.idx)"
 
 """Coherent state in defined Fock basis."""
-@withmetadata struct CoherentState <: SpecialKet
+@withmetadata struct CoherentState <: AbstractSingleBosonState
     alpha::Number # TODO parameterize
-    basis::FockBasis
 end
-CoherentState(alpha::Number) = CoherentState(alpha, inf_fock_basis)
 symbollabel(x::CoherentState) = "$(x.alpha)"
 
 """Squeezed vacuum state in defined Fock basis."""
-@withmetadata struct SqueezedState <: SpecialKet
+@withmetadata struct SqueezedState <: AbstractSingleBosonState
     z::Number
-    basis::FockBasis
 end
-SqueezedState(z::Number) = SqueezedState(z, inf_fock_basis)
 symbollabel(x::SqueezedState) = "0,$(x.z)"
 
-const inf_fock_basis = FockBasis(Inf,0.)
-"""Vacuum basis state of n"""
+"""Two-mode squeezed vacuum state, or EPR state, in defined Fock basis."""
+@withmetadata struct TwoSqueezedState <: AbstractTwoBosonState
+    z::Number
+end
+symbollabel(x::TwoSqueezedState) = "0,$(x.z)"
+
+"""Single-mode vacuum state"""
 const vac = const F₀ = const F0 = FockState(0)
-"""Single photon basis state of n"""
+"""Single photon state"""
 const F₁ = const F1 = FockState(1)
 
 ##
@@ -37,9 +43,13 @@ const F₁ = const F1 = FockState(1)
 ##
 
 abstract type AbstractSingleBosonOp <: Symbolic{AbstractOperator} end
+abstract type AbstractTwoBosonOp <: Symbolic{AbstractOperator} end
 abstract type AbstractSingleBosonGate <: AbstractSingleBosonOp end # TODO maybe an IsUnitaryTrait is a better choice
+abstract type AbstractTwoBosonGate <: AbstractTwoBosonOp end
 isexpr(::AbstractSingleBosonGate) = false
 basis(x::AbstractSingleBosonOp) = inf_fock_basis
+isexpr(::AbstractTwoBosonGate) = false
+basis(x::AbstractTwoBosonOp) = inf_fock_basis^2
 
 """Number operator.
 
@@ -54,10 +64,7 @@ julia> qsimplify(num*f, rewriter=qsimplify_fock)
 2|2⟩
 ```
 """
-@withmetadata struct NumberOp <: AbstractSingleBosonOp 
-    basis::FockBasis
-end
-NumberOp() = NumberOp(inf_fock_basis)
+@withmetadata struct NumberOp <: AbstractSingleBosonOp end
 symbollabel(::NumberOp) = "n"
 
 """Creation (raising) operator.
@@ -73,10 +80,7 @@ julia> qsimplify(create*f, rewriter=qsimplify_fock)
 (sqrt(3))|3⟩
 ```
 """
-@withmetadata struct CreateOp <: AbstractSingleBosonOp
-    basis::FockBasis
-end
-CreateOp() = CreateOp(inf_fock_basis)
+@withmetadata struct CreateOp <: AbstractSingleBosonOp end
 symbollabel(::CreateOp) = "a†"
 
 """Annihilation (lowering or destroy) operator in defined Fock basis.
@@ -92,10 +96,7 @@ julia> qsimplify(destroy*f, rewriter=qsimplify_fock)
 (sqrt(2))|1⟩
 ```
 """
-@withmetadata struct DestroyOp <: AbstractSingleBosonOp
-    basis::FockBasis
-end
-DestroyOp() = DestroyOp(inf_fock_basis)
+@withmetadata struct DestroyOp <: AbstractSingleBosonOp end
 symbollabel(::DestroyOp) = "a"
 
 """Phase-shift operator in defined Fock basis.
@@ -111,11 +112,9 @@ julia> qsimplify(phase*c, rewriter=qsimplify_fock)
 |1.2246467991473532e-16 - 1.0im⟩
 ```
 """
-@withmetadata struct PhaseShiftOp <: AbstractSingleBosonOp
+@withmetadata struct PhaseShiftOp <: AbstractSingleBosonGate
     phase::Number
-    basis::FockBasis
 end
-PhaseShiftOp(phase::Number) = PhaseShiftOp(phase, inf_fock_basis)
 symbollabel(x::PhaseShiftOp) = "U($(x.phase))"
 
 """Displacement operator in defined Fock basis.
@@ -131,11 +130,9 @@ julia> qsimplify(displace*f, rewriter=qsimplify_fock)
 |im⟩
 ```
 """
-@withmetadata struct DisplaceOp <: AbstractSingleBosonOp
+@withmetadata struct DisplaceOp <: AbstractSingleBosonGate
     alpha::Number
-    basis::FockBasis
 end
-DisplaceOp(alpha::Number) = DisplaceOp(alpha, inf_fock_basis)
 symbollabel(x::DisplaceOp) = "D($(x.alpha))"
 
 """Number operator, also available as the constant `n̂`, in an infinite dimension Fock basis."""
@@ -156,9 +153,25 @@ julia> qsimplify(S*vac, rewriter=qsimplify_fock)
 |0,π⟩
 ```
 """
-@withmetadata struct SqueezeOp <: AbstractSingleBosonOp
+@withmetadata struct SqueezeOp <: AbstractSingleBosonGate
     z::Number
-    basis::FockBasis
 end
-SqueezeOp(z::Number) = SqueezeOp(z, inf_fock_basis)
-symbollabel(x::SqueezeOp) = "S($(x.z))"
\ No newline at end of file
+symbollabel(x::SqueezeOp) = "S($(x.z))"
+
+"""Thermal bosonic state in defined Fock basis."""
+@withmetadata struct BosonicThermalState <: AbstractSingleBosonOp
+    photons::Number
+end
+symbollabel(x::BosonicThermalState) = "ρₜₕ($(x.photons))"
+
+"""Two-mode squeezing operator in defined Fock basis."""
+@withmetadata struct TwoSqueezeOp <: AbstractTwoBosonGate
+    z::Number
+end
+symbollabel(x::TwoSqueezeOp) = "S₂($(x.z))"
+
+"""Two-mode beamsplitter operator in defined Fock basis."""
+@withmetadata struct BeamSplitterOp <: AbstractTwoBosonGate
+    transmit::Number
+end
+symbollabel(x::BeamSplitterOp) = "B($(x.transmit))"
\ No newline at end of file
diff --git a/src/QSymbolicsBase/rules.jl b/src/QSymbolicsBase/rules.jl
index 9ae6fda..1664ffe 100644
--- a/src/QSymbolicsBase/rules.jl
+++ b/src/QSymbolicsBase/rules.jl
@@ -11,7 +11,6 @@ function hasscalings(xs)
     end
 end
 _isa(T) = x->isa(x,T)
-_isequal(obj) = x->(x==obj)
 _vecisa(T) = x->all(_isa(T), x)
 
 ##
@@ -33,36 +32,36 @@ RULES_PAULI = [
     @rule(~o1::_isa(HGate)*~o2::_isa(YGate)*~o3::_isa(HGate) => -Y),
     @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(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(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(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(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(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(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)
+    @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
@@ -86,8 +85,8 @@ RULES_ANTICOMMUTATOR = [
     @rule(anticommutator(~o1::_isa(YGate), ~o2::_isa(XGate)) => 0),
     @rule(anticommutator(~o1::_isa(ZGate), ~o2::_isa(YGate)) => 0),
     @rule(anticommutator(~o1::_isa(XGate), ~o2::_isa(ZGate)) => 0),
-    @rule(commutator(~o1::_isa(DestroyOp), ~o2::_isa(CreateOp)) => IdentityOp((~o1).basis)),
-    @rule(commutator(~o1::_isa(CreateOp), ~o2::_isa(DestroyOp)) => -IdentityOp((~o1).basis)),
+    @rule(commutator(~o1::_isa(DestroyOp), ~o2::_isa(CreateOp)) => IdentityOp(basis(~o1))),
+    @rule(commutator(~o1::_isa(CreateOp), ~o2::_isa(DestroyOp)) => -IdentityOp(basis(~o1))),
     @rule(commutator(~o1::_isa(NumberOp), ~o2::_isa(DestroyOp)) => -(~o2)),
     @rule(commutator(~o1::_isa(DestroyOp), ~o2::_isa(NumberOp)) => (~o1)),
     @rule(commutator(~o1::_isa(NumberOp), ~o2::_isa(CreateOp)) => (~o2)),
@@ -95,7 +94,7 @@ RULES_ANTICOMMUTATOR = [
 ]
 
 RULES_FOCK = [
-    @rule(~o::_isa(DestroyOp) * ~k::_isequal(vac) => SZeroKet()),
+    @rule(~o::_isa(DestroyOp) * ~k::isequal(vac) => SZeroKet()),
     @rule(~o::_isa(CreateOp) * ~k::_isa(FockState) => Term(sqrt,[(~k).idx+1])*FockState((~k).idx+1)),
     @rule(~o::_isa(DestroyOp) * ~k::_isa(FockState) => Term(sqrt,[(~k).idx])*FockState((~k).idx-1)),
     @rule(~o::_isa(NumberOp) * ~k::_isa(FockState) => (~k).idx*(~k)),
@@ -103,10 +102,11 @@ RULES_FOCK = [
     @rule(~o::_isa(PhaseShiftOp) * ~k::_isa(CoherentState) => CoherentState((~k).alpha * exp(-im*(~o).phase))),
     @rule(dagger(~o1::_isa(PhaseShiftOp)) * ~o2::_isa(DestroyOp) * ~o1 => ~o2*exp(-im*((~o1).phase))),
     @rule(~o1::_isa(PhaseShiftOp) * ~o2::_isa(DestroyOp) * dagger(~o1) => ~o2*exp(im*((~o1).phase))),
-    @rule(dagger(~o1::_isa(DisplaceOp)) * ~o2::_isa(DestroyOp) * ~o1 => (~o2) + (~o1).alpha*IdentityOp((~o2).basis)),
-    @rule(dagger(~o1::_isa(DisplaceOp)) * ~o2::_isa(CreateOp) * ~o1 => (~o2) + conj((~o1).alpha)*IdentityOp((~o2).basis)),
+    @rule(dagger(~o1::_isa(DisplaceOp)) * ~o2::_isa(DestroyOp) * ~o1 => (~o2) + (~o1).alpha*IdentityOp(basis(~o2))),
+    @rule(dagger(~o1::_isa(DisplaceOp)) * ~o2::_isa(CreateOp) * ~o1 => (~o2) + conj((~o1).alpha)*IdentityOp(basis(~o2))),
     @rule(~o::_isa(DisplaceOp) * ~k::((x->(isa(x,FockState) && x.idx == 0))) => CoherentState((~o).alpha)),
-    @rule(~o::_isa(SqueezeOp) * ~k::_isequal(vac) => SqueezedState((~o).z, (~o).basis))
+    @rule(~o::_isa(SqueezeOp) * ~k::isequal(vac) => SqueezedState((~o).z)),
+    @rule(~o::_isa(TwoSqueezeOp) * ~k::isequal(vac ⊗ vac) => TwoSqueezedState((~o).z))
 ]
 
 RULES_SIMPLIFY = [RULES_PAULI; RULES_COMMUTATOR; RULES_ANTICOMMUTATOR; RULES_FOCK]
diff --git a/test/test_express_opt.jl b/test/test_express_opt.jl
index 287ff2a..e4e8f91 100644
--- a/test/test_express_opt.jl
+++ b/test/test_express_opt.jl
@@ -36,9 +36,9 @@
 
     state = (3im*(2*dagger(Z1)+dagger(Y1))) * (3im*(2*X1+X2))
 
-    cstate = CoherentState(im, inf_fock_basis)
-    displace = DisplaceOp(im,inf_fock_basis)
-    phase = PhaseShiftOp(im, inf_fock_basis)
+    cstate = CoherentState(im)
+    displace = DisplaceOp(im)
+    phase = PhaseShiftOp(im)
     @test express(N*F1) ≈ express(N)*express(F1)
     @test express(Create*F1) ≈ express(Create)*express(F1)
     @test express(Destroy*F1) ≈ express(Destroy)*express(F1)
diff --git a/test/test_fock.jl b/test/test_fock.jl
index 33a7495..d4db1ef 100644
--- a/test/test_fock.jl
+++ b/test/test_fock.jl
@@ -7,7 +7,9 @@
     phase2 = PhaseShiftOp(pi)
     displace = DisplaceOp(im)
     squeezeop = SqueezeOp(pi)
+    twosqueezeop = TwoSqueezeOp(pi)
     sstate = SqueezedState(pi)
+    tsstate = TwoSqueezedState(pi)
 
     @testset "ladder and number operators" begin
         @test isequal(qsimplify(Destroy*vac, rewriter=qsimplify_fock), SZeroKet())
@@ -28,5 +30,6 @@
 
     @testset "Squeeze operators" begin
         @test isequal(qsimplify(squeezeop*vac, rewriter=qsimplify_fock), sstate)
+        @test isequal(qsimplify(twosqueezeop*(vac ⊗ vac), rewriter=qsimplify_fock), tsstate)
     end
 end
\ No newline at end of file