qojulia / quantumcumulants.jl Goto Github PK

Generalized mean-field equations in open quantum systems

License: Other

Julia 100.00%

julia quantum quantum-mechanics quantum-optics cumulants

quantumcumulants.jl's Introduction

QuantumCumulants.jl

QuantumCumulants.jl is a package for the symbolic derivation of mean-field equations for quantum mechanical operators in Julia. The equations are derived using fundamental commutation relations of operators. When averaging these equations they can be automatically expanded in terms of cumulants to an arbitrary order (generalized mean-field approximation). This results in a closed set of symbolic differential equations, which can also be solved numerically.

For the application of commutation relations QuantumCumulants.jl implements a simple noncommutative algebra, where any commutation relations are applied immediately. All other symbolic simplification and rewriting is done using the Symbolics.jl package.

To obtain a numerical solution, equations derived with QuantumCumulants.jl can be converted to ModelingToolkit.jl and subsequently solved with DifferentialEquations.jl.

Development status

Note that QuantumCumulants.jl is still at an early stage of development.

Installation

The package can be installed with

|pkg> add QuantumCumulants

Documentation

Please refer to the latest Documentation for more details and examples.

Short example

To briefly illustrate how QuantumCumulants.jl works, here's how you can implement a first-order mean-field model of a laser with a single atom as a gain medium:

using QuantumCumulants

h_cav = FockSpace(:cavity)
h_atom = NLevelSpace(:atom, (:g,:e))
h = tensor(h_cav, h_atom)

@cnumbers Δ g κ γ ν
@qnumbers a::Destroy(h) σ::Transition(h)

H = Δ*a'*a + g*(a'*σ(:g,:e) + a*σ(:e,:g))
J = [a,σ(:g,:e),σ(:e,:g)]
rates = [κ,γ,ν]

eqs = meanfield([a,σ(:g,:e),σ(:e,:e)], H, J; rates=rates, order=1)

using ModelingToolkit, OrdinaryDiffEq
@named sys = ODESystem(eqs)
p0 = (Δ=>0, g=>1.5, κ=>1, γ=>0.25, ν=>4)
u0 = ComplexF64[1e-2, 0, 0]
prob = ODEProblem(sys,u0,(0.0,50.0),p0)
sol = solve(prob,RK4())

using Plots
n = abs2.(sol[a])
plot(sol.t, n, xlabel="t", label="n")

The above code implements the Jaynes-Cummings Hamiltonian describing an optical cavity mode that couples to a two-level atom. Additionally, the decay processes are specified. Then, mean-field equations for the average values of the operators [a,σ(:g,:e),σ(:e,:e)] are derived and expanded to first order (average values of products are factorized). For the numerical solution an ODESystem (from ModelingToolkit.jl) is created and solved with the OrdinaryDiffEq.jl library. Finally, the time dynamics of the photon number n is plotted.

Citing

If you find QuantumCumulants.jl useful in your research, please consider citing this paper:

@article{plankensteiner2022quantumcumulants,
  doi = {10.22331/q-2022-01-04-617},
  url = {https://doi.org/10.22331/q-2022-01-04-617},
  title = {Quantum{C}umulants.jl: {A} {J}ulia framework for generalized mean-field equations in open quantum systems},
  author = {Plankensteiner, David and Hotter, Christoph and Ritsch, Helmut},
  journal = {{Quantum}},
  issn = {2521-327X},
  publisher = {{Verein zur F{\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
  volume = {6},
  pages = {617},
  month = jan,
  year = {2022}
}

quantumcumulants.jl's People

Contributors

Stargazers

Watchers

Forkers

christophhotter ritschh kllrak benneti j-moser tureci-rc jkosata tnadolny dywang12 nico0oblr karolpezet lixingli1994 oameye

quantumcumulants.jl's Issues

metadata function for QAdd

The package defines a function to extract the metadata of the QAdd type but QAdd does not have a metadata field.

QuantumCumulants.jl/src/qnumber.jl

Line 232 in 1d144b1

SymbolicUtils.metadata(q::QAdd) = q.metadata

The function is used when you try to use Symbolics.expand on a QAdd expression:

using QuantumCumulant
using Symbolics

h = FockSpace(:cavity)

@qnumbers a::Destroy(h)

expr = (a'+a) |> simplify
expr |> expand

How to represent the number $N$ of atom in an atomic ensemble with collective spin operators?

Model

I want to study the dynamic evolution of an atomic ensemble with $N$ atoms, and the master equation of the system is:
Screenshot from 2024-04-20 11-40-24
where the $\hat{H}$ is:

So I try to solve it with the QuantumCumulants v0.3.1.

Code

using QuantumCumulants, OrdinaryDiffEq
using ModelingToolkit, DelimitedFiles
using Symbolics

@cnumbers   λ N ξ η;

h = SpinSpace("spinspace")
Jx = Spin(h,"J1",1)
Jy = Spin(h,"J1",2)
Jz = Spin(h,"J1",3)
Jm = Jx - 1im*Jy
Jp = Jx + 1im*Jy

H = -ξ*λ^2/N*(N^2/4-Jz^2)
J = [Jm]
rates = [2*η*λ^2/N]
ops = [Jx]
eqs = meanfield(ops,H,J,rates = rates,order=2)

And the outcome is:

How to simplify the number $N$ of atom.

Adjoint jump operators get ignored

We're being inconsistent with the jump operators. While we allow explicitly specifying Jdagger in meanfield, this doesn't get stored in the MeanFieldEquations type. This means that if there is a case where Jdagger != adjoint.(J), you'd get different results depending on whether you specify a list of operators and Jdagger in meanfield or if you use complete. The latter will just infer Jdagger=adjoint.(J).

We should either store Jdagger and use it where appropriate, or disallow setting Jdagger altogether.

Superradiant Example: Correlation between atoms

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.

Including noise terms in the symbolic derivation

Hello, I was having a look at your interesting paper and code and I was wondering if it would be difficult to include the noise terms of the quantum Langevin equations into the derivation? It would be interesting since the inclusion of these terms in the (stochastic) evolution have been shown to yield good results for the description of the dynamics/ steady states of several models.

It seems to me that all building blocks should be present here, i.e. the symbolic derivation, transformation to an ODESystem and then solving it with an ODEProblem. If the symbolic derivation could output a new ODESystem describing the stochastic noise contribution to the dynamics then I suppose this ODESystem can be used, in addition to the already existing functionality for the deterministic ODEProblem, to solve an according SDESystem.

I am checking with you if you believe this approach would be feasible, as I am not yet very familiar with the inner workings of the package. Any remarks and suggestions would be much appreciated, as well as possibly pointing me in the direction of the functions I would need to take a closer look at to implement this.

Thanks in advance!

timedependent parameters

Is it possible (or even in scope) to explicitly time dependent, to for example study the interaction of a pulse and a cavity it would be important to have a time dependent parameter to control the pulse slope.

Sums with IndexedVariable simplifies incorrectly

Related to #188 when using IndexedVariables some resulting sums get simplified incorrectly

Example:

@cnumbers N
ha = NLevelSpace(:atoms, 2)
i = Index(ha,:i,N,ha)
gi = IndexedVariable(:g, i)
s=∑(gi,i)+∑(gi,i)
SymbolicUtils.simplify(s)

outputs

2N*gi

which is not correct. Is the missing of an operator an issue?

Typed parameters

The update to MTK v9 #197 required an overload of concrete_symtype to work:

MTK.concrete_symtype(::Symbolics.BasicSymbolic{T}) where T <: CNumber = ComplexF64

That isn't exactly great but should work for most use cases.

However, this now enables us to actually change the implementation of Parameter such that we can represent actual types. The whole parallel implementation for RealParameter could go and we can represent also e.g. 32-bit floats.

We'd just need to add another type parameter, something like

struct Parameter{T<:Number} <: CNumber
...
end

Then, changing the overload to something like

MTK.concrete_symtype(::Symbolics.BasicSymbolic{T}) where {S<:Number, T <: Parameter{S}} = S

could work.

This could be complemented by some syntactic sugar for the @cnumbers macro, basically just stealing the syntax from Symbolics, @cnumbers p::Float32...

UndefVarError: avg not defined

Hi,

I have been trying to do some calculations with QuantumCumulants.jl but when I try to solve the resulting ODEProblem, I keep getting this UndefVarError saying that avg is not defined. Here is a short script that reproduces the issue. I am using v0.2.7.

using QuantumCumulants
using OrdinaryDiffEq, ModelingToolkit
using Plots

@syms t::Real
@register r(t)
@register i(t)

hq = NLevelSpace(:qubit,(:g,:e))
a = Transition(hq, :annihilate, :g, :e)
σ = Transition(hq, :σ)
σz = σ(:g, :g) - σ(:e, :e)
σy = 1im*(σ(:e, :g)-σ(:g, :e))
σx = σ(:g, :e) + σ(:e, :g)

#Hamiltonian 

H = r(t) * (a' + a) + i(t)*1im*(a' - a)

eqns = meanfield([σx, σy, σz], H)

@named sys = ODESystem(eqns)


function r(t)
    return 2π*19e-3
end

function i(t)
    return 0.
end

u0 = [1.0 + 0im, 0. + 0im, 0. + 0im]

prob = ODEProblem(sys, u0, (0.0, 1/(4*19e-3)))

sol = solve(prob, Tsit5()) #where error happens

t = sol.t
s22 = real.(sol[σz])
plot(t, s22, xlabel="t", ylabel="⟨σz⟩", legend=false, size=(600,300))

And here is the resulting stacktrace.

Is this some kind of versioning issue? I can't seem to install v0.2.8.

Thank you so much,

Aditya

How can I define the transition operator with two subscripts？

for example σim,
i indicates different atoms and m indicates different clusters.

`change_index` not working for a summed term

I have successfully used the following expression:
change_index(σ(1,2,j),j,i) =$\sigma_i^{12}$, It works as expected!

However, when attempting to calculate the product of collective operators, such as:

$\sum_i\sigma^{12}_i \sum_i \sigma^{23}_i$;
$\sum_i\sigma^{12}_i \sum_i \sigma^{12}_i$;

Using the code:

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

h1 = Σ(σ(1,2,j),j);

h1*h1

I encountered an error: "Specification of an extra Index is needed!"

Yet, if I utilize:

change_index(h1, j, i)

It still returns $\sum_{j}^N\sigma_j^{12}$. Although I understand that h1 is a sum that is not dependent on the index j, the h1 * h1 operation requires a different index label.

Any insights on resolving this issue would be greatly appreciated!

Periodic boundary conditions

Hello,

I was reading the Symbolic-Sums-and-Indices section of the docs, but I could not find there whether I can define a model on a lattice with periodic boundary conditions.
For example, for a simple tight-binding model, I wanted to to do something like that
ham = ∑(a(j)*a(j+1)',j)+∑(a(j)'*a(j+1),j)
This does not work. Is there another way of achieving this?

CorrelationFunction with steady_state = false

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?

Simplify expressions

Hello,

I want to find a simplify method for this simple test code

using QuantumCumulants

# Define hilbert space
hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = hf ⊗ ha

# Define the fundamental operators
@qnumbers a::Destroy(h) σ::Transition(h,:g,:e)

σm = σ
σp = σ'
sz = σp * σm - σm * σp

sz

which returns (-1+σee+σee) instead of (-1+2*σee).
Is there a way to simplify it?

`complete` and `scale` fail when using linear combinations of operators

Below is the Short Example from Symbolic Sums and Indices.

using QuantumCumulants

ha = NLevelSpace(:atoms,2)
hc = FockSpace(:cavity)
h = hc ⊗ ha

@cnumbers N Δ κ γ ν

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

@qnumbers b::Destroy(h)
σ(x,y,z) = IndexedOperator(Transition(h,:σ,x,y),z)
gi = IndexedVariable(:g,i)

H = Δ*b'*b + ∑(gi*(b*σ(2,1,i) + b'*σ(1,2,i)),i)
ops = [b'b, σ(2,2,j)]
J = [b, σ(1,2,i), σ(2,1,i)]
rates = [κ, γ, ν]

eqs = meanfield(ops,H,J;rates=rates,order=2)
eqs_complete=complete(eqs)

This produces a repeated set of equations of motion, as expected. However, if I replace b'*b in ops with 2*b'*b or (it appears) any linear combination of operators, the call to complete returns an error

ERROR: LoadError: MethodError: no method matching _inconj(::SymbolicUtils.Mul{CNumber, Int64, Dict{Any, Number}, Nothing})
Closest candidates are:
  _inconj(::SymbolicUtils.Term{<:QuantumCumulants.AvgSym}) at ~/.julia/packages/QuantumCumulants/MXRNL/src/index_utils.jl:59
Stacktrace:
  [1] iterate
    @ ./generator.jl:47 [inlined]
  [2] _collect(c::Vector{SymbolicUtils.Symbolic}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic}, typeof(QuantumCumulants._inconj)}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:807
  [3] collect_similar(cont::Vector{SymbolicUtils.Symbolic}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic}, typeof(QuantumCumulants._inconj)})
    @ Base ./array.jl:716
  [4] map(f::Function, A::Vector{SymbolicUtils.Symbolic})
    @ Base ./abstractarray.jl:2933
  [5] indexed_complete!(de::QuantumCumulants.IndexedMeanfieldEquations; order::Int64, multithread::Bool, filter_func::Nothing, mix_choice::Function, simplify::Bool, extra_indices::Vector{Symbol}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ QuantumCumulants ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:321
  [6] indexed_complete!(de::QuantumCumulants.IndexedMeanfieldEquations)
    @ QuantumCumulants ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:203
  [7] #indexed_complete#405
    @ ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:194 [inlined]
  [8] indexed_complete
    @ ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:192 [inlined]
  [9] #complete#450
    @ ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:644 [inlined]
 [10] complete(eqs::QuantumCumulants.IndexedMeanfieldEquations)
    @ QuantumCumulants ~/.julia/packages/QuantumCumulants/MXRNL/src/index_meanfield.jl:644
 [11] top-level scope

I similarly find scale(eqs) does not work, this time with MethodError: no method matching hilbert(::Int64). Is this expected behaviour? I am new to using Julia and the package so apologies if I am missing something obvious. My use case is taht I would like to directly determine equations of motion for operators such as sigma^z=\sigma(2,2)-\sigma(1,1) and so on, if possible.

to_numeric yields wrong results for larger Hilbert space

When taking a larger Fock cutoff, to_numeric yields wrong results in some cases.
Example:

using QuantumCumulants
using QuantumOptics

hfock = FockSpace(:fock)
@qnumbers a::Destroy(hfock)
bfock = FockBasis(100)

diff = (2*create(bfock)+2*destroy(bfock)) - to_numeric((2*(a)+2*(a')), bfock)

The difference diff between these operators is non-zero.

The conversion from symbolic to numeric operators in QuantumOptics causes solvers to fail since I think it breaks Hermiticity.
Any idea what causes this strange behavior?

scale sums: prefactor N or (N-1)

Hi,
I noticed that some terms in the mean field equation scale like N, while others scale like N-1.
For example:

#v0.2.19
@cnumbers N
h= NLevelSpace(:s,2)
i1 = Index(h,:i1,N,h)
i2 = Index(h,:i2,N,h)
s(α,β,i) = IndexedOperator(Transition(h, :σ, α, β),i)

scale(average(Σ(s(1,2,i1)*s(2,2,i2),i1,[i2]))) # (−1+𝑁)⟨𝜎121𝜎222⟩
scale(average(Σ(s(1,2,i1),i1,[i2]))) #𝑁⟨𝜎121⟩

Is this intended? Importantly, I get

scale(average(Σ(s(1,2,i1)*s(1,1,i2),i1,[i2]))) # 𝑁⟨𝜎121⟩+(1−1𝑁)⟨𝜎121𝜎222⟩

and because of such terms, I sometimes cannot divide some other parameter by N to eliminate N from the mean field equations.

Negative photon number with a coherent driving

Dear QuantumCumulants developer,

I attempted to add a coherent driving term \eta*(a'-a) in the Hamiltonian, and the resultant time evolution of the photon number <a'a> contains huge negative values when \eta was simply chosen to be 1 (ode solver failed when \eta becomes larger, e.g. 10^5). When \eta is 0, the time evolution of the photon number looks correct. Please see the attached figure.

May I ask whether QuantumCumulants is suitable for simulating the coupling between many atoms and a single cavity mode under a coherent driving of the cavity mode?

Thank you!

Best wishes,
Eric

How to choose the right algorithm to compute the differential equation

How to understand rigid and non-rigid problems and Why the code can't run when i use rigid algorithms to solve differential equations.

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

I am attempting to solve the Tavis-Cummings model using [email protected].
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 :

Decay rate matrices are broken

There are currently two issues concerning correlated decay.

Decay matrices get turned into vectors (where only a single row is considered) because of the call to _expand_clusters here:

QuantumCumulants.jl/src/meanfield.jl

Line 43 in 3cd511e

J_, Jdagger_, rates_ = _expand_clusters(J,Jdagger,rates)
The rates field in MeanFieldEquations is mistyped (should be AbstractVecOrMat) here:

QuantumCumulants.jl/src/equations.jl

Line 35 in 53cab7b

rates::Vector

@ChristophHotter can you look into this? We also need tests.

Compilation Error

Thank you for the work you've done to make this code and documentation available.

I believe this problem is related to Symbolics.jl, but when precompiling, I get the following error:

Failed to precompile Qumulants [35bcea6d-e19f-57db-af74-8011de6c7255] to /home/ubuntu/.julia/compiled/v1.4/Qumulants/td6ql_AH2mt.ji.
error(::String) at error.jl:33
compilecache(::Base.PkgId, ::String) at loading.jl:1272
_require(::Base.PkgId) at loading.jl:1029
require(::Base.PkgId) at loading.jl:927
require(::Module, ::Symbol) at loading.jl:922

Thank you for any feedback.

Best regards,
Larry

Adjoint of CNumber not defined

I would like to avoid the following error:

# 1st Example

h = FockSpace(:osc)
@cnumbers c
adjoint(c)

# MethodError: no method matching adjoint(::SymbolicUtils.Sym{Parameter, Base.ImmutableDict{DataType, Any}})

This can also happen via

# 2nd Example

@qnumbers a::Destroy(h) 
q = (c*a+a')^2  # = (c^2*(a*a)+c*(a′*a)+(a′*a′)+c*(a′*a)+c)
adjoint(q)

# MethodError: no method matching adjoint(::SymbolicUtils.Sym{Parameter, Base.ImmutableDict{DataType, Any}})

or by calling meanfield with one of the loss terms containing q.

Is it a solution to add Base.adjoint(c::SymbolicUtils.Sym) = conj(c)?

However, if the error in the first example is desired, maybe one can avoid the error in the second example by changing
Base.adjoint(q::QAdd) = QAdd(map(adjoint, q.arguments)) to
Base.adjoint(q::QAdd) = QAdd([q_arg isa QNumber ? adjoint(q_arg) : conj(q_arg) for q_arg in q.arguments])

define real number parameter?

There is a macro @cnumber.
However sometimes, the expression can be further simplified when the parameters are real.
Is there a way to define real number parameters?

Defining Pauli Matrices

Hello, I want to get the cumulant expansion equations with real coefficients so I can evolve them using an ordinary differential equation solver for real-valued ODEs. This is doable when no jump operators if the cumulants can be recast in terms of Pauli matrices. Can you provide me with an example of how to code the Pauli-X matrix? Something like this: PauliX(n) is the Pauli-x matrix acting on site n where Pauli_x is the matrix |g><e|+|e><g|.

Thanks!

How to represent sigma_z ?

When I use sum(σ(2,2) - σ(1,1)) the system reports an error.

Unknown rates type in macOS

Dear contributors of QuantumCumulants,

First, thanks for the super useful package!

I recently encountered an issue that when 'rates' contains operations, e.g. '+' and '*', the code can't run in macOS but works well in Windows. The error message is "Unknown rates type!".

I wonder is there any solution to this?

My macOS is Catalina 10.15.7 installed in MacBook Pro (13 inches, 2015 early).

Here is the code (reproduced from Zhang et al., arXiv:2203.04102v1, 2022).

using QuantumCumulants
kbT = 1.380649e-23293; # J/KK =J, room temperature
kb = 1.380649e-23; # J/K
hh = 6.62607015e-34; # Planck's constant,unit J.s
hbar = 6.62607015e-34/(2pi);
@cnumbers N Δ31 ξ ksp k47 k57 k67 k73 k72 k71 k31 k13 k21 k12 χ2 χ3
@cnumbers Δm g31 κ nmth ωd
hm = FockSpace(:resonator)
hNV_ = NLevelSpace(:NV,7)
hNV = ClusterSpace(hNV_,N,2)
h = hm ⊗ hNV
@qnumbers a::Destroy(h)
σ(i,j) = Transition(h,:σ,i,j,2)
H = Δ31sum(σ(3,3)) + Δma'a + g31(a'sum(σ(1,3)) + asum(σ(3,1))) ;
J = [σ(4,1),σ(5,2),σ(6,3),σ(1,4),σ(2,5),σ(3,6),σ(7,4),σ(7,5),
σ(7,6),σ(2,7),σ(3,7),σ(1,7),σ(1,2),σ(1,3),σ(2,1),σ(3,1),σ(2,2),σ(3,3),a,a']
rates = [ξ,ξ,ξ,ξ+ksp,ξ+ksp,ξ+ksp,k47,k57,k67,k73,k72,k71,k21,k31,k12,k13,2χ2,2χ3,(1.0+nmth)κ,nmthκ];
ops = [a'*a,σ(2,2)[1],σ(3,3)[1]]
eqs = meanfield(ops,H,J;rates=rates,order=2)
eqs_c = complete(eqs,order=2);

and error message.

Thank you!

Best wishes,
Eric

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Sums with different sum indices don't get simplified

@cnumbers N
h = NLevelSpace(:atoms, 2)

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

sigma(i,j,k) = IndexedOperator(Transition(h,:σ,i,j),k)
sigmam(k) = sigma(1, 2, k)

SymbolicUtils.simplify(∑(sigmam(i), i)+∑(sigmam(j), j))

outputs

(Σ(i=1:N)σ12i+Σ(j=1:N)σ12j)

Even though it should

2(Σ(i=1:N)σ12i)

Correlation functions with time dependent Hamiltonians fails to solve

I am trying to calculate the correlation functions for a time-dependent Hamiltonian:

# Hilbert space
h1 = FockSpace(:cavity1)
h2 = FockSpace(:cavity2)
h3 = FockSpace(:cavity3)

h = h1 ⊗ h2 ⊗ h3 

# Operators
a1 = Destroy(h, :(a1), 1)
a2 = Destroy(h, :(a2), 2)
a3 = Destroy(h, :(a3), 3)

# Symbolically define the cavity parameters
@cnumbers η_1 η_3 θ_1 θ_2 g_2
@syms t::Real

@register f_1(t) 
@register f_2(t) 
@register f_3(t) 
@register f_4(t) 
@register f_5(t)

# Define the Hamiltonian
H = -f_1(t)*η_3*(a3 + a3') + f_2(t)*θ_2*(a2'*a3 + a2*a3') + f_3(t)*g_2*(a2*a2 + a2'*a2') +
    f_4(t)*θ_1*(a1'*a2 + a1*a2') + f_5(t)*η_1*(a1 + a1')

# Define the dissipators
J = [a1, a2, a3]
rates = γ

Solving the meanfield equations works as expected,

@named sys = ODESystem(me)

# Parameter values
ps = [η_3; θ_2; g_2; θ_1; η_1]

# Initial state
u0 = zeros(ComplexF64, length(me))

prob = ODEProblem(sys, u0, (start_time, stop_time), ps.=>pulse_amps)
sol = solve(prob, Tsit5(), maxiters=1e7, saveat=t_list)

but solving for the correlation functions give an error which I was not able to isolate:


@named csys = ODESystem(c)
u0_c = correlation_u0(c, sol.u[end])
p0_c = correlation_p0(c, sol.u[end], ps.=>pulse_amps)

prob_c = ODEProblem(csys, u0_c, (start_time, 5*stop_time), p0_c)

The ODEProblem object is created as intended, but trying to solve the problem returns an error:

sol_c = solve(prob_c, RK4(), save_idxs=1)

>> UndefVarError: `t` not defined

Stacktrace:
  [1] macro expansion
    @ [~/.julia/packages/SymbolicUtils/H684H/src/code.jl:395](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/SymbolicUtils/H684H/src/code.jl:395) [inlined]
  [2] macro expansion
    @ [~/.julia/packages/Symbolics/3jLt1/src/build_function.jl:520](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/Symbolics/3jLt1/src/build_function.jl:520) [inlined]
  [3] macro expansion
    @ [~/.julia/packages/SymbolicUtils/H684H/src/code.jl:352](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/SymbolicUtils/H684H/src/code.jl:352) [inlined]
  [4] macro expansion
    @ [~/.julia/packages/RuntimeGeneratedFunctions/TAGuS/src/RuntimeGeneratedFunctions.jl:139](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/RuntimeGeneratedFunctions/TAGuS/src/RuntimeGeneratedFunctions.jl:139) [inlined]
  [5] macro expansion
    @ [./none:0](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/none:0) [inlined]
  [6] generated_callfunc
    @ [./none:0](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/none:0) [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :τ), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x56476eeb, 0x038659b5, 0xa85c227f, 0xc226b7b7, 0xcc501726), Expr})(::Vector{ComplexF64}, ::Vector{ComplexF64}, ::Tuple{Float64, Float64, ComplexF64, Float64, Float64, Float64}, ::Float64)
    @ RuntimeGeneratedFunctions [~/.julia/packages/RuntimeGeneratedFunctions/TAGuS/src/RuntimeGeneratedFunctions.jl:126](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/RuntimeGeneratedFunctions/TAGuS/src/RuntimeGeneratedFunctions.jl:126)
  [8] f
    @ [~/.julia/packages/ModelingToolkit/xecyK/src/systems/diffeqs/abstractodesystem.jl:286](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/ModelingToolkit/xecyK/src/systems/diffeqs/abstractodesystem.jl:286) [inlined]
  [9] Void
    @ [~/.julia/packages/SciMLBase/qp2gL/src/utils.jl:468](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/SciMLBase/qp2gL/src/utils.jl:468) [inlined]
 [10] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{ModelingToolkit.var"#f#522"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :τ), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x2592f7ee, 0x0935354f, 0xa063172c, 0x5b7ceb2b, 0xe473f033), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :τ), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x56476eeb, 0x038659b5, 0xa85c227f, 0xc226b7b7, 0xcc501726), Expr}}}, arg1::Vector{ComplexF64}, arg2::Vector{ComplexF64}, arg3::Tuple{Float64, Float64, ComplexF64, Float64, Float64, Float64}, arg4::Float64)
    @ FunctionWrappers [~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65)
 [11] macro expansion
    @ [~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137) [inlined]
...
    @ [~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:945](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:945) [inlined]
 [25] solve(prob::ODEProblem{Vector{ComplexF64}, Tuple{Int64, Int64}, true, Tuple{Float64, Float64, ComplexF64, Float64, Float64, Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#522"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :τ), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x2592f7ee, 0x0935354f, 0xa063172c, 0x5b7ceb2b, 0xe473f033), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :τ), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x56476eeb, 0x038659b5, 0xa85c227f, 0xc226b7b7, 0xcc501726), Expr}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#570#generated_observed#530"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:save_idxs,), Tuple{Int64}}})
    @ DiffEqBase [~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:882](https://file+.vscode-resource.vscode-cdn.net/Users/leonbello/embedded-amplifiers-control/_research/notebooks/~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:882)
 [26] top-level scope

Would appreciate any help on what I may be doing wrong.

Implementation example in the documentation doesn't get fully simplified

Currently, the example in the documentation show-casing the implementation of position and momentum operators produces equations that aren't fully simplified:

This is probably due to the Div type introduced in SymbolicUtils. We may need to call out to different simplification routines, e.g. simplify_fractions when simplifying in the meanfield function. I couldn't find any recommendation in the docs on how to best simplify fractions though.

_new_operator(nOp::NumberedOperator)

Should line 7 of index_correlation.jl read

_new_operator(nOp::NumberedOperator,h,aon=acts_on(nOp)) = NumberedOperator(_new_operator(nOp.op,h,aon),nOp.numb)

instead of

_new_operator(nOp::NumberedOperator,h,aon=acts_on(nOp)) = NumberedOperator(_new_operator(op.op,h,aon),nOp.numb)

(op.op -> nOp.op)?

`IndexedVariable` with a time-dependence

I'm trying to solve a problem with all-to-all coupling between different two-level systems, with the following parameters,

# Parameters
ha = NLevelSpace(:atoms, 2)
hb = FockSpace(:bus)
h = hb ⊗ ha

b = Destroy(h, :b)

@cnumbers L
j = Index(h, :j, L, ha)
k = Index(h, :k, L, ha)

@cnumbers Δ_b ϵ_r κ γ ν
χ(j) = IndexedVariable(:χ, j)
Δ(j) = IndexedVariable(:Δ, j)

σz(i) = IndexedOperator(Transition(h, :σ, 2, 2), i)
σp(i) = IndexedOperator(Transition(h, :σ, 1, 2), i)
σm(i) = IndexedOperator(Transition(h, :σ, 2, 1), i)

Some of my parameters are double-indexed,

u(j, k) = IndexedVariable(:u, j, k)
J(j, k) = IndexedVariable(:J, j, k)
δ(j,k) = IndexedVariable(:δ, j, k)

I'm able to run all the definitions, including carrying out the meanfield equations,

l = Index(h, :l, L, ha)
ops = [b'b, σz(l)]
Jops = [b, σp(j), σm(j)]
rates = [κ, γ, ν]

eqs = meanfield(ops, H, Jops; rates=rates, order=2)
eqs = indexed_complete(eqs)

After that, I run into two problems. First, when using evaluate,

evaled = evaluate(eqs; limits=(L=>2))

L is not completely removed from the equations and I have to include it in the variable map.

Second, and more severe, I get and UndefVarErrror when I try to run solve the resulting ODESystem,

@named sys = ODESystem(evaled)
u0 = zeros(ComplexF64, length(evaled))

p = [
    L,
    Δ_b,
    ϵ_r,
    χ(j),
    Δ(j),
    u(j,k),
    J(j,k),
    δ(j,k), 
    γ,
    κ,
    ν
]

p0 = [
        2,
        0, 
        0.1,
        [-2.5, -1.9],
        [1, 1.1], 
        [0 2; 1 0], 
        [0 0.1; 0.12 0], 
        [0 10; 15 0], 
        0,
        0,
        0
    ] 

p_ = value_map(p, p0; limits=(L=>2))
prob = ODEProblem(sys, u0, (0.0, 10.0), p_)

When I try to solve prob I get an error,

sol = solve(prob,RK4())
>> UndefVarError: `δjk` not defined

I suspect I may be doing something wrong in the assignment of the values in value_map, but I'm not sure what I'm doing wrong.

Note that I've only used L=2 for testing purposes, so I can't just define my parameters explicitly and have to use the IndexedOperator scheme.

Full error trace:

{
	"name": "UndefVarError",
	"message": "UndefVarError: `δjk` not defined",
	"stack": "UndefVarError: `δjk` not defined

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:418 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/gBKZv/src/build_function.jl:537 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:375 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc
    @ ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0x1c9db295, 0xf2731d10, 0x92b1b07e, 0x8515d200, 0x08e29352), Nothing})(::Vector{ComplexF64}, ::Vector{ComplexF64}, ::Vector{ComplexF64}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
  [8] k
    @ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/diffeqs/abstractodesystem.jl:371 [inlined]
  [9] Void
    @ ~/.julia/packages/SciMLBase/8XHkk/src/utils.jl:481 [inlined]
 [10] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{ModelingToolkit.var\"#k#551\"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0xce8523c1, 0xd478be1f, 0x696dbb94, 0xd4a4fef1, 0x7d9d6b9a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0x1c9db295, 0xf2731d10, 0x92b1b07e, 0x8515d200, 0x08e29352), Nothing}}}, arg1::Vector{ComplexF64}, arg2::Vector{ComplexF64}, arg3::Vector{ComplexF64}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [11] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [12] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [13] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [14] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
 [15] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [16] ODEFunction
    @ ~/.julia/packages/SciMLBase/8XHkk/src/scimlfunctions.jl:2407 [inlined]
 [17] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{ComplexF64}, Nothing, Float64, Vector{ComplexF64}, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.RK4Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, SciMLBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, OrdinaryDiffEq.RK4Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.RK4Cache{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/perform_step/fixed_timestep_perform_step.jl:336
 [18] __init(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ComplexF64}, Vector{ComplexF64}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{Complex{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}, Vector{ComplexF64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, ComplexF64}, ComplexF64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:512
 [19] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:10 [inlined]
 [20] #__solve#746
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:5 [inlined]
 [21] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:1 [inlined]
 [22] #solve_call#34
    @ ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:608 [inlined]
 [23] solve_call
    @ ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:566 [inlined]
 [24] solve_up(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var\"#k#551\"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0xce8523c1, 0xd478be1f, 0x696dbb94, 0xd4a4fef1, 0x7d9d6b9a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0x1c9db295, 0xf2731d10, 0x92b1b07e, 0x8515d200, 0x08e29352), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{ComplexF64}, p::Vector{ComplexF64}, args::RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:1057
 [25] solve_up
    @ ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:1043 [inlined]
 [26] #solve#40
    @ ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:980 [inlined]
 [27] solve(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, Vector{ComplexF64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var\"#k#551\"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0xce8523c1, 0xd478be1f, 0x696dbb94, 0xd4a4fef1, 0x7d9d6b9a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0x1c9db295, 0xf2731d10, 0x92b1b07e, 0x8515d200, 0x08e29352), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var\"#661#generated_observed#561\"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Parameter}}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::RK4{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/7s9cb/src/solve.jl:970
 [28] top-level scope
    @ ~/qbus-rc/basic_qbus.ipynb:1"
}

ModelingToolkit?

It doesn't really make sense that this is using SymbolicUtils.jl instead of ModelingToolkit.jl, which is the higher level of SymbolicUtils.jl that does the derivations for lots of the things in this package. I'm curious why that choice was made so we can make sure to fix the messaging.

MTKParameterCtx ist not defined

@cnumbers stopped working for me, probably because of ModelingToolkit commit 194064ea462c0573e042c771cd8ae875335533e3, which removes the struct MTKParameterCtx.

Workaround: using a version of ModelingToolkit before this commit seems to work: e.g.:
pkg> add ModelingToolkit#af025524ba90cf51f8e5354eb1d0d940c6b06618

Change how averages are represented

Currently, averages are symbolic functions with a QNumber argument, i.e. avg(::QNumber)::CNumber. Once the metadata system of SymbolicUtils is fixed, we can change that so that averages are directly represented in the format required by ModelingToolkit. We can write them as time-dependent variables x(t) and store the QNumber argument as metadata. This saves us the step of converting to ODESystems with the varmap.

`test_correlation.jl` fails for newer `Symbolics.jl` versions

I am trying to use QuantumCumulants.jl functionality in HarmonicBalance.jl, https://github.com/NonlinearOscillations/HarmonicBalance.jl . The latter relies on newer version of Symbolics.jl, but the former is capped at < 4.3.1

Updating to 4.11.0 gives

Use ModelingToolkit's parameters macro

Currently we have the @cnumbers macro for this, but we can just get rid of this once MTK properly supports complex numbers. As of yet, using @parameters from MTK will lead to Complex{Num} types which don't play nice with simplification yet.