Superradiant Example: Correlation between atoms

[tnadolny]

I am trying to obtain the correlation functions between indexed transitions. My issue can be seen in the superradiant example, set up by

using QuantumCumulants
using OrdinaryDiffEq, SteadyStateDiffEq, ModelingToolkit

# Hilbertspace
hc = FockSpace(:cavity)
ha = NLevelSpace(:atom,2)
h = hc ⊗ ha

# operators
@qnumbers a::Destroy(h)
σ(α,β,i) = IndexedOperator(Transition(h, :σ, α, β),i)

@cnumbers N Δ κ Γ R ν
g(i) = IndexedVariable(:g, i)

i = Index(h,:i,N,ha)
j = Index(h,:j,N,ha)

# Hamiltonian
H = -Δ*a'a + Σ(g(i)*( a'*σ(1,2,i) + a*σ(2,1,i) ),i)

# Jump operators wth corresponding rates
J = [a, σ(1,2,i), σ(2,1,i), σ(2,2,i)]
rates = [κ, Γ, R, ν]


# Derive equations
ops = [a'*a, σ(2,2,j)]
eqs = meanfield(ops,H,J;rates=rates,order=2)

eqs_c = complete(eqs)#; filter_func=phase_invariant);
eqs_sc = scale(eqs_c);

@named sys = ODESystem(eqs_sc)

# Initial state
u0 = zeros(ComplexF64, length(eqs_sc))
# System parameters
N_ = 2e5
Γ_ = 1.0 #Γ=1mHz
Δ_ = 2500Γ_ #Δ=2.5Hz
g_ = 1000Γ_ #g=1Hz
κ_ = 5e6*Γ_ #κ=5kHz
R_ = 1000Γ_ #R=1Hz
ν_ = 1000Γ_ #ν=1Hz

ps = [N, Δ, g(1), κ, Γ, R, ν]
p0 = [N_, Δ_, g_, κ_, Γ_, R_, ν_]

prob = ODEProblem(sys,u0,(0.0, 1.0/50Γ_), ps.=>p0)

# Solve the numeric problem
sol = solve(prob,Tsit5(),maxiters=1e7);

One can now proceed to obtain the cavity spectrum via

corr = CorrelationFunction(a',a, eqs_c; steady_state=true) # works 

I obtain an error trying to evaluate the correlations between the atoms, e.g.,

corr = CorrelationFunction(σ(1,2,j), σ(2,1,j), eqs_c; steady_state=true) 
# could not calculate meanfield-equations for operator (σ12j*σ21j)
# AssertionError: isequal(ind.hilb, hilbert(op))

The error occurs in line 59 of index_meanfield.jl, when evaluating rhs_ = commutator(imH,a[i]).

I summarize the operators and their Hilbert spaces at this point:

op_ = op1_*op2_
hilbert(op_ )  is   ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0)
hilbert(op1_)  is   ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0)
hilbert(op2_)  is   ℋ(cavity) ⊗ ℋ(atom)
hilbert(op2_0) is   ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0)
a[i] is (σ12j*σ21j) # = op_
hilbert(a[i]) is ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0)
imH is (Σ(i=1:N)(im)*gi*(a′*σ12i)+Σ(i=1:N)(im)*gi*(a*σ21i)+(0 - 1im)*Δ*(a′*a))
hilbert(imH) is ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0)

I guess, since one of the terms in a[i] contains an operator with Hilbert space ℋ(cavity) ⊗ ℋ(atom), it is not possible to compute the commutator with imH. hilbert(op2_) should also be ℋ(cavity) ⊗ ℋ(atom) ⊗ ℋ(atom0), right?

Proposed fix:
I found that the line 127 in index_correlation.jl: op2_ = _new_operator(op2, h, length(h.spaces); add_subscript=add_subscript) calls the function _new_operator(t, h, aon=nothing; kwargs...) in correlation.jl, which is probably not intended.
I guess one has to alter the function function _new_operator(op::IndexedOperator,h,aon=acts_on(op)) in the beginning of index_correlation.jl to
function _new_operator(op::IndexedOperator,h,aon=acts_on(op); add_subscript=nothing) ... , allowing for the newly generated IndexedOperators to have a Symbol with a subscript (and the properly extended Hilbert space).

The same issue occurs for NumberedOperators.