Question about Delayed Operators

[phantomlsh]

Hi,

Sorry to disturb you again. I might misunderstand the principle behind the delayed operators. In the following program, I tried to simulate a system with two qubits (2-level cavity FockBasis(1)) connected by a waveguide, where the waveguide should have a finite length, causing the traveling time between the qubits to be $\Delta$.

For simplicity, let's assume that the couplings are all constant. With the only initial excitation in qubit 1 (Q1), I wrote the following simulation code:

using QuantumOptics
using WaveguideQED
using LinearAlgebra
using Plots

times = 0:0.1:10
dt = times[2] - times[1]

bw = WaveguideBasis(1, times)
bc = [FockBasis(1), FockBasis(1)] # 2 qubits
Iw = identityoperator(bw)
Ic = identityoperator.(bc)

Δ = 5
wd = create(bw)
wdΔ = create(bw; delay=Δ/dt)
w = destroy(bw)
wΔ = destroy(bw; delay=Δ/dt)
ad = create.(bc)
a = destroy.(bc)

γ = 5
H = im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗w - a[1]⊗Ic[2]⊗wd) +
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗w - Ic[1]⊗a[2]⊗wd) +
    im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗wΔ - a[1]⊗Ic[2]⊗wdΔ) + 
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗wΔ - Ic[1]⊗a[2]⊗wdΔ)  

ncc = (ad[1]*a[1])⊗(ad[2]*a[2])⊗Iw
nc1 = (ad[1]*a[1])⊗Ic[2]⊗Iw
nc2 = Ic[1]⊗(ad[2]*a[2])⊗Iw
nw = Ic[1]⊗Ic[2]⊗(wd*w)
function expval(time, psi)
  (expect(ncc, psi), expect(nc1, psi), expect(nc2, psi), expect_waveguide(nw, psi))
end

ψ = fockstate(bc[1], 1) ⊗ fockstate(bc[2], 0) ⊗ zerophoton(bw)
ψ, ecc, ec1, ec2, ew = waveguide_evolution(times, ψ, H, fout=expval)
plot!(times, abs.(ec1); label="Qubit 1", xlabel="Time", ylabel="Expectation Value")
plot!(times, abs.(ec2); label="Qubit 2")
plot!(times, abs.(ew); label="Waveguide")
savefig("./plot.pdf")

This program yields the following plot:

which is confusing because qubit 2 (Q2) rises at t=0. The photon released by Q1 should travel for time $\Delta$ before getting to Q2. I suspect that the Hamiltonian I wrote is not correct. Would you offer some clues about simulating such a system (2 qubits connected by a waveguide with some length)?

Thank you very much!

[mabuni1998]

Hi,

Very exciting that you are testing out delayed operators! I realize that perhaps the documentation is misleading in the way one should think about the delayed operators. It is kind of counterintuitive, but when you have a delayed operator, it acts in the future, meaning, let's say, 5 timesteps AHEAD. Thus, in your example, you actually want to start the initial state in qubit 2 since this acts in the future (more forward on the conveyor belt), and qubit 1 will then 5 timesteps later see this spot on the conveyor belt. This is not really intuitive and not obvious from the documentation, but I'm currently thinking about how it could be illustrated better.

Another comment is the Hamiltonian. You have too many terms much? The way you write it, both qubits talk with the waveguide at both times. Instead, you should only have qubit one talk with the waveguide at a "not delayed time" and qubit two at the delayed time. Fixing this, I have:

using QuantumOptics
using WaveguideQED
using LinearAlgebra
using Plots

times = 0:0.1:10
dt = times[2] - times[1]

bw = WaveguideBasis(1, times)
bc = [FockBasis(1), FockBasis(1)] # 2 qubits
Iw = identityoperator(bw)
Ic = identityoperator.(bc)

Δ = 5
wd = create(bw)
wdΔ = create(bw; delay=Δ/dt)
w = destroy(bw)
wΔ = destroy(bw; delay=Δ/dt)
ad = create.(bc)
a = destroy.(bc)

γ = 5
H = im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗w - a[1]⊗Ic[2]⊗wd) +
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗wΔ - Ic[1]⊗a[2]⊗wdΔ)  
ncc = (ad[1]*a[1])⊗(ad[2]*a[2])⊗Iw
nc1 = (ad[1]*a[1])⊗Ic[2]⊗Iw
nc2 = Ic[1]⊗(ad[2]*a[2])⊗Iw
nw = Ic[1]⊗Ic[2]⊗(wd*w)
function expval(time, psi)
  (expect(ncc, psi), expect(nc1, psi), expect(nc2, psi), expect_waveguide(nw, psi))
end

ψ = fockstate(bc[1], 0) ⊗ fockstate(bc[2], 1) ⊗ zerophoton(bw)
ψ, ecc, ec1, ec2, ew = waveguide_evolution(times, ψ, H, fout=expval)

Plots.CURRENT_PLOT.nullableplot = nothing
plot!(times, abs.(ec1); label="Qubit 1", xlabel="Time", ylabel="Expectation Value")
plot!(times, abs.(ec2); label="Qubit 2")
plot!(times, abs.(ew); label="Waveguide")

savefig("./plot.png")

which gives me:

Obviously, the above configuration of the operators makes it so that after qubit 1 reemits into the waveguide qubit 2 will not pick up the photon again. For this to happen, one should add another term with an operator that is two times delayed: wd2Δ = create(bw; delay=2 * Δ/dt). But this quickly grows impractical when you then want to consider a third bounce. Instead, it should be possible to code the waveguide to loop around on itself. Would this be interesting?

[phantomlsh]

Thank you for your reply! That makes perfect sense. Yes, I am very interested in all kinds of looping/connecting systems. One further question on the system mentioned above is how to simulate the system with initial excitation in both qubits? Namely, to simulate a symmetric two-photon system (Q1 - W - Q2).

We are particularly interested in this because our ultimate goal is to simulate the full HOM dip setup, i.e. a system like

Q1 - W1 - BS - W2 - Q2

where Q is qubit, W is waveguide, BS is beamsplitter. Ideally, the waveguides should have a finite length. I find it particularly tricky with the delay operators to simulate two photons in one system. Is there any clever method to build a system like this?

Thank you very much!

[mabuni1998]

Alright, super cool! I will look into looping. I almost have it done now, but I will do some testing and see if I can make the time ordering and so on more intuitive. I expect to have an update at the start of the following week.

For your question about two initially excited qubits, I would simply just have a waveguide with a maximum of two photons, e.g., WaveguideBasis(2,times). I just tested this, and it should work. Please let me know if problems arise.

For the final configuration Q1 - W1 - BS - W2 - Q2, I'm guessing that you have a detector after the beamsplitter? So that Q1 and Q2 are inputs on one side of the BS and then you detect? In that case, have you tried something like a waveguide basis with two waveguides and then Q1 coupled to W1 and Q2 to W2, and then via delay operators, you place the beamsplitter Hamiltonian $H_bs = V(w_1^\dagger w_2 + w_1 w_2^\dagger)$ a certain distance from the two qubits? I don't think this is an easy exercise, and I would suggest you draw a conveyor belt and try to visualize how the system evolves as the conveyor belt moves forward. I think in this instance, your solution is simply having H_bs being delayed by $\Delta \tau$ where $\Delta \tau$ is the time it takes the light to travel from Q1 and Q2 to the BS.

[phantomlsh]

Regarding the simple system (Q1 - W - Q2) with two photons, how to write the Hamiltonian? The question is that if we write a symmetric Hamiltonian regarding both qubits to release their photons at t=0, we might write

H = im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗w - a[1]⊗Ic[2]⊗wd) +
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗w - Ic[1]⊗a[2]⊗wd) +
    im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗wΔ - a[1]⊗Ic[2]⊗wdΔ) + 
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗wΔ - Ic[1]⊗a[2]⊗wdΔ)

then the released photon from Q1 might interact with Q2 immediately without a delay?

Regarding the beamsplitter settings, our physical experiment would capture the photons after the beamsplitter. The resulting qubit states serve as the "detection". I understand that this scenario is more complicated, and maybe I should first understand the simple system discussed above.

Thank you very much!

[mabuni1998]

Ahh, now I see the confusion. The Hamiltonian should not change, but instead, the waveguides should loop. I will try to address this case of two photons in an update sometime next week, where I introduce the looping (which is necessary). Hopefully, this will help.

As for the beamsplitter experiment, I think this could already be implemented now using a waveguide basis with two waveguides, by having the BS Hamiltonian be delayed and the two qubits interacting with W1 and W2, respectively.

[mabuni1998]

I have just released an update 0.2.6, where I have implemented the above mentioned looping. I suggest you take a look at the example I provide in the documentation: https://qojulia.github.io/WaveguideQED.jl/dev/time_delay/#emitters
Please let me know if you have any questions.

[phantomlsh]

That's great! The QWQ model works completely fine. Thank you very much!

Regarding the Q1 - W1 - BS - W2 - Q2 model, there might be a glitch(?). The following code will raise error when solving the waveguide evolution:

dt = 0.1
times = 0:dt:10+dt
times_sim = 0:dt:20
bw = WaveguideBasis(2, 2, times)
bc = [FockBasis(2), FockBasis(2)]

Δ1 = 5
Δ2 = 5
Iw, Ic = identityoperator(bw), identityoperator.(bc)
wd = [create(bw, 1), create(bw, 2)]
wdΔ = [create(bw, 1; delay=Δ1/dt+1), create(bw, 2; delay=Δ2/dt+1)]
w = [destroy(bw, 1), destroy(bw, 2)]
wΔ = [destroy(bw, 1; delay=Δ1/dt+1), destroy(bw, 2; delay=Δ2/dt+1)]
ad = create.(bc)
a = destroy.(bc)

n = 1/2
γ = 5
H = im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗wΔ[1] - a[1]⊗Ic[2]⊗wdΔ[1]) +
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗wΔ[2] - Ic[1]⊗a[2]⊗wdΔ[2]) +
    asin(√n)/dt*(Ic[1]⊗Ic[2]⊗(wd[1]*w[2]) + Ic[1]⊗Ic[2]⊗(wd[2]*w[1]))

ncc = (ad[1]*a[1])⊗(ad[2]*a[2])⊗Iw
nc1 = (ad[1]*a[1])⊗Ic[2]⊗Iw
nc2 = Ic[1]⊗(ad[2]*a[2])⊗Iw
nw1 = Ic[1]⊗Ic[2]⊗(wd[1]*w[1])
nw2 = Ic[1]⊗Ic[2]⊗(wd[2]*w[2])
function expval(time, psi)
  expect(ncc, psi), expect(nc1, psi), expect(nc2, psi), expect_waveguide(nw1, psi), expect_waveguide(nw2, psi)
end

ψ = fockstate(bc[1], 1) ⊗ fockstate(bc[2], 0) ⊗ zerophoton(bw)
ψ, ecc, ec1, ec2, ew1, ew2 = waveguide_evolution(times_sim, ψ, H, fout=expval)

This raises error:

ERROR: LoadError: BoundsError: attempt to access 21115-element view(::UnsafeArrays.UnsafeArray{ComplexF64, 1}, 1:9:190027) with eltype ComplexF64 at index [21116]

This error only appears when times_sim exceeds the range of times and the last line in Hamiltonian (interaction between waveguides) exists, so I believe that it is related to waveguide loop. I am not sure where the problem is.

Thank you very much!

[mabuni1998]

Hi,

Yes, I'm sorry, I didn't include looping in waveguide interaction operators. I have now fixed this in the update 0.2.7, and the above should work. Below, I also plot the evolution of the above system, where I included the coincidence counts in the two waveguides. It is clear that after the beamsplitter, there are no coincidence counts; thus, there is a full HOM dip. Note that calculating the coincidence counts is quite expensive to compute, and thus, calculating it at each step makes the code take quite a while to run. In reality, just computing the coincidence counts in the end is fine, and this should not be too expensive.

using WaveguideQED
using QuantumOptics

dt = 0.1
times = 0:dt:10+dt
times_sim = 0:dt:10
bw = WaveguideBasis(2, 2, times)
bc = [FockBasis(2), FockBasis(2)]

Δ1 = 2
Δ2 = 2
Iw, Ic = identityoperator(bw), identityoperator.(bc)
wd = [create(bw, 1), create(bw, 2)]
wdΔ = [create(bw, 1; delay=Δ1/dt+1), create(bw, 2; delay=Δ2/dt+1)]
w = [destroy(bw, 1), destroy(bw, 2)]
wΔ = [destroy(bw, 1; delay=Δ1/dt+1), destroy(bw, 2; delay=Δ2/dt+1)]
ad = create.(bc)
a = destroy.(bc)

n = 1/2
γ = 5
H = im*√(γ/2/dt)*(ad[1]⊗Ic[2]⊗wΔ[1] - a[1]⊗Ic[2]⊗wdΔ[1]) +
    im*√(γ/2/dt)*(Ic[1]⊗ad[2]⊗wΔ[2] - Ic[1]⊗a[2]⊗wdΔ[2]) +
    pi/4/dt*(Ic[1]⊗Ic[2]⊗(wd[1]*w[2]) + Ic[1]⊗Ic[2]⊗(wd[2]*w[1]))

    

ψ = fockstate(bc[1], 1) ⊗ fockstate(bc[2], 1) ⊗ zerophoton(bw)

using LinearAlgebra
psi_c = copy(ψ)
function expect_coincidence(n1, n2, psi, times)
  expval = 0
  for i in eachindex(times)
    set_waveguidetimeindex!(n1, i)
    for j in eachindex(times)
      set_waveguidetimeindex!(n2, j)
      QuantumOpticsBase.mul!(psi_c, n1*n2, psi)
      expval += dot(psi_c.data, psi.data)
    end
  end
  return expval
end

ncc = (ad[1]*a[1])⊗(ad[2]*a[2])⊗Iw
nc1 = (ad[1]*a[1])⊗Ic[2]⊗Iw
nc2 = Ic[1]⊗(ad[2]*a[2])⊗Iw
nw1 = Ic[1]⊗Ic[2]⊗(wd[1]*w[1])
nw2 = Ic[1]⊗Ic[2]⊗(wd[2]*w[2])
function expval(time, psi)
  expect(ncc, psi), expect(nc1, psi), expect(nc2, psi), expect_waveguide(nw1, psi), expect_waveguide(nw2, psi),expect_coincidence(nw1,nw2,psi,times)
end
#ψ, ecc, ec1, ec2, ew1, ew2,coincedence = waveguide_evolution(times_sim, ψ, H, fout=expval)
psi = waveguide_evolution(times_sim, ψ, H, fout=1)

out = zeros(ComplexF64,6, length(psi))
for (i,p) in enumerate(psi)
  out[1,i] = expect(ncc, p)
  out[2,i] = expect(nc1, p)
  out[3,i] = expect(nc2, p)
  out[4,i] = expect_waveguide(nw1, p)
  out[5,i] = expect_waveguide(nw2, p)
  out[6,i] = expect_coincidence(nw1,nw2,p,times)
end

using PyPlot
pygui(true)
fig,ax = subplots(1,2,figsize=(9,4))
ax[1].plot(times, out[2,:], label="Q1")
ax[1].plot(times, out[3,:], label="Q2")

ax[2].plot(times, out[4,:], label="W1")
ax[2].plot(times, out[5,:], label="W2")
ax[2].plot(times, out[6,:], label="Coincidence")
ax[2].legend()

Giving the plot:

Note that it might be possible to just simulate the two qubits emitting their photons in one simulation and then take the output state of the first simulation and apply the beamsplitter operation in a subsequent simulation.