timedependent parameters

[benneti]

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.

[david-pl]

That's already possible, just not documented well. We should add an example on that to the documentation.

You only have to make sure that you use the same time parameter when deriving the equations of motion as you do in your time-dependent Hamiltonian. You can do this by passing in the iv kwarg when calling the meanfield function. For example, here's how you could simulate a Gaussian pulse on a cavity:

using QuantumCumulants
using Symbolics
using ModelingToolkit, OrdinaryDiffEq

@syms t::Real  # time parameter
@cnumbers Ω0 κ t0  # other parameters

# Gaussian pulse centred at t0
Ω(t,t0) = Ω0*exp(-(t - t0)^2)

# Hilbert space and operator
h = FockSpace(:cavity)
@qnumbers a::Destroy(h)

# Hamiltonian
H = Ω(t,t0)*(a + a')

# Derive equations -- set iv here
eqs = meanfield(a,H,[a];rates=[κ],iv=t,order=1)

# Proceed as usual
@named sys = ODESystem(eqs)
p = (κ=>1, Ω0=>1, t0=>5)
u0 = zeros(ComplexF64,length(eqs))
prob = ODEProblem(sys,u0,(0.0,15.0),p)
sol = solve(prob,RK4())

For context: ModelingToolkit expects time-dependent variables of the form x(t). Therefore, averages are internally converted to this form when using meanfield to derive the equations, so they always require setting a time-parameter. If you don't explicitly set iv it just defaults to Sym{Real}(:t).

Note that in the example above the symbolic variables t and t0 actually trace through the function Ω(t,t0), so your Hamiltonian and the equations will explicitly contain the exponential. This is usually fine, but the simplification will attempt to simplify the exponential as well. So if you have a problem with a very complex time-dependence and you want to avoid spending time on the simplification you can also use Symbolic.jl's function registration to avoid tracing into the function. This would also allow you to change the definition of Ω(t,t0) without having to re-derive the symbolic equations.

[benneti]

Ah cool thanks, if I get it working like I want to i'll write some documentation!

[benneti]

I think there is a slight improvement possible in the Constructor of CorrelationFunction.
Right now to calculate a correlation function one needs to set iv to the same variable as for eqs,
this could probably be solved by substituting iv of eqs (t)->iv of CorrelationFunction (τ) in the Jump operators/rates and the Hamiltonian.

eqs = meanfield([a'*a,a],H,[a];rates=[κ],iv=t,order=2)
c = CorrelationFunction(a', a, eqs, iv=t)
@named sys_corr = ODESystem(c)
_t = 5
p0_corr = correlation_p0(c,sol(_t),(p...,t=>_t))
u0_corr = correlation_u0(c,sol(_t))
prob_corr = ODEProblem(sys_corr,u0_corr,(0.0,15.0),(p0_corr...,t=>sol.t[1]))
sol_corr = solve(prob_corr,Tsit5())
plot(t->sol_corr(t,idxs=1), sol.t[[1,end]]...)

[david-pl]

I'm not 100% sure that's what you want. If you set the time parameter in the correlation function to the same one as in the standard time evolution you use as input, you're effectively resetting your time-dependent Hamiltonian to time 0. However, what you really want is to start at time sol.t[end] and integrating up to sol.t[end] + sol_corr.t[end]. Though I'm not really sure, as you seem to be setting the independent variable as parameter in your example, so I don't really know what the result here is.

In the example of the pulse, that would mean that another pulse occurs when calculating the correlation function. That won't match the state of the variables which have already been evolved for a longer time. You can actually see the pulse being applied in your solution.

The problem is that it actually doesn't seem to work properly. Besides having to use a workaround to update the time-dependent parts in the original equations, I still get errors. Here's how I would do it:

@syms τ::Real
@cnumbers t1
# Workaround to get time-dependent parts right
eqs_corr = substitute(eqs, Dict(t=>τ+t1))

c = CorrelationFunction(a', a, eqs_corr, iv=τ)
@named sys_corr = ODESystem(c)
_t = 5
p0_corr = correlation_p0(c,sol(_t),p)
u0_corr = correlation_u0(c,sol(_t))
prob_corr = ODEProblem(sys_corr,u0_corr,(0.0,15.0),(p0_corr...,t1=>sol.t[end]))
sol_corr = solve(prob_corr,Tsit5())

Unfortunately, that results in an UndefVarError: t not defined. I've never tried to calculate the correlation function with a time-dependent Hamiltonian before, so there is some undetected bug there. When I find time I'll do some digging.

[benneti]

I think I uderstand the problem, i.e. substitude(de) does only substitute in the equation, when it should also substitude in H, J and kappa. I'll try to write some code that tries to decide on which fieldnames substitution makes sense

[benneti]

Ok I think it is a bit more complicated as QuantumCumulants.QAdd and so on do not already support substitude.

Maybe for this part it would be sufficient to define
Base.length(q::QuantumCumulants.QAdd) = length(q.arguments)
and get substitude(x::QuantumCumulants.QMul{Nothing}) to behave.

[david-pl]

I think I uderstand the problem, i.e. substitude(de) does only substitute in the equation, when it should also substitude in H, J and kappa.

Ah you're right of course. You need to create a full copy of eqs with updated Hamiltonian and equations. Here's a quick fix for that:

H(t,t0) = Ω(t,t0)*(a + a')

@syms τ::Real
@cnumbers t1
# Workaround to get time-dependent parts right
eqs_corr = meanfield([a'*a,a],H(t1+τ,t0),[a];rates=[κ],iv=t,order=2)

I.e., just re-derive the equations accordingly. That's ugly though and annoying if the symbolic derivation of equations takes a long time. We should probably expose t0 (by which now I mean the time at which the time evolution of the correlation function starts) as a kwarg taking a symbolic when constructing the CorrelationFunction. Internally, we can then substitute eqs.iv => t0 wherever necessary.

Ok I think it is a bit more complicated as QuantumCumulants.QAdd and so on do not already support substitude.

That's a bug. E.g. substitute(t*(a+a'), Dict(t=>0)) works, but substitute(t*(a+a'), Dict(t => t0)) errors.

[benneti]

ok if the following assumptions are correct this should hopefully work not need rederivation of the equations and be only slightly confusing:
Assumptions:

• de0 is used to replace expectations values of the a_0 operators, i.e. there the iv of the correlator is 0
• the de0 equations seem to not be integrated and solely be given via the used initial condition
• then it makes sense to reuse t to derive "de" of the correlator and simply use another cnumber parameter to fix the wrong driving

using QuantumCumulants
using Symbolics
using ModelingToolkit, OrdinaryDiffEq
using Plots

@syms t::Real # τ::Real  # time parameter
# We define τ as a cnumber because else it is not detected as a parameter
@cnumbers Ω0 κ t0 τ # other parameters

# Gaussian pulse centred at t0
Ω(t,t0) = Ω0*exp(-(t - t0)^2)

# Hilbert space and operator
h = FockSpace(:cavity)
@qnumbers a::Destroy(h)

# Hamiltonian
H(t) = Ω(t,t0)*(a + a')

# Derive equations -- set iv here
# Workaround to get time-dependent parts right
eqs_corr = meanfield([a'*a,a],H(τ+t),[a];rates=[κ],iv=t,order=2)

@named sys = ODESystem(eqs_corr)
p = (κ=>1, Ω0=>1, t0=>5)
u0 = zeros(ComplexF64,length(eqs))
tmax=30.0
# for the meanfieldequations tau should be 0
prob = ODEProblem(sys,u0,(0.0,tmax),(p..., τ=>0))
sol = solve(prob,Tsit5())
# due to the limitations of cnumbers and Sym{Real} we now need to use t as tau and vice versa
# now t has the meaning of tau and vice versa
c = CorrelationFunction(a', a, eqs_corr, iv=t)
@named sys_corr = ODESystem(c)
function corr(t)
    p0_corr = correlation_p0(c,sol(t),(p...,τ=>t))
    u0_corr = correlation_u0(c,sol(t))
    prob_corr = ODEProblem(sys_corr,u0_corr,(0.0,tmax),p0_corr)
    sol_corr = solve(prob_corr,Tsit5())
end
τs = range(0,tmax,length=100)
hm = hcat([corr(t)(τs, idxs=1) for t in τs]...)
heatmap(τs,τs,abs.(hm))

Also I would be glad to help out a bit, it just might take some time for me to figure out the internals better first.

[david-pl]

If I remember correctly those assumptions should be fine. So in the end, this should work. It's a more efficient solution, but yeah, it's a bit confusing 😅

I still think an actual solution to this would be to expose the original time evolution stopping point, something like

@syms τ t0
c = CorrelationFunction(a', a, eqs, iv=τ, iv0=t0)

and internally substitute (which should be reasonably fast) eqs.iv => iv0 in the Hamiltonian and jumps.

Also I would be glad to help out a bit, it just might take some time for me to figure out the internals better first.

That would be great! The CorrelationFunction may be a bit of a tough place to start as the internal logic is a bit of a mess. There might be substantial changes required to make the construction and everything easier to understand. Feel free to give it a go if you want and ping me if you get stuck somewhere.

[benneti]

That's a bug. E.g. substitute(t*(a+a'), Dict(t=>0)) works, but substitute(t*(a+a'), Dict(t => t0)) errors.

Do you think this is a bug in QuantumCumulants or Symbolic{s,Utils}, if the former any pointer on where to fix it?
The error MethodError: no method matching +(::typeof(*), ::typeof(*)) is not too helpful for me.

[david-pl]

I'm not sure. It's probably something we missed when updating to the latest versions of the symbolics packages. Maybe there is some functionality missing that we need to implement in order for susbtitute to work properly. I didn't have time to look into it yet.

[ChristophHotter]

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.

There is now an example with a time dependent pulse on the documentaion.
With the @register macro you can use every function you like, also interpolated data points or something like this.
Maybe you need this on some point.

[benneti]

thanks for the heads up and the nice example!
Do you want to leave this open as a reminder for the correlation function?