CorrelationFunction with steady_state = false

[tnadolny]

Hi,
I am trying to work with correlation functions, when no steady state reached. My issue can be seen using the Superradiant Laser example. For completeness, here is the basic setup, copied from the example (without the custom filter function):

using QuantumCumulants
using OrdinaryDiffEq, SteadyStateDiffEq, ModelingToolkit
using Plots

# 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)

# Plot time evolution
t = sol.t
n = real.(sol[a'a])
s22 = real.(sol[σ(2,2,1)])
# Plot
p1 = plot(t, n, xlabel="tΓ", ylabel="⟨a⁺a⟩", legend=false)
p2 = plot(t, s22, xlabel="tΓ", ylabel="⟨σ22⟩", legend=false)
plot(p1, p2, layout=(1,2), size=(700,300))

The problem comes when trying to solve the correlation equations:

# corr = CorrelationFunction(a', a, eqs_c; steady_state=true) # works as expected
corr = CorrelationFunction(a', a, eqs_c; steady_state=false) # with steady_state=false, we encounter an error below
corr_sc = scale(corr)

p0_c = correlation_p0(corr_sc, sol.u[end],ps.=>p0)
u0_c = correlation_u0(corr_sc, sol.u[end])

@named csys = ODESystem(corr_sc) # 5 equations, 8 parameters

println(p0_c)
# (N => 200000.0, Δ => 2500.0, g_1 => 1000.0, κ => 5.0e6, Γ => 1.0, R => 1000.0, ν => 1000.0, var"⟨a⟩" => 0.0 + 0.0im)
println(csys.eqs[2])
# Differential(τ)(var"⟨σ211*a_0⟩"(τ)) ~ (0 + 1im)*g_1*var"⟨a′*a_0⟩"(τ) + (0 - 2im)*g_1*(⟨σ221⟩*var"⟨a′*a_0⟩"(τ) + var"⟨σ221*a_0⟩"(τ)*⟨var"⟨a⟩"⟩ + var"⟨a⟩"*⟨a′*σ221⟩ - 2var"⟨a⟩"*⟨σ221⟩*⟨var"⟨a⟩"⟩) - 0.5R*var"⟨σ211*a_0⟩"(τ) - 0.5Γ*var"⟨σ211*a_0⟩"(τ) - 0.5ν*var"⟨σ211*a_0⟩"(τ)

prob_c = ODEProblem(csys,u0_c,(0.0,1.0),p0_c)

sol_c = solve(prob_c,Tsit5(),save_idxs=1); #UndefVarError: avg not defined

To me, it seems to be an issue that csys contains objects like ⟨σ221⟩ that are neither states of the ODE, nor parameters which are replaced by p0_c. Is there something that I am missing?

[ChristophHotter]

@tnadolny thank you for reporting this. @j-moser fixed this in #158. It should work now in version 0.2.15.

[tnadolny]

Thank you!
The new test works well for me. However, I also tried the test without the filter_function, expecting it to work as well.
I get two different errors:
1st error:

u0_c = correlation_u0(corr_sc, sol.u[end])
# Could not find initial value for ⟨σ211*σ222⟩ !

2nd error:

@named csys = ODESystem(corr_sc)
# MethodError: no method matching has_cluster(::Sym{Parameter, Base.ImmutableDict{DataType, Any}})

Regarding the first error, I checked, that ⟨σ211*σ222⟩ is not in the ODESystem:

for s in sys.states
    println(s)
end

var"⟨a′*a⟩"(t)
var"⟨σ221⟩"(t)
var"⟨a′*σ121⟩"(t)
var"⟨a′⟩"(t)
var"⟨a*σ221⟩"(t)
var"⟨σ211*σ122⟩"(t)
var"⟨σ211⟩"(t)
var"⟨a*a⟩"(t)
var"⟨a*σ121⟩"(t)
var"⟨σ221*σ122⟩"(t)
var"⟨σ121*σ122⟩"(t)
var"⟨σ221*σ222⟩"(t)

Is the missing state ⟨σ211*σ222⟩ the same as the conjugate of the state var"⟨σ221*σ122⟩"(t)? If so, could one extend correlation_u0() to check for such permutations?