Multiple Waveguides

In the tutorial, we considered one waveguide containing one or two waveguides. In WaveguideQED, we can also consider problems where we have multiple waveguides or equivalently multiple waveguide modes. We can create a waveguide basis containing two photons and two waveguides with:

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

julia> times = 0:0.1:100.0:0.1:10.0
julia> bw = WaveguideBasis(2,times)┌ Warning: Note that for only one waveguide the lengths keyword is not needed since the length of the waveguidebins will be set to the length of the times vector.
└ @ WaveguideQED ~/work/WaveguideQED.jl/WaveguideQED.jl/src/basis.jl:35
WaveguideBasis{2, 1}([5253], 5252, 0, 101, 0.1, [101])
julia> Btotal = bw ⊗ bw[WaveguideBasis{2, 1}([5253], 5252, 0, 101, 0.1, [101]) ⊗ WaveguideBasis{2, 1}([5253], 5252, 0, 101, 0.1, [101])]

julia> bw = WaveguideBasis(2,2,times)WaveguideBasis{2, 2}([20706], 20705, 0, 101, 0.1, [101, 101])

We can then create the annihilation and creation operators for each waveguide by providing an index to destroy and create.

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

Here wd1 creates a photon in waveguide mode $1$: $w_{k,1}^\dagger \ket{\emptyset}_1 \ket{\emptyset}_2 = \ket{1_k}_1 \ket{\emptyset}_2$, and equivalently wd2: $w_{k,2}^\dagger \ket{\emptyset}_1 \ket{\emptyset}_2 = \ket{\emptyset}_1 \ket{1_k}_2$. $\ket{1_k}_i$ here denotes an excitation in time-bin $k$ in waveguide $i$.

Similarly, if we want to create a one-photon state in waveguide 1 or 2 we give onephoton an extra waveguide index:

ξ(t,σ,t0) = sqrt(2/σ)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2)
w_idx = 1
ψ_single_1 = onephoton(bw,w_idx,ξ,2,5)

w_idx = 2
ψ_single_2 = onephoton(bw,w_idx,ξ,2,5)

where ψ_single_1 here denotes the state: $\ket{\psi}_1 = \sum_k \sqrt{\Delta t} \xi(t_k) w_{k,1}^\dagger \ket{\emptyset} \otimes \ket{\emptyset}$ and ψ_single_2 equivalently $\ket{\psi}_2 = \ket{\emptyset} \otimes \sum_k \sqrt{\Delta t} \xi(t_k) w_{k,2}^\dagger \ket{\emptyset}$. We can also define two-photons states with:

ξ2(t1,t2,σ,t0) = ξ(t1,σ,t0)*ξ(t2,σ,t0)
w_idx = 1
ψ_double_1 = twophoton(bw,w_idx,ξ2,2,5)/sqrt(2)

w_idx = 2
ψ_double_2 = twophoton(bw,w_idx,ξ2,2,5)/sqrt(2)

We can even describe a simultaneous excitation in both waveguide states like $\ket{1_i}_1\ket{1_j }_2$ by providing two indices:

ψ_single_1_and_2 = twophoton(bw,[1,2],ξ2,2,5)

Note the absence of $1/\sqrt{2}$ when defining a state with one photon in each waveguide.

Scattering on two-level system

As an example of a system with multiple waveguide modes, we consider a two-level system coupled with two directional waveguide modes; a left and a right propagating mode. A sketch of the system can be seen below here:

Such a system is also considered experimentally in Le Jeannic, et al. Nat. Phys. 18, 1191–1195 (2022) and we here reproduce the theoretical results shown in the paper. The emitter here couples equally to the two waveguide modes with a total rate of $\gamma$, the Hamiltonian is here:

\[H = i \sqrt{\gamma/2/\Delta t} ( a^\dagger w_1 - a w_1^\dagger) + i \sqrt{\gamma/2/\Delta t} ( a^\dagger w_2 - a w_2^\dagger)\]

We can define this Hamiltonian, by creating the basis of the emitter and combining it with the waveguide operators:

be = FockBasis(1)
a = destroy(be)
ad = create(be)
γ = 1
H = im*sqrt(γ/2/dt)*( ad ⊗ w1 - a ⊗ wd1  ) + im*sqrt(γ/2/dt)*( ad ⊗ w2 - a ⊗ wd2  )

In the above, we denote waveguide $1$ as the right propagating mode and waveguide $2$ as the left propagating mode.

We can now study how single or two-photon states scatter on the emitter. We define the initial one-photon or two-photon Gaussian state and solve it using the defined Hamiltonian:

ξ₁(t1,σ,t0) = sqrt(2/σ)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t1-t0)^2/σ^2)
ξ₂(t1,t2,σ,t0) = ξ₁(t1,σ,t0) * ξ₁(t2,σ,t0)
w = 1
t0 = 5
ψ1 = fockstate(be,0) ⊗ onephoton(bw,1,ξ₁,w,t0)
ψ2 = fockstate(be,0) ⊗ twophoton(bw,1,ξ₂,w,t0)/sqrt(2)
ψScat1 = waveguide_evolution(times,ψ1,H)
ψScat2 = waveguide_evolution(times,ψ2,H)

Viewing the scattered states is then done using TwoPhotonView and the index for the corresponding waveguide. Giving two indices returns the combined single photon state in both waveguides $\sum_{j,k} \ket{1_j}_1 \ket{1_k}_2$:

ψ2RightScat = TwoPhotonView(ψScat2,1,[1,:])
ψ2LeftScat = TwoPhotonView(ψScat2,2,[1,:])
ψ2LeftRightScat = TwoPhotonView(ψScat2,2,1,[1,:])

We now want to compare the difference between a two-photon state scattering on the emitter and a single-photon state scattering. The scattered two-photon wavefunction can be expressed as $\xi^{(2)}(t_1,t_2) = \xi^{(1)}(t_1)\xi^{(1)}(t_2) + N(t_1,t_2)$, where $\xi^{(1)}(t)$ is the scattered single-photon state and $N(t_1,t_2)$ is the bound-state or non-linear contribution of the interaction with the emitter (stimulated emission). Thus, we compare the scattered two-photon state with the product of the scattered single-photon state. We thus calculate the two-time scattered distribution as:

ψ1RightScat = zeros(ComplexF64,(length(times),length(times)))
ψ1LeftScat = zeros(ComplexF64,(length(times),length(times)))
ψ1LeftRightScat = zeros(ComplexF64,(length(times),length(times)))
ψ1Right = OnePhotonView(ψScat1,1,[1,:])
ψ1Left = OnePhotonView(ψScat1,2,[1,:])

for i in eachindex(times)
    for j in eachindex(times)
        ψ1RightScat[i,j] = ψ1Right[i]*ψ1Right[j]
        ψ1LeftScat[i,j] = ψ1Left[i]*ψ1Left[j]
        ψ1LeftRightScat[i,j] = ψ1Left[i]*ψ1Right[j]
    end
end

Finally, we can plot the scattered wavefunctions, and we note that this matches Fig. 3 in Ref.^[1]:

fig,axs = subplots(3,2,figsize=(6,9))
plot_list = [ψ2RightScat,ψ2LeftScat,ψ2LeftRightScat,ψ1RightScat,ψ1LeftScat,ψ1LeftRightScat]
for (i,ax) in enumerate(axs)
    plot_twophoton!(ax,plot_list[i],times)
end
axs[1].set_ylabel("\$C^{RR}\$ \n t2 [a.u]")
axs[2].set_ylabel("\$C^{LL}\$ \n t2 [a.u]")
axs[3].set_ylabel("\$C^{LR}\$ \n t2 [a.u]")
axs[3].set_xlabel("t1 [a.u]")
axs[6].set_xlabel("t1 [a.u]")
plt.tight_layout()

If we consider the single photon state, we can also visualize the temporal evolution as:

Different lengths

It is also possible to consider waveguides of different effective lengths. This can be useful when modeling systems with feedback loops or delays, where photons travel different distances before interacting again with the system.

The waveguide length can be specified through the optional keyword argument lengths when constructing the WaveguideBasis. Each entry corresponds to the number of time bins stored for that waveguide mode.

As an example, consider a system consisting of two cavities connected through a circular waveguide loop. Photons emitted from one emitter propagate through the waveguide and return after a delay determined by the waveguide length. A sketch of the system is shown below:

We begin by defining the waveguide basis with two waveguides of different lengths. Here the first waveguide contains the full number of time bins, while the second waveguide contains only half as many:

using QuantumOptics
using WaveguideQED

N_ROUND_TRIPS = 5
FULL_WAVEGUIDE = 200

dt = 0.05
times = range(0, step=dt, length=FULL_WAVEGUIDE)
times_sim = range(0, step=dt, length=FULL_WAVEGUIDE*N_ROUND_TRIPS)

bw = WaveguideBasis(2,2,times,lengths=[FULL_WAVEGUIDE, FULL_WAVEGUIDE/2])

Here the argument lengths=[FULL_WAVEGUIDE, FULL_WAVEGUIDE/2] specifies that the first waveguide contains FULL_WAVEGUIDE time bins while the second contains only half as many. Physically this corresponds to photons in the second waveguide returning to the system after a shorter delay.

We can now construct the creation and annihilation operators for the two waveguides.

wd1 = create(bw,1)
w1 = destroy(bw,1)

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

Next we define two emitter subspaces which couple to the two waveguides. The emitter operators are defined using a tensor product basis.

bc = FockBasis(1)

c1  = destroy(bc) ⊗ identityoperator(bc)
c1d = create(bc)  ⊗ identityoperator(bc)

c2  = identityoperator(bc) ⊗ destroy(bc)
c2d = identityoperator(bc) ⊗ create(bc)

I_c = identityoperator(bc)

The interaction Hamiltonian describing the coupling between the cavities and their respective waveguides is

\[\begin{equation*} H = i \sqrt{\frac{\gamma}{\Delta t}} \left( \sigma _1 ^ \dagger w_1 - \sigma _1 w_1^\dagger\right) + i \sqrt{\frac{\gamma}{\Delta t}} \left(\sigma _2^\dagger w_2 - \sigma _2 w_2^\dagger\right) \end{equation*}\]

which can be implemented as:

γ = 5

H = im*sqrt(γ/dt)*((c1d ⊗ w1) - (c1 ⊗ wd1)) +
    im*sqrt(γ/dt)*((c2d ⊗ w2) - (c2 ⊗ wd2))

To monitor the emitter populations during the evolution we define the number operators.

n_c1 = (c1d*c1) ⊗ identityoperator(bw)
n_c2 = (c2d*c2) ⊗ identityoperator(bw)

function f_out(times, ψ)
    return [expect(n_c1, ψ), expect(n_c2, ψ)]
end

We initialize the system with one excitation in each emitter and no photons in the waveguide.

Ψ_in = fockstate(bc,1) ⊗ fockstate(bc,1) ⊗ zerophoton(bw)

Finally, we evolve the state using waveguide_evolution.

Ψ_out, exp_c1, exp_c2 = waveguide_evolution(times_sim, Ψ_in, H, fout=f_out)

We can visualize the emitter populations during the evolution.

fig,ax = subplots(1,1,figsize=(6,4))
ax.plot(times_sim, real.(exp_c1), label="emitter 1")
ax.plot(times_sim, real.(exp_c2), label="emitter 2")
ax.set_xlabel(L"time [$1/\gamma$]")
ax.set_ylabel("⟨n⟩")
legend()
plt.tight_layout()

Because the two waveguides have different lengths, photons emitted from the cavities return to the system after different delay times. This produces oscillatory dynamics as photons repeatedly circulate through the waveguide loop and re-enter the cavities.