qojulia · ChristophHotter · Jun 1, 2026
diff --git a/Project.toml b/Project.toml
@@ -17,6 +17,9 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
+[sources]
+SecondQuantizedAlgebra = {path = "/home/christoph/git/SecondQuantizedAlgebra.jl"}
+
 [compat]
 Aqua = "0.8"
 CheckConcreteStructs = "0.1"
@@ -29,14 +32,14 @@ LinearAlgebra = "1.10"
 ModelingToolkit = "10"
 NumericalIntegration = "0.3"
 OrdinaryDiffEq = "6"
-QuantumCumulants = "0.4"
+QuantumCumulants = "0.5"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.4, 0.5"
-SecondQuantizedAlgebra = "0.4.5"
+SecondQuantizedAlgebra = "0.5"
 SpecialFunctions = "2"
 StaticArrays = "1"
-Symbolics = "6"
-SymbolicUtils = "3.6"
+Symbolics = "7"
+SymbolicUtils = "4.30"
 Test = "1.10"
 julia = "1.10"
 

diff --git a/docs/src/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.md b/docs/src/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.md
@@ -30,8 +30,8 @@ c = Destroy(h, :c, 2)
 av = Destroy(h, :a_v, 3)
 
 # symbolic parameters
-gu, Δ, γ = rnumbers("g_u Δ γ")
-gv = cnumber("g_v")
+gu, Δ, γ = real_vars("g_u Δ γ")
+gv = complex_var("g_v")
 nothing # hide
 ````
 

diff --git a/docs/src/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.md b/docs/src/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.md
@@ -29,8 +29,8 @@ au = Destroy(h, :a_u, 1)
 av = Destroy(h, :a_v, 3)
 
 # symbolic parameters
-@rnumbers γ Δ Γ
-gv = rnumber("g_v")
+@variables γ::Real Δ::Real Γ::Real
+gv = real_var("g_v")
 nothing # hide
 ````
 

diff --git a/docs/src/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.md b/docs/src/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.md
@@ -28,8 +28,8 @@ au = Destroy(h, :a_u, 1)
 c = Destroy(h, :c, 2)
 
 # symbolic parameters
-@rnumbers γ Δ γ_p
-gu, gv = rnumbers("g_u g_v")
+@variables γ::Real Δ::Real γ_p::Real
+gu, gv = real_vars("g_u g_v")
 nothing # hide
 ````
 

diff --git a/docs/src/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.md b/docs/src/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.md
@@ -36,8 +36,8 @@ av1 = Destroy(h, :a_v1, 3)
 av2 = Destroy(h, :a_v2, 4)
 
 # symbolic parameters
-@rnumbers γ g ω12
-gv1, gv2 = cnumbers("g_v1 g_v2")
+@variables γ::Real g::Real ω12::Real
+gv1, gv2 = complex_vars("g_v1 g_v2")
 nothing # hide
 ````
 

diff --git a/docs/src/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.md b/docs/src/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.md
@@ -35,8 +35,8 @@ au1 = Destroy(h, :au_1, 2)
 av1 = Destroy(h, :av_2, 4)
 
 # symbolic parameters
-@rnumbers γ Δ
-gu1, gu2, gv1 = cnumbers("gu_1 gu_2 gv_1");
+@variables γ::Real Δ::Real
+gu1, gu2, gv1 = complex_vars("gu_1 gu_2 gv_1");
 nothing #hide
 ````
 

diff --git a/docs/src/examples/04-1_two-sided-cavity_with-atom_coh-drive.md b/docs/src/examples/04-1_two-sided-cavity_with-atom_coh-drive.md
@@ -17,7 +17,7 @@ using Plots
 ````
 
 ````@example 04-1_two-sided-cavity_with-atom_coh-drive
-@rnumbers E κ_L κ_R Δ g γ
+@variables E::Real κ_L::Real κ_R::Real Δ::Real g::Real γ::Real
 Natoms = 2
 
 hc = FockSpace(:cavity)

diff --git a/docs/src/examples/04-2_two-sided-cavity_with-atom_coh-drive__cumulants.md b/docs/src/examples/04-2_two-sided-cavity_with-atom_coh-drive__cumulants.md
@@ -10,14 +10,14 @@ Here we show how to solve the dynamics of the example `Two-sided Cavity with Ato
 using QuantumInputOutput
 using SecondQuantizedAlgebra
 using QuantumCumulants
-using ModelingToolkit
+using ModelingToolkitBase: @named, unknowns
 using OrdinaryDiffEq
 using QuantumOpticsBase
 using Plots
 ````
 
 ````@example 04-2_two-sided-cavity_with-atom_coh-drive__cumulants
-@rnumbers E κ_L κ_R Δ g γ
+@variables E::Real κ_L::Real κ_R::Real Δ::Real g::Real γ::Real
 Natoms = 2
 
 hc = FockSpace(:cavity)

diff --git a/docs/src/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.md b/docs/src/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.md
@@ -27,11 +27,11 @@ h = tensor([ha(i) for i = 1:N]...)
 σ(α, i, j) = Transition(h, "σ_$(α)", i, j, α)
 
 # symbolic parameters
-γR(i) = rnumber("γ^{($(i))}_R") # right-moving decay rate
-γL(i) = rnumber("γ^{($(i))}_L") # left-moving decay rate
-Δ(i) = rnumber("Δ_{$(i)}") # detuning
-ϕ(i, j) = rnumber("ϕ_{$(i)$(j)}") # phase between QD-i and QD-j
-Ein = rnumber("E_{in}") # coherent drive in the right-moving input
+γR(i) = real_var("γ^{($(i))}_R") # right-moving decay rate
+γL(i) = real_var("γ^{($(i))}_L") # left-moving decay rate
+Δ(i) = real_var("Δ_{$(i)}") # detuning
+ϕ(i, j) = real_var("ϕ_{$(i)$(j)}") # phase between QD-i and QD-j
+Ein = real_var("E_{in}") # coherent drive in the right-moving input
 nothing # hide
 ````
 

diff --git a/docs/src/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.md b/docs/src/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.md
@@ -28,11 +28,11 @@ h = tensor([ha(i) for i = 1:N]...)
 σ(α, i, j) = Transition(h, "σ_$(α)", i, j, α)
 
 # symbolic parameters
-γR(i) = rnumber("γ^{($(i))}_R")
-γL(i) = rnumber("γ^{($(i))}_L")
-Δ(i) = rnumber("Δ_{$(i)}")
-ϕ(i, j) = rnumber("ϕ_{$(i)$(j)}")
-Ein = rnumber("E_{in}")
+γR(i) = real_var("γ^{($(i))}_R")
+γL(i) = real_var("γ^{($(i))}_L")
+Δ(i) = real_var("Δ_{$(i)}")
+ϕ(i, j) = real_var("ϕ_{$(i)$(j)}")
+Ein = real_var("E_{in}")
 nothing # hide
 ````
 

diff --git a/docs/src/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.md b/docs/src/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.md
@@ -30,7 +30,7 @@ av_sym = Destroy(h, :a_v, 3)
 σ12_sym = Transition(h, :σ, 1, 2, 2)
 
 # symbolic parameters
-gu, γ, gv = rnumbers("gu γ gv")
+gu, γ, gv = real_vars("gu γ gv")
 
 # cascade the SLH elements
 G_u = SLH(1, gu' * au_sym, 0)
@@ -60,7 +60,7 @@ To do so, we first subtract $H_{uv}$ from $H$ and then replace the virtual cavit
 H_int_sym_ = simplify(H - H_uv)
 
 # symbolic coefficient matrix $M(t)$
-M(i, j) = cnumber("M_{$(i)$(j)}")
+M(i, j) = complex_var("M_{$(i)$(j)}")
 a0_ls = [au_sym, av_sym]
 la = length(a0_ls)
 a_int_ls = [sum(M(i, j)*a0_ls[j] for j = 1:la) for i = 1:la]

diff --git a/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md b/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md
@@ -26,7 +26,7 @@ au = Destroy(h, :a_u, 1)
 av = Destroy(h, :a_v, 2)
 
 # symbolic parameters
-@rnumbers gu gv r t
+@variables gu::Real gv::Real r::Real t::Real
 nothing # hide
 ````
 

diff --git a/docs/src/examples/07-2_hong-ou-mandel__quantum-pulse.md b/docs/src/examples/07-2_hong-ou-mandel__quantum-pulse.md
@@ -31,7 +31,7 @@ av1 = Destroy(h, :a_v1, 3)
 av2 = Destroy(h, :a_v2, 4)
 
 # symbolic parameters
-@rnumbers gu1 gu2 gv1 gv2 t r
+@variables gu1::Real gu2::Real gv1::Real gv2::Real t::Real r::Real
 nothing # hide
 ````
 

diff --git a/docs/src/examples/08-1_pulse-delay__simple.md b/docs/src/examples/08-1_pulse-delay__simple.md
@@ -32,7 +32,7 @@ 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")
+gu, gin, gout, gv = complex_vars("g_u g_in g_out g_v")
 nothing # hide
 ````
 
@@ -141,7 +141,7 @@ G_d_in = SLH(S2, [gin*ad, 0], 0)
 H_ud = hamiltonian(cascade(G_u2, G_d_in))
 H_int_sym_ = simplify(H - H_ud)
 
-M(i, j) = cnumber("M_{$(i)$(j)}")
+M(i, j) = complex_var("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]

diff --git a/docs/src/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.md b/docs/src/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.md
@@ -22,12 +22,12 @@ using Plots
 hc = FockSpace(:c)
 
 # symbolic operator
-a = Destroy(hc, :a, 1)
+a = Destroy(hc, :a)
 
 # symbolic parameters
-κ = rnumber("κ")
-ϵ = rnumber("ϵ")
-η = rnumber("η")
+κ = real_var("κ")
+ϵ = real_var("ϵ")
+η = real_var("η")
 nothing # hide
 ````
 

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
@@ -5,6 +5,7 @@
 # We start by loading the needed packages and specifying the model. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -25,8 +26,8 @@ c = Destroy(h, :c, 2)
 av = Destroy(h, :a_v, 3)
 
 ## symbolic parameters
-gu, Δ, γ = rnumbers("g_u Δ γ")
-gv = cnumber("g_v")
+gu, Δ, γ = real_vars("g_u Δ γ")
+gv = complex_var("g_v")
 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. 

diff --git a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl
@@ -5,6 +5,7 @@
 # We start by loading the packages and defining the symbolic operators and parameters. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -24,8 +25,8 @@ au = Destroy(h, :a_u, 1)
 av = Destroy(h, :a_v, 3)
 
 ## symbolic parameters
-@rnumbers γ Δ Γ
-gv = rnumber("g_v")
+@variables γ::Real Δ::Real Γ::Real
+gv = real_var("g_v")
 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. 

diff --git a/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl b/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl
@@ -5,6 +5,7 @@
 # We start by loading the packages and defining the symbolic operators and parameters. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -23,8 +24,8 @@ au = Destroy(h, :a_u, 1)
 c = Destroy(h, :c, 2)
 
 ## symbolic parameters
-@rnumbers γ Δ γ_p
-gu, gv = rnumbers("g_u g_v")
+@variables γ::Real Δ::Real γ_p::Real
+gu, gv = real_vars("g_u g_v")
 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. 

diff --git a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl
@@ -9,6 +9,7 @@
 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -32,8 +33,8 @@ av1 = Destroy(h, :a_v1, 3)
 av2 = Destroy(h, :a_v2, 4)
 
 ## symbolic parameters
-@rnumbers γ g ω12
-gv1, gv2 = cnumbers("g_v1 g_v2")
+@variables γ::Real g::Real ω12::Real
+gv1, gv2 = complex_vars("g_v1 g_v2")
 nothing # hide
 
 # The localized system consists of a cavity mode coupled to a three-level

diff --git a/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl b/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl
@@ -7,6 +7,7 @@
 # As usual, we start by loading the packages and defining the symbolic operators and parameters. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -30,8 +31,8 @@ au1 = Destroy(h, :au_1, 2)
 av1 = Destroy(h, :av_2, 4)
 
 ## symbolic parameters
-@rnumbers γ Δ
-gu1, gu2, gv1 = cnumbers("gu_1 gu_2 gv_1");
+@variables γ::Real Δ::Real
+gu1, gu2, gv1 = complex_vars("gu_1 gu_2 gv_1");
 
 # We use the symbolic operators and parameters to define the SLH triples and cascade them to obtain the Hamiltonian and Lindblad for the system. 
 

diff --git a/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl b/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl
@@ -6,13 +6,14 @@
 # However, for the numerical simulation of the empty cavity we provide a dictionary of the actual QuantumOptics.jl operators we want to use. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
 
 #
 
-@rnumbers E κ_L κ_R Δ g γ
+@variables E::Real κ_L::Real κ_R::Real Δ::Real g::Real γ::Real
 Natoms = 2
 
 hc = FockSpace(:cavity)

diff --git a/examples/04-2_two-sided-cavity_with-atom_coh-drive__cumulants.jl b/examples/04-2_two-sided-cavity_with-atom_coh-drive__cumulants.jl
@@ -3,16 +3,17 @@
 # Here we show how to solve the dynamics of the example `Two-sided Cavity with Atoms` in the Heisenberg picture with a higher-order mean-field approach (cumulant expansion), which is done with the package [QuantumCumulants.jl](https://github.com/qojulia/QuantumCumulants.jl).    
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumCumulants
-using ModelingToolkit
+using ModelingToolkitBase: @named, unknowns
 using OrdinaryDiffEq
 using QuantumOpticsBase
 using Plots
 
 #
 
-@rnumbers E κ_L κ_R Δ g γ
+@variables E::Real κ_L::Real κ_R::Real Δ::Real g::Real γ::Real
 Natoms = 2
 
 hc = FockSpace(:cavity)

diff --git a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl
@@ -5,6 +5,7 @@
 # mode), and we compute the time evolution of the transmitted and reflected intensities. 
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using SecondQuantizedAlgebra
 using QuantumOptics
 using Plots
@@ -22,11 +23,11 @@ h = tensor([ha(i) for i = 1:N]...)
 σ(α, i, j) = Transition(h, "σ_$(α)", i, j, α)
 
 ## symbolic parameters
-γR(i) = rnumber("γ^{($(i))}_R") # right-moving decay rate
-γL(i) = rnumber("γ^{($(i))}_L") # left-moving decay rate
-Δ(i) = rnumber("Δ_{$(i)}") # detuning
-ϕ(i, j) = rnumber("ϕ_{$(i)$(j)}") # phase between QD-i and QD-j
-Ein = rnumber("E_{in}") # coherent drive in the right-moving input
+γR(i) = real_var("γ^{($(i))}_R") # right-moving decay rate
+γL(i) = real_var("γ^{($(i))}_L") # left-moving decay rate
+Δ(i) = real_var("Δ_{$(i)}") # detuning
+ϕ(i, j) = real_var("ϕ_{$(i)$(j)}") # phase between QD-i and QD-j
+Ein = real_var("E_{in}") # coherent drive in the right-moving input
 nothing # hide
 
 # We use the symbolic operators and parameters to define the SLH triples, cascade the left and right moving channels, and concatenate them to obtain the Hamiltonian and Lindblad for the system. 

diff --git a/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.jl b/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.jl
@@ -6,6 +6,7 @@
 # As usual, we start by loading the packages and defining the operators and parameters of the system.
 
 using QuantumInputOutput
+using QuantumInputOutput: dagger
 using QuantumOptics
 using Plots