diff --git a/docs/src/illustrations/emitters_delay_box.svg b/docs/src/illustrations/emitters_delay_box.svg
index 66a3a38..f2b7ace 100644
--- a/docs/src/illustrations/emitters_delay_box.svg
+++ b/docs/src/illustrations/emitters_delay_box.svg
@@ -408,7 +408,7 @@
     <g
        id="g52680"
        transform="translate(-93.313559,-99.066128)"
-       inkscape:export-filename="..\scripts\open_quantum_optics\WaveguideQED\docs\src\illustrations\horseshoe.svg"
+       inkscape:export-filename="..\scripts\open_quantum_optics\WaveguideQED\docs\src\illustrations\emitters_delay_box.svg"
        inkscape:export-xdpi="1000.59"
        inkscape:export-ydpi="1000.59">
       <rect
@@ -2124,7 +2124,7 @@
        inkscape:export-xdpi="300"
        inkscape:export-ydpi="300">
       <ellipse
-         style="mix-blend-mode:normal;fill:#000000;fill-opacity:1;stroke:#3da458;stroke-width:1.2577;stroke-opacity:1;filter:url(#filter121)"
+         style="mix-blend-mode:normal;fill:#000000;fill-opacity:1;stroke:#bd6021;stroke-width:1.2577;stroke-opacity:1;filter:url(#filter121)"
          id="ellipse44"
          cx="142.65512"
          cy="131.52335"
@@ -2162,7 +2162,7 @@
        id="g86"
        transform="matrix(-1,0,0,1,170.95569,-68.90602)">
       <ellipse
-         style="mix-blend-mode:normal;fill:#000000;fill-opacity:1;stroke:#3da458;stroke-width:1.2577;stroke-opacity:1;filter:url(#filter121-2)"
+         style="mix-blend-mode:normal;fill:#000000;fill-opacity:1;stroke:#c36322;stroke-width:1.2577;stroke-opacity:1;filter:url(#filter121-2)"
          id="ellipse85"
          cx="142.65512"
          cy="131.52335"
diff --git a/docs/src/time_delay.md b/docs/src/time_delay.md
index 19552ff..13d3763 100644
--- a/docs/src/time_delay.md
+++ b/docs/src/time_delay.md
@@ -141,8 +141,8 @@ using WaveguideQED #hide
 using LinearAlgebra #hide
 using PyPlot #hide
 
-tau = 5
-dt = 0.05
+tau = 10
+dt = 0.1
 times = 0:dt:2*tau+dt
 bw = WaveguideBasis(1, times) #waveguide basis
 be = FockBasis(1) #emitter basis
@@ -162,9 +162,9 @@ 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
-γ = 5
-H_e1 = im*sqrt(γ/dt/2)*(sd⊗Ie⊗w - s⊗Ie⊗wd) 
-H_e2 = im*sqrt(γ/dt/2)*(Ie⊗sd⊗w_tau - Ie⊗s⊗wd_tau) 
+γ = 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
@@ -187,7 +187,7 @@ 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("Time")
+ax.set_xlabel(L"Time [$1/\gamma$]")
 ax.set_ylabel("Population")
 tight_layout()
 ax.legend()
@@ -206,7 +206,7 @@ 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("Time")
+ax.set_xlabel(L"Time [$1/\gamma$]")
 ax.set_ylabel("Expectation Value")
 tight_layout()
 ax.legend()
@@ -234,9 +234,9 @@ function expval(time, psi) #hide
   expect(ne_a, psi), expect(ne_b, psi), expect_waveguide(nw, psi) #hide
 end #hide
 
-γ = 2
-H_e1 = im*sqrt(γ/dt/2)*(sd⊗Ie⊗w - s⊗Ie⊗wd)
-H_e2 = im*sqrt(γ/dt/4)*(Ie⊗sd⊗w_tau - Ie⊗s⊗wd_tau)
+γ = 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
 
 
@@ -248,7 +248,7 @@ 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("Time")
+ax.set_xlabel(L"Time [$1/\gamma$]")
 ax.set_ylabel("Expectation Value")
 tight_layout()
 ax.legend()
diff --git a/src/WaveguideOperator.jl b/src/WaveguideOperator.jl
index 6513e6a..b9656c9 100644
--- a/src/WaveguideOperator.jl
+++ b/src/WaveguideOperator.jl
@@ -341,9 +341,7 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,1,idx},b,alpha,beta) wher
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1
-    end
+    
     @inbounds result[1] += alpha*a.factor*b[timeindex+(idx-1)*nsteps+1]
     return
 end
@@ -356,9 +354,7 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,2,1},b,alpha,beta) where
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1
-    end
+    
     @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))
@@ -374,9 +370,7 @@ function waveguide_mul!(result,a::WaveguideDestroy{B,B,2,idx},b,alpha,beta) wher
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1 
-    end
+    
     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)
@@ -402,9 +396,7 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,1,idx},b,alpha,beta) where
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1 
-    end
+    
     @inbounds result[timeindex+(idx-1)*nsteps+1] += alpha*a.factor*b[1]
     return
 end
@@ -417,9 +409,7 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,2,1},b,alpha,beta) where {
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1 
-    end
+    
     @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)
@@ -435,9 +425,7 @@ function waveguide_mul!(result,a::WaveguideCreate{B,B,2,idx},b,alpha,beta) where
     end
     nsteps = a.basis_l.nsteps
     timeindex = (a.timeindex + a.delay -1) % (nsteps) +1
-    if timeindex == 0 || timeindex > nsteps
-        timeindex = 1 
-    end
+    
     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)