build: move to SQA v0.10

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"

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

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"

docs/src/api.md

@@ -50,6 +50,10 @@ to_numeric
 Gaussian
 ```
 
+```@docs
+PulseCoupling
+```
+
 ```@docs
 coupling_input
 ```

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)

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)

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
 
 #

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

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
 

src/QuantumInputOutput.jl

@@ -30,6 +30,7 @@ export SLH,
     to_numeric,
     substitute,
     # Pulse coupling
+    PulseCoupling,
     coupling_input,
     coupling_output,
     coupling_delay_in,

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)