WIP: Symmetric ordering

src/QuantumCumulants.jl

@@ -24,6 +24,7 @@ export HilbertSpace, ProductSpace, ⊗, tensor,
         find_missing, complete, complete!, find_operators, fundamental_operators,
             unique_ops, unique_ops!,
         CorrelationFunction, Spectrum, correlation_u0, correlation_p0,
+        Symmetrized, symmetrize,
         ClusterSpace,
         scale,
         transition_superscript
@@ -39,6 +40,7 @@ include("average.jl")
 include("utils.jl")
 include("diffeq.jl")
 include("correlation.jl")
+include("symmetrize.jl")
 include("cluster.jl")
 include("scale.jl")
 include("latexify_recipes.jl")

src/printing.jl

@@ -14,6 +14,11 @@ 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::Transition) = write(io, Symbol(x.name,x.i,x.j))
+function Base.show(io::IO,x::Symmetrized)
+    write(io, "S(")
+    show(io, x.operator)
+    write(io, ")")
+end
 
 show_brackets = Ref(true)
 function Base.show(io::IO,x::QTerm)

src/symmetrize.jl

@@ -0,0 +1,71 @@
+"""
+    Symmetrized <: QSym
+    Symmetrized(op::QNumber)
+
+A [`QSym`](@ref) representing the symmetrized expression of an operator given by
+`S(x) = 0.5*(x + dagger(x))`, where `S` is the symmetrization operator.
+"""
+struct Symmetrized{T} <: QSym
+    operator::T
+    function Symmetrized(operator::T) where T<:QNumber
+        new{T}(operator)
+    end
+end
+Symmetrized(x::Average) = _average(Symmetrized(x.arguments[1]))
+Base.adjoint(s::Symmetrized) = s
+Base.isequal(s1::Symmetrized, s2::Symmetrized) = isequal(s1.operator, s2.operator)
+
+for f ∈ [:hilbert, :acts_on]
+    @eval $(f)(s::Symmetrized) = $(f)(s.operator)
+end
+
+function Base.getproperty(s::Symmetrized, field::Symbol)
+    if field === :name
+        _get_name(s.operator)
+    else
+        return getfield(s, field)
+    end
+end
+
+# Generate a name when necessary
+_get_name(x::QSym) = x.name
+_get_name(x::QMul) = Symbol(map(_get_name, x.args_nc)...)
+
+
+"""
+    symmetrize(eqs)
+
+Compute the equations that follow symmetric ordering of operators out of a set
+of equations that uses normal ordering. This is done by taking the equation of
+an operator `x` adding `adjoint(x)` to it and dividing by 2. Operators are wrapped
+as [`Symmetrized`](@ref) to avoid repeated application of normal-ordered commutation
+relations.
+"""
+function symmetrize(eqs::MeanfieldEquations)
+    lhs = eqs.states
+    rhs = getfield.(eqs.equations, :rhs)
+
+    # Substitute normal ordering by symmetric one
+    lhs_sym = Symmetrized.(lhs)
+    subs = Dict(lhs .=> lhs_sym)
+    rhs_sym = [substitute(r, subs) for r ∈ rhs]
+
+    # Compute correct form of equations by adding the adjoint
+    rhs_sym = [Symbolics.simplify(0.5*r + 0.5*_conj(r)) for r ∈ rhs_sym]
+
+    eqs_sym = lhs_sym .~ rhs_sym
+
+    varmap = make_varmap(lhs_sym, eqs.iv)
+
+    return MeanfieldEquations(eqs_sym,
+                                eqs.operator_equations,
+                                lhs_sym,
+                                eqs.operators,
+                                eqs.hamiltonian,
+                                eqs.jumps,
+                                eqs.rates,
+                                eqs.iv,
+                                varmap,
+                                eqs.order
+                                )
+end

test/runtests.jl

@@ -9,6 +9,7 @@ names = [
     "test_mixed-order.jl"
     "test_correlation.jl"
     "test_two-level-laser.jl"
+    "test_symmetrize.jl"
     "test_cluster.jl"
     "test_scaling.jl"
     "test_higher-order.jl"

test/test_symmetrize.jl

@@ -0,0 +1,41 @@
+using QuantumCumulants
+using ModelingToolkit, OrdinaryDiffEq
+using Test
+
+@testset "symmetrize" begin
+
+# Implement example using x and p
+struct Position <: QSym
+    hilbert
+    name
+    aon
+end
+
+struct Momentum <: QSym
+    hilbert
+    name
+    aon
+end
+
+for T ∈ [:Position,:Momentum]
+    @eval Base.adjoint(op::($T)) = op
+end
+
+QuantumCumulants.ismergeable(::Position,::Momentum) = true
+Base.:*(x::Position,p::Momentum) = p*x + im
+
+h = FockSpace(:oscillator)
+x = Position(h,:x,1)
+p = Momentum(h,:p,1)
+
+@cnumbers ω m
+H = p^2/(2m) + 0.5m*ω^2*x^2
+
+eqs = meanfield([x,p,x^2,p^2,p*x],H)
+
+s_eqs = symmetrize(eqs)
+
+@test all(st.arguments[1] isa Symmetrized for st ∈ s_eqs.states)
+@test isempty(find_missing(s_eqs))
+
+end # testset