diff --git a/docs/src/beamsplitter.md b/docs/src/beamsplitter.md index 6e018e1..3d337e4 100644 --- a/docs/src/beamsplitter.md +++ b/docs/src/beamsplitter.md @@ -2,7 +2,7 @@ Having introduced multiple waveguides in [`Multiple Waveguides`](@ref multiple), it is natural to implement the beamsplitter operation and consider some of the common measurements one would make on states having undergone a beamsplitter transformation. -We start by introducing how a beamsplitter can be implemented using two waveguides that interact. For simplicity, we start by considering only a single photon in one waveguide. First we create the basis, operators of the multiple waveguides and an initial single photon state with a gaussian wavefunction residing in waveguide 1: +We start by introducing how a beamsplitter can be implemented using two waveguides that interact. For simplicity, we consider only a single photon in one waveguide. First, we create the basis, operators of the multiple waveguides, and an initial single photon state with a gaussian wavefunction residing in waveguide 1: ```@example bs using WaveguideQED #hide @@ -23,7 +23,7 @@ psi = onephoton(bw,1, ξ) nothing #hide ``` -We now want to consider the following beamsplitter transformation: $w_{k,1} \rightarrow \cos(V) w_{k,1} - i \sin(V) w_{k,2}$ and equivalently for waveguide 2: $w_{k,2} \rightarrow - i \sin(V) w_{k,1} + \cos(V) w_{k,2}$. Having access to our initial Gaussian wavefunction we could just create the transformed state as: +We now want to consider the following beamsplitter transformation: $w_{k,1} \rightarrow \cos(V) w_{k,1} - i \sin(V) w_{k,2}$ and equivalently for waveguide 2: $w_{k,2} \rightarrow - i \sin(V) w_{k,1} + \cos(V) w_{k,2}$. Having access to our initial Gaussian wavefunction, we could just create the transformed state as: ```@example bs V = pi/4 @@ -46,12 +46,12 @@ n1 = wd1 *w1 n2 = wd2 *w2 psi = onephoton(bw,1, ξ) for (i,V) in enumerate(Vs) - H = V/dt*(wd1 * w2 + wd2*w1) - - psi_trans = waveguide_evolution(times,psi,H) + H  = V/dt*(wd1 * w2 + wd2*w1) +    + psi_trans = waveguide_evolution(times,psi,H) - transmission[i] = expect_waveguide(n1,psi_trans) - reflection[i] = expect_waveguide(n2,psi_trans) + transmission[i] = expect_waveguide(n1,psi_trans) + reflection[i] = expect_waveguide(n2,psi_trans) end fig,ax = subplots(1,1,figsize=(9,4.5)) ax.plot(Vs/pi,reflection,"b-",label="Waveguide a") @@ -65,16 +65,16 @@ nothing #hide ![beamsplitter](beamsplitter_trans.svg) -Here we see this population of waveguide 1 and 2 after the transformation vary as cosines and sines as we change the interaction parameter V. Thus, we confirm that we are applying the desired transformation. For an even beamsplitter, we thus choose $V=\pi/4$ +Here, we see this population of waveguides 1 and 2 after the transformation varies as cosines and sines as we change the interaction parameter V. Thus, we confirm that we are applying the desired transformation. For an even beamsplitter, we thus choose $V=\pi/4$ ## [Hong Ou Mandel with twophotons](@id hom) -As a more advanced example, we now consider a Hong Ou Mandel setup, where we have one photon in each waveguide impinging on a beamsplitter. If the two photons equivalent, we will see the Hong Ou Mandel effect and thus expect no photons in both waveguide simultanouesly after the transformation. As a measure of this, we calcuate the chance of having a coincidences count where one photon is in waveguide 1 while the other is in waveguide 2. This calculated using the two projection operators: +As a more advanced example, we now consider a Hong Ou Mandel setup, where we have one photon in each waveguide impinging on a beamsplitter. If the two photons equivalent, we will see the Hong Ou Mandel effect and thus expect no photons in both waveguide simultanouesly after the transformation. As a measure of this, we thus calculate the chance of having a coincidence count where one photon is in waveguide one while the other is in waveguide 2. This coincendence is calculated using the two projection operators: $$P_1 = \int_0^T dt w_1^\dagger(t) |0\rangle\langle0| w_1(t) \qquad P_2 = \int_0^T dt w_2^\dagger(t) |0\rangle\langle0| w_2(t)$$ where $w_1(t)$ and $w_2(t)$ are lowering operators for two waveguides. The chance of coincidence count is computed by $\langle\psi|P_1 P_2 |\psi\rangle$. -To compute the coincidence count expectation we create our own custom expectation value function: +To compute the coincidence count expectation, we create our own custom expectation value function: ```@example bs n1 = wd1*w1 @@ -82,35 +82,35 @@ 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 + 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 nothing #hide ``` -Here we evaluate $w_1^\dagger(t) w_1(t)$ at one timeidex `i`, while we evaluate $w_2^\dagger(t) w_2(t)$ at another timeidex `j`. Together this gives us the total coincedence count chance. In the following, we use the Hamiltonian from the previous section with $V=\pi/4$ and consider two Gaussian photons in each their waveguide with different centers of time $t_0$. By changing the difference in $t_0$, we can see the transition from a perfect overlap meaning no coincedence count, to no overlap meaning that the two photons never interact. In this case, the two photons will split up randomly and $1/4$ of the time they will end up in waveguide 1, similarly $1/4$ of the time they will end up in waveguide 2, and the remaining $1/2$ time they will end up in each of their waveguides. Thus, we expect a coincedence count of $1/2$ when the two pulses are fully seperated. Note that in the above function, we can just use the waveguide operators as projectors as we never have twophotons in both waveguides. +Here we evaluate $w_1^\dagger(t) w_1(t)$ at one timeidex `i`, while we evaluate $w_2^\dagger(t) w_2(t)$ at another timeidex `j`. Together, this gives us the total coincidence count chance. In the following, we use the Hamiltonian from the previous section with $V=\pi/4$ and consider two Gaussian photons in each their waveguide with different centers of time $t_0$. By changing the difference in $t_0$, we can see the transition from a perfect overlap, meaning no coincidence count, to no overlap, meaning that the two photons never interact. In this case, the two photons will split up randomly, and $1/4$ of the time, they will end up in waveguide 1. Similarly, $1/4$ of the time they will end up in waveguide 2, and the remaining $1/2$ time they will end up in each of their waveguides. Thus, we expect a coincidence count of $1/2$ when the two pulses are fully separated. Note that in the above function, we can just use the waveguide operators as projectors as we never have twophotons in both waveguides. ```@example bs taus = -3:0.1:3 ξ_twophoton(t1, t2, t01, t02) = ξ(t1, t01) * ξ(t2, t02) V = pi/4 -H = V/dt*(wd1*w2 + wd2*w1) +H  = V/dt*(wd1*w2 + wd2*w1) coincidences = zeros(length(taus)) t01 = 5 for (i, τ) in enumerate(taus) - t02 = 5 + τ - psi_pre = twophoton(bw, [1,2], ξ_twophoton, t01,t02) - ψ = waveguide_evolution(times,psi_pre,H) - coincidences[i] = expect_waveguide2(expval_op, ψ,times) + t02 = 5 + τ   + psi_pre = twophoton(bw, [1,2], ξ_twophoton, t01,t02) + ψ = waveguide_evolution(times,psi_pre,H) + coincidences[i] = expect_waveguide2(expval_op, ψ,times) end fig, ax = subplots(1,1, figsize=(9,4.5)) @@ -123,7 +123,7 @@ nothing #hide ``` ![hom_plot](hom.svg) -We could also have created this plot by performing the the beamsplitter operation by hand and instead initializing the state directly in this state. The initial state we consider is a single photon in each waveguide: $$|\psi \rangle = \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_1^\dagger(t_1) w_2^\dagger(t_2) \ket{\emptyset}$$, where $$\xi_1(t_1)$$ and $$\xi_2(t_2)$$ denote the wavefunction of the photon in waveguide 1 and 2, respectively. Notice that there is not factor of $$1/\sqrt(2)$$ in front of the initial state as the two photons occupy each their waveguide. If they initially occupied the same waveguide, we would need a factor of $$1/\sqrt(2)$$ for the state to be normalized. Performing the beamsplitter operation $w_1(t) \rightarrow 1/\sqrt(2) ( w_1(t) - i w_2(t))$ and $w_2(t) \rightarrow 1/\sqrt(2) ( - i w_1(t) + w_2(t))$, we arrive at the transformed state: +We could also have created this plot by performing the beamsplitter operation by hand and instead initializing the state directly in this state. The initial state we consider is a single photon in each waveguide: $$|\psi \rangle = \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_1^\dagger(t_1) w_2^\dagger(t_2) \ket{\emptyset}$$, where $$\xi_1(t_1)$$ and $$\xi_2(t_2)$$ denote the wavefunction of the photon in waveguide 1 and 2, respectively. Notice that there is no factor of $$1/\sqrt(2)$$ in front of the initial state as the two photons occupy different waveguides. If they initially occupied the same waveguide, we would need a factor of $$1/\sqrt(2)$$ for the state to be normalized. Performing the beamsplitter operation $w_1(t) \rightarrow 1/\sqrt(2) ( w_1(t) - i w_2(t))$ and $w_2(t) \rightarrow 1/\sqrt(2) ( - i w_1(t) + w_2(t))$, we arrive at the transformed state: $$\begin{align*} @@ -136,10 +136,10 @@ taus = -3:0.1:3 coincidences_manual = zeros(length(taus)) t01 = 5 for (i, τ) in enumerate(taus) - t02 = 5 + τ - psi_trans = 1/2*( im*twophoton(bw, 1, ξ_twophoton, t01,t02) - im * twophoton(bw, 2, ξ_twophoton, t01,t02) - + twophoton(bw, [1,2], ξ_twophoton, t01,t02) - twophoton(bw, [2,1], ξ_twophoton, t01,t02)) - coincidences_manual[i] = expect_waveguide2(expval_op, psi_trans,times) + t02 = 5 + τ   + psi_trans = 1/2*( im*twophoton(bw, 1, ξ_twophoton, t01,t02) - im * twophoton(bw, 2, ξ_twophoton, t01,t02) + + twophoton(bw, [1,2], ξ_twophoton, t01,t02) - twophoton(bw, [2,1], ξ_twophoton, t01,t02)) + coincidences_manual[i] = expect_waveguide2(expval_op, psi_trans,times) end fig, ax = subplots(1,1, figsize=(9,4.5)) diff --git a/docs/src/illustrations/emitters_circle.svg b/docs/src/illustrations/emitters_circle.svg new file mode 100644 index 0000000..228df2f --- /dev/null +++ b/docs/src/illustrations/emitters_circle.svg @@ -0,0 +1,1119 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 2 3 4 5 6 7 8 9 10 + + + 1 2 3 4 5 6 7 8 9 10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + A + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + B + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/illustrations/emitters_delay.svg b/docs/src/illustrations/emitters_delay.svg new file mode 100644 index 0000000..cd608eb --- /dev/null +++ b/docs/src/illustrations/emitters_delay.svg @@ -0,0 +1,481 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 2 3 4 5 6 7 8 9 10 + A + B + + diff --git a/docs/src/illustrations/emitters_delay_box.svg b/docs/src/illustrations/emitters_delay_box.svg new file mode 100644 index 0000000..f2b7ace --- /dev/null +++ b/docs/src/illustrations/emitters_delay_box.svg @@ -0,0 +1,2192 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + + + + + + + + + + + + + + + 5 + + + + + + + + + + 1 2 3 4 5 6 7 8 9 10 + A + B + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/index.md b/docs/src/index.md index bc55ec2..75290b1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,13 +1,13 @@ # WaveguideQED.jl -WaveguideQED.jl is a package for simulating continuous fockstates in waveguides. It expands on [`QuantumOptics.jl`](https://qojulia.org/) by adding a custom basis, operators, and routines for doing detection. +WaveguideQED.jl is a package for simulating continuous fockstates in waveguides. It expands on [`QuantumOptics.jl`](https://qojulia.org/) by adding a custom basis and operators for efficiently representing time-binned photon states. ## Dev docs Added functionalities: * [`WaveguideBasis`](@ref) for representing the waveguide Hilbert space and the related functions for generating states in this Hilbert space: [`zerophoton`](@ref), [`onephoton`](@ref), and [`twophoton`](@ref). Also see [`OnePhotonView`](@ref), [`TwoPhotonView`](@ref), and [`plot_twophoton!`](@ref) for viewing the waveguide states and plotting them. Note that [`WaveguideBasis`](@ref) can contain multiple waveguides. * [`WaveguideOperator`](@ref) are specialized operators allowing efficient annihilation and creation operators at each time-bin in the waveguide. They are created by giving a basis to [`WaveguideQED.destroy`](@ref) and [`WaveguideQED.create`](@ref) * Since the interaction between the waveguide time-bin mode $k$ and cavity/emitter is given as: $a^\dagger w_k - a w_k^\dagger$ there are specially optimized functions for doing these operations called [`CavityWaveguideOperator`](@ref) which are created using a fockbasis and a waveguide basis and the functions [`emission`](@ref) and [`absorption`](@ref). -* [`Detector`](@ref), [`LazyTensorKet`](@ref), and [`LazySumKet`](@ref), together with [`detect_single_click`](@ref) and [`detect_double_click`](@ref) allow one to do a beamsplitter interference and subsequent detection of photons coming from two waveguides. +* (OBSOLETE. SEE [Beamsplitter](@ref) INSTEAD). [`Detector`](@ref), [`LazyTensorKet`](@ref), and [`LazySumKet`](@ref), together with [`detect_single_click`](@ref) and [`detect_double_click`](@ref) allow one to do a beamsplitter interference and subsequent detection of photons coming from two waveguides. ![alt text](./animations/firstgif.gif) diff --git a/docs/src/time_delay.md b/docs/src/time_delay.md index 20314a5..66922f9 100644 --- a/docs/src/time_delay.md +++ b/docs/src/time_delay.md @@ -117,5 +117,143 @@ nothing #hide ``` ![alt text](feedback_0.svg) -For comparison, we here also plotted an exponential decay with the rate $\gamma$ which corresponds to an infinite delay time. We see that for $\phi=0$ we have constructive interference for the excitation going to the right and the population of the emitter decays faster than had there been no feedback. The photon thus experiences stimulated emission with itself! +For comparison, we here also plotteds an exponential decay with the rate $\gamma$ which corresponds to an infinite delay time. We see that for $\phi=0$ we have constructive interference for the excitation going to the right and the population of the emitter decays faster than had there been no feedback. +## [Spatially separated emitters](@id emitters) +Another interesting configuration to investigate is two spatially seperated emitters. If the separation between the emitters $\tau$ is comparable to the emission time $1/\gamma$, then a full description of the waveguide state is necessary. A sketch of the considered system can be seen here: + +![alt text](./illustrations/emitters_delay.svg) + + +An emission from emitter A thus arrives at emitter B via the waveguide after $\tau$ time and vice versa. We can represent this in the `Waveguide.jl` formalism as a looped conveyor belt, where the two emitters interact with the waveguide or conveyor belt at two different times/places. This is illustrated below: + +![alt text](./illustrations/emitters_circle.svg) + +We can create such a loop by using delayed operators. This way, the two emitters interact with the same waveguide in different time bins. The waveguide operators, per default, loop around themselves, and thus, by placing the delayed emitter exactly in the middle of the waveguide state, we recreate the loop above. This is also illustrated here: + +![alt text](./illustrations/emitters_delay_box.svg) + +Here, the looping mechanism is illustrated as "portals" that move the waveguide state's last box to the waveguide's beginning whenever the end is reached due to the progression of time. We thus create a waveguide basis with a length of $2 \tau$ and a set of waveguide operators which are delayed by $\tau$: + +```@example twoemitters +using QuantumOptics #hide +using WaveguideQED #hide +using LinearAlgebra #hide +using PyPlot #hide + +tau = 10 +dt = 0.1 +times = 0:dt:2*tau+dt +bw = WaveguideBasis(1, times) #waveguide basis +be = FockBasis(1) #emitter basis + +#Defining operators +w = destroy(bw) +wd = create(bw) +w_tau = destroy(bw; delay=tau/dt+1)#Delay of tau +wd_tau = create(bw; delay=tau/dt+1)#Delay of tau +sd = create(be) #Emitter operator +s = destroy(be) #Emitter operator +Iw = identityoperator(bw) #Identityoperator +Ie = identityoperator(be) #Identityoperator +nothing #hide +``` + +We can then create the Hamiltonian, where emitter A interacts with the not-delayed operators and emitter B with the delayed operators. We also create the expectation value operators here: + +```@example twoemitters +γ = 1 +H_e1 = im*sqrt(γ/dt)*(sd⊗Ie⊗w - s⊗Ie⊗wd) +H_e2 = im*sqrt(γ/dt)*(Ie⊗sd⊗w_tau - Ie⊗s⊗wd_tau) +H = H_e1 + H_e2 + +ne_a = (sd*s)⊗Ie⊗Iw +ne_b = Ie⊗(sd*s)⊗Iw +nw = Ie⊗Ie⊗(wd*w) +function expval(time, psi) +  expect(ne_a, psi), expect(ne_b, psi), expect_waveguide(nw, psi) +end +nothing #hide +``` + +Finally, we can simulate the emission of a photon from emitter A, which, after a delay, impinges on emitter B: + +```@example twoemitters +times_sim = 0:dt:2*tau +ψ = fockstate(be, 1) ⊗ fockstate(be, 0) ⊗ zerophoton(bw) +ψ, e1, e2, ew = waveguide_evolution(times_sim, ψ, H, fout=expval) + +fig,ax = subplots(1,1,figsize=(9,4.5)) +ax.plot(times_sim, abs.(e1), label="Emitter A") +ax.plot(times_sim, abs.(e2), label="Emitter B") +ax.plot(times_sim, abs.(ew), label="Waveguide") +ax.set_xlabel(L"Time [$1/\gamma$]") +ax.set_ylabel("Population") +tight_layout() +ax.legend() +savefig("e1_excited.svg") #hide +``` +![alt text](e1_excited.svg) + +In the above, we see emitter A decay and emitter B absorps the emitted photon, albeit not entirely. In reality, part of the incoming pulse from emitter A will reflect, and part of it will be absorbed. This means that if we run the simulation for a longer period, we will see emitter A being excited before emitter B is finished emitting. We can show this by extending the simulation time to $6 \tau$: + +```@example twoemitters +times_sim = 0:dt:8*tau +ψ = fockstate(be, 1) ⊗ fockstate(be, 0) ⊗ zerophoton(bw) +ψ, e1, e2, ew = waveguide_evolution(times_sim, ψ, H, fout=expval) + +fig,ax = subplots(1,1,figsize=(9,4.5)) +ax.plot(times_sim, abs.(e1), label="Emitter A") +ax.plot(times_sim, abs.(e2), label="Emitter B") +ax.plot(times_sim, abs.(ew), label="Waveguide") +ax.set_xlabel(L"Time [$1/\gamma$]") +ax.set_ylabel("Expectation Value") +tight_layout() +ax.legend() +savefig("e1_excited_slushing.svg") #hide +``` +![alt text](e1_excited_slushing.svg) + +Here, we see the pulse doing multiple roundtrips, and the excitation of the emitters gets more and more complicated as reflections and reflections of reflections interfere. We see that as more round trips are completed, there is a simultaneous excitation of emitter A and B, indicating the non-trivial nature of the feedback loop. This gets even more complicated if one considers a lifetime of the emitters $1/\gamma$ that is longer than the delay $1/\tau$. + + +One could also consider both emitters being excited in the beginning. This is done by allowing the waveguide to have two photons in it and changing the initial state. However, if there is no difference between the emitters, we will just see a completely symmetric excitation of them both, as the photons bounce back and forth. In the following, we thus consider two different decay rates of the emitters: + +```@example twoemitters +bw = WaveguideBasis(2, times) #waveguide basis with twophotons +be = FockBasis(1) #hide +wd = create(bw) #hide +wd_tau = create(bw; delay=tau/dt+1) #hide +w = destroy(bw) #hide +w_tau = destroy(bw; delay=tau/dt+1) #hide +Iw = identityoperator(bw) #hide +ne_a = (sd*s)⊗Ie⊗Iw #hide +ne_b = Ie⊗(sd*s)⊗Iw #hide +nw = Ie⊗Ie⊗(wd*w) #hide +function expval(time, psi) #hide +  expect(ne_a, psi), expect(ne_b, psi), expect_waveguide(nw, psi) #hide +end #hide + +γ = 1 +H_e1 = im*sqrt(γ/dt)*(sd⊗Ie⊗w - s⊗Ie⊗wd) +H_e2 = im*sqrt(γ/dt/2)*(Ie⊗sd⊗w_tau - Ie⊗s⊗wd_tau) +H = H_e1 + H_e2 + + +times_sim = 0:dt:6*tau +ψ = fockstate(be, 1) ⊗ fockstate(be, 1) ⊗ zerophoton(bw) +ψ, e1, e2, ew = waveguide_evolution(times_sim, ψ, H, fout=expval) + +fig,ax = subplots(1,1,figsize=(9,4.5)) +ax.plot(times_sim, abs.(e1), label="Emitter A") +ax.plot(times_sim, abs.(e2), label="Emitter B") +ax.plot(times_sim, abs.(ew), label="Waveguide") +ax.set_xlabel(L"Time [$1/\gamma$]") +ax.set_ylabel("Expectation Value") +tight_layout() +ax.legend() +savefig("both_excited.svg") #hide +``` +![alt text](both_excited.svg) + +Here, we see that emitter A, which has a larger decay rate than emitter B, quickly decays so that emitter B is reexcited before it is fully de-excited. Afterward, the evolution quickly grows complicated due to the multiple scattering events. \ No newline at end of file diff --git a/src/WaveguideOperator.jl b/src/WaveguideOperator.jl index 9b48f01..b9656c9 100644 --- a/src/WaveguideOperator.jl +++ b/src/WaveguideOperator.jl @@ -38,7 +38,7 @@ function Base.:eltype(x::WaveguideOperator) typeof(x.factor) end #Base.:*(x::WaveguideOperator{B1,B2},y::WaveguideOperator{B1,B2}) where {B1,B2} = LazyProduct((x,y),x.factor*y.factor) @inline function set_time!(o::WaveguideOperator, t::Number) - o.timeindex = min(round(Int,t/o.basis_l.dt,RoundDown) + 1,o.basis_l.nsteps) + o.timeindex = round(Int,t/o.basis_l.dt,RoundDown) + 1 return o end @@ -340,10 +340,8 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,1,idx},b,alpha,beta) wher rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + @inbounds result[1] += alpha*a.factor*b[timeindex+(idx-1)*nsteps+1] return end @@ -355,10 +353,8 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,2,1},b,alpha,beta) where rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + @inbounds result[1] += alpha*a.factor*b[timeindex+1] #twophotonview = TwoPhotonTimestepView(b,timeindex,nsteps,nsteps+1) #axpy!(alpha*a.factor,twophotonview,view(result,2:1:nsteps+1)) @@ -373,10 +369,8 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,2,idx},b,alpha,beta) wher rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + Nw = get_number_of_waveguides(a.basis_l) @inbounds result[1] += alpha*a.factor*b[timeindex+(idx-1)*a.basis_l.nsteps+1] #twophotonview = TwoPhotonTimestepView(b,timeindex,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) @@ -401,10 +395,8 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,1,idx},b,alpha,beta) where rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + @inbounds result[timeindex+(idx-1)*nsteps+1] += alpha*a.factor*b[1] return end @@ -416,10 +408,8 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,2,1},b,alpha,beta) where { rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + @inbounds result[timeindex+1] += alpha*a.factor*b[1] #twophotonview = TwoPhotonTimestepView(result,timeindex,nsteps,nsteps+1) #axpy!(alpha*a.factor,view(b,2:1:nsteps+1),twophotonview) @@ -434,10 +424,8 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,2,idx},b,alpha,beta) where rmul!(result, beta) end nsteps = a.basis_l.nsteps - timeindex = a.timeindex + a.delay - if timeindex < 0 || timeindex > nsteps - return - end + timeindex = (a.timeindex + a.delay -1) % (nsteps) +1 + Nw = get_number_of_waveguides(a.basis_l) @inbounds result[timeindex+(idx-1)*a.basis_l.nsteps+1] += alpha*a.factor*b[1] #twophotonview = TwoPhotonTimestepView(result,timeindex,nsteps,Nw*nsteps+1+(idx-1)*(nsteps*(nsteps+1))÷2) diff --git a/src/solver.jl b/src/solver.jl index 82ebfeb..1a69360 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -26,13 +26,13 @@ function waveguide_evolution(times,psi,H;fout=nothing,kwargs...) dt = get_dt(H.basis_l) nsteps = get_nsteps(H.basis_l) #Note that dt is added to the last input time to ensure that - tend = min(times[end],(nsteps)*dt)+dt + tend = times[end]+dt times_sim = 0:dt:tend isapprox(norm(psi),1,rtol=10^(-6)) || @warn "Initial waveguidestate is not normalized. Consider passing norm=true to the state generation function." function get_hamiltonian(time,psi) - tidx = min(round(Int,time/dt,RoundDown) + 1,nsteps) + tidx = round(Int,time/dt,RoundDown) + 1 set_waveguidetimeindex!(ops,tidx) return H end