Liouvillian superoperator way cannot run by `timeevolution.master_dynamic` of QuantumOptics.jl

Specific Domains Quantum
, , , ,
1

When i use the package QuantumOptics.jl to simulate the evolution of Liouvillian superoperator, i find it cannot run by the command timeevolution.master_dynamic. Because after obtaining the superoperator, i ready to solve the propagator of this superoperator. I believe it cannot use the command liouvillian, and i cannot find out which type of program my program belongs to. Thus, i need to ask if it can help me solve my difficulties.
My program is

using QuantumOptics
using PyPlot
using SpecialFunctions
using BenchmarkTools
using LinearAlgebra
elapsed_time_1 = @elapsed begin

b = NLevelBasis(4)
s0 = nlevelstate(b,1)
s1 = nlevelstate(b,2)
sr = nlevelstate(b,3)
se = nlevelstate(b,4)
plus = (s0 + s1)/sqrt(2)
minus = (s0 - s1)/sqrt(2)

eye4 = identityoperator(b)
eye64 = identityoperator(b⊗b⊗b)


sigma_x = s0⊗dagger(s1) + s1⊗dagger(s0)
sigma_y = -1im*s0⊗dagger(s1) + 1im*s1⊗dagger(s0)
sigma_z = s0⊗dagger(s0) - s1⊗dagger(s1)

GHZ_00 = ( tensor(plus, s0, plus) + tensor(minus, s1, minus) )/sqrt(2)
GHZ_01 = ( tensor(plus, s0, plus) - tensor(minus, s1, minus) )/sqrt(2)
GHZ_10 = ( tensor(plus, s0, minus) + tensor(minus, s1, plus) )/sqrt(2)
GHZ_11 = ( tensor(plus, s0, minus) - tensor(minus, s1, plus) )/sqrt(2)
GHZ_20 = ( tensor(plus, s1, plus) + tensor(minus, s0, minus) )/sqrt(2)
GHZ_21 = ( tensor(plus, s1, plus) - tensor(minus, s0, minus) )/sqrt(2)
GHZ_30 = ( tensor(plus, s1, minus) + tensor(minus, s0, plus) )/sqrt(2)
GHZ_31 = ( tensor(plus, s1, minus) - tensor(minus, s0, plus) )/sqrt(2)

# parameters
Omega_r = 0.2*gamma_P1;   Omega_a = 0.2*gamma_P2;

Omega = 2*pi*1*10^6;   Urr = 2*pi*150*10^6;
alpha = 43.184769223265;
omega_0 = Urr/14;
Delta_0 = 0*omega_0;
delta = alpha*omega_0;

m0 = Delta_0/omega_0;
m1 = (Delta_0-Urr)/omega_0;
m2 = (Delta_0-2*Urr)/omega_0;
phase = m0*pi/2-alpha-pi/2 + pi*sign(besselj(m0,alpha));

t1 = abs(pi/(Omega*besselj(m0,alpha)))
maxtime = t1

# initial and final states
psi_0 = (GHZ_00 + GHZ_01 + GHZ_10 + GHZ_11 + GHZ_20 + GHZ_21 + GHZ_30 + GHZ_31)/sqrt(8);
psi_f = GHZ_00;
psi0_l = tensor(psi_0,psi_0)
psif_l = tensor(psi_f,psi_f)

propag_0 = identityoperator(b⊗b⊗b⊗b⊗b⊗b)

tlist = range(0, stop=maxtime, length=1000)

# Hamiltonian
Hv = Urr*tensor(sr⊗dagger(sr), sr⊗dagger(sr), eye4) + Urr*tensor(sr⊗dagger(sr), eye4, sr⊗dagger(sr)) + Urr*tensor(eye4, sr⊗dagger(sr), sr⊗dagger(sr))

function Liouvillian(t,rho)
        
    detuning = Delta_0*t + (2*delta*sin((omega_0*t)/2)^2)/omega_0;

    HP = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(plus) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(plus));
    H0 = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s0) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s0));
    H1 = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s1) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s1));

    SEP_123 = tensor(H0, eye4, eye4) + tensor(eye4, HP, eye4) + tensor(eye4, eye4, H1) + Hv;
    Hami = SEP_123;
    Liouvillian = -1im*(tensor(Hami, eye64) - tensor(eye64, transpose(Hami)))

    return Liouvillian
end

tout, rho = timeevolution.master_dynamic(tlist, propag_0, Liouvillian)

exp_psi0_l = real(expect(psif_l ⊗ dagger(psif_l), rho*(psi0_l ⊗ dagger(psi0_l))))


figure(figsize=(9,3))
# subplot(1,2,1)
# ylim([0, 2])
plot(tlist, exp_psi0_l);
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()


end

println("run time：", elapsed_time_1)

gcf()

And it gives the Error:

ERROR: MethodError: no method matching OrdinaryDiffEq.DP5Cache(::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::typeof(OrdinaryDiffEq.trivial_limiter!), ::typeof(OrdinaryDiffEq.trivial_limiter!), ::Static.False)

Closest candidates are:
  OrdinaryDiffEq.DP5Cache(::var"#2018#uType", ::var"#2018#uType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2018#uType", ::var"#2018#uType", ::var"#2020#uNoUnitsType", ::var"#2021#StageLimiter", ::var"#2022#StepLimiter", ::var"#2023#Thread") where {var"#2018#uType", var"#2019#rateType", var"#2020#uNoUnitsType", var"#2021#StageLimiter", var"#2022#StepLimiter", var"#2023#Thread"}
   @ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\caches\low_order_rk_caches.jl:552

Stacktrace:
  [1] alg_cache(alg::OrdinaryDiffEq.DP5{…}, u::Diagonal{…}, rate_prototype::Diagonal{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Diagonal{…}, uprev2::Diagonal{…}, f::SciMLBase.ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, ::Val{…})
    @ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\caches\low_order_rk_caches.jl:600
  [2] __init(prob::SciMLBase.ODEProblem{…}, alg::OrdinaryDiffEq.DP5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::SciMLBase.CallbackSet{…}, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Float64, reltol::Float64, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(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::@Kwargs{})
    @ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:342
  [3] __solve(::SciMLBase.ODEProblem{…}, ::OrdinaryDiffEq.DP5{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:5
  [4] solve_call(_prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEq.DP5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:609
  [5] solve_up(prob::SciMLBase.ODEProblem{…}, sensealg::Nothing, u0::Diagonal{…}, p::SciMLBase.NullParameters, args::OrdinaryDiffEq.DP5{…}; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:1058
  [6] solve(prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEq.DP5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:981
  [7] integrate(tspan::StepRangeLen{…}, df::QuantumOptics.timeevolution.var"#dmaster_#41"{…}, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::QuantumOptics.timeevolution.var"#fout#7"; alg::OrdinaryDiffEq.DP5{…}, steady_state::Bool, tol::Float64, save_everystep::Bool, saveat::StepRangeLen{…}, callback::Nothing, kwargs::@Kwargs{})
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:59
  [8] integrate(tspan::StepRangeLen{…}, df::Function, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:14
  [9] #integrate#6
    @ C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:75 [inlined]
 [10] integrate(tspan::StepRangeLen{…}, df::Function, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:70 [inlined]   
 [11] integrate_master(tspan::StepRangeLen{…}, df::Function, rho0::Operator{…}, fout::Nothing; kwargs::@Kwargs{})
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:264
 [12] integrate_master
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:260 [inlined]
 [13] #master_dynamic#40
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:217 [inlined]
 [14] master_dynamic(tspan::StepRangeLen{…}, rho0::Operator{…}, f::Function)
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:211
 [15] macro expansion
    @ d:\OneDrive(重重重重重重重重)\OneDrive\博士期间------工作\第二篇\程序\Liouvillian方法\GHZ\Liouville 传播子\stabilzer\request_help.jl:78 [inlined]
 [16] top-level scope
    @ .\timing.jl:395
Some type information was truncated. Use `show(err)` to see complete types.

Thanks for your help!

2

I have not looked in detail in the code, but I suspect the output of the function Liouvillian is not what master_dynamic expects. If I run ?timeevolution.master_dynamic I see among the output the following:

  Arguments
  ≡≡≡≡≡≡≡≡≡

    •  tspan: Vector specifying the points of time for which output
       should be displayed.

    •  rho0: Initial density operator. Can also be a state vector which
       is automatically converted into a density operator.

    •  f: Function f(t, rho) -> (H, J, Jdagger) or f(t, rho) -> (H, J,
       Jdagger, rates)

It seems the callable is supposed to return a tuple containing a Hamiltonian, a collapse operator and the conjugate of the collapse operator. There are certainly fancier options (e.g. if you want to prepare in advance some superoperator or if you have multiple collapse operators, etc), but this seems to be the basic one.

Is it easy to restate your problem in terms of “collapse operator as expected in Lindblad master equation”? That is what is expected here. But if that is not easy, we can figure out a different way to set up your master equation.

3

I review the example in package QuantumOptics.jl, and find your comment is right. However, i want to solve \frac{d}{dt}U = \mathcal{L}U, where U can be regarded as “propagator”. I think it is more similar to master equation. Furthermore, it doesn’t need the collapse operator. I cannot find the counterpart in those commands of package QuantumOptics.jl.

4

Just skimming through the documentation, it does not seem like there is a built-in way in QuantumOptics.jl to evolve an arbitrary time-dependent superoperator. It is easy only for Liouvillian’s that can be expressed in the Lindblad form.

Can you rephrase your problem to be in the Lindblad form? That would be the easiest solution.

If that is not doable, we can discuss how to set up evolution with an arbitrary time-dependent superoperator, but that would require looking into some of the internals of QuantumOptics.jl (which can also be an opportunity to contribute new functionality to the package).

5

Separate question: looking at the actual code in your example, this seems like purely unitary evolution. Why do you need to even use a master equation or Liouvillians in such a situation? Why not just a unitary evolution equation?

6

For such unitary evolution timeevolution.schroedinger_dynamic should be enough (even if the starting state is a density matrix instead of a ket)

7

As you said, the program can use unitary evolution. However, the program is simplified by me , in fact, it includes collapse operator. For the sake of simplicity, I didn’t give it specifically, which is why I use Liouvillian superoperator. To find a way to describe the dynamical evolution including dissipation, I use the liouvillian method \dot{U} = \mathcal{L}U. The propagator U can describe the evolution of system, independent with initial state.

8

I am not sure I completely understand. Here are a few comments around this though:

  • Generally, in this library, you are not supposed to construct a time-dependent Liouvillian. There are simply no convenient tools in this library for this approach. If you consider the explicit construction of a Liouvillian important for solving your problem, we can discuss how this can be implemented.
  • If you want to do time-dependent unitary evolution, you are supposed to use schroedinger_dynamic and provide a Hamiltonian (not a Liouvillian).
  • If you want to do time-dependent non-unitary evolution, you are supposed to use master_dynamic and provide a Hamiltonian and Lindblad collapse operators (not a Liouvillian).
  • Lastly, if you want to construct superoperators in this library, this is usually not done with tensor products. In this library we try to distinguish “operator on a larger hilbert space” from “superoperator”, as this distinction makes programming a bit easier, even if these objects can mathematically be represented by the same big matrix. Check out functions mentioned at the bottom of this page: Super-operators · QuantumOptics.jl
9

Before, I wrote this code in MATLAB. However, my scheme is multi-step floquet scheme (this program here shows the most important one of the steps, which will also be repeated multiple times in the final scheme), numerical calculations have high requirements (the difference between the maximum and minimum values is too large), my Matlab program needs to run for 4-7 days and is running on 3-atoms GHZ state. Besides, my later goal require calculating 4-atoms cluster state. So I want to use the Liouvillian superoperator method to obtain the “propagator” of the dynamical evolution of a dissipative system, in order to easily repeat multiple operational steps. The difficulty discovered again is that Liouvillian way cannot be executed or is very slow in Matlab due to its large dimensions (4096*4096). I heard that Julia has a very deep optimization of programs related to quantum optics (and in my own practical experience, I have indeed found that its running speed exceeds that of MATLAB when executes the same scheme). Therefore, I turned to Julia, and learn Julia to achieve my programming goals. So, I need to solve my simplified programming solution and then expand to complete my final code

10

1

Here is one plausible way to do it. I am not sure it is the best one, it is also not very ergonomic (we probably should provide more convenient interface for this), etc, but it can be worthwhile to try it:

rho_basis = [...] # a vector of basis density matrices (they do not really need to be valid density matrices, just a set of matrices that form a basis for the space of matrices)
rho_end_states = [] # empty vector
for rho in rho_basis
    end_state = timeevolution.master_dynamic(rho, ...)... # evolve this particular basis density matrix
    push!(rho_end_states, end_state)

Now from rho_basis and rho_end_state you can construct your superoperator.

2

Here is another approach:

If you are permitted to discretize time you can write your final superoperator as S = S_1 * S_2 * S_3 ... where S_i = exp(L_i * delta_t) where L_i is the Liouvillian at time-step i.

3

FYI, you might also want to consider:

  • maybe it is better to just run Monte Carlo for noisy Clifford Circuits given that you are interested only in stabilizer states. But maybe it would be difficult to represent the type of noise you care about. The QuantumClifford.jl package is pretty good for this type of stuff and it has good interoperability with QuantumOptics.jl
  • maybe it is worthwhile to investigate the use of tensor networks. PastaQ and ITensor might be valuable.
11

Thanks for your suggestion very much. For the approach 1, I may not understand your thoughts. My final goal is to solve the propagator U including dissipative dynamic, but it seems that your idea only provides the final state. For the approach 2, I believe it is feasible, and I will try it this way. For approach 3, I am more concerned about the dynamical evolution process of my scheme in actual physical systems (Rydberg atom), therefore, I’m not sure if they (Clifford Circuits, PastaQ and ITensor) can meet my coding expectations (More likely, it’s because I don’t understand these tools). In the future, I can delve deeper into learning them.

Nonetheless, I still want to consult with you: is it possible to obtain U according to the coding idea I provided initially? Or, can we use a universal method for solving ODE to solve this problem?

12

For 1: It does not get just one final state, it gets a mapping from all basis states to all corresponding final states, which is one way to define a superoperator (just like how a unitary operator is equivalently defined by how it maps all basis states).

For 3: The Clifford approach might not work, but I suspect the tensor network approach will still be pretty fruitful.

Concerning your question about setting it up as an ODE: Yes, internally everything done by QuantumOptics.jl is to convert the nice complicated unitaries or states (with all their basis and tensor product structure) into simple matrices and feed these to solvers from DifferentialEquations.jl. If you want to, you can use QuantumOptics.jl as a convenient way to prepare these matrices and then feed everything to DifferentialEquations.jl yourself.

E.g. after you have created whatever QuantumOptics.jl you like, you can extract the underlying matrix from the data attribute. Here is an example where the underlying matrix is a sparse matrix:

image
image765×87 7.34 KB

You can then use that however you wish, with some of the tools in DifferentialEquations.jl.

You probably will face performance issues in your first implementation: I see from your current code that you create a lot of temporary objects, where the cost of computation stops being the actual numerics and starts being dominated from reserving and freeing memory (i.e. you have a lot of “allocation” overhead in your code). That is something we can address as well, probably making the code 10x or 100x faster, but we can discuss this after you have your first prototype working.

13

Thanks for your patient explanation. Due to other reasons, this matter has been delayed. I will immediately engage in resolving this issue as soon as I am at ease, and I hope to continue seeking your advice at that time. Thanks again for your help.

1 Like