Simple Pulse Delay with a Virtual Cavity

In this example, a single-photon pulse is emitted from an input cavity, delayed by a virtual delay cavity, and finally captured by an output cavity. The delay cavity is driven by an incoming pulse u(t) and simultaneously emits a delayed pulse u(t-τ) using the pulse-shaping couplings introduced in V. R. Christiansen and K. Mølmer, Phys. Rev. A 113, 013730 (2026).

using QuantumInputOutput
using SecondQuantizedAlgebra
using QuantumOptics
using SymbolicUtils
using LinearAlgebra
using Plots
using LaTeXStrings

# symbolic Hilbert space
hu = FockSpace(:u)
hd = FockSpace(:d)
hv = FockSpace(:v)
h = hu ⊗ hd ⊗ hv

# symbolic operators
au = Destroy(h, :a_u, 1)
ad = Destroy(h, :a_d, 2)
av = Destroy(h, :a_v, 3)

# symbolic parameters
gu, gin, gout, gv = cnumbers("g_u g_in g_out g_v")

The input cavity couples ($g_u(t)'*a_u$) into the input port of the delay cavity ($g_{in}(t)*a_d$) and the delay cavity couples the photons via the output port ($g_{out}(t)*a_d$) into the the output cavity ($g_v(t)'*a_v$). This leads to the following cascade of SLH elements.

G_u = SLH(1, gu'*au, 0)
G_u2 = concatenate(G_u, SLH(1, 0, 0))

S2 = Matrix(I, 2, 2)
G_d = SLH(S2, [gin*ad, gout*ad], 0)

G_v = SLH(1, gv'*av, 0)
G_v2 = concatenate(SLH(1, 0, 0), G_v)

G_cas = cascade(G_u2, G_d, G_v2)
H = get_hamiltonian(G_cas)
L = get_lindblad(G_cas)

# short pulse delay
σ = 1.0
tp = 6*σ
τ = 0.5σ # pulse delay
u(t) = 1/(√(σ)*π^(1/4)) * exp(-(t - tp)^2 / (2*σ^2))
u_del(t_) = u(t_ - τ)

Tend = 2tp + τ
dt = Tend/5e2
T = [0:dt:Tend;]

gu_ = u_to_gu(u, T)
gout_ = uv_to_gout(u_del, u, T)
gin_ = uv_to_gin(u_del, u, T)
gv_ = v_to_gv(u_del, T)

dict_p_t = Dict([gu, gout, gin, gv] .=> [gu_, gout_, gin_, gv_])

# numeric bases
n = 3
bu = FockBasis(n)
bd = FockBasis(n)
bv = FockBasis(n)
b = bu ⊗ bd ⊗ bv

H_QO = translate_qo(H, b; time_parameter = dict_p_t)
L_QO = [translate_qo(L[i], b; time_parameter = dict_p_t) for i = 1:length(L)]

function input_output(t, ρ)
    Ht = H_QO(t)
    Jt = [L_QO[i](t) for i = 1:length(L_QO)]
    return Ht, Jt, dagger.(Jt)
end

# time evolution
ψ0 = fockstate(bu, n) ⊗ fockstate(bd, 0) ⊗ fockstate(bv, 0)
t_, ρt = timeevolution.master_dynamic(T, ψ0, input_output)

au_qo = translate_qo(au, b)
ad_qo = translate_qo(ad, b)
av_qo = translate_qo(av, b)

nu = real.(expect(dagger(au_qo)*au_qo, ρt))
nd = real.(expect(dagger(ad_qo)*ad_qo, ρt))
nv = real.(expect(dagger(av_qo)*av_qo, ρt))

p = plot(T, nu; label = L"\langle a_u^\dagger a_u \rangle")
plot!(p, T, nd; label = L"\langle a_d^\dagger a_d \rangle")
plot!(p, T, nv; label = L"\langle a_v^\dagger a_v \rangle")
plot!(
    p;
    xlabel = "time",
    ylabel = "mean photon number",
    grid = true,
    legend = :best,
    size = (500, 300),
)
p

We can see that the pulse is perfectly absorbed by the delayed output mode $v(t) = u(t-\tau)$

Interaction picture for the input and delay cavities

Introducing a separate cavity to delay the pulse can be a big disadvantage if one has multiple modes. To eliminate the delay cavity, we can transform into the interaction picture of the output and delay cavity coupling. This is, however, only possible if the delay is larger than the pulse, because only then we have $g_{in}(t) \approx g_{v=u}(t)$. Nevertheless, since in most cases only the relative delay between different modes is crucial, e.g. for two arms of an interferometer, we can simply add a constant delay $T_c \gg \sigma$ to all modes.

G_d_in = SLH(S2, [gin*ad, 0], 0)
H_ud = get_hamiltonian(cascade(G_u2, G_d_in))
H_int_sym_ = simplify(H - H_ud)

M(i, j) = cnumber("M_{$(i)$(j)}")
a0_ls = [au, ad]
la = length(a0_ls)
a_int_ls = [sum(M(i, j)*a0_ls[j] for j = 1:la) for i = 1:la]

int_dict = Dict([a0_ls; adjoint.(a0_ls)] .=> [a_int_ls; adjoint.(a_int_ls)])

H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict))

\[0.5 i {g_out{^{*}}} {g_v{^{*}}} {M_{21}{^{*}}} a_u^\dagger a_{v} + 0.5 i {g_out{^{*}}} {M_{22}{^{*}}} {g_v{^{*}}} a_d^\dagger a_{v} -0.5 i {{g_v{^{*}}}{^{*}}} g_{out} {{M_{21}{^{*}}}{^{*}}} a_{u} a_v^\dagger -0.5 i {{g_v{^{*}}}{^{*}}} g_{out} {{M_{22}{^{*}}}{^{*}}} a_{d} a_v^\dagger\]

L_int_sym = simplify.(substitute_operators.(L, Ref(int_dict)))
L_int_sym[1]

\[{g_u{^{*}}} M_{{11}} a_{u} + g_{in} {{M_{21}{^{*}}}{^{*}}} a_{u} + M_{{12}} {g_u{^{*}}} a_{d} + g_{in} {{M_{22}{^{*}}}{^{*}}} a_{d}\]

L_int_sym[2]

\[g_{out} {{M_{22}{^{*}}}{^{*}}} a_{d} + {g_v{^{*}}} a_{v} + g_{out} {{M_{21}{^{*}}}{^{*}}} a_{u}\]

# long pulse delay
σ = 1.0
tp = 6*σ
τ = 6σ # pulse delay
u(t) = 1/(√(σ)*π^(1/4)) * exp(-(t - tp)^2 / (2*σ^2))
u_del(t_) = u(t_ - τ)

Tend = 2tp + τ
dt = Tend/5e2
T = [0:dt:Tend;]

gu_ = u_to_gu(u, T)
gout_ = uv_to_gout(u_del, u, T)
gin_ = uv_to_gin(u_del, u, T)
gv_ = v_to_gv(u_del, T)

# interaction-picture coefficient matrix M(t) for u ↔ d
A_ud = interaction_picture_A_2modes(gu_, gin_)
M_t = interaction_picture_M(A_ud, T)

M_ls = [M(i, j) for i = 1:la for j = 1:la]
M_t_ls = [t -> M_t(t)[i, j] for i = 1:la for j = 1:la]

p_t_sym = [gu, gin, gout, gv, M_ls...]
p_t_num = [gu_, gin_, gout_, gv_, M_t_ls...]
dict_p_t_int = Dict(p_t_sym .=> p_t_num)

au_int = destroy(bu) ⊗ one(bv)
ad_int = one(bu ⊗ bv)
av_int = one(bu) ⊗ destroy(bv)
operators = Dict(
    [au, au', ad, ad', av, av'] .=> [au_int, au_int', ad_int, ad_int', av_int, av_int'],
)

H_int_QO = translate_qo(H_int_sym, b; time_parameter = dict_p_t_int, operators)
L_int_QO = [
    translate_qo(L_int_sym[i], b; time_parameter = dict_p_t_int, operators) for
    i = 1:length(L_int_sym)
]

function input_output_int(t, ρ)
    Ht = H_int_QO(t)
    Jt = [L_int_QO[i](t) for i = 1:length(L_int_QO)]
    return Ht, Jt, dagger.(Jt)
end

ψ0_int = fockstate(bu, 3) ⊗ fockstate(bv, 0)
t_int, ρt_int = timeevolution.master_dynamic(T, ψ0_int, input_output_int)

nu_int = real.(expect(dagger(au_int)*au_int, ρt_int))
nv_int = real.(expect(dagger(av_int)*av_int, ρt_int))

Above we introduced the unity matrix for the delay cavity operators which guarantees that it does not have an effect. We can see that the delayed pulse is perfectly absorbed.

p = plot(T, nu_int; label = L"\langle a_u^\dagger a_u \rangle_{IP}")
plot!(p, T, nv_int; label = L"\langle a_v^\dagger a_v \rangle_{IP}")
plot!(
    p;
    xlabel = "time",
    ylabel = "mean photon number",
    grid = true,
    legend = :best,
    size = (500, 300),
)
p

Package versions

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(
    [
        "QuantumInputOutput",
        "SecondQuantizedAlgebra",
        "QuantumOptics",
        "Plots",
        "LaTeXStrings",
    ],
    mode = PKGMODE_MANIFEST,
)

Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 9V74 80-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEBUG = Documenter
  JULIA_NUM_THREADS = auto
Status `~/work/QuantumInputOutput.jl/QuantumInputOutput.jl/docs/Manifest.toml`
  [b964fa9f] LaTeXStrings v1.4.0
  [91a5bcdd] Plots v1.41.6
  [18f9eda6] QuantumInputOutput v0.1.0 `~/work/QuantumInputOutput.jl/QuantumInputOutput.jl`
  [6e0679c1] QuantumOptics v1.2.4
  [f7aa4685] SecondQuantizedAlgebra v0.4.4

This page was generated using Literate.jl.