From e9cb961b4d41458a03edae7f22a2c33c11ee24ab Mon Sep 17 00:00:00 2001
From: Benedikt Tissot <benedikt@tissot.fun>
Date: Wed, 11 Dec 2024 10:57:42 +0100
Subject: [PATCH] add GroundStateProjection operator for FockSpace

---
 docs/src/api.md         |  6 +++-
 src/QuantumCumulants.jl |  2 +-
 src/fock.jl             | 67 ++++++++++++++++++++++++++++++++++++-----
 src/latexify_recipes.jl |  4 +++
 src/printing.jl         |  1 +
 test/test_fock.jl       | 12 ++++++++
 6 files changed, 83 insertions(+), 9 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index f185b56d..545c72c3 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -68,6 +68,10 @@ Destroy
 Create
 ```
 
+```@docs
+GroundStateProjection
+```
+
 ```@docs
 Transition
 ```
@@ -330,4 +334,4 @@ MeanfieldNoiseEquations
 
 ```@docs
 IndexedMeanfieldNoiseEquations
-```
\ No newline at end of file
+```
diff --git a/src/QuantumCumulants.jl b/src/QuantumCumulants.jl
index a117c6d2..7b79c1e6 100644
--- a/src/QuantumCumulants.jl
+++ b/src/QuantumCumulants.jl
@@ -20,7 +20,7 @@ import QuantumOpticsBase: ⊗, tensor
 
 export HilbertSpace, ProductSpace, ⊗, tensor,
         QSym, QTerm, @qnumbers,
-        FockSpace, Destroy, Create,
+        FockSpace, Destroy, Create, GroundStateProjection,
         NLevelSpace, Transition,
         PauliSpace, Pauli, SpinSpace, Spin,
         MeanfieldEquations,
diff --git a/src/fock.jl b/src/fock.jl
index 794c3c90..c8da03ff 100644
--- a/src/fock.jl
+++ b/src/fock.jl
@@ -43,11 +43,31 @@ struct Create{H<:HilbertSpace,S,A,M} <: QSym
     end
 end
 
-for T in (:Create, :Destroy)
+"""
+    GroundStateProjection <: QSym
+
+Bosonic operator on a [`FockSpace`](@ref) representing the quantum harmonic
+oscillator ground state projection.
+
+Warning: The ground state projection is not meant to be used in the cumulant
+ expansion.
+"""
+struct GroundStateProjection{H<:HilbertSpace,S,A,M} <: QSym
+    hilbert::H
+    name::S
+    aon::A
+    metadata::M
+    function GroundStateProjection{H,S,A,M}(hilbert::H, name::S, aon::A, metadata::M) where {H,S,A,M}
+        @assert has_hilbert(FockSpace,hilbert,aon)
+        new(hilbert,name,aon,metadata)
+    end
+end
+
+for T in (:Create, :Destroy, :GroundStateProjection)
     @eval Base.isequal(a::$T, b::$T) = isequal(a.hilbert, b.hilbert) && isequal(a.name, b.name) && isequal(a.aon, b.aon)
 end
 
-for f in [:Destroy,:Create]
+for f in [:Destroy,:Create, :GroundStateProjection]
     @eval $(f)(hilbert::H, name::S, aon::A; metadata::M=NO_METADATA) where {H,S,A,M} = $(f){H,S,A,M}(hilbert,name,aon,metadata)
     @eval $(f)(hilbert::FockSpace, name; metadata=NO_METADATA) = $(f)(hilbert,name,1; metadata)
     @eval function $(f)(hilbert::H, name::S, aon::A; metadata::M=NO_METADATA) where {H<:ProductSpace,S,A<:Int,M}
@@ -75,6 +95,7 @@ end
 
 Base.adjoint(op::Destroy) = Create(op.hilbert,op.name,acts_on(op); op.metadata)
 Base.adjoint(op::Create) = Destroy(op.hilbert,op.name,acts_on(op); op.metadata)
+Base.adjoint(op::GroundStateProjection) = op
 
 # Commutation relation in simplification
 function *(a::Destroy,b::Create)
@@ -89,11 +110,43 @@ function *(a::Destroy,b::Create)
         return QMul(1, [b,a])
     end
 end
-# ismergeable(::Destroy,::Create) = true
 
-# TODO: test if faster; delete if and elseif in *-function above?
-function ismergeable(a::Destroy,b::Create)
+function *(a::GroundStateProjection,b::GroundStateProjection)
+    check_hilbert(a,b)
     aon_a = acts_on(a)
     aon_b = acts_on(b)
-    return aon_a == aon_b
-end
\ No newline at end of file
+    if aon_a == aon_b
+        return a
+    elseif aon_a < aon_b
+        return QMul(1, [a,b])
+    else
+        return QMul(1, [b,a])
+    end
+end
+
+for (T1,T2) in ((:Destroy, :GroundStateProjection), (:GroundStateProjection, :Create))
+    @eval function *(a::$T1,b::$T2)
+        check_hilbert(a,b)
+        aon_a = acts_on(a)
+        aon_b = acts_on(b)
+        if aon_a == aon_b
+            return 0
+        elseif aon_a < aon_b
+            return QMul(1, [a,b])
+        else
+            return QMul(1, [b,a])
+        end
+    end
+end
+
+for T1 in (:Destroy, :GroundStateProjection)
+    for T2 in (:Create, :GroundStateProjection)
+        # ismergeable(::$T1,::$T2) = true
+        # TODO: test if faster; delete if and elseif in *-function above?
+        @eval function ismergeable(a::$T1,b::$T2)
+            aon_a = acts_on(a)
+            aon_b = acts_on(b)
+            return aon_a == aon_b
+        end
+    end
+end
diff --git a/src/latexify_recipes.jl b/src/latexify_recipes.jl
index 2464fa8e..f3c46ba5 100644
--- a/src/latexify_recipes.jl
+++ b/src/latexify_recipes.jl
@@ -28,6 +28,9 @@ function _postwalk_func(x)
     elseif MacroTools.@capture(x, dagger(arg_))
         s = "$(arg)^\\dagger"
         return s
+    elseif MacroTools.@capture(x, GroundStateProjection(arg_))
+        s = "\\vert \\emptyset_$(arg) \\rangle \\langle \\emptyset_$(arg) \\vert"
+        return s
     elseif MacroTools.@capture(x, Transition(arg_,i_,j_))
         s = "{$(arg)}$(transition_idx_script[]){{$(i)$(j)}}"
         return s
@@ -202,6 +205,7 @@ function _to_expression(x::Complex{Symbolics.Num}) # forward complex Nums to Sym
 end
 _to_expression(op::QSym) = op.name
 _to_expression(op::Create) = :(dagger($(op.name)))
+_to_expression(op::GroundStateProjection) = :(GroundStateProjection($(op.name)))
 _to_expression(op::Transition) = :(Transition($(op.name),$(op.i),$(op.j)) )
 xyz_sym=[:x,:y,:z]
 _to_expression(op::Pauli) = :(Pauli($(op.name),$(xyz_sym[op.axis])))
diff --git a/src/printing.jl b/src/printing.jl
index c14cf712..7d64b050 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -13,6 +13,7 @@ end
 
 Base.show(io::IO,x::QSym) = write(io, x.name)
 Base.show(io::IO,x::Create) = write(io, string(x.name, "′"))
+Base.show(io::IO,x::GroundStateProjection) = write(io, string("∅_", x.name))
 Base.show(io::IO,x::Transition) = write(io, Symbol(x.name,x.i,x.j))
 Base.show(io::IO,x::Pauli) = write(io, Symbol(x.name,xyz_sym[x.axis]))
 Base.show(io::IO,x::Spin) = write(io, Symbol(x.name,xyz_sym[x.axis]))
diff --git a/test/test_fock.jl b/test/test_fock.jl
index b0c2a7a3..ac29dc17 100644
--- a/test/test_fock.jl
+++ b/test/test_fock.jl
@@ -8,20 +8,26 @@ using Test
 hf = FockSpace(:c)
 
 a = Destroy(hf,:a)
+P0 = GroundStateProjection(hf, :a)
 ad = a'
 
 @test !isequal(hash(a), hash(ad))
 b = Destroy(hf,:b)
 @test !isequal(hash(a), hash(b))
+@test !isequal(hash(a), hash(P0))
 
 @test isequal(a,ad')
 @test isequal(simplify(a+a),2*a)
 @test isequal(simplify(a/2 + 0.5*a),a)
 @test isequal(a*a' , 1+a'*a)
 @test isequal(simplify(a*a' + 1) , 2 + a'*a)
+@test isequal(simplify(P0 + P0) , 2*P0)
+@test isequal(simplify(P0^2) , P0)
+@test isequal(simplify(a'*a*(a')^3*P0^2) , 3*(a')^3*P0)
 
 @test isequal(simplify(-1*(a'*a + 1)*a + a) , -1*a'*a^2)
 @test isequal(simplify(a'*a*a - a*a'*a) , -1*a)
+@test isequal(simplify(P0 - P0) , 0)
 
 # Single mode
 ωc = 0.1313
@@ -34,11 +40,17 @@ da = simplify(1.0im*(H*a - a*H))
 @test iszero(substitute(x*a, Dict(x=>0)))
 @test isequal(substitute(x*a, Dict(x=>1)), a)
 @test iszero(substitute(x*a, Dict(a=>0)))
+@test iszero(substitute(x*P0, Dict(x=>0)))
+@test isequal(substitute(x*P0, Dict(x=>1)), P0)
+@test iszero(substitute(x*P0, Dict(P0=>0)))
 
 # Test substitute by syms
 @syms y
 @test isequal(substitute(x*a, Dict(x=>y)), y*a)
 @test isequal(substitute(x*a, Dict(a=>y)), x*y)
 @test isequal(substitute(x*(a+a'), Dict(x => y)), y*(a + a'))
+@test isequal(substitute(x*P0, Dict(x=>y)), y*P0)
+@test isequal(substitute(x*P0, Dict(P0=>y)), x*y)
+@test isequal(substitute(x*(P0+a'*P0), Dict(x => y)), y*(P0 + a'*P0))
 
 end # testset