I have another question on `timeevolution.master_dynamic`

and `timeevolution.schroedinger_dynamic`

. I run the same code on time-evolution system by `timeevolution.master_dynamic`

with density operator and `timeevolution.schroedinger_dynamic`

and ket state. However, i can obtain correct answer with `timeevolution.schroedinger_dynamic`

, while show an error with `timeevolution.master_dynamic`

.

Below is my program:

```
using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin
dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms
b = NLevelBasis(dim) # Basis for a single atom
s0 = nlevelstate(b,1) # |0> state
s1 = nlevelstate(b,2) # |1> state
sr = nlevelstate(b,3) # |r> state
plus = s0 # |+> state
minus = s1 # |-> state
# plus = (s0 + s1)/sqrt(2) # |+> state
# minus = (s0 - s1)/sqrt(2) # |-> state
Identity4 = identityoperator(b)
Hzeros = 0*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)
XZ = tensor(sigma_x, sigma_z, Identity4)
ZXZ = tensor(sigma_z, sigma_x, sigma_z)
ZX = tensor(Identity4, sigma_z, sigma_x)
Omega_1 = 2*pi*3*10^6; Omega_2 = 2*pi*3*10^6;
Delta = 20*Omega_1; U_rr = Delta;
NN = 6;
MM = 16;
t2 = pi/NN/Omega_2;
t1 = 7*t2;
t_pump = t1 + t2;
t_decay = 3.02*10^-6;
t_cycle = 3*NN*t_pump;
gamma_r = 2*pi*0*10^3;
gammar_P = 1/(26.2*10^-9);
Omega_pr = 0.2*gammar_P;
gammar_eff = Omega_pr^2/gammar_P;
gamma_0r = 2/3*gammar_eff;
gamma_1r = 1/3*gammar_eff;
Hv = U_rr*( tensor(srâdagger(sr), srâdagger(sr), Identity4) + tensor(srâdagger(sr), Identity4, srâdagger(sr)) + tensor(Identity4, srâdagger(sr), srâdagger(sr)) );
function URP123(t,args=nothing)
H0_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(s0) + Omega_1/2*exp(1im*Delta*t)*s0âdagger(sr) - (Omega_1/2)^2/Delta*(s0âdagger(s0) - srâdagger(sr));
H0_ome2 = 1im*Omega_2/2*srâdagger(s0) - 1im*Omega_2/2*s0âdagger(sr);
HP_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(plus) + Omega_1/2*exp(-1im*Delta*t)*plusâdagger(sr) - (Omega_1/2)^2/Delta*(plusâdagger(plus) - srâdagger(sr));
HP_ome2 = 1im*Omega_2/2*srâdagger(plus) - 1im*Omega_2/2*plusâdagger(sr);
H1_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(s1) + Omega_1/2*exp(-1im*Delta*t)*s1âdagger(sr) - (Omega_1/2)^2/Delta*(s1âdagger(s1) - srâdagger(sr));
H1_ome2 = 1im*Omega_2/2*srâdagger(s1) - 1im*Omega_2/2*s1âdagger(sr);
HM_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(minus) + Omega_1/2*exp(-1im*Delta*t)*minusâdagger(sr) - (Omega_1/2)^2/Delta*(minusâdagger(minus) - srâdagger(sr));
HM_ome2 = 1im*Omega_2/2*srâdagger(minus) - 1im*Omega_2/2*minusâdagger(sr);
URP123_ome1 = tensor(H0_ome1, Identity4, Identity4) + tensor(Identity4, HP_ome1, Identity4) + tensor(Identity4, Identity4, H1_ome1) + Hv;
URP123_ome2 = tensor(H0_ome2, Identity4, Identity4) + tensor(Identity4, HP_ome2, Identity4) + tensor(Identity4, Identity4, H1_ome2) + Hv;
return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1)
end
psi_0 = tensor(s0,minus,s0);
psi_f = tensor(sr,minus,s0);
dm0 = psi_0âdagger(psi_0);
rho_final = psi_fâdagger(psi_f);
tlist_pump = range(0, stop=t_pump, length=1000);
exp_psi0 = [];
for jj = 1:NN
global tout, rho_XZ, dm0, psi_0
tout, rho_XZ = timeevolution.master_dynamic(tlist_pump, dm0, URP123, []);
append!(exp_psi0, real(expect(rho_final, rho_XZ)))
dm0 = rho_XZ[end];
end
plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(exp_psi0)
xlabel(L"T")
ylabel(L"Fidelity")
tight_layout()
gcf()
```

which show the error

```
ERROR: MethodError: no method matching master_dynamic(::StepRangeLen{âŠ}, ::Operator{âŠ}, ::typeof(URP123), ::Vector{âŠ})
Closest candidates are:
master_dynamic(::Any, ::Ket, ::Any...; kwargs...)
@ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:228
master_dynamic(::Any, ::Operator, ::AbstractTimeDependentOperator, ::Any; kwargs...)
@ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:220
master_dynamic(::Any, ::Operator, ::Any; rates, fout, kwargs...)
@ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:211
Stacktrace:
[1] top-level scope
@ d:\fqguo\software\Onedrive_education\OneDrive - nenu.edu.cn\Açș§éèŠ---çŹŹäșçŻèźșæçšćș\éȘèŻURP-NæŹĄæŒć\URP_unitary.jl:100
Some type information was truncated. Use `show(err)` to see complete types.
```

and another program

```
using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin
dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms
b = NLevelBasis(dim) # Basis for a single atom
s0 = nlevelstate(b,1) # |0> state
s1 = nlevelstate(b,2) # |1> state
sr = nlevelstate(b,3) # |r> state
plus = s0 # |+> state
minus = s1 # |-> state
# plus = (s0 + s1)/sqrt(2) # |+> state
# minus = (s0 - s1)/sqrt(2) # |-> state
Identity4 = identityoperator(b)
Hzeros = 0*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)
XZ = tensor(sigma_x, sigma_z, Identity4)
ZXZ = tensor(sigma_z, sigma_x, sigma_z)
ZX = tensor(Identity4, sigma_z, sigma_x)
Omega_1 = 2*pi*3*10^6; Omega_2 = 2*pi*3*10^6;
Delta = 20*Omega_1; U_rr = Delta;
NN = 6;
MM = 16;
t2 = pi/NN/Omega_2;
t1 = 7*t2;
t_pump = t1 + t2;
t_decay = 3.02*10^-6;
t_cycle = 3*NN*t_pump;
gamma_r = 2*pi*0*10^3;
gammar_P = 1/(26.2*10^-9);
Omega_pr = 0.2*gammar_P;
gammar_eff = Omega_pr^2/gammar_P;
gamma_0r = 2/3*gammar_eff;
gamma_1r = 1/3*gammar_eff;
Hv = U_rr*( tensor(srâdagger(sr), srâdagger(sr), Identity4) + tensor(srâdagger(sr), Identity4, srâdagger(sr)) + tensor(Identity4, srâdagger(sr), srâdagger(sr)) );
function URP123(t,args=nothing)
H0_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(s0) + Omega_1/2*exp(1im*Delta*t)*s0âdagger(sr) - (Omega_1/2)^2/Delta*(s0âdagger(s0) - srâdagger(sr));
H0_ome2 = 1im*Omega_2/2*srâdagger(s0) - 1im*Omega_2/2*s0âdagger(sr);
HP_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(plus) + Omega_1/2*exp(-1im*Delta*t)*plusâdagger(sr) - (Omega_1/2)^2/Delta*(plusâdagger(plus) - srâdagger(sr));
HP_ome2 = 1im*Omega_2/2*srâdagger(plus) - 1im*Omega_2/2*plusâdagger(sr);
H1_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(s1) + Omega_1/2*exp(-1im*Delta*t)*s1âdagger(sr) - (Omega_1/2)^2/Delta*(s1âdagger(s1) - srâdagger(sr));
H1_ome2 = 1im*Omega_2/2*srâdagger(s1) - 1im*Omega_2/2*s1âdagger(sr);
HM_ome1 = Omega_1/2*exp(-1im*Delta*t)*srâdagger(minus) + Omega_1/2*exp(-1im*Delta*t)*minusâdagger(sr) - (Omega_1/2)^2/Delta*(minusâdagger(minus) - srâdagger(sr));
HM_ome2 = 1im*Omega_2/2*srâdagger(minus) - 1im*Omega_2/2*minusâdagger(sr);
URP123_ome1 = tensor(H0_ome1, Identity4, Identity4) + tensor(Identity4, HP_ome1, Identity4) + tensor(Identity4, Identity4, H1_ome1) + Hv;
URP123_ome2 = tensor(H0_ome2, Identity4, Identity4) + tensor(Identity4, HP_ome2, Identity4) + tensor(Identity4, Identity4, H1_ome2) + Hv;
return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1)
end
psi_0 = tensor(s0,minus,s0);
psi_f = tensor(sr,minus,s0);
dm0 = psi_0âdagger(psi_0);
rho_final = psi_fâdagger(psi_f);
tlist_pump = range(0, stop=t_pump, length=1000);
exp_psi0 = [];
for jj = 1:NN
global tout, rho_XZ, dm0, psi_0
tout, rho_XZ = timeevolution.schroedinger_dynamic(tlist_pump, psi_0, URP123);
append!(exp_psi0, real(expect(rho_final, rho_XZ)))
psi_0 = rho_XZ[end];
end
plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(exp_psi0)
xlabel(L"T")
ylabel(L"Fidelity")
tight_layout()
gcf()
```