diff --git a/Project.toml b/Project.toml
index 22b3834..e0b5594 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,10 +29,10 @@ LinearAlgebra = "1.10"
 ModelingToolkitBase = "1.42.2"
 NumericalIntegration = "0.3"
 OrdinaryDiffEq = "6"
-QuantumCumulants = "0.5"
+QuantumCumulants = "0.7"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.5"
-SecondQuantizedAlgebra = "0.9"
+SecondQuantizedAlgebra = "0.10"
 SpecialFunctions = "2"
 StaticArrays = "1"
 SymbolicUtils = "4"
diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl
index 3e31a3b..8a6c597 100644
--- a/benchmarks/correlations.jl
+++ b/benchmarks/correlations.jl
@@ -42,24 +42,22 @@ function benchmark_correlations!(SUITE)
     H_QO = to_numeric(H_sym, b; parameter = dict_p, time_parameter = dict_p_t)
     L_QO = to_numeric(L_sym, b; parameter = dict_p, time_parameter = dict_p_t)
 
-    function input_output(t, ρ)
-        Ht = H_QO(t)
-        J = [L_QO(t)]
-        return Ht, J, QuantumOpticsBase.dagger.(J)
-    end
-
     ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0)
-    _, ρt = timeevolution.master_dynamic(T, ψ0, input_output)
+    _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [L_QO])
 
     au_qo = to_numeric(au, b)
     c_qo = to_numeric(c, b)
     Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo
 
     ## --- Two-time correlation ---
+    # Pass the time-dependent operators (`H_QO`, `[L_QO]`) directly to the solver rather
+    # than wrapping them in a `(t, ρ)` closure: since v0.10 `to_numeric` returns a lazy
+    # `TimeDependentSum`, the direct path lets the solver build its integrator once and is
+    # markedly faster than re-reading the operator from a function at every step.
 
     SUITE["Correlations"]["two-time"] = BenchmarkGroup()
 
     SUITE["Correlations"]["two-time"]["single photon cavity"] =
-        @benchmarkable correlation_matrix($T, $ρt, $input_output, $Ls)
+        @benchmarkable correlation_matrix($T, $ρt, $H_QO, [$L_QO], $Ls)
     return nothing
 end
diff --git a/docs/Project.toml b/docs/Project.toml
index 99b5057..b9f6089 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -24,8 +24,9 @@ MacroTools = "0.5"
 ModelingToolkitBase = "1"
 OrdinaryDiffEq = "6"
 Plots = "1"
-QuantumCumulants = "0.5"
+QuantumCumulants = "0.7"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.5"
+SecondQuantizedAlgebra = "0.10"
 SymbolicUtils = "4"
 Symbolics = "7"
diff --git a/docs/src/api.md b/docs/src/api.md
index 85650a1..d647501 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -50,6 +50,10 @@ to_numeric
 Gaussian
 ```
 
+```@docs
+PulseCoupling
+```
+
 ```@docs
 coupling_input
 ```
diff --git a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl
index 5edf9b4..982c46a 100644
--- a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl
+++ b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl
@@ -1,8 +1,8 @@
 # # Cavity Scattering of a Single Photon
 
-# In this example, we simulate the scattering of a resonant single photon on an empty one-sided cavity. The temporal mode of the light pulse is a Gaussian with width $\sigma$ and the cavity has a decay rate $\gamma$. This system has been studied in [A. Kiilerich, et al., Phys. Rev. Lett. 123, 123604 (2019)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.123604). 
+# In this example, we simulate the scattering of a resonant single photon on an empty one-sided cavity. The temporal mode of the light pulse is a Gaussian with width $\sigma$ and the cavity has a decay rate $\gamma$. This system has been studied in [A. Kiilerich, et al., Phys. Rev. Lett. 123, 123604 (2019)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.123604).
 
-# We start by loading the needed packages and specifying the model. 
+# We start by loading the needed packages and specifying the model.
 
 using QuantumInputOutput
 using SecondQuantizedAlgebra
@@ -29,9 +29,9 @@ av = Destroy(h, :a_v, 3)
 @variables g_u::Number Δ::Real γ::Real g_v::Number
 nothing # hide
 
-# We use the symbolic operators and parameters to define the SLH triples and cascade them to obtain the Hamiltonian and Lindblad for the system. 
+# We use the symbolic operators and parameters to define the SLH triples and cascade them to obtain the Hamiltonian and Lindblad for the system.
 
-G_u = SLH(1, g_u'*au, 0) # input cavity 
+G_u = SLH(1, g_u'*au, 0) # input cavity
 G_c = SLH(1, √(γ)*c, Δ*c'c) # system cavity
 G_v = SLH(1, g_v'*av, 0) # output cavity
 
@@ -46,9 +46,9 @@ H = hamiltonian(G_cas)
 
 L = lindblad(G_cas)[1] # only one Lindblad term in this example
 
-# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of QuantumOptics.jl. To do so, we define the numerical parameters and operator basis. 
+# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of QuantumOptics.jl. To do so, we define the numerical parameters and operator basis.
 
-## numerical parameters 
+## numerical parameters
 γ_ = 1.0
 Δ_ = 0.0
 
@@ -67,21 +67,21 @@ T = [0:0.002:1;]*T_end
 gu_t = coupling_input(u1, T)
 dict_p_t = Dict(g_u => gu_t)
 
-## numeric bases 
+## numeric bases
 bu1 = FockBasis(1)
 bc1 = FockBasis(1)
 bv1 = FockBasis(1)
 b = bu1 ⊗ bc1 ⊗ bv1
 nothing # hide
 
-# We now use the function [`to_numeric`](@ref) to create the numeric operators. If the kwarg `time_parameter` is provided the created operator is a time-dependent function. 
+# We now use the function [`to_numeric`](@ref) to create the numeric operators. If the kwarg `time_parameter` is provided the created operator is a time-dependent function.
 
 H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t)
 L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t)
 nothing # hide
 
 
-# To solve the dynamics we use the QuantumOptics.jl function `timeevolution.master_dynamic`. 
+# To solve the dynamics we use the QuantumOptics.jl function `timeevolution.master_dynamic`.
 
 ## time-dependent function for timeevolution.master_dynamic that returns H(t), J(t) and Jd(t)
 function input_output_1(t, ρ)
@@ -107,10 +107,10 @@ n_c_t = real.(expect(c_qo'*c_qo, ρt))
 n_u1_t = real.(expect(au_qo'*au_qo, ρt))
 nothing # hide
 
-# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$.  
+# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$.
 
 Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo
-g1_m = correlation_matrix(T, ρt, input_output_1, Ls)
+g1_m = correlation_matrix(T, ρt, H_QO, [L_QO], Ls)
 nothing # hide
 
 #
@@ -127,7 +127,7 @@ p = heatmap(
 )
 p
 
-# The eigenvalues and corresponding eigenvectors are sorted in ascending order, which means the last eigenvalue corresponds to the highest populated temporal mode. 
+# The eigenvalues and corresponding eigenvectors are sorted in ascending order, which means the last eigenvalue corresponds to the highest populated temporal mode.
 
 F = eigen(g1_m)
 n_avg = round.(real.(F.values)*ΔT; digits = 3)
@@ -137,7 +137,7 @@ v_mode = (modes[:, end]) / sqrt(ΔT)
 @show n_avg[(end-1):end]
 nothing # hide
 
-# We use now the mode with the highest mean photon number as our out-mode to determine its quantum state. 
+# We use now the mode with the highest mean photon number as our out-mode to determine its quantum state.
 p_sym_2 = [γ, Δ]
 p_num_2 = [γ_, Δ_]
 dict_p_2 = Dict(p_sym_2 .=> p_num_2)
@@ -198,7 +198,7 @@ plot(p1, p2; layout = (2, 1), size = (600, 550))
 
 # ## Cavity with phase noise
 
-# We slightly adapt the above example by assuming the initial pulse to be in a coherent state and adding phase noise to the cavity. This results in scattering into multiple modes. 
+# We slightly adapt the above example by assuming the initial pulse to be in a coherent state and adding phase noise to the cavity. This results in scattering into multiple modes.
 
 ## new basis of the system
 bu1_3 = FockBasis(12)
@@ -212,7 +212,7 @@ c_3 = embed(b_3, 2, destroy(bc1_3))
 av_3 = embed(b_3, 3, destroy(bv1_3))
 cdc_3 = c_3'c_3
 
-## we use the same Hamiltonian as before but add a depasing term to the dissipation 
+## we use the same Hamiltonian as before but add a depasing term to the dissipation
 H_QO_3 = to_numeric(H, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2)
 L_QO_3 = to_numeric(L, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2)
 function input_output_3(t, ρ)
@@ -223,7 +223,7 @@ end;
 
 #
 
-# Due to the larger Hilbert space the time evolution takes a few seconds. 
+# Due to the larger Hilbert space the time evolution takes a few seconds.
 
 ψ0_3 = coherentstate(bu1_3, 2) ⊗ fockstate(bc1_3, 0) ⊗ fockstate(bv1_3, 0)
 t_3, ρt_3 = timeevolution.master_dynamic(T, ψ0_3, input_output_3)
diff --git a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl
index ab91ea5..910fc10 100644
--- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl
+++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl
@@ -8,7 +8,7 @@ using QuantumInputOutput
 using SecondQuantizedAlgebra
 using Symbolics: Symbolics
 using QuantumOptics
-using QuantumOpticsBase: dagger
+using QuantumOpticsBase: dagger, static_operator
 using Plots
 using LaTeXStrings
 
@@ -159,9 +159,11 @@ lT = length(T)
 G2 = zeros(lT, lT) # transmission
 G2_ref = zeros(lT, lT) # reflection
 
-L0(t) = L_R_QO(t)
+# Materialize the lazy `TimeDependentSum` to a concrete operator at each time, so the
+# quantum-regression products below give a plain operator usable as the solver's initial state.
+L0(t) = dense(static_operator(L_R_QO(t)))
 L0_dag(t) = dagger(L0(t))
-L0_ref(t) = L_L_QO(t)
+L0_ref(t) = dense(static_operator(L_L_QO(t)))
 L0_ref_dag(t) = dagger(L0_ref(t))
 
 for it1 = 1:(lT-1)
diff --git a/examples/08-1_pulse-delay__simple.jl b/examples/08-1_pulse-delay__simple.jl
index 78a1a59..16715be 100644
--- a/examples/08-1_pulse-delay__simple.jl
+++ b/examples/08-1_pulse-delay__simple.jl
@@ -190,6 +190,9 @@ p_t_sym = [g_u, g_in, g_out, g_v, M_ls...]
 p_t_num = [gu_, gin_, gout_, gv_, M_t_ls...]
 dict_p_t_int = Dict(p_t_sym .=> p_t_num)
 
+# The interaction picture eliminates the delay cavity `d`, so the numeric operators live
+# on the two-mode basis `bu ⊗ bv` (with `ad` mapped to the identity). Build on that basis.
+b_int = bu ⊗ bv
 au_int = destroy(bu) ⊗ one(bv)
 ## This is amazing! # hide
 ad_int = one(bu ⊗ bv)
@@ -198,13 +201,13 @@ operators = Dict(
     [au, au', ad, ad', av, av'] .=> [au_int, au_int', ad_int, ad_int', av_int, av_int'],
 )
 
-H_int_QO = to_numeric(H_int_sym, b; time_parameter = dict_p_t_int, operators)
+H_int_QO = to_numeric(H_int_sym, b_int; time_parameter = dict_p_t_int, operators)
 L_int_QO = [
-    to_numeric(L_int_sym[i], b; time_parameter = dict_p_t_int, operators) for
+    to_numeric(L_int_sym[i], b_int; time_parameter = dict_p_t_int, operators) for
     i = 1:length(L_int_sym)
 ]
-## H_int_QO = to_numeric(H_int_sym, b; time_parameter=dict_p_t_int) # hide 
-## L_int_QO = [to_numeric(L_int_sym[i], b; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide
+## H_int_QO = to_numeric(H_int_sym, b_int; time_parameter=dict_p_t_int) # hide
+## L_int_QO = [to_numeric(L_int_sym[i], b_int; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide
 nothing # hide
 
 #
diff --git a/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl b/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl
index 6d43b19..f16f70b 100644
--- a/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl
+++ b/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl
@@ -13,6 +13,7 @@
 using QuantumInputOutput
 using SecondQuantizedAlgebra
 using QuantumOptics
+using QuantumOpticsBase: dagger
 using Plots
 using LinearAlgebra
 using DataInterpolations
@@ -23,7 +24,7 @@ using DataInterpolations
 hc = FockSpace(:c)
 
 ## symbolic operator
-a = Destroy(hc, :a, 1)
+a = Destroy(hc, :a)
 
 ## symbolic parameters
 @variables γ::Real K::Real Δ::Real p::Number
diff --git a/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl b/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl
index 38ddc01..40dcca0 100644
--- a/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl
+++ b/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl
@@ -11,6 +11,7 @@ using QuantumInputOutput
 using SecondQuantizedAlgebra
 using Symbolics: Symbolics
 using QuantumOptics
+using QuantumOpticsBase: dagger
 using SymbolicUtils
 using LinearAlgebra
 using Plots
@@ -181,9 +182,9 @@ t_, ρt = timeevolution.master_dynamic(T, ψ0_fock, input_output)
 
 expect(au_QO'au_QO, ρt)
 
-p = plot(T, expect(au_QO' * au_QO, ρt); label = "")
-plot!(p, T, expect(ad1_QO' * ad1_QO, ρt); label = "")
-plot!(p, T, expect(ad2_QO' * ad2_QO, ρt); label = "")
+p = plot(T, real.(expect(au_QO' * au_QO, ρt)); label = "")
+plot!(p, T, real.(expect(ad1_QO' * ad1_QO, ρt)); label = "")
+plot!(p, T, real.(expect(ad2_QO' * ad2_QO, ρt)); label = "")
 p
 
 
@@ -197,8 +198,6 @@ expect(σ22_QO, ρt[end])
 pop_fock = zeros(length(Δ_ls))
 pop_coh = zeros(length(Δ_ls))
 idx_t1 = findmin(abs.(T .- t1))[2]
-using ProgressMeter
-prog = Progress(length(Δ_ls))
 
 for (it, Δn) in enumerate(Δ_ls)
     Δn = Δ_ls[it]
@@ -219,7 +218,6 @@ for (it, Δn) in enumerate(Δ_ls)
 
     # t_, ρt = timeevolution.master_dynamic(T, ψ0_coh, input_output)
     # pop_coh[it] = real(expect(σ22_QO, ρt[idx_t1]))
-    next!(prog)
 end
 nothing # hide
 
diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl
index 600b1c0..3f6e7b0 100644
--- a/src/QuantumInputOutput.jl
+++ b/src/QuantumInputOutput.jl
@@ -30,6 +30,7 @@ export SLH,
     to_numeric,
     substitute,
     # Pulse coupling
+    PulseCoupling,
     coupling_input,
     coupling_output,
     coupling_delay_in,
diff --git a/src/correlations.jl b/src/correlations.jl
index 177dbc5..3ea313a 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -5,35 +5,61 @@
 Compute the two-time correlation matrix
 ``g^{(1)}(t_1, t_2) = \\langle L_s^\\dagger(t_1) L_s(t_2) \\rangle``
 on the time grid `T`. Writes directly into output matrix.
+
+Supply the dynamics either as a `master_dynamic`-style function `f(t, ρ)`, or as operators
+passed straight to the solver: a time-dependent `H` (e.g. the `TimeDependentSum` from
+[`to_numeric`](@ref)) with jump operators `J`, or a constant `H` with constant `J`. The
+operator form is much faster for time-dependent problems (the integrator is built once).
+`Ls` is either a constant operator or a function `Ls(t)` returning the operator at `t`.
 """
-function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls::Function; kwargs...)
-    Ls_vec = Ls.(T)
-    Ls_dag_vec = dagger.(Ls_vec)
+function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...)
+    Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls)
     _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0
         timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...)
     end
 end
 
-function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...)
-    Ls_dag = dagger(Ls)
-    l_T = length(T)
-    Ls_vec = fill(Ls, l_T)
-    Ls_dag_vec = fill(Ls_dag, l_T)
+function correlation_matrix(
+    T::Vector,
+    ρt::Vector,
+    H::QuantumOpticsBase.AbstractTimeDependentOperator,
+    J::AbstractVector,
+    Ls;
+    kwargs...,
+)
+    Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls)
     _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0
-        timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...)
+        timeevolution.master_dynamic(T_slice, ρ0, copy(H), [copy(j) for j in J]; kwargs...)
     end
 end
 
-function correlation_matrix(T::Vector, ρt::Vector, H, J::Vector, Ls; kwargs...)
-    Ls_dag = dagger(Ls)
-    l_T = length(T)
-    Ls_vec = fill(Ls, l_T)
-    Ls_dag_vec = fill(Ls_dag, l_T)
+function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Ls; kwargs...)
+    Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls)
     _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0
         timeevolution.master(T_slice, ρ0, H, J; kwargs...)
     end
 end
 
+function _sample_operator_and_adjoint(T::Vector, op::Function)
+    vals = op.(T)
+    return vals, dagger.(vals)
+end
+function _sample_operator_and_adjoint(
+    ::Vector,
+    op::QuantumOpticsBase.AbstractTimeDependentOperator,
+)
+    throw(
+        ArgumentError(
+            "`Ls` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " *
+            "function `Ls(t)` returning the concrete operator at time `t`.",
+        ),
+    )
+end
+function _sample_operator_and_adjoint(T::Vector, op)
+    nt = length(T)
+    return fill(op, nt), fill(dagger(op), nt)
+end
+
 function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_dag_vec)
     l_T = length(T)
     @assert l_T == length(ρt)
diff --git a/src/pulses.jl b/src/pulses.jl
index 0b9b616..c5efc87 100644
--- a/src/pulses.jl
+++ b/src/pulses.jl
@@ -6,9 +6,23 @@ const _tol_div = 1e-10
 const _extrapolate = ExtrapolationType.Extension
 const _ϵ = 1e-10
 
-_mode_interp(mode::AbstractVector, T::AbstractVector) =
+"""
+    PulseCoupling(f)
+
+A single-argument time function `g(t)` used as a `to_numeric` `time_parameter` value.
+Wraps any evaluator `f` (a sampled interpolation, an analytic closure, ...) behind exactly
+one call method, so it satisfies SecondQuantizedAlgebra's single-arity `time_parameter`
+contract by construction (a raw interpolation is rejected for having methods of conflicting
+arity). Callable as `g(t)` and broadcastable as `g.(T)`.
+"""
+struct PulseCoupling{F}
+    f::F
+end
+(g::PulseCoupling)(t::Real) = g.f(t)
+
+_interpolate_mode(mode::AbstractVector, T::AbstractVector) =
     LinearInterpolation(mode, T; extrapolation = _extrapolate)
-_mode_interp(mode, T::AbstractVector) = mode
+_interpolate_mode(mode, T::AbstractVector) = mode
 
 # ──────────────────────────────────────────────
 # Gaussian pulse type
@@ -42,19 +56,19 @@ end
 # Shared coupling core
 # ──────────────────────────────────────────────
 
-function _compute_coupling(mode::Vector, T::Vector, denom_fn)
-    l_T = length(T)
-    buf = Vector{Float64}(undef, l_T)
-    map!(abs2, buf, mode)
-    ∫m2 = cumul_integrate(T, buf)
-    g = zeros(ComplexF64, l_T)
-    @inbounds for i = 2:l_T
-        d = denom_fn(∫m2[i])
-        if sqrt(abs(d)) > _tol_div
-            g[i] = mode[i]' / sqrt(d)
+function _coupling_from_mode(mode::Vector, T::Vector, denominator)
+    nt = length(T)
+    mode_sq = Vector{Float64}(undef, nt)
+    map!(abs2, mode_sq, mode)
+    ∫mode2 = cumul_integrate(T, mode_sq)
+    coupling = zeros(ComplexF64, nt)
+    @inbounds for i = 2:nt
+        denom = denominator(∫mode2[i])
+        if sqrt(abs(denom)) > _tol_div
+            coupling[i] = mode[i]' / sqrt(denom)
         end
     end
-    return LinearInterpolation(g, T; extrapolation = _extrapolate)
+    return PulseCoupling(LinearInterpolation(coupling, T; extrapolation = _extrapolate))
 end
 
 # ──────────────────────────────────────────────
@@ -65,9 +79,10 @@ end
     coupling_input(u, T)
 
 Compute the virtual-cavity input coupling ``g_u(t)`` from an input mode `u(t)`
-sampled on time grid `T`. Returns the `LinearInterpolation` directly (callable as `g(t)`).
+sampled on time grid `T`. Returns a [`PulseCoupling`](@ref) (callable as `g(t)`,
+broadcastable as `g.(T)`) that plugs directly into a `to_numeric` `time_parameter`.
 """
-coupling_input(u::Vector, T::Vector) = _compute_coupling(u, T, x -> abs(1 - x) + _ϵ)
+coupling_input(u::Vector, T::Vector) = _coupling_from_mode(u, T, x -> abs(1 - x) + _ϵ)
 coupling_input(u::Function, T::Vector) = coupling_input(u.(T), T)
 coupling_input(u::LinearInterpolation, T::Vector) = coupling_input(u.(T), T)
 
@@ -78,16 +93,17 @@ Analytical input coupling ``g_u(t)`` for a Gaussian pulse.
 """
 function coupling_input(g::Gaussian)
     mode, ∫m2 = _gaussian_mode(g)
-    return t -> mode(t)' / √(abs(1 - ∫m2(t)) + _ϵ)
+    return PulseCoupling(t -> mode(t)' / √(abs(1 - ∫m2(t)) + _ϵ))
 end
 
 """
     coupling_output(v, T)
 
 Compute the virtual-cavity output coupling ``g_v(t)`` from an output mode `v(t)`
-sampled on time grid `T`. Returns the `LinearInterpolation` directly (callable as `g(t)`).
+sampled on time grid `T`. Returns a [`PulseCoupling`](@ref) (callable as `g(t)`,
+broadcastable as `g.(T)`) that plugs directly into a `to_numeric` `time_parameter`.
 """
-coupling_output(v::Vector, T::Vector) = _compute_coupling(-v, T, x -> x + _ϵ)
+coupling_output(v::Vector, T::Vector) = _coupling_from_mode(-v, T, x -> x + _ϵ)
 coupling_output(v::Function, T::Vector) = coupling_output(v.(T), T)
 coupling_output(v::LinearInterpolation, T::Vector) = coupling_output(v.(T), T)
 
@@ -98,7 +114,7 @@ Analytical output coupling ``g_v(t)`` for a Gaussian pulse.
 """
 function coupling_output(g::Gaussian)
     mode, ∫m2 = _gaussian_mode(g)
-    return t -> -mode(t)' / √(∫m2(t) + _ϵ)
+    return PulseCoupling(t -> -mode(t)' / √(∫m2(t) + _ϵ))
 end
 
 # ──────────────────────────────────────────────
@@ -118,14 +134,13 @@ All kwargs are passed on to the ODE solver.
 function effective_output_mode(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...)
     @assert i > 1
     n = i - 1
-    # Capture as tuples for concrete closure types
-    gv_t = ntuple(k -> gv_fcts[k], n)
+    # Capture as a tuple for concrete closure types
+    gv = ntuple(k -> gv_fcts[k], n)
     v_i = v_fcts[i]
-    gv_all = ntuple(k -> gv_fcts[k], n)
     function multiple_outputs_α!(dα, α, p, t)
         gv_buf = p
         @inbounds for k = 1:n
-            gv_buf[k] = gv_t[k](t)
+            gv_buf[k] = gv[k](t)
         end
         vi_t = v_i(t)
         @inbounds for j = 1:n
@@ -145,7 +160,7 @@ function effective_output_mode(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...)
         α_t = sol_α(t)
         result = v_i(t)
         @inbounds for k = 1:n
-            result += gv_all[k](t)' * α_t[k]
+            result += gv[k](t)' * α_t[k]
         end
         return result
     end
@@ -168,8 +183,8 @@ function effective_output_mode(
     alg = Tsit5(),
     kwargs...,
 )
-    v_fcts = [_mode_interp(v_, T) for v_ in v_data]
-    gv_fcts = [_mode_interp(gv_, T) for gv_ in gv_data]
+    v_fcts = [_interpolate_mode(v_, T) for v_ in v_data]
+    gv_fcts = [_interpolate_mode(gv_, T) for gv_ in gv_data]
     return effective_output_mode(v_fcts, gv_fcts, T, i; alg, kwargs...)
 end
 
@@ -205,15 +220,14 @@ All kwargs are passed on to the ODE solver.
 function effective_input_mode(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...)
     @assert i > 1
     n = i - 1
-    # Capture as tuples for concrete closure types
-    gu_t = ntuple(k -> gu_fcts[k], n)
+    # Capture as a tuple for concrete closure types
+    gu = ntuple(k -> gu_fcts[k], n)
     u_i = u_fcts[i]
     gu_i = gu_fcts[i]
-    gu_all = ntuple(k -> gu_fcts[k], n)
     function multiple_inputs_α!(dα, α, p, t)
         gu_buf = p
         @inbounds for k = 1:n
-            gu_buf[k] = gu_t[k](t)
+            gu_buf[k] = gu[k](t)
         end
         ui_t = u_i(t)
         gui_t = gu_i(t)
@@ -234,7 +248,7 @@ function effective_input_mode(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...)
         α_t = sol_α(t)
         result = u_i(t)
         @inbounds for k = 1:n
-            result -= gu_all[k](t)' * α_t[k]
+            result -= gu[k](t)' * α_t[k]
         end
         return result
     end
@@ -252,8 +266,8 @@ function effective_input_mode(
     alg = Tsit5(),
     kwargs...,
 )
-    u_fcts = [_mode_interp(u_, T) for u_ in u_data]
-    gu_fcts = [_mode_interp(gu_, T) for gu_ in gu_data]
+    u_fcts = [_interpolate_mode(u_, T) for u_ in u_data]
+    gu_fcts = [_interpolate_mode(gu_, T) for gu_ in gu_data]
     return effective_input_mode(u_fcts, gu_fcts, T, i; alg, kwargs...)
 end
 
@@ -276,30 +290,30 @@ effective_input_mode(
 # coupling_delay_out / coupling_delay_in
 # ──────────────────────────────────────────────
 
-function _compute_coupling_delay(num_mode::Vector, u::Vector, v::Vector, T::Vector)
-    l_T = length(T)
-    buf = Vector{Float64}(undef, l_T)
-    map!(abs2, buf, u)
-    ∫u2 = cumul_integrate(T, buf)
-    map!(abs2, buf, v)
-    ∫v2 = cumul_integrate(T, buf)
-    g = zeros(ComplexF64, l_T)
-    @inbounds for i = 2:l_T
-        d = abs(∫v2[i] - ∫u2[i])
-        if sqrt(abs(d)) > _tol_div
-            g[i] = num_mode[i]' / sqrt(d + _ϵ)
+function _delay_coupling_from_modes(num_mode::Vector, u::Vector, v::Vector, T::Vector)
+    nt = length(T)
+    mode_sq = Vector{Float64}(undef, nt)
+    map!(abs2, mode_sq, u)
+    ∫u2 = cumul_integrate(T, mode_sq)
+    map!(abs2, mode_sq, v)
+    ∫v2 = cumul_integrate(T, mode_sq)
+    coupling = zeros(ComplexF64, nt)
+    @inbounds for i = 2:nt
+        denom = abs(∫v2[i] - ∫u2[i])
+        if sqrt(abs(denom)) > _tol_div
+            coupling[i] = num_mode[i]' / sqrt(denom + _ϵ)
         end
     end
-    return LinearInterpolation(g, T; extrapolation = _extrapolate)
+    return PulseCoupling(LinearInterpolation(coupling, T; extrapolation = _extrapolate))
 end
 
 """
     coupling_delay_out(u, v, T)
 
 Compute the out-coupling strength for a delay cavity.
-Returns `LinearInterpolation` directly.
+Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, broadcastable as `g.(T)`).
 """
-coupling_delay_out(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(u, u, v, T)
+coupling_delay_out(u::Vector, v::Vector, T::Vector) = _delay_coupling_from_modes(u, u, v, T)
 coupling_delay_out(u::Function, v::Function, T::Vector) =
     coupling_delay_out(u.(T), v.(T), T)
 coupling_delay_out(u::LinearInterpolation, v::LinearInterpolation, T::Vector) =
@@ -309,9 +323,9 @@ coupling_delay_out(u::LinearInterpolation, v::LinearInterpolation, T::Vector) =
     coupling_delay_in(u, v, T)
 
 Compute the in-coupling strength for a delay cavity.
-Returns `LinearInterpolation` directly.
+Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, broadcastable as `g.(T)`).
 """
-coupling_delay_in(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(-v, u, v, T)
+coupling_delay_in(u::Vector, v::Vector, T::Vector) = _delay_coupling_from_modes(-v, u, v, T)
 coupling_delay_in(u::Function, v::Function, T::Vector) = coupling_delay_in(u.(T), v.(T), T)
 coupling_delay_in(u::LinearInterpolation, v::LinearInterpolation, T::Vector) =
     coupling_delay_in(u.(T), v.(T), T)
diff --git a/test/test_code_quality.jl b/test/test_code_quality.jl
index 7e8e89c..8c8b30a 100644
--- a/test/test_code_quality.jl
+++ b/test/test_code_quality.jl
@@ -38,4 +38,5 @@ end
 
     all_concrete(typeof(QIO.SLH(1, 0, 0))) # TODO
     all_concrete(QIO.Gaussian)
+    all_concrete(QIO.PulseCoupling)
 end
diff --git a/test/test_correlations.jl b/test/test_correlations.jl
index 718aa86..eb4ca20 100644
--- a/test/test_correlations.jl
+++ b/test/test_correlations.jl
@@ -1,6 +1,6 @@
 using QuantumInputOutput
 using QuantumOptics
-using QuantumOpticsBase: dagger
+using QuantumOpticsBase: dagger, TimeDependentSum
 using LinearAlgebra
 using Test
 
@@ -45,4 +45,17 @@ using Test
     mode_num_aligned = mode_num * exp(-1im * angle(overlap))
     overlap_aligned = abs(sum(conj.(mode_exp) .* mode_num_aligned) * ΔT)
     @test overlap_aligned > 0.99
+
+    @testset "time-dependent Ls rejected" begin
+        td_Ls = TimeDependentSum([t -> 1.0 + 0im], J)
+        @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Ls)
+    end
+
+    @testset "time-dependent operator path matches static" begin
+        H_td = TimeDependentSum([t -> 0.0 + 0im], [one(ba)])
+        J_td = [TimeDependentSum([t -> sqrt(γ) + 0im], [σm])]
+        _, ρt_td = timeevolution.master_dynamic(T, ψ0, H_td, J_td)
+        g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Ls)
+        @test maximum(abs.(g1_td .- g1_static)) < 1e-8
+    end
 end
diff --git a/test/test_translate.jl b/test/test_translate.jl
index 6b16ddc..5fe1168 100644
--- a/test/test_translate.jl
+++ b/test/test_translate.jl
@@ -2,9 +2,11 @@ using QuantumInputOutput
 using SecondQuantizedAlgebra
 using SymbolicUtils
 using QuantumOpticsBase
-using QuantumOpticsBase: dagger
+using QuantumOpticsBase: dagger, static_operator, TimeDependentSum
 using Test
 
+mat(F, t) = dense(static_operator(F(t)))
+
 @testset "translate" begin
     @variables κ_L::Real κ_R::Real Δ::Real g::Real γ::Real
     @variables E::Complex E1::Complex E2::Complex
@@ -56,8 +58,8 @@ using Test
             time_parameter = dict_p_t2,
             operators = ops_dict,
         )
-        @test isa(F1, Function)
-        @test isequal(F1(0.1), a_QO*3)
+        @test F1 isa TimeDependentSum
+        @test sum(abs.((mat(F1, 0.1) - dense(a_QO*3)).data)) < 1e-12
         F2 = to_numeric(
             a*E,
             bc1;
@@ -65,7 +67,7 @@ using Test
             time_parameter = dict_p_t2,
             operators = ops_dict,
         )
-        @test isequal(F2(0.1), a_QO*E_t(0.1))
+        @test sum(abs.((mat(F2, 0.1) - dense(a_QO*E_t(0.1))).data)) < 1e-12
         # association of the scalar prefactor differs from `Δn*A'*A`, compare numerically
         @test sum(
             abs.(
@@ -101,38 +103,37 @@ using Test
         abs.((dense(to_numeric(Δ, b; parameter = dict_p1)) - dense(one(b)*Δn)).data),
     ) < 1e-12
     F3 = to_numeric(Δ, b; parameter = dict_p1, time_parameter = dict_p_t2)
-    @test sum(abs.((dense(F3(4)) - dense(one(b)*Δn)).data)) < 1e-8
+    @test sum(abs.((mat(F3, 4) - dense(one(b)*Δn)).data)) < 1e-8
     F4 = to_numeric(
         a*3*conj(E) + Δ*σ(2, 2),
         b;
         parameter = dict_p1,
         time_parameter = dict_p_t2,
     )
-    @test sum(abs.((F4(0.2) - dense(a_QO2*3*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8
+    @test sum(abs.((mat(F4, 0.2) - dense(a_QO2*3*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8
     F5 = to_numeric(
         a*conj(E) + Δ*σ(2, 2),
         b;
         parameter = dict_p1,
         time_parameter = dict_p_t2,
     )
-    @test sum(abs.((F5(0.2) - dense(a_QO2*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8
-    @inferred F5(0.2)  # QAdd time-dependent path should return concrete type
+    @test sum(abs.((mat(F5, 0.2) - dense(a_QO2*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8
     F5_ = to_numeric(a*conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2)
-    @test sum(abs.((dense(F5_(0.2)) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8
+    @test sum(abs.((mat(F5_, 0.2) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8
     F6 = to_numeric(conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2)
-    @test sum(abs.((dense(F6(0.2)) - dense(E_t_c(0.2)*one(b))).data)) < 1e-8
+    @test sum(abs.((mat(F6, 0.2) - dense(E_t_c(0.2)*one(b))).data)) < 1e-8
     F7 = to_numeric(conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, time_parameter = dict_p_t2)
-    @test sum(abs.((F7(0.2) - dense(E_t_c(0.2)*one(b) + Δn*σ_QO(2, 2))).data)) < 1e-8
+    @test sum(abs.((mat(F7, 0.2) - dense(E_t_c(0.2)*one(b) + Δn*σ_QO(2, 2))).data)) < 1e-8
     # a bare symbolic scalar without a numeric/time value cannot be translated
     @test_throws ArgumentError to_numeric(conj(E), b; parameter = dict_p1)
     F8 = to_numeric(E^2, b; parameter = dict_p1, time_parameter = dict_p_t2)
-    @test sum(abs.((dense(F8(0.2)) - dense(E_t(0.2)^2*one(b))).data)) < 1e-8
+    @test sum(abs.((mat(F8, 0.2) - dense(E_t(0.2)^2*one(b))).data)) < 1e-8
 
 
     @testset "time_parameter_normalization" begin
         dict_p_t_num = Dict([E] .=> [2.5])
         F_num = to_numeric(a*E, b; parameter = dict_p1, time_parameter = dict_p_t_num)
-        @test sum(abs.((dense(F_num(0.2)) - dense(a_QO2 * 2.5)).data)) < 1e-8
+        @test sum(abs.((mat(F_num, 0.2) - dense(a_QO2 * 2.5)).data)) < 1e-8
     end
 
     @testset "multiple_time_prefactors" begin
@@ -147,7 +148,7 @@ using Test
             time_parameter = dict_p_t_multi,
         )
         expected = dense(a_QO2 * E1_t(0.4) * E2_t_c(0.4))
-        @test sum(abs.((dense(F_multi(0.4)) - expected).data)) < 1e-8
+        @test sum(abs.((mat(F_multi, 0.4) - expected).data)) < 1e-8
     end
 
     @testset "to_numeric vector overload" begin
@@ -160,10 +161,10 @@ using Test
             operators = ops_dict,
         )
         @test length(translated_ops) == length(ops)
-        @test all(op -> op isa Function, translated_ops)
-        @test isequal(translated_ops[1](0.3), a_QO)
-        @test isequal(translated_ops[2](0.3), one(bc1) * Δn)
-        @test isequal(translated_ops[3](0.3), a_QO * E_t(0.3))
+        @test all(op -> op isa TimeDependentSum, translated_ops)
+        @test sum(abs.((mat(translated_ops[1], 0.3) - dense(a_QO)).data)) < 1e-12
+        @test sum(abs.((mat(translated_ops[2], 0.3) - dense(one(bc1) * Δn)).data)) < 1e-12
+        @test sum(abs.((mat(translated_ops[3], 0.3) - dense(a_QO * E_t(0.3))).data)) < 1e-12
     end
 
     @testset "to_numeric SLH overload" begin
@@ -176,11 +177,12 @@ using Test
             operators = ops_dict,
         )
 
-        @test H_QO isa Function
+        @test H_QO isa TimeDependentSum
         @test length(L_QO) == 1
-        @test L_QO[1] isa Function
-        @test sum(abs.((H_QO(0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8
-        @test sum(abs.((L_QO[1](0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8
+        @test L_QO[1] isa TimeDependentSum
+        @test sum(abs.((mat(H_QO, 0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8
+        @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) <
+              1e-8
     end
 
     @testset "substitute operators qmul" begin
@@ -231,7 +233,29 @@ using Test
         @variables gR::Real
         gR_t(t) = 1.0 + 2.0im
         F = to_numeric(im * gR * a3, b3; time_parameter = Dict(gR => gR_t))
-        @test F isa Function
-        @test sum(abs.((F(0.0) - dense((im * (1.0 + 2.0im)) * a3_QO)).data)) < 1e-8
+        @test F isa TimeDependentSum
+        @test sum(abs.((mat(F, 0.0) - dense((im * (1.0 + 2.0im)) * a3_QO)).data)) < 1e-8
+    end
+
+    @testset "PulseCoupling as time_parameter" begin
+        # A `PulseCoupling` from a coupling constructor must plug into `to_numeric`'s
+        # `time_parameter` without the conflicting-arity error a raw interpolation triggers.
+        hcp = FockSpace(:cp)
+        acp = Destroy(hcp, :a)
+        bcp = FockBasis(3)
+        acp_QO = destroy(bcp)
+        @variables gu::Complex
+        Tg = collect(0.0:0.05:1.0)
+        umode(t) = exp(-(t - 0.5)^2)
+        gu_t = coupling_input(umode, Tg)
+        @test gu_t isa PulseCoupling
+        Fpc = to_numeric(
+            acp * gu,
+            bcp;
+            time_parameter = Dict(gu => gu_t),
+            operators = Dict([acp, acp'] .=> [acp_QO, dagger(acp_QO)]),
+        )
+        @test Fpc isa TimeDependentSum
+        @test sum(abs.((mat(Fpc, 0.3) - dense(gu_t(0.3) * acp_QO)).data)) < 1e-10
     end
 end
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 1ee3b63..e92fadf 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -1,4 +1,5 @@
 using QuantumInputOutput
+using DataInterpolations: LinearInterpolation
 using Test
 
 @testset "utils" begin
@@ -23,6 +24,14 @@ using Test
     gu_ana = coupling_input(Gaussian(τ, σ; δ = δ))
     gv_ana = coupling_output(Gaussian(τ, σ; δ = δ))
 
+    # Every coupling constructor returns a `PulseCoupling`: callable and broadcastable.
+    @test gu_num isa PulseCoupling
+    @test gv_num isa PulseCoupling
+    @test gu_ana isa PulseCoupling
+    @test gv_ana isa PulseCoupling
+    @test gu_num(T[5]) isa ComplexF64
+    @test gu_num.(T[1:3]) isa AbstractVector{ComplexF64}
+
     gv_err = maximum(abs.(gv_num.(T[2:end]) .- gv_ana.(T[2:end])))
     gu_err = maximum(abs.(gu_num.(T[2:end]) .- gu_ana.(T[2:end])))
 
@@ -53,4 +62,39 @@ using Test
     u_eff_data2 = effective_input_mode([u1.(T), u2.(T)], [gu1.(T), gu2.(T)], T, 2)
     @test maximum(abs.(u_eff_f.(T) .- u_eff_data2.(T))) /
           maximum(abs.(u_eff_f.(T) .+ u_eff_data2.(T))) < 1e-2
+
+    # zero coupling ⇒ no distortion ⇒ effective mode equals the bare mode (exercises i = 3)
+    zero_c(t) = 0.0 + 0.0im
+    v3(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2)
+    v_eff3 = effective_output_mode([v1, v2, v3], [zero_c, zero_c, zero_c], T, 3)
+    u_eff3 = effective_input_mode([u1, u2, v3], [zero_c, zero_c, zero_c], T, 3)
+    @test maximum(abs.(v_eff3.(T) .- v3.(T))) < 1e-8
+    @test maximum(abs.(u_eff3.(T) .- v3.(T))) < 1e-8
+
+    # constant modes give the closed form g_out(tᵢ) = u / sqrt((|v|² - |u|²) i Δt)
+    a_c = 0.5
+    b_c = 0.8
+    Td = collect(0.0:0.01:2.0)
+    dt_c = Td[2] - Td[1]
+    uc(t) = a_c
+    vc(t) = b_c
+    g_out = coupling_delay_out(uc, vc, Td)
+    g_in = coupling_delay_in(uc, vc, Td)
+    exp_out = [a_c / sqrt((b_c^2 - a_c^2) * (i * dt_c)) for i = 2:length(Td)]
+    exp_in = [-b_c / sqrt((b_c^2 - a_c^2) * (i * dt_c)) for i = 2:length(Td)]
+    @test g_out isa PulseCoupling
+    @test maximum(abs.(g_out.(Td[2:end]) .- exp_out)) < 1e-6
+    @test maximum(abs.(g_in.(Td[2:end]) .- exp_in)) < 1e-6
+    # the shared denominator cancels, so in/out = -v/u exactly
+    @test maximum(abs.(g_in.(Td[2:end]) ./ g_out.(Td[2:end]) .- (-b_c / a_c))) < 1e-12
+
+    # LinearInterpolation dispatch must reproduce the analytic references
+    u_itp = LinearInterpolation(u.(T), T)
+    @test maximum(abs.(coupling_input(u_itp, T).(T[2:end]) .- gu_ana.(T[2:end]))) < 5e-4
+    @test maximum(abs.(coupling_output(u_itp, T).(T[2:end]) .- gv_ana.(T[2:end]))) < 5e-4
+    uc_itp = LinearInterpolation(uc.(Td), Td)
+    vc_itp = LinearInterpolation(vc.(Td), Td)
+    @test maximum(abs.(coupling_delay_out(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_out)) <
+          1e-6
+    @test maximum(abs.(coupling_delay_in(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_in)) < 1e-6
 end