Correlation function for indexed objects?

[zhouqingyi616]

Hi there,
I'm working with a large number of TLSs and am trying to construct equations using indexed operators.

The equations can be derived without problem. However I encounter the following error when trying to construct correlation function between two operators:
"could not calculate meanfield-equations for operator (σ21k*σ_012l)
ERROR: AssertionError: isequal(acts_on(op), ind.aon)"

Here's the code I'm using:

using QuantumCumulants
using ModelingToolkit
using OrdinaryDiffEq
using Plots

# Parameters
@cnumbers N Δ

# Hilbert space
ha = NLevelSpace(:atom, 2)

# Operators
σ(x, y, z) = IndexedOperator(Transition(ha, :σ, x, y), z)
Sp(i) = σ(2, 1, i)
Sm(i) = σ(1, 2, i)
S_ee(i) = σ(2, 2, i)

# Indexed variables for coupling and rates
J(i, j) = IndexedVariable(:J, i, j)  # Coupling constant
Γ(i, j) = IndexedVariable(:Γ, i, j)  # Dissipation rate

# Indices
i = Index(ha, :i, N, ha)
j = Index(ha, :j, N, ha)
k = Index(ha, :k, N, ha)
l = Index(ha, :l, N, ha)

# Hamiltonian
H_free = Σ( Δ * S_ee(i), i)  # Non-interacting Hamiltonian
H_coupling = Σ((i ≠ j) * J(i, j) * Sp(i) * Sm(j), i, j)  # Coupling term
H = H_free + H_coupling  # Total Hamiltonian

# Jump operators with rates
Jump_op = [] #[Sm(i)]
rates = [] #[Γ(i, j)]

# Derive equations
ops = [Sm(k), S_ee(k)]
eqs = meanfield(ops, H, Jump_op; rates=rates, order=2)
println(eqs)

C = CorrelationFunction(Sp(k), Sm(l), eqs; steady_state=true)

I'm not sure whether it's correct to use C = CorrelationFunction(Sp(k), Sm(l), eqs; steady_state=true) to find correlation function for abstract indexed operators?

Thanks so much for the help!

[zhouqingyi616]

Version information for the packages I'm using:
QuantumCumulants v0.3.7
QuantumOptics v1.2.1
SymbolicUtils v3.7.2
Symbolics v6.22.0