How to measure the coincidence expectation value

[phantomlsh]

When doing simulations involving HOM dip, I am wondering if we can compute the expectation value for coincidence:

expect(wd1*w1 ⊗ wd2*w2, ψ)

Where ψ is a state with only 2 waveguides, basis defined by WaveguideBasis(2, 2, times) with proper ladder operators.

Currently, Julia throws the error for this expression:

ERROR: LoadError: BoundsError: attempt to access Int64 at index [2]
Stacktrace:
 [1] getindex(x::Int64, i::Int64)
   @ Base ./number.jl:98
 [2] tensor(a::WaveguideInteraction{WaveguideBasis{2, 2}, WaveguideBasis{2, 2}}, b::WaveguideCreate{WaveguideBasis{2, 2}, WaveguideBasis{2, 2}, 2, 2})
   @ WaveguideQED ~/.julia/packages/WaveguideQED/yljRX/src/WaveguideInteraction.jl:104

Is there a way to properly compute the expectation value?

Thank you very much!

[mabuni1998]

Hi again and sorry for the delay!

How exactly did you define wd1,w1, w2, and wd2? If you are using the 2 waveguide basis: WaveguideBasis(2,2,times), then the following should work:

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

bw = WaveguideBasis(2,2,times)
w1 = destroy(bw,1)
wd1 = create(bw,1)
w2 = destroy(bw,2)
wd2 = create(bw,2)
n1 = wd1*w1
n2 = wd2*w2
expval_op = n2*n1

ξfun(t1,t2,σ1,σ2,t0) = sqrt(2/σ1)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t1-t0)^2/σ1^2)*sqrt(2/σ2)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t2-t0)^2/σ2^2)
psi_0220 = 1/sqrt(2)*(twophoton(bw,1,ξfun,1,1,2.5) + twophoton(bw,2,ξfun,1,1,2.5)) 
psi_11 = twophoton(bw,1,2,ξfun,1,1,5) 

function expect_waveguide(expval_op, ρ0,times)
    ρ1 = copy(ρ0)
    expval = 0
    for i in eachindex(times)
        set_time!(expval_op, times[i])
        mul!(ρ1, expval_op, ρ0)
        expval += sum(ρ1.data)
    end
    return expval
end

expect_waveguide(expval_op, psi_0220,times) #Equals 0
expect_waveguide(expval_op, psi_11,times) #Equals 1

I am here defining a twophoton state that either resides as a |02>+|20> photon state or in a |11> photon state. As expected, your operator only gives 1 for the later state, as the first state has no simultaneous photons in both waveguides.

From our interactions, I learned that using expect() to define the expectation value of the whole waveguide is not consistent when the user has defined many different operators. I am, therefore, working on an update to introduce the function expect_waveguide(), which I defined for you in the above. From now on, this operator should be used when wanting expectation values of waveguide states / operators. I also noticed some inconsistencies when trying to combine operators the way you did through (WaveguideInteraction times WaveguideOperator), which I am also working on fixing. For now, the above should solve the problem for you.

best Matias

[mabuni1998]

A correction to the above is that probably what you are interested in when you write: expect(wd1*w1 ⊗ wd2*w2, ψ) is something like: $\sum_k \sum_l \langle \psi | w_1(t_k)^\dagger w_1(t_k) w_2(t_l)^\dagger w_2(t_l) |\psi \rangle$, basically asking given I have a photon in waveguide 1 at time t_k what is the probability of also having it in waveguide 2 at time t_l and then summing over all of these. Again, expectation values are tricky, and I think they are best defined on a case-to-case basis. For the above, one would do:

function expect_waveguide2(O,psi,times)
    psi_c = copy(psi)
    expval = 0
    for i in eachindex(times)
        set_waveguidetimeindex!(n1,i)
        for j in eachindex(times)
            set_waveguidetimeindex!(n2,j)
            mul!(psi_c, O, psi)
            expval += dot(psi_c.data,psi.data)
        end
    end
    return expval
end

[phantomlsh]

This is very useful. Thank you very much! I realized that the sum in expect_waveguide might be a dot product, as in expect_waveguide2. And after the fix with the dot product, expect_waveguide behaviors similar to expect.

I do have a further question about constructing an arbitrary state. Is it possible to construct a waveguide state as the following: We have 2 waveguides defined by their ladder operators $a^\dagger(t), b^\dagger(t)$ and we have two photons in the state output from a 50-50 beamsplitter, the state should be

$$|\psi\rangle = \int\int dt_1 dt_2 \frac{1}{2}\xi(t_1,t_2)[a^\dagger(t_1)a^\dagger(t_2) - b^\dagger(t_1)b^\dagger(t_2) + b^\dagger(t_1)a^\dagger(t_2) - a^\dagger(t_1)b^\dagger(t_2)]|0\rangle$$

I tried to use twophoton functions to construct this state but I am not very confident about my expression, as the normalization seems to be broken when the wavefunction for the two incident photons are the same $\xi(t_1,t_2) = \xi(t_1)\xi(t_2)$

bw = WaveguideBasis(2, 2, times)
ξ(t, t0=5, σ=1) = complex(√(2/σ)*(log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ξ1(t1, t2, t01, t02) = ξ(t1, t01) * ξ(t2, t02)
ξ2(t1, t2, t01, t02) = ξ(t2, t01) * ξ(t1, t02)
t01 = 5
t02 = 5 + τ
ψ = twophoton(bw, 1, ξ1, t01, t02)/2 - twophoton(bw, 2, ξ1, t01, t02)/2 + twophoton(bw, [1, 2], ξ2, t01, t02)/2 - twophoton(bw, [1, 2], ξ1, t01, t02)/2
expval1 = expect(wd1*w1, ψ) # This is 1 for large τ, but drops to 0.5 for τ = 0, which is confusing.

Would you have any immediate suggestions? Thank you very much!

[mabuni1998]

I see... The problem is indeed with normalization in the above example. I'm not sure I immediately understand what goes wrong, but a quick fix is the following:

ψ = twophoton(bw, 1, ξ1, t01, t02)/2 - twophoton(bw, 2, ξ1, t01, t02)/2 + twophoton(bw, [1, 2], ξ2, t01, t02)/2 - twophoton(bw, [1, 2], ξ1, t01, t02)/2
ψ = normalize(ψ)

I will have to look into how multiple twophoton creation objects are combined correctly, to give it a more permanent fix.

[phantomlsh]

Hi,

I am wondering if there is any update on this. I am thinking about doing projective measurements on such 2-waveguide superposition state, for example, evaluating expectation value for the operator:

$$P = \int_0^T a_1^\dagger |0\rangle\langle0| a_1$$

Would you have a chance to provide an example of a beamsplitter simulation? I guess that would be super useful for other users as well.

Thank you very much!

[mabuni1998]

Okay, so there are many details of this that need to be addressed, and I'm not entirely sure how to deal with all of them.

First of all, can you provide more details on the operator P? I guess there should be some time dependence on $a_1^\dagger$ and $a_1$ that depends on dt in the integral? This could help me create the example you suggested, which I agree with you would be very helpful.

EDIT: Disregard the following, I didn't account for the integral in the below which solves the problem. I will get back with an update once I have time.

Secondarily, the normalization problem is not trivial to solve since the code assumes the wavefunction (of two photons in the same waveguide) to be symmetric under exchange of two times: $\xi^{(2)}(t_1,t_2)=\xi^{(2)}(t_2,t_1)$. As far as I know this is a standard requirement for indistinguishable particles. However, I am actually now unsure whether this is necessarily a true assumption. For example, take the example that you are describing: We have two photons that start in waveguides 1 and 2, respectively. They each have an associated wavefunction $\xi_1(t)$ and $\xi_2(t)$, respectively, which does not have to be the same. If we combine them on a beamsplitter we get the state:

$$|\psi\rangle = \int\int dt_1 dt_2 \frac{1}{2}\xi_1(t_1)\xi_2(t_2)[a^\dagger(t_1)a^\dagger(t_2) - b^\dagger(t_1)b^\dagger(t_2) + b^\dagger(t_1)a^\dagger(t_2) - a^\dagger(t_1)b^\dagger(t_2)]|0\rangle$$

The wavefunction for the two-photon state in, e.g., waveguide 1 does now only obey the symmetry relation if $\xi_1(t) = \xi_2(t)$. But as the code is written right now, it is assumed that the wavefunctions of two-photon states in the same waveguide always obey the symmetry relation. I think this might be where the strange results you are seeing are coming from. If you have any inputs on the above reasoning I would also be interested to hear them. I guess one can interpret $\xi_1(t) \neq \xi_2(t)$ as having to be DISTINQUISHABLE particles. If the above reasoning is sound, then I might have to change some parts of the code to allow for non-symmetric wavefunctions. I am unsure how big of a task this will be and will look into this. After this I will work towards providing the above example.

Thank you for your patience.

[phantomlsh]

Yes, sorry for the delay. The projective measurement operator should be

$$P_a = \int_0^T dt a^\dagger(t) |0\rangle\langle0| a(t) \qquad P_b = \int_0^T dt b^\dagger(t) |0\rangle\langle0| b(t) $$

where $a$ and $b$ are lowering operators for two waveguides. We are interested in the coincidence whose expected value is computed by $\langle\psi_{ab}|P_aP_b|\psi_{ab}\rangle$.

Interestingly, I recently evaluated the expectation value $n_1 * n_2$ on the beamsplitter state as described above. If I call the normalize function, the expectation value of the photon number in one waveguide ($n_1$) will be constantly 1 (good), but the expectation value $n_1 * n_2$ will not agree with the theoretical value of coincidence, especially when two photons overlap. However, if I do not call the normalize function, the expectation value of $n_1$ is not correct with overlapping photons, but the expectation value $n_1 * n_2$ perfectly agrees with the theoretical result.

Thank you very much!

[mabuni1998]

Okay, perfect, thank you! Please see the example I provided in the documentation under the tab Beamsplitter, hopefully this is what you requested.

I fixed a few bugs in this update, so be sure to update waveguideQED. The biggest was how I initialized the two-photon waveguide state. Something is still funny about normalization, but for now, your fix with normalizing also works. See this:

using WaveguideQED
ξ(t, t0=5, σ=1) = complex(√(2/σ)*(log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ξ1(t1, t2, t01, t02) = ξ(t1, t01) * ξ(t2, t02)
times = 0:0.1:10
using PyPlot
pygui(true)
taus = 0:1:4
expvals = zeros(length(taus))
bw = WaveguideBasis(2, 2, times)
w1 = destroy(bw,1)
wd1 = create(bw,1)

w2 = destroy(bw,2)
wd2 = create(bw,2)


n1 = wd1*w1
n2 = wd2*w2
expval_op = n1*n2
using LinearAlgebra
function expect_waveguide2(O,psi,times)
  psi_c = copy(psi)
  expval = 0
  for i in eachindex(times)
      set_waveguidetimeindex!(n1,i)
      for j in eachindex(times)
          set_waveguidetimeindex!(n2,j)
          mul!(psi_c, O, psi)
          expval += dot(psi_c.data,psi.data)
      end
  end
  return expval
end

t01=5
for (i, τ) in enumerate(taus)
    t02 = 5 + τ
    ψ = twophoton(bw,1,ξ1, t01, t02,norm=true)/2 - twophoton(bw, 2, ξ1, t01, t02,norm=true)/2 + twophoton(bw, [2,1], ξ1, t01, t02,norm=true)/2 - twophoton(bw, [1,2], ξ1, t01, t02,norm=true)/2
    ψ = normalize(ψ)
    expvals[i] = expect_waveguide2(expval_op,ψ,times)
    println("Norm: ", norm(ψ))
end

fig, ax = subplots(1,1, figsize=(6,4))
ax.plot(taus, expvals)

This works now, but as I show in the documentation, one could avoid this problem totally by performing the beamsplitter operation by letting the waveguide state evolve under the Hamiltonian: $H = V(w_1^\dagger(t) w_2(t) + w_2^\dagger(t) w_1(t) )$. Please let me know if you have further questions.

[phantomlsh]

Hi,

Thank you for the helpful example! The simulation with Hamiltonian fixes our problem. But we do want to mention a potential problem here. As demonstrated in the following code, when using the constructed state, the unnormalized state has the same coincidence as the evolved state (but the wrong expectation value of the photon number in a single waveguide), while the normalized state fixes the expectation value but breaks the coincidence. The plot is attached for reference.

Maybe there's a problem caused by the state construction? Thank you for your effort!

using QuantumOptics
using WaveguideQED
using LinearAlgebra
using LaTeXStrings
using Plots

times = 0:0.1:10
dt = times[2] - times[1]
bw = WaveguideBasis(2, 2, times)

wd1 = create(bw, 1)
w1 = destroy(bw, 1)
wd2 = create(bw, 2)
w2 = destroy(bw, 2)
n1 = wd1*w1
n2 = wd2*w2

ξ(t, t0=5, σ=1) = complex(√(2/σ)*(log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ξ1(t1, t2, t01, t02) = ξ(t1, t01) * ξ(t2, t02)
ξ2(t1, t2, t01, t02) = ξ(t2, t01) * ξ(t1, t02)

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

H = π/4/dt*(wd1*w2 + wd2*w1)

values = [[], [], []]

τs = -3:0.1:3
for τ ∈ τs
  t01 = 5
  t02 = 5 + τ
  ψ1 = waveguide_evolution(times, twophoton(bw, [1, 2], ξ1, t01, t02), H)
  expval1 = expect_coincidence(n1, n2, ψ1, times) # by evolution, agree with analytical result
  ψ2 = twophoton(bw, 1, ξ1, t01, t02, norm=true)/2 - twophoton(bw, 2, ξ1, t01, t02, norm=true)/2 + twophoton(bw, [2, 1], ξ1, t01, t02, norm=true)/2 - twophoton(bw, [1, 2], ξ1, t01, t02, norm=true)/2
  expval2 = expect_coincidence(n1, n2, ψ2, times) # by construction, unnormalized, agree with analytical result
  ψ3 = twophoton(bw, 1, ξ1, t01, t02, norm=true)/2 - twophoton(bw, 2, ξ1, t01, t02, norm=true)/2 + twophoton(bw, [2, 1], ξ1, t01, t02, norm=true)/2 - twophoton(bw, [1, 2], ξ1, t01, t02, norm=true)/2
  ψ3 = normalize(ψ3)
  expval3 = expect_coincidence(n1, n2, ψ3, times) # by construction, normalized, DISAGREE with analytical result
  append!(values[1], expval1)
  append!(values[2], expval2)
  append!(values[3], expval3)
end

plot(τs, values[1], label="Evolution", lw=3, ylim=[0, 1.02], title=latexstring("Coincidence (\$\\eta = 1/2\$)"), xlabel="Pulse Delay Time (τ)", ylabel="Expectation Value")
plot!(τs, values[2], label="Unnormalized")
plot!(τs, values[3], label="Normalized")
savefig("../plots/coincidence.pdf")

plot.pdf

[mabuni1998]

There was indeed a problem with state construction. Forcing normalization on construction (norm=true by default) was causing the combined state ψ2 and ψ3 to not be normalized. Now, two-photon states are initialized without a prefactor of 1/sqrt(2) and are not, per default, normalized. This makes the above example work (see the updated documentation: https://qojulia.github.io/WaveguideQED.jl/dev/beamsplitter/) because now ψ2 and ψ3 do not need to be normalized as they are normalized from construction (remember to not enforce normalization on each individual term, that is: norm=false). I hope this helps!

[phantomlsh]

Thank you very much for your work! With [email protected], all issues regarding this are fixed!