⟨σ_112⟩(t) exceed 1 when setting order=3

[Papageno2]

I am attempting to solve the Tavis-Cummings model using QuantumCumulants.jl@v0.2.25.
When I use order=2, everything appears reasonable.
However, when I set order=3, 4, or 5, I obtain results that seem unphysical.

Is it acceptable for ⟨σ_112⟩(t) to exceed 1? (I believe there might be an error.)
The results appear significantly different from those obtained with order=2.

using QuantumCumulants, OrdinaryDiffEq
using ModelingToolkit, DelimitedFiles
using Symbolics
using Plots
using CSV, DataFrames

@cnumbers   ω_a κ_a g_a ω_p Ω_p ω21 Δ_a Δ_c N Γ21;
@syms t::Real

hc = FockSpace(:cav1); 
ha_ = NLevelSpace(:atom, 2);
ha = ClusterSpace(ha_, N, 2);
h = hc⊗ha

@qnumbers a::Destroy(h,1)
σ(i,j) = Transition(h, :σ, i, j, 2)

Δ_a = ω21-ω_a;

H = -Δ_a*(a'*a)+g_a*(a'*sum(σ(1,2))+a*sum(σ(2,1)));

J = [a, σ(1,2)];
rates = [κ_a, Γ21];
ops = [a'*a, σ(1,2)[1], σ(2,2)[1], σ(1,2)[1]*σ(1,2)[2], σ(2,2)[1]*σ(2,2)[2]];

eqs = meanfield(ops, H, J, rates=rates, order=3, iv=t);
eqs_c = complete(eqs)

u0=zeros(ComplexF64, length(eqs_c));

u0[2]=0.50; #σ12
u0[3]=0.50; # σ22
u0[4]=u0[2]*u0[2]  #σ12 σ12
u0[5]=u0[3]*u0[3]  # σ22 σ22
u0[10]=u0[3]*u0[2]  #σ22 σ12
u0[11]=conj(u0[2])*u0[2]  #σ21 σ12

@named  sys=ODESystem(eqs_c);
tmax = 10e0
prob0=ODEProblem(sys, u0, (0.0,tmax), (ps.=>ps_));
sol0=solve(prob0, RK4(), saveat=100e-6, maxiters=2e10, reltol=1e-12, abstol=1e-12);

df=DataFrame(sol0)

The results for order=2:

The results for order=3 and order=4 are the same :

[Papageno2]

For collective spin calculations, if either $Re(\sigma_1^{12})$ or $Im(\sigma_1^{12})$ exceeds 0.5, it implies that the projection of the collective spin exceeds $N/2$ for a system of $N$ spin-1/2 atoms.
$$
\begin{align}
\langle \hat{J}_x\rangle &= (N/2)(\langle \hat{\sigma}_1^{12}\rangle+\langle \hat{\sigma}_1^{21}\rangle),\
\langle \hat{J}_y\rangle &= (iN/2)(\langle \hat{\sigma}_1^{12}\rangle - \langle \hat{\sigma}_1^{21}\rangle),\
\langle \hat{J}_z\rangle &= (N/2)(2\langle \hat{\sigma}_1^{22}\rangle - 1).
\end{align}$$

[ChristophHotter]

Hi @Papageno2,
You used the ClusterSpace which is outdated. Please use the Symbolic Sums and Indices instead. This should resolve your problem.
The superradiant laser example might be useful for you: https://qojulia.github.io/QuantumCumulants.jl/stable/examples/superradiant_laser_indexed/
Here's some documentation about the symbolic sums: https://qojulia.github.io/QuantumCumulants.jl/stable/symbolic_sums/

[Papageno2]

Hi @Papageno2, You used the ClusterSpace which is outdated. Please use the Symbolic Sums and Indices instead. This should resolve your problem. The superradiant laser example might be useful for you: https://qojulia.github.io/QuantumCumulants.jl/stable/examples/superradiant_laser_indexed/ Here's some documentation about the symbolic sums: https://qojulia.github.io/QuantumCumulants.jl/stable/symbolic_sums/

Thank you for reminding me about the new API: Sumbolic Sums and Indices, now I get the same results for order=2 and order=3. The inconsistency in results that I previously encountered may be attributed to the instability of the ODE due to the sizable detuned parameter $\Delta$.