waveguide_evolution for time-dependent Hamiltonian

[phantomlsh]

Hi,

(This is just a minor suggestion)

It seems that waveguide_evolution does not work for the time-dependent Hamiltonian. I figured out that we can use timeevolution.schroedinger_dynamic to do the simulation.

Is it possible to make waveguide_evolution supporting the time-dependent Hamiltonian, so that the interface could be unified?

Thank you very much!

[mabuni1998]

Hi,

Thanks for the suggestion and pointing this out. Perhaps this should be addressed in the documentation. Schroedinger_dynamic is the way one would add time-dependency, and in fact, waveguide_evolution is just a wrapper for calling the solver, so its perfectly possible to do so.

I also want to add, that one can utilize TimeDependentSum from the QuantumOptics.jl package which is a very neat and efficient way to include time dependency with waveguide operators. If you are interested in an immediate solution, then you can see how it could be implemented in the following (here I use time dependency to simulate the system considered in: https://www.nature.com/articles/s41534-022-00604-5 )

using WaveguideQED
using QuantumOptics
using PyPlot

#parameters are the same as in fig. 4
Ωg =1
τg = 4 * log(2)/Ωg
g = 0.4 * Ωg
κc = 6*Ωg
Ω0 = 15*g
Tin = 4.3/g

times = 0:0.1:16/g
dt = times[2] - times[1]

#Create waveguide basis
bw = WaveguideBasis(1,times)

#Deine basis of cavity
bc = FockBasis(1)
#Two level system is just fockbasis cutoff at 1. One could also use spinbasis for three level system. 
be = FockBasis(1)

#Define operators
#Cavity A operators
a = identityoperator(be) ⊗ identityoperator(bc) ⊗ destroy(bc) ⊗ identityoperator(bw)
ad = identityoperator(be) ⊗ identityoperator(bc) ⊗ create(bc) ⊗ identityoperator(bw)
#Cavity B operators
b = identityoperator(be) ⊗ destroy(bc) ⊗ identityoperator(bc) ⊗ identityoperator(bw)
bd = identityoperator(be) ⊗ create(bc) ⊗ identityoperator(bc) ⊗ identityoperator(bw)
#Emitter operators
s =  destroy(be) ⊗ identityoperator(bc) ⊗ identityoperator(bc) ⊗ identityoperator(bw)
sm =  create(be) ⊗ identityoperator(bc) ⊗ identityoperator(bc) ⊗ identityoperator(bw)
#Waveguide operators
w = destroy(bw)
wd = create(bw)
#Waveguide-cavity operators
adw = identityoperator(be) ⊗ identityoperator(bc) ⊗ create(bc) ⊗ w
wda = identityoperator(be) ⊗ identityoperator(bc) ⊗ destroy(bc) ⊗ wd


#Define time independent part of Hamiltonian (from line 3 in npj)
H0 = im*sqrt(κc/dt)*( adw - wda ) + g*( b*sm + bd*s )


#Define input state function (gaussian)
function xi_in(t, t0, τg)
    sqrt(2/τg)*(log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/τg^2)
end

#Define input state with the onephoton state in waveguide and zero excitations in everything else.
psi_in = fockstate(be, 0) ⊗ fockstate(bc,0) ⊗ fockstate(bc,0) ⊗ onephoton(bw, xi_in,Tin,τg)


#Now we define the time-dependent part of the Hamiltonian. Which is the control field Λ(t) and the detuning Ω(t).

#Used to define control functions, not part of original parameters
Tcut1 = 6/g
Tcut2 = 10/g
t_control = 3/g


#Turn detuning off during middle of pulse
function Ω(t)
    if t < Tcut1 || t > Tcut2
        return Ω0
    else
        return 0
    end
end

#Define control field Λ(t) (not optimum just chosen as a gaussian) before and after loading.
function Λ(t)
    if t < Tcut1
        return κc/10*xi_in(t, 3/g, τg)/xi_in(τg, 3/g, τg*3)
    elseif t>Tcut2
        return κc/10*xi_in(t, 12/g, τg)/xi_in(τg, 3/g, τg*3)
    else
        return 0
    end
end

#define time-dependent part of Hamiltonian
H1 = TimeDependentSum(Ω => sm*s,Λ=>a*bd,Λ=>ad*b)

H = H1+H0

tout, ψ = timeevolution.schroedinger_dynamic(times, psi_in, H;d_discontinuities=times)
#Get the final state
psi_end = OnePhotonView(ψ[end])

psi_in_plot = OnePhotonView(psi_in)

#Plot results
fig, ax = subplots(1,1,figsize=(9,4.5))
ax.plot(times*g,τg * psi_in_plot,label="Input τg ξin(t)","k-")
ax.plot(times*g,τg * psi_end,label="Output τg ξout(t)","k--")
ax.plot(times*g,Λ.(times)/κc,label="Control field Λ(t)/κ","r-")
ax.plot(times*g,Ω.(times)/max(Ω.(times)...),label="Control field Ω(t)/Ω0","b-")
ax.set_xlabel("t [1/g]")
ax.set_ylabel("Normalized units")
ax.legend(fontsize=12)
tight_layout()

[phantomlsh]

Sorry for the delay. Yes I realized how to use the quantumoptics function, but it would be great if you can make them unified as one interface! Thank you very much!