diff --git a/.other/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl b/.other/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl index a9047c1..e6de942 100644 --- a/.other/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl +++ b/.other/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl @@ -59,8 +59,8 @@ L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] G_atom = SLH(I4, L_atom, Δ*σ(2, 2)) G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom -H = get_hamiltonian(G_u_bs_d1_d2_atom) -L = get_lindblad(G_u_bs_d1_d2_atom) +H = hamiltonian(G_u_bs_d1_d2_atom) +L = lindblad(G_u_bs_d1_d2_atom) # diff --git a/.other/backups/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4__bu.jl b/.other/backups/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4__bu.jl index 9c14d09..67e6fda 100644 --- a/.other/backups/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4__bu.jl +++ b/.other/backups/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4__bu.jl @@ -60,8 +60,8 @@ L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] G_atom = SLH(I4, L_atom, Δ*σ(2, 2)) G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom -H = get_hamiltonian(G_u_bs_d1_d2_atom) -L = get_lindblad(G_u_bs_d1_d2_atom) +H = hamiltonian(G_u_bs_d1_d2_atom) +L = lindblad(G_u_bs_d1_d2_atom) # diff --git a/Project.toml b/Project.toml index c00639d..8123ba3 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0" [deps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503" +FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" @@ -12,6 +13,7 @@ QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" @@ -19,6 +21,7 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" # Aqua = "0.8" # CheckConcreteStructs = "0.1" DataInterpolations = "8" +FunctionWrappers = "1" DiffEqNoiseProcess = "5.27.0" # ExplicitImports = "1.12" LinearAlgebra = "1.10" @@ -30,6 +33,7 @@ QuantumOptics = "1" QuantumOpticsBase = "0.4, 0.5" SecondQuantizedAlgebra = "0.4.4" SpecialFunctions = "2" +StaticArrays = "1" Symbolics = "6" SymbolicUtils = "3.6" Test = "1.10" diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml index 59d1325..0d95222 100644 --- a/benchmarks/Project.toml +++ b/benchmarks/Project.toml @@ -1,5 +1,6 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuantumInputOutput = "18f9eda6-924c-47c7-a881-996695cfd7c6" QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl index a6e93d6..e0141e1 100644 --- a/benchmarks/correlations.jl +++ b/benchmarks/correlations.jl @@ -23,14 +23,14 @@ function benchmark_correlations!(SUITE) G_v = SLH(1, gv * av, 0) G_cas = ▷(G_u, G_c, G_v) - H_sym = get_hamiltonian(G_cas) - L_sym = get_lindblad(G_cas)[1] + H_sym = hamiltonian(G_cas) + L_sym = lindblad(G_cas)[1] γ_ = 1.0 σ_pulse = 1 / γ_ T = [0:0.002:1;] * 12σ_pulse u1(t) = 1 / (sqrt(σ_pulse) * π^(1 / 4)) * exp(-(t - 4σ_pulse)^2 / (2 * σ_pulse^2)) - gu_t = u_to_gu(u1, T) + gu_t = coupling_input(u1, T) bu1 = FockBasis(1) bc1 = FockBasis(1) @@ -61,7 +61,7 @@ function benchmark_correlations!(SUITE) SUITE["Correlations"]["two-time"] = BenchmarkGroup() SUITE["Correlations"]["two-time"]["single photon cavity"] = - @benchmarkable two_time_corr_matrix($T, $ρt, $input_output, $Ls) + @benchmarkable correlation_matrix($T, $ρt, $input_output, $Ls) return nothing end diff --git a/benchmarks/interaction_picture.jl b/benchmarks/interaction_picture.jl index f449c37..1d5c12e 100644 --- a/benchmarks/interaction_picture.jl +++ b/benchmarks/interaction_picture.jl @@ -1,6 +1,5 @@ """ Interaction picture benchmarks — computing coefficient matrices for mode transformations. -Based on Christiansen et al., PRA 107, 013706 (2023). """ function benchmark_interaction_picture!(SUITE) SUITE["Interaction Picture"] = BenchmarkGroup() @@ -14,9 +13,9 @@ function benchmark_interaction_picture!(SUITE) T_end = 12.0 T = [0:0.005:1;] * T_end - gu_t = u_to_gu(u, T) - gv_t = v_to_gv(u, T) - A_uv = interaction_picture_A_2modes(gu_t, gv_t) + gu_t = coupling_input(u, T) + gv_t = coupling_output(u, T) + A_uv = coupling_matrix((gu_t, gv_t)) ## --- Coupling matrix A(t) evaluation (called every ODE step) --- @@ -28,11 +27,11 @@ function benchmark_interaction_picture!(SUITE) SUITE["Interaction Picture"]["coupling matrix evaluation"]["2 modes"] = @benchmarkable $A_uv($t_mid) - # 4-mode A(t) — more realistic for multi-port systems + # 4-mode A(t) u2(t) = 1 / (sqrt(τ) * π^(1 / 4)) * exp(-0.5 * ((t - t_p * 1.5) / τ)^2) - g3_t = u_to_gu(u2, T) - g4_t = v_to_gv(u2, T) - A_4m = interaction_picture_A_4modes(gu_t, gv_t, g3_t, g4_t) + g3_t = coupling_input(u2, T) + g4_t = coupling_output(u2, T) + A_4m = coupling_matrix((gu_t, gv_t, g3_t, g4_t)) SUITE["Interaction Picture"]["coupling matrix evaluation"]["4 modes"] = @benchmarkable $A_4m($t_mid) @@ -42,10 +41,10 @@ function benchmark_interaction_picture!(SUITE) SUITE["Interaction Picture"]["coefficient matrix M"] = BenchmarkGroup() SUITE["Interaction Picture"]["coefficient matrix M"]["numerical (ODE)"] = - @benchmarkable interaction_picture_M($A_uv, $T) + @benchmarkable solve_mode_evolution($A_uv, $T) SUITE["Interaction Picture"]["coefficient matrix M"]["analytical (2 equal modes)"] = - @benchmarkable interaction_picture_M_2modes_equal($u, $T) + @benchmarkable solve_mode_evolution_symmetric($u, $T) ## --- Symbolic operator substitution --- @@ -65,9 +64,9 @@ function benchmark_interaction_picture!(SUITE) G_v = SLH(1, gv_sym' * av_sym, 0) G_cas = ▷(G_u, G_s, G_v) - H = get_hamiltonian(G_cas) - L = get_lindblad(G_cas)[1] - H_uv = get_hamiltonian(▷(G_u, G_v)) + H = hamiltonian(G_cas) + L = lindblad(G_cas)[1] + H_uv = hamiltonian(▷(G_u, G_v)) H_int_ = simplify(H - H_uv) M_sym(i, j) = cnumber("M_{$(i)$(j)}") diff --git a/benchmarks/pulse_couplings.jl b/benchmarks/pulse_couplings.jl index a0c639d..4cb20bc 100644 --- a/benchmarks/pulse_couplings.jl +++ b/benchmarks/pulse_couplings.jl @@ -17,9 +17,18 @@ function benchmark_pulse_couplings!(SUITE) SUITE["Pulse Couplings"]["single-pulse"] = BenchmarkGroup() - SUITE["Pulse Couplings"]["single-pulse"]["input"] = @benchmarkable u_to_gu($u, $T) + SUITE["Pulse Couplings"]["single-pulse"]["input"] = + @benchmarkable coupling_input($u, $T) - SUITE["Pulse Couplings"]["single-pulse"]["output"] = @benchmarkable v_to_gv($v, $T) + SUITE["Pulse Couplings"]["single-pulse"]["output"] = + @benchmarkable coupling_output($v, $T) + + # Gaussian analytical + SUITE["Pulse Couplings"]["single-pulse"]["input Gaussian"] = + @benchmarkable coupling_input(Gaussian($τ, $σ; δ = $δ)) + + SUITE["Pulse Couplings"]["single-pulse"]["output Gaussian"] = + @benchmarkable coupling_output(Gaussian($τ, $σ; δ = $δ)) ## --- Effective multi-pulse couplings --- @@ -31,10 +40,10 @@ function benchmark_pulse_couplings!(SUITE) u2(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - (τ + 2σ)) / σ)^2) SUITE["Pulse Couplings"]["multi-pulse"]["output 2 modes"] = - @benchmarkable v_eff([$v1, $v2], $T, 2) + @benchmarkable effective_output_mode([$v1, $v2], $T, 2) SUITE["Pulse Couplings"]["multi-pulse"]["input 2 modes"] = - @benchmarkable u_eff([$u1, $u2], $T, 2) + @benchmarkable effective_input_mode([$u1, $u2], $T, 2) return nothing end diff --git a/benchmarks/runbenchmarks.jl b/benchmarks/runbenchmarks.jl index f8abb40..f6c2d37 100644 --- a/benchmarks/runbenchmarks.jl +++ b/benchmarks/runbenchmarks.jl @@ -5,6 +5,7 @@ using QuantumOptics using QuantumOpticsBase using SymbolicUtils using LinearAlgebra +using FunctionWrappers const SUITE = BenchmarkGroup() diff --git a/benchmarks/slh_algebra.jl b/benchmarks/slh_algebra.jl index 551b327..a8617ea 100644 --- a/benchmarks/slh_algebra.jl +++ b/benchmarks/slh_algebra.jl @@ -8,7 +8,6 @@ function benchmark_slh_algebra!(SUITE) SUITE["SLH Algebra"]["symbolic"] = BenchmarkGroup() - # Setup: 3-cavity cascade (example 01-1) hu1 = FockSpace(:u1) hc1 = FockSpace(:c1) hv1 = FockSpace(:v1) @@ -27,8 +26,8 @@ function benchmark_slh_algebra!(SUITE) SUITE["SLH Algebra"]["symbolic"]["3-cavity cascade"] = @benchmarkable begin G_cas = ▷($G_u, $G_c, $G_v) - get_hamiltonian(G_cas) - get_lindblad(G_cas) + hamiltonian(G_cas) + lindblad(G_cas) end SUITE["SLH Algebra"]["symbolic"]["concatenate + cascade"] = @benchmarkable begin @@ -37,7 +36,7 @@ function benchmark_slh_algebra!(SUITE) G_R ⊞ G_L end - # Feedback: coherent-feedback OPO loop (from test_feedback.jl) + # Feedback: coherent-feedback OPO loop hs = FockSpace(:s) a_fb = Destroy(hs, :a, 1) κ_fb = 0.7 @@ -52,11 +51,10 @@ function benchmark_slh_algebra!(SUITE) SUITE["SLH Algebra"]["symbolic"]["feedback OPO loop"] = @benchmarkable feedback($G_fb_unconnected, 1 => 2, 2 => 1) - ## --- Numeric SLH (SLHqo with QuantumOptics operators) --- + ## --- Numeric SLH (unified SLH with QuantumOptics operators) --- - SUITE["SLH Algebra"]["numeric (SLHqo)"] = BenchmarkGroup() + SUITE["SLH Algebra"]["numeric"] = BenchmarkGroup() - # Setup: 2-QD bidirectional waveguide (example 05-2, realistic basis) N = 2 bu = FockBasis(20) ba = NLevelBasis(2) @@ -76,43 +74,45 @@ function benchmark_slh_algebra!(SUITE) t0 = 4σt T = [0:0.005:1;] * 3t0 u1(t) = 1 / (sqrt(σt) * π^(1 / 4)) * exp(-(t - t0)^2 / (2 * σt^2)) - gu_t = u_to_gu(u1, T) - - G_u_qo = SLHqo(1, t -> gu_t(t) * a_u, 0 * one(b)) - G_ϕ_12 = SLHqo(exp(1im * ϕn[1]), 0 * one(b), 0 * one(b)) - G_R1 = SLHqo(1, √(γRn[1]) * σ(1, 1, 2), 0 * one(b)) - G_R2 = SLHqo(1, √(γRn[2]) * σ(2, 1, 2), 0 * one(b)) - G_L1 = SLHqo(1, √(γLn[1]) * σ(1, 1, 2), 0 * one(b)) - G_L2 = SLHqo(1, √(γLn[2]) * σ(2, 1, 2), 0 * one(b)) - - SUITE["SLH Algebra"]["numeric (SLHqo)"]["2-QD waveguide composition"] = - @benchmarkable begin - G_R_t = $G_u_qo ▷ $G_R1 ▷ $G_ϕ_12 ▷ $G_R2 - G_L_t = $G_L2 ▷ $G_ϕ_12 ▷ $G_L1 - G_t = G_R_t ⊞ G_L_t - get_hamiltonian(G_t) - get_lindblad(G_t) - end + gu_t = coupling_input(u1, T) + + G_u_qo = SLH(1, t -> gu_t(t) * a_u, 0 * one(b)) + G_ϕ_12 = SLH(exp(1im * ϕn[1]), 0 * one(b), 0 * one(b)) + G_R1 = SLH(1, √(γRn[1]) * σ(1, 1, 2), 0 * one(b)) + G_R2 = SLH(1, √(γRn[2]) * σ(2, 1, 2), 0 * one(b)) + G_L1 = SLH(1, √(γLn[1]) * σ(1, 1, 2), 0 * one(b)) + G_L2 = SLH(1, √(γLn[2]) * σ(2, 1, 2), 0 * one(b)) + + SUITE["SLH Algebra"]["numeric"]["2-QD waveguide composition"] = @benchmarkable begin + G_R_t = $G_u_qo ▷ $G_R1 ▷ $G_ϕ_12 ▷ $G_R2 + G_L_t = $G_L2 ▷ $G_ϕ_12 ▷ $G_L1 + G_t = G_R_t ⊞ G_L_t + hamiltonian(G_t) + lindblad(G_t) + end ## --- Time-dependent closure evaluation (the ODE hot loop) --- SUITE["SLH Algebra"]["closure evaluation"] = BenchmarkGroup() - # Build the composed network once, then benchmark calling H(t) and L(t) G_R_t = G_u_qo ▷ G_R1 ▷ G_ϕ_12 ▷ G_R2 G_L_t = G_L2 ▷ G_ϕ_12 ▷ G_L1 G_t = G_R_t ⊞ G_L_t - H_f = get_hamiltonian(G_t) - L_f = get_lindblad(G_t) + H_f = hamiltonian(G_t) + L_f = lindblad(G_t) t_mid = T[length(T)÷2] + _callable(x) = x isa Union{Function,FunctionWrappers.FunctionWrapper} + Hf = _callable(H_f) ? H_f : (t -> H_f) + L_callables = [_callable(Li) ? Li : (t -> Li) for Li in L_f] + SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide H(t)"] = - @benchmarkable $H_f($t_mid) + @benchmarkable $Hf($t_mid) SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide L(t)"] = @benchmarkable begin - for Li in $L_f + for Li in $L_callables Li($t_mid) end end diff --git a/benchmarks/translation.jl b/benchmarks/translation.jl index c48f3c5..4a4b5c1 100644 --- a/benchmarks/translation.jl +++ b/benchmarks/translation.jl @@ -29,8 +29,8 @@ function benchmark_translation!(SUITE) G_c = SLH(1, √(γ_sym) * c_, Δ_sym * c_'c_) G_v = SLH(1, gv_sym * av_, 0) G_cas = ▷(G_u, G_c, G_v) - H_sym = get_hamiltonian(G_cas) - L_sym = get_lindblad(G_cas)[1] + H_sym = hamiltonian(G_cas) + L_sym = lindblad(G_cas)[1] ## --- Static translation (no time dependence) --- @@ -42,7 +42,6 @@ function benchmark_translation!(SUITE) dict_p_static = Dict([κ_R, κ_L, Δ_sym] .=> [1.5, 1.0, 0.2]) - # Composite expression: operator + parameter + transition expr_composite = a * 3 + Δ_sym * σ(2, 2) SUITE["Translation"]["static"]["atom-cavity"] = @@ -55,7 +54,6 @@ function benchmark_translation!(SUITE) E_t(t) = 2 * t + 1im dict_p_t = Dict(E => E_t) - # Expression with time-dependent + static parameters and operators expr_td = a * 3 * conj(E) + Δ_sym * σ(2, 2) SUITE["Translation"]["time-dependent"]["atom-cavity"] = @benchmarkable translate_qo( @@ -65,15 +63,14 @@ function benchmark_translation!(SUITE) time_parameter = $dict_p_t, ) - # Full cascade H and L translation (realistic workflow from example 06-1) + # Full cascade H and L translation γ_ = 1.0 σ_pulse = 1 / γ_ T = [0:0.002:1;] * 12σ_pulse u_pulse(t) = 1 / (sqrt(σ_pulse) * π^(1 / 4)) * exp(-(t - 4σ_pulse)^2 / (2 * σ_pulse^2)) - gu_t = u_to_gu(u_pulse, T) - gv_t = v_to_gv(u_pulse, T) + gu_t = coupling_input(u_pulse, T) + gv_t = coupling_output(u_pulse, T) - # Realistic basis sizes from interaction picture example (n_ph=20) bu = FockBasis(20) bc3 = FockBasis(6) bv = FockBasis(6) diff --git a/docs/Project.toml b/docs/Project.toml index 57b6048..bf715ca 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,6 +9,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255" +QuantumInputOutput = "18f9eda6-924c-47c7-a881-996695cfd7c6" QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907" diff --git a/docs/src/api.md b/docs/src/api.md index ca17600..1b8eca1 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -7,19 +7,15 @@ SLH ``` ```@docs -SLHqo +scattering ``` ```@docs -get_scattering +lindblad ``` ```@docs -get_lindblad -``` - -```@docs -get_hamiltonian +hamiltonian ``` ```@docs @@ -51,63 +47,51 @@ translate_qo ## [Pulses](@id API: Pulses) ```@docs -u_to_gu -``` - -```@docs -v_to_gv +Gaussian ``` ```@docs -u_to_gu_Gauss +coupling_input ``` ```@docs -v_to_gv_Gauss +coupling_output ``` ```@docs -u_eff +effective_input_mode ``` ```@docs -v_eff +effective_output_mode ``` ```@docs -uv_to_gin +coupling_delay_in ``` ```@docs -uv_to_gout +coupling_delay_out ``` ## [Interaction Picture](@id API: Interaction Picture) ```@docs -interaction_picture_M -``` - -```@docs -interaction_picture_M_2modes_equal -``` - -```@docs -interaction_picture_A_2modes +coupling_matrix ``` ```@docs -interaction_picture_A_3modes +solve_mode_evolution ``` ```@docs -interaction_picture_A_4modes +solve_mode_evolution_symmetric ``` ## [Correlations](@id API: Correlations) ```@docs -two_time_corr_matrix +correlation_matrix ``` ## [Utilities](@id API: Utilities) 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 index 5e1190e..ffe0a0b 100644 --- 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 @@ -38,20 +38,20 @@ 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. ````@example 01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3 -G_u = SLH(1, gu*au, 0) # input cavity +G_u = SLH(1, gu'*au, 0) # input cavity G_c = SLH(1, √(γ)*c, Δ*c'c) # system cavity -G_v = SLH(1, gv*av, 0) # output cavity +G_v = SLH(1, gv'*av, 0) # output cavity G_cas = ▷(G_u, G_c, G_v) nothing # hide ```` ````@example 01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3 -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) ```` ````@example 01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3 -L = get_lindblad(G_cas)[1] # only one Lindblad term in this example +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. @@ -63,7 +63,7 @@ To solve the dynamics of the system we translate the symbolic expressions into n p_sym = [γ, Δ, gv] p_num = [γ_, Δ_, 0] # gv=0 -dict_p = Dict(p_sym .=> p_num); +dict_p = Dict(p_sym .=> p_num) # Gaussian input mode σ = 1/γ_ @@ -73,7 +73,7 @@ u1(t) = 1/(sqrt(σ)*π^(1/4)) * exp(-(t - 4σ)^2 / (2*σ^2)) T = [0:0.002:1;]*T_end ΔT = T[2] - T[1] -gu_t = u_to_gu(u1, T) +gu_t = coupling_input(u1, T) dict_p_t = Dict(gu => gu_t) # numeric bases @@ -126,7 +126,7 @@ In order to determine suitable temporal output modes we calculate the two-time a ````@example 01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3 Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo -g1_m = two_time_corr_matrix(T, ρt, input_output_1, Ls) +g1_m = correlation_matrix(T, ρt, input_output_1, Ls) nothing # hide ```` @@ -164,7 +164,7 @@ p_num_2 = [γ_, Δ_] dict_p_2 = Dict(p_sym_2 .=> p_num_2) # time-dependent coupling for the output mode $v(t)$ -gv_t = v_to_gv(v_mode, T) +gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]) 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 index f3db2d1..e8723fc 100644 --- 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 @@ -39,7 +39,7 @@ We use the symbolic operators and parameters to define the SLH triples and casca ````@example 01-2_stimulated-emission__PRL2019_123-123604_fig4 G_u = SLH(1, √(Γ)*au, 0) # input cavity G_c = SLH(1, √(γ)*σ(1, 2), Δ*σ(2, 2)) # two-level atom -G_v = SLH(1, gv*av, 0) # output cavity +G_v = SLH(1, gv'*av, 0) # output cavity G_cas = G_u ▷ G_c ▷ 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 index 51f7302..1e6a213 100644 --- 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 @@ -36,7 +36,7 @@ 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. ````@example 02-1_cavity-phase-noise__PRA2020_102-_023717_fig2 -G_u = SLH(1, gu*au, 0) # input cavity +G_u = SLH(1, gu'*au, 0) # input cavity G_c = SLH(1, √(γ)*c, Δ*c'c) # scattering cavity G_cas = ▷(G_u, G_c) @@ -44,11 +44,11 @@ nothing # hide ```` ````@example 02-1_cavity-phase-noise__PRA2020_102-_023717_fig2 -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) ```` ````@example 02-1_cavity-phase-noise__PRA2020_102-_023717_fig2 -L = get_lindblad(G_cas)[1] # only one Lindblad in this example +L = lindblad(G_cas)[1] # only one Lindblad in this example ```` ````@example 02-1_cavity-phase-noise__PRA2020_102-_023717_fig2 @@ -68,7 +68,7 @@ u(t) = 1/(sqrt(τ)*π^(1/4)) * exp(-(t - tp)^2 / (2*τ^2)) T = [0:0.002:1;]*14τ ΔT = T[2] - T[1] -gu_ = u_to_gu(u, T) +gu_ = coupling_input(u, T) dict_p_t = Dict(gu => gu_) # numeric bases @@ -105,7 +105,7 @@ We calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L ````@example 02-1_cavity-phase-noise__PRA2020_102-_023717_fig2 Ls(t) = (gu_(t))'*au_qo + √(γ_)*c_qo -g1_m = two_time_corr_matrix(T, ρt, input_output, Ls); +g1_m = correlation_matrix(T, ρt, input_output, Ls); 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 new file mode 100644 index 0000000..9824683 --- /dev/null +++ b/docs/src/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.md @@ -0,0 +1,221 @@ +```@meta +EditURL = "../../../examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl" +``` + +# Photon Number and Mode Entanglement with a Quantum Emitter + +We simulate the decay of a three-level Λ emitter through a cavity into two dominant temporal +output modes. We identify these modes from the field autocorrelation +function, capture them with two virtual output cavities, and visualize the +resulting mode populations and the final atom-mode density matrix. +This example, reproduces Fig. 4 of +[A. Kiilerich and K. Molmer, Phys. Rev. A 102, 023717 (2020)](https://doi.org/10.1103/PhysRevA.102.023717). + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumOptics +using Plots +using LaTeXStrings +using LinearAlgebra +using DataInterpolations +```` + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +# symbolic Hilbert space +hc = FockSpace(:c) +hs = NLevelSpace(:atom, 3) +hv1 = FockSpace(:v1) +hv2 = FockSpace(:v2) +h = hc ⊗ hs ⊗ hv1 ⊗ hv2 + +# symbolic operators +a = Destroy(h, :a, 1) +σ(i, j) = Transition(h, :σ, i, j, 2) +av1 = Destroy(h, :a_v1, 3) +av2 = Destroy(h, :a_v2, 4) + +# symbolic parameters +@rnumbers γ g ω12 +gv1, gv2 = cnumbers("g_v1 g_v2") +nothing # hide +```` + +The localized system consists of a cavity mode coupled to a three-level +Λ emitter. The transition |g₁⟩ ↔ |e⟩ is resonant with the cavity and +|g₂⟩ ↔ |e⟩ is detuned by ω₁₂. + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +H_s = g * (a' * σ(1, 3) + a * σ(3, 1) + a' * σ(2, 3) + a * σ(3, 2)) + ω12 * σ(2, 2) +G_s = SLH(1, √(γ) * a, H_s) + +G_v1 = SLH(1, gv1' * av1, 0) +G_v2 = SLH(1, gv2' * av2, 0) +G = cascade(G_s, G_v1, G_v2) + +H = hamiltonian(G) +L = lindblad(G)[1] +nothing # hide +```` + +We use the parameters quoted in the paper and initialize the emitter in the excited state. + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +γ_ = 1.0 +g_ = 0.1γ_ +ω12_ = 0.5γ_ + +T = [0:0.005:1;]*100/γ_ +ΔT = T[2] - T[1] + +dict_p_1 = Dict([γ, g, ω12, gv1, gv2] .=> [γ_, g_, ω12_, 0.0, 0.0]) +nothing # hide + +# numeric bases and operators +bc = FockBasis(1) +bs = NLevelBasis(3) +bv1 = FockBasis(1) +bv2 = FockBasis(1) +b = bc ⊗ bs ⊗ bv1 ⊗ bv2 + +a_qo = destroy(bc) ⊗ one(bs) ⊗ one(bv1) ⊗ one(bv2) +σ_qo(i, j) = one(bc) ⊗ transition(bs, i, j) ⊗ one(bv1) ⊗ one(bv2) +av1_qo = one(bc) ⊗ one(bs) ⊗ destroy(bv1) ⊗ one(bv2) +av2_qo = one(bc) ⊗ one(bs) ⊗ one(bv1) ⊗ destroy(bv2) +nothing # hide +```` + +We first solve the decay dynamics without explicit output cavities and +determine the two dominant temporal output modes from the first-order +correlation matrix g⁽¹⁾(t₁, t₂). + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +H_QO_1 = translate_qo(H, b; parameter = dict_p_1) +L_QO_1 = translate_qo(L, b; parameter = dict_p_1) + +function input_output_1(t, ρ) + J = [L_QO_1] + return H_QO_1, J, dagger.(J) +end + +ψ0 = fockstate(bc, 0) ⊗ nlevelstate(bs, 3) ⊗ fockstate(bv1, 0) ⊗ fockstate(bv2, 0) +t_1, ρt_1 = timeevolution.master_dynamic(T, ψ0, input_output_1) +nothing # hide + +Ls = √(γ_) * a_qo +g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) + +F = eigen(g1_m) +n_avg = real.(F.values) * ΔT + +v1_mode = F.vectors[:, end] / √(ΔT) +v2_mode = F.vectors[:, end-1] / √(ΔT) +n1 = n_avg[end] +n2 = n_avg[end-1] +nothing # hide +```` + +After identifying the two modes, we add two cascaded virtual output cavities. +The coupling of the second cavity must be corrected for the reshaping caused by +the first output cavity, which is done with [`effective_output_mode`](@ref). + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +gv1_t = coupling_output(v1_mode, T) +v2_eff = effective_output_mode([v1_mode, v2_mode], T, 2) +gv2_t = coupling_output(v2_eff, T) + +dict_p_2 = Dict([γ, g, ω12] .=> [γ_, g_, ω12_]) +dict_p_t_2 = Dict(gv1 => gv1_t, gv2 => gv2_t) + +H_QO_2 = translate_qo(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +L_QO_2 = translate_qo(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) + +function input_output_2(t, ρ) + J = [L_QO_2(t)] + return H_QO_2(t), J, dagger.(J) +end + +t_2, ρt_2 = timeevolution.master_dynamic(T, ψ0, input_output_2) +nothing # hide +```` + +We monitor the excited-state population, the cavity population, and the +populations transferred to the two output modes. + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +P_e_t = real.(expect(σ_qo(3, 3), ρt_2)) +n_cavity_t = real.(expect(a_qo' * a_qo, ρt_2)) +n1_t = real.(expect(av1_qo' * av1_qo, ρt_2)) +n2_t = real.(expect(av2_qo' * av2_qo, ρt_2)) +nothing # hide +```` + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +p_a = heatmap( + T, + T, + real.(g1_m); + c = :inferno, + xlabel = L"\gamma t_2", + ylabel = L"\gamma t_1", + colorbar_title = L"g^{(1)}(t_1,t_2)", + xlims = (0, 100), + ylims = (0, 100), + aspect_ratio = 1, +) + +p_b = plot( + T, + real.(v1_mode); + color = :blue, + ls = :dash, + lw = 2, + label = "n₁ = $(round(n1; digits = 2))", +) +plot!(p_b, T, real.(v2_mode); color = :red, lw = 2, label = "n₂ = $(round(n2; digits = 2))") +plot!(p_b; xlabel = L"\gamma t", ylabel = L"\Re[v_i(t)]", legend = :topright) + +p_c = plot(T, n_cavity_t; color = :green, ls = :dot, lw = 2, label = L"n_\mathrm{cavity}") +plot!(p_c, T, P_e_t; color = :black, ls = :dashdot, lw = 2, label = L"P(|e\rangle)") +plot!(p_c, T, n1_t; color = :blue, lw = 2, label = L"n_1") +plot!(p_c, T, n2_t; color = :red, ls = :dash, lw = 2, label = L"n_2") +plot!( + p_c; + xlabel = L"\gamma t", + ylabel = "mean excitation", + ylims = (0, 1), + legend = :topright, +) + +plot(p_a, p_b, p_c, layout = (1, 3), size = (1200, 360)) +```` + +Note that the plotted real part of the output-mode functions are not the same +as in the paper due the arbitrary global phase. + +## Package versions + +These results were obtained using the following versions: + +````@example 02-3_mode-entanglement__PRA2020_102-023717_fig4 +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status( + [ + "QuantumInputOutput", + "SecondQuantizedAlgebra", + "QuantumOptics", + "Plots", + "DataInterpolations", + "LaTeXStrings", + ], + mode = PKGMODE_MANIFEST, +) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + 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 index 948a45b..a7b9c00 100644 --- 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 @@ -43,21 +43,21 @@ 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. ````@example 03-1_beam-combiner__PRA2023_107-023715_fig2-fig3 -G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 -G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 +G_u2 = SLH(1, gu2'*au2, 0) # input cavity 2 +G_u1 = SLH(1, gu1'*au1, 0) # input cavity 1 G_a = SLH(1, √(γ)*σ(1, 2), Δ*σ(2, 2)) # scattering atom -G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 +G_v1 = SLH(1, gv1'*av1, 0) # output cavity 1 G_cas = cascade(G_u2, G_u1, G_a, G_v1) nothing # hide ```` ````@example 03-1_beam-combiner__PRA2023_107-023715_fig2-fig3 -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) ```` ````@example 03-1_beam-combiner__PRA2023_107-023715_fig2-fig3 -L = get_lindblad(G_cas)[1] # only one Lindblad in this example +L = lindblad(G_cas)[1] # only one Lindblad in this example ```` Next, the numerical parameters and functions of the system are defined. @@ -77,7 +77,7 @@ u(t) = 1/(sqrt(τ)*π^(1/4)) * exp(-(t - tp)^2 / (2*τ^2)) T = [0:0.002:1;]*20 ΔT = T[2] - T[1] -gu_ = u_to_gu(u, T) +gu_ = coupling_input(u, T) dict_p_t = Dict(gu1 => gu_) nothing # hide ```` @@ -121,7 +121,7 @@ au1_qo = translate_qo(au1, b) σ_qo(i, j) = translate_qo(σ(i, j), b) Ls(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) -g1_m = two_time_corr_matrix(T, ρt, input_output, Ls); +g1_m = correlation_matrix(T, ρt, input_output, Ls); p = heatmap( T, @@ -173,24 +173,21 @@ nothing # hide ```` The pulse from input cavity $u_2$ is scattered on the cavity $u_1$. This distortion needs to be taken into account for the coupling of $u_2$, -which is done with the function [`u_eff`](@ref). +which is done with the function [`effective_input_mode`](@ref). The coupling of the $u_1$ cavity needs no adaptation, since it directly couples to the two-level system. ````@example 03-1_beam-combiner__PRA2023_107-023715_fig2-fig3 -gu1_ = u_to_gu(u1_new, T) +gu1_ = coupling_input(u1_new, T) u_new_data = [u1_new, u2_new] u_new_fct = [LinearInterpolation(u, T) for u in u_new_data] -```` - -effective u2 mode and corresponding coupling -````@example 03-1_beam-combiner__PRA2023_107-023715_fig2-fig3 -u2_for_gu2 = u_eff(u_new_fct, T, 2) -gu2_ = u_to_gu(u2_for_gu2, T) +# effective u2 mode and corresponding coupling +u2_for_gu2 = effective_input_mode(u_new_fct, T, 2) +gu2_ = coupling_input(u2_for_gu2, T) # coupling of the output mode -gv1_ = v_to_gv(v1_new, T) +gv1_ = coupling_output(v1_new, T) # dictionary for the time-dependent functions g_sym = [gu1, gu2, gv1] 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 index 16a3a8b..10ebd51 100644 --- 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 @@ -13,7 +13,6 @@ However, for the numerical simulation of the empty cavity we provide a dictionar using QuantumInputOutput using SecondQuantizedAlgebra using QuantumOptics -using QuantumCumulants using Plots ```` @@ -51,28 +50,18 @@ nothing # hide Note that one needs to be careful to not double-count the Hamiltonian terms with the concatenation rule. ````@example 04-1_two-sided-cavity_with-atom_coh-drive -H1 = get_hamiltonian(G_cav_L_R_drive) +H1 = hamiltonian(G_cav_L_R_drive) ```` ````@example 04-1_two-sided-cavity_with-atom_coh-drive -L1_L = get_lindblad(G_cav_L_R_drive)[1] +L1_L = lindblad(G_cav_L_R_drive)[1] ```` ````@example 04-1_two-sided-cavity_with-atom_coh-drive -L1_R = get_lindblad(G_cav_L_R_drive)[2] +L1_R = lindblad(G_cav_L_R_drive)[2] ```` -Here, the usual classical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ appears as a combination of Hamiltonian and Lindblad term. To show the meanfield equation for the intra-cavity field we use the function `meanfield` of [QuantumCumulants.jl](https://github.com/qojulia/QuantumCumulants.jl). We could, in principle, also proceed by solving this equation, see e.g. the example `Mean-field Two-sided Cavity`. - -````@example 04-1_two-sided-cavity_with-atom_coh-drive -eqs_a = meanfield([a], H1, [L1_L, L1_R]) -nothing # hide -```` - -```math -\begin{align} \frac{d}{dt} \langle a\rangle &= 1 i \Delta \langle a\rangle -1.0 \sqrt{\kappa{L}} E -0.5 \left( \left( \sqrt{\kappa{L}} \right)^{2} + \left( \sqrt{\kappa_{R}} \right)^{2} \right) \langle a\rangle \end{align} -``` - +Here, the usual classical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ appears as a combination of Hamiltonian and Lindblad term. To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl). Since we do not want to include the basis of the atoms, we provide a dictionary of operators with the kwarg `operators` in the function [`translate_qo`](@ref). ````@example 04-1_two-sided-cavity_with-atom_coh-drive 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 index 5e94c59..bf4cff1 100644 --- 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 @@ -49,15 +49,15 @@ nothing # hide Note that one needs to be careful to not double-count the Hamiltonian terms with the concatenation rule. ````@example 04-2_two-sided-cavity_with-atom_coh-drive__cumulants -H1 = get_hamiltonian(G_cav_L_R_drive) +H1 = hamiltonian(G_cav_L_R_drive) ```` ````@example 04-2_two-sided-cavity_with-atom_coh-drive__cumulants -L1_L = get_lindblad(G_cav_L_R_drive)[1] +L1_L = lindblad(G_cav_L_R_drive)[1] ```` ````@example 04-2_two-sided-cavity_with-atom_coh-drive__cumulants -L1_R = get_lindblad(G_cav_L_R_drive)[2] +L1_R = lindblad(G_cav_L_R_drive)[2] ```` The typical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ is a combination of Hamiltonian term and Lindblad. @@ -178,6 +178,7 @@ We derive the equations of motion for system with a second-order mean-field appr ````@example 04-2_two-sided-cavity_with-atom_coh-drive__cumulants J_add = [√(γ)*σ(α, 1, 2) for α = 1:Natoms] eqs2 = meanfield([a'a, σ(1, 2, 2)], H2, [L2_L, L2_R, J_add...]; order = 2) +nothing # hide ```` ```math 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 index ecc2185..bb6e222 100644 --- 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 @@ -56,11 +56,11 @@ G_L_t = G_L(2) ▷ G_ϕ(1, 2) ▷ G_L(1) G_t = G_R_t ⊞ G_L_t nothing # hide -H = get_hamiltonian(G_t) +H = hamiltonian(G_t) ```` ````@example 05-1_N-QDs_bidirectional-waveguide_coherent-pulse -L = get_lindblad(G_t) +L = lindblad(G_t) L_R = L[1] ```` diff --git a/docs/src/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.md b/docs/src/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.md index ded95d7..8763f0b 100644 --- a/docs/src/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.md +++ b/docs/src/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.md @@ -4,7 +4,7 @@ EditURL = "../../../examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo # Quantum Pulse Bi-Directional Waveguide -This example mirrors the example `Bi-Directional Waveguide` but uses the numeric SLH struct [`SLHqo`](@ref), to circumvent the symbolic derivation part. +This example mirrors the example `Bi-Directional Waveguide` but passes numeric operators directly to [`SLH`](@ref), circumventing the symbolic derivation part. Furthermore it drives the system with a *quantum* single-photon pulse via a virtual cavity. As usual, we start by loading the packages and defining the operators and parameters of the system. @@ -48,17 +48,17 @@ T = [0:0.005:1;]*Tend # u1(t) = sqrt(1 / (σt * √(2π)) * exp(-0.5 * (t - t0)^2 / σt^2)) # hide u1(t) = 1/(sqrt(σt)*π^(1/4)) * exp(-(t - t0)^2 / (2*σt^2)) -gu_t = u_to_gu(u1, T) +gu_t = coupling_input(u1, T) nothing # hide ```` -We use the [`SLHqo`](@ref) to directly use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions to the model the system. +We use the [`SLH`](@ref) to directly use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions to the model the system. ````@example 05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo -G_u = SLHqo(1, t -> gu_t(t) * a_u, 0*one(b)) -G_ϕ(i, j) = SLHqo(exp(1im * ϕn[i]), 0*one(b), 0*one(b)) -G_R(i) = SLHqo(1, √(γRn[i]) * σ(i, 1, 2), -Δn[i] * σ(i, 2, 2)) -G_L(i) = SLHqo(1, √(γLn[i]) * σ(i, 1, 2), 0*one(b)) +G_u = SLH(1, t -> gu_t(t) * a_u, 0*one(b)) +G_ϕ(i, j) = SLH(exp(1im * ϕn[i]), 0*one(b), 0*one(b)) +G_R(i) = SLH(1, √(γRn[i]) * σ(i, 1, 2), -Δn[i] * σ(i, 2, 2)) +G_L(i) = SLH(1, √(γLn[i]) * σ(i, 1, 2), 0*one(b)) # Cascade right-moving channel # G_R_t = G_d ▷ cascade([G_R(i) ▷ G_ϕ(i, i + 1) for i=1:N-1]...) ▷ G_R(N) # for N > 1 # hide @@ -75,8 +75,8 @@ nothing # hide The full Hamiltonian and Lindblad terms are extracted from the final SLH element. Note that as soon as one time-dependent function is involved in a cascade or concatenate, the returned $H$ and $L$ will also be time-dependent. ````@example 05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo -H = get_hamiltonian(G_t) -L = get_lindblad(G_t) +H = hamiltonian(G_t) +L = lindblad(G_t) L_R = L[1] L_L = L[2] 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 new file mode 100644 index 0000000..b4d6f47 --- /dev/null +++ b/docs/src/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.md @@ -0,0 +1,167 @@ +```@meta +EditURL = "../../../examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl" +``` + +# Bi-Directional Waveguide via Feedback Reduction + +This example reconstructs the `N=2` bidirectional-waveguide model using the +SLH feedback reduction rule. The resulting Hamiltonian and Lindblad operators +are then used to simulate the transmitted and reflected intensities for a +coherent input pulse. + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumOptics +using Plots +using LinearAlgebra +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +N = 2 + +# symbolic Hilbert space +ha(i) = NLevelSpace("a$(i)", 2) +h = tensor([ha(i) for i = 1:N]...) + +# symbolic operators +σ(α, 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}") +nothing # hide +```` + +Each quantum dot is treated as a two-port component: port 1 couples to the +right-moving field and port 2 couples to the left-moving field. We concatenate +the source, the two dots, and the phase shifter, and then eliminate the +internal waveguide links with six feedback reductions. + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +I2 = Matrix{Int}(I, 2, 2) +G_in = SLH(I2, [Ein, 0], 0) +G_qd(i) = SLH(I2, [√(γR(i)) * σ(i, 1, 2), √(γL(i)) * σ(i, 1, 2)], -Δ(i) * σ(i, 2, 2)) +G_phase(i, j) = SLH([exp(1im * ϕ(i, j)) 0; 0 exp(1im * ϕ(i, j))], [0, 0], 0) + +G_unc = G_in ⊞ G_qd(1) ⊞ G_phase(1, 2) ⊞ G_qd(2) +G_t = feedback(G_unc, 1 => 3, 3 => 5, 5 => 7, 8 => 6, 6 => 4, 4 => 2) +nothing # hide +```` + +The reduced model has two external ports. In the remaining port order, the +first Lindblad operator corresponds to the reflected left-moving output and the +second one to the transmitted right-moving output. + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +H = hamiltonian(G_t) +L = lindblad(G_t) +L_L = L[1] +L_R = L[2] +nothing # hide +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +γ_ = 1.0 +β = 0.9 +γRn = fill(γ_ * β / 2, N) +γLn = fill(γ_ * β / 2, N) +γ_add = fill(γ_ * (1 - β), N) +Δn = fill(0.0, N) +ϕn = fill(π / 10, max(N - 1, 0)) + +σt = 0.8 +α0 = √(0.1) +t0 = 4σt +Tend = 3t0 +u1(t) = 1 / (sqrt(σt) * π^(1 / 4)) * exp(-(t - t0)^2 / (2 * σt^2)) +Ein_t(t) = α0 * u1(t) + +p_sym = [ + [γR(i) for i = 1:N]; + [γL(i) for i = 1:N]; + [Δ(i) for i = 1:N]; + [ϕ(i, i + 1) for i = 1:(N-1)] +] +p_num = [γRn; γLn; Δn; ϕn] +dict_p = Dict(p_sym .=> p_num) +dict_p_t = Dict(Ein => Ein_t) +nothing # hide +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +ba = NLevelBasis(2) +b = tensor([ba for _ = 1:N]...) + +H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) +L_R_QO = translate_qo(L_R, b; parameter = dict_p, time_parameter = dict_p_t) +L_L_QO = translate_qo(L_L, b; parameter = dict_p, time_parameter = dict_p_t) + +σ_qo(α, i, j) = translate_qo(σ(α, i, j), b) +J_add = [√(γ_add[i]) * σ_qo(i, 1, 2) for i = 1:N] + +function input_output(t, ρ) + Ht = H_QO(t) + J = [L_R_QO(t), L_L_QO(t), J_add...] + return Ht, J, dagger.(J) +end +nothing # hide +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +T = collect(0.0:0.005:1.0) .* Tend +ψ0 = tensor([nlevelstate(ba, 1) for _ = 1:N]...) +t, ρt = timeevolution.master_dynamic(T, ψ0, input_output) +nothing # hide +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +I_R = zeros(length(t)) +I_L = zeros(length(t)) + +for (i, ti) in enumerate(t) + LR = L_R_QO(ti) + LL = L_L_QO(ti) + I_R[i] = real(expect(LR' * LR, ρt[i])) + I_L[i] = real(expect(LL' * LL, ρt[i])) +end +nothing # hide +```` + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +p = plot(t, I_R; label = "Transmission") +plot!(p, t, I_L; label = "Reflection") +plot!(p, t, abs2.(Ein_t.(t)); color = :grey, ls = :dash, label = "Input") +plot!( + p; + xlabel = "time", + ylabel = "intensity", + legend = :best, + grid = true, + size = (500, 320), +) +p +```` + +## Package versions + +These results were obtained using the following versions: + +````@example 05-3_N-QDs_bidirectional-waveguide_feedback-reduction +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status( + ["QuantumInputOutput", "SecondQuantizedAlgebra", "QuantumOptics", "Plots"], + mode = PKGMODE_MANIFEST, +) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + 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 index f25ae6c..2ba1b03 100644 --- 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 @@ -33,25 +33,25 @@ av_sym = Destroy(h, :a_v, 3) gu, γ, gv = rnumbers("gu γ gv") # cascade the SLH elements -G_u = SLH(1, gu * au_sym, 0) +G_u = SLH(1, gu' * au_sym, 0) G_s = SLH(1, sqrt(γ) * σ12_sym, 0) -G_v = SLH(1, gv * av_sym, 0) +G_v = SLH(1, gv' * av_sym, 0) G_cas = ▷(G_u, G_s, G_v) nothing # hide ```` ````@example 06-1_interaction-picture__PRA2023_107-013706_fig2 -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) ```` ````@example 06-1_interaction-picture__PRA2023_107-013706_fig2 -L = get_lindblad(G_cas)[1] +L = lindblad(G_cas)[1] ```` Usually we deal with the above derived Hamiltonian and Lindblad. In this example, however, we transform the system into the interaction picture of the virtual cavity-cavity interaction Hamiltonian $H_{uv}$. ````@example 06-1_interaction-picture__PRA2023_107-013706_fig2 -H_uv = get_hamiltonian(▷(G_u, G_v)) +H_uv = hamiltonian(▷(G_u, G_v)) ```` To do so, we first subtract $H_{uv}$ from $H$ and then replace the virtual cavity operators $a_v$ and $a_u$ as described in the [Theory](@ref) section. @@ -94,12 +94,12 @@ T_end = 12.0 T = [0:0.005:1;] * T_end # virtual-cavity couplings for input/output modes -gu_t = u_to_gu(u, T) -gv_t = v_to_gv(u, T) # identical output mode v(t) = u(t) +gu_t = coupling_input(u, T) +gv_t = coupling_output(u, T) # identical output mode v(t) = u(t) # interaction-picture coefficient matrix M(t) for u ↔ v -A_uv = interaction_picture_A_2modes(gu_t, gv_t) -M_t = interaction_picture_M(A_uv, T) +A_uv = coupling_matrix((gu_t, gv_t)) +M_t = solve_mode_evolution(A_uv, T) # constant and time-dependent parameters dict_p = Dict(γ => γ_) diff --git a/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md b/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md index 4bfb42e..e5df7d7 100644 --- a/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md +++ b/docs/src/examples/07-1_beamsplitter_loss__quantum-pulse.md @@ -41,10 +41,10 @@ G_p = SLH(1, 0, 0) S_bs = [r t; t -r] # input cavity, beam splitter, and output cavity -G_u = SLH(1, gu * au, 0) +G_u = SLH(1, gu' * au, 0) G_in = G_u ⊞ G_p G_bs = SLH(S_bs, [0, 0], 0) -G_v = SLH(1, gv * av, 0) +G_v = SLH(1, gv' * av, 0) # G_out = G_v ⊞ G_p # reflection is tracked G_out = G_p ⊞ G_v # transmission is tracked @@ -53,11 +53,16 @@ nothing # hide ```` ````@example 07-1_beamsplitter_loss__quantum-pulse -H = get_hamiltonian(G) +H = hamiltonian(G) ```` ````@example 07-1_beamsplitter_loss__quantum-pulse -L = get_lindblad(G) +L = lindblad(G) +L[1] +```` + +````@example 07-1_beamsplitter_loss__quantum-pulse +L[2] ```` ````@example 07-1_beamsplitter_loss__quantum-pulse @@ -71,8 +76,8 @@ T = [0:0.004:1;] * T_end ΔT = T[2] - T[1] # time-dependent coupling for the virtual cavities -gu_t = u_to_gu(u, T) -gv_t = v_to_gv(u, T) # v(t) = u(t) +gu_t = coupling_input(u, T) +gv_t = coupling_output(u, T) # v(t) = u(t) dict_p_t = Dict(gu => gu_t, gv => gv_t) # beam splitter parameters 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 index aa49f80..fa8908e 100644 --- a/docs/src/examples/07-2_hong-ou-mandel__quantum-pulse.md +++ b/docs/src/examples/07-2_hong-ou-mandel__quantum-pulse.md @@ -38,23 +38,28 @@ nothing # hide ````@example 07-2_hong-ou-mandel__quantum-pulse # input cavities, beam splitter, and output cavities S_bs = [r t; t -r] -G_u1 = SLH(1, gu1 * au1, 0) -G_u2 = SLH(1, gu2 * au2, 0) +G_u1 = SLH(1, gu1' * au1, 0) +G_u2 = SLH(1, gu2' * au2, 0) G_in = G_u1 ⊞ G_u2 G_bs = SLH(S_bs, [0, 0], 0) -G_v1 = SLH(1, gv1 * av1, 0) -G_v2 = SLH(1, gv2 * av2, 0) +G_v1 = SLH(1, gv1' * av1, 0) +G_v2 = SLH(1, gv2' * av2, 0) G_out = G_v1 ⊞ G_v2 G = G_in ▷ G_bs ▷ G_out nothing # hide ```` ````@example 07-2_hong-ou-mandel__quantum-pulse -H = get_hamiltonian(G) +H = hamiltonian(G) ```` ````@example 07-2_hong-ou-mandel__quantum-pulse -L = get_lindblad(G) +L = lindblad(G) +L[1] +```` + +````@example 07-2_hong-ou-mandel__quantum-pulse +L[2] ```` ````@example 07-2_hong-ou-mandel__quantum-pulse @@ -71,13 +76,13 @@ T = [0:0.002:1;] * T_end # time-dependent couplings for input and output modes u1 = u u2 = u -gu1_t = u_to_gu(u1, T) -gu2_t = u_to_gu(u2, T) +gu1_t = coupling_input(u1, T) +gu2_t = coupling_input(u2, T) v1(t) = u(t) v2(t) = u(t) -gv1_t = v_to_gv(v1, T) -gv2_t = v_to_gv(v2, T) +gv1_t = coupling_output(v1, T) +gv2_t = coupling_output(v2, T) dict_p_t = Dict(gu1 => gu1_t, gu2 => gu2_t, gv1 => gv1_t, gv2 => gv2_t) diff --git a/docs/src/examples/08-1_pulse-delay__simple.md b/docs/src/examples/08-1_pulse-delay__simple.md index 4545576..5335391 100644 --- a/docs/src/examples/08-1_pulse-delay__simple.md +++ b/docs/src/examples/08-1_pulse-delay__simple.md @@ -36,23 +36,23 @@ gu, gin, gout, gv = cnumbers("g_u g_in g_out g_v") nothing # hide ```` -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$). +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. ````@example 08-1_pulse-delay__simple -G_u = SLH(1, gu*au, 0) +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_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) +H = hamiltonian(G_cas) +L = lindblad(G_cas) nothing # hide ```` @@ -68,10 +68,10 @@ 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) +gu_ = coupling_input(u, T) +gout_ = coupling_delay_out(u_del, u, T) +gin_ = coupling_delay_in(u_del, u, T) +gv_ = coupling_output(u_del, T) dict_p_t = Dict([gu, gout, gin, gv] .=> [gu_, gout_, gin_, gv_]) nothing # hide @@ -138,7 +138,7 @@ delay between different modes is crucial, e.g. for two arms of an interferometer ````@example 08-1_pulse-delay__simple G_d_in = SLH(S2, [gin*ad, 0], 0) -H_ud = get_hamiltonian(cascade(G_u2, G_d_in)) +H_ud = hamiltonian(cascade(G_u2, G_d_in)) H_int_sym_ = simplify(H - H_ud) M(i, j) = cnumber("M_{$(i)$(j)}") @@ -156,6 +156,11 @@ H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict)) ````@example 08-1_pulse-delay__simple L_int_sym = simplify.(substitute_operators.(L, Ref(int_dict))) +L_int_sym[1] +```` + +````@example 08-1_pulse-delay__simple +L_int_sym[2] ```` ````@example 08-1_pulse-delay__simple @@ -170,17 +175,17 @@ 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) +gu_ = coupling_input(u, T) +gout_ = coupling_delay_out(u_del, u, T) +gin_ = coupling_delay_in(u_del, u, T) +gv_ = coupling_output(u_del, T) nothing # hide ```` ````@example 08-1_pulse-delay__simple # 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) +A_ud = coupling_matrix((gu_, gin_)) +M_t = solve_mode_evolution(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] 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 new file mode 100644 index 0000000..05ab768 --- /dev/null +++ b/docs/src/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.md @@ -0,0 +1,96 @@ +```@meta +EditURL = "../../../examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl" +``` + +# Coherent-Feedback Squeezing with a Beam-Splitter Loop + +This example implements the feedback loop to enhance field squeezing using coherent feedback described in +[J. Gough and S. Wildfeuer, PRA 80, 042107 (2009)](http://dx.doi.org/10.1103/PhysRevA.80.042107), +see also Example VI.1 of [J. Combes, et al. Advances in Physics: X, 2:3, 784-888 (2017)](https://doi.org/10.1080/23746149.2017.1343097). +A degenerate parametric oscillator is concatenated with a beam splitter, +and the internal wires are eliminated with the SLH feedback reduction rule. + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +using QuantumInputOutput +using SecondQuantizedAlgebra +using SymbolicUtils +using Plots +```` + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +# symbolic Hilbert space +hc = FockSpace(:c) + +# symbolic operator +a = Destroy(hc, :a, 1) + +# symbolic parameters +κ = rnumber("κ") +ϵ = rnumber("ϵ") +η = rnumber("η") +nothing # hide +```` + +The open-loop OPO has one port, while the beam splitter has two ports. We first +form the unconnected network and then apply the two feedback reductions. + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +r = √(1 - η^2) +G_opo = SLH(1, √(κ) * a, 1im * ϵ * (a'^2 - a^2)) +G_bs = SLH([-r η; η r], [0, 0], 0) +G_unconnected = G_opo ⊞ G_bs +G_loop = feedback(G_unconnected, 1 => 2, 2 => 1) +nothing # hide +```` + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +S_loop = scattering(G_loop) +L_loop = lindblad(G_loop)[1] +```` + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +H_loop = hamiltonian(G_loop) +```` + +The feedback loop leaves the OPO Hamiltonian unchanged but rescales the +coupling operator by $l = \eta / (1 + \sqrt{1-\eta^2})$. +For a numerical illustration, we evaluate the effective damping factor $l^2$ as +a function of the beam-splitter transmission coefficient. + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +η_grid = collect(0.0:0.005:0.999) +l2_grid = @. (η_grid / (1 + sqrt(1 - η_grid^2)))^2 + +p = plot( + η_grid, + l2_grid; + lw = 2, + label = "", + xlabel = "η", + ylabel = "effective damping fraction", + grid = true, + size = (520, 320), +) +hline!(p, [1.0]; color = :grey, ls = :dash, label = "open loop") +p +```` + +## Package versions + +These results were obtained using the following versions: + +````@example 09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009 +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status( + ["QuantumInputOutput", "SecondQuantizedAlgebra", "Plots"], + mode = PKGMODE_MANIFEST, +) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/implementation.md b/docs/src/implementation.md index 91f1582..df5b446 100644 --- a/docs/src/implementation.md +++ b/docs/src/implementation.md @@ -21,18 +21,18 @@ av = Destroy(h, :a_v, 3) An SLH component is represented as `(S, L, H)` by the [`SLH`](@ref) type. The cascade [`▷`](@ref), concatenation [`⊞`](@ref), and feedback reduction [`feedback`](@ref) rules implement the standard network composition from the SLH framework. -The resulting effective operators are accessed by [`get_hamiltonian`](@ref) and [`get_lindblad`](@ref) and remain symbolic until translation. This is especially useful when you want to further manipulate the expressions, e.g. to transform into the interaction picture. +The resulting effective operators are accessed by [`hamiltonian`](@ref) and [`lindblad`](@ref) and remain symbolic until translation. This is especially useful when you want to further manipulate the expressions, e.g. to transform into the interaction picture. ```julia G_cas = ▷(G_u, G_s, G_v) -H = get_hamiltonian(G_cas) -L = get_lindblad(G_cas) +H = hamiltonian(G_cas) +L = lindblad(G_cas) ``` For networks with internal loops, the symbolic model can be reduced directly with [`feedback`](@ref), which applies the SLH feedback reduction rule before translation. This keeps the symbolic workflow consistent: build a network from cascades and concatenations, eliminate internal connections symbolically, and only then translate the reduced Hamiltonian and Lindblad operators to numerics. -If you directly want to use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions, the [`SLHqo`](@ref) type skips the symbolic layer and allows time-dependent `L` or `H` as callables while still using the same cascade, concatenate, and feedback rules. This can be much faster. +If you directly want to use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions, you can pass numeric operators and callables directly to [`SLH`](@ref), which supports both symbolic and numeric operator types as well as time-dependent `L` or `H` while still using the same cascade, concatenate, and feedback rules. This can be much faster. ## Translate to numerics @@ -51,7 +51,7 @@ bv = FockBasis(2) b = bu ⊗ bs ⊗ bv dict_p = Dict(γ => 1.0, gv => 0.0) -gu_t = u_to_gu(t -> exp(-t^2), 0:0.01:5) +gu_t = coupling_input(t -> exp(-t^2), 0:0.01:5) dict_p_t = Dict(gu => gu_t) H_QO = translate_qo(H, b; parameter=dict_p, time_parameter=dict_p_t) @@ -71,8 +71,8 @@ Quantum pulses are encoded through virtual cavities. Given a normalized input mo The implementation uses cumulative numerical integration on a time grid `T` and a linear interpolation. -- [`u_to_gu`](@ref) and [`v_to_gv`](@ref) build interpolated couplings from sampled modes -- [`u_to_gu_Gauss`](@ref) and [`v_to_gv_Gauss`](@ref) provide analytic expressions for Gaussian pulses +- [`coupling_input`](@ref) and [`coupling_output`](@ref) build interpolated couplings from sampled modes +- For Gaussian pulses, pass a [`Gaussian`](@ref) descriptor to [`coupling_input`](@ref) or [`coupling_output`](@ref) for analytic expressions ```julia T = [0:0.002:1;]*12 @@ -80,11 +80,11 @@ T = [0:0.002:1;]*12 τ = 4.0 u(t) = 1/(√(σ)*π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2) -gu_t = u_to_gu(u, T) -gv_t = v_to_gv(u, T) +gu_t = coupling_input(u, T) +gv_t = coupling_output(u, T) ``` -For multiple input/output modes the distortion of the pulse due to the subsequent/preceding virtual cavities needs to be taken into account. The effective input mode ``u_i^{\mathrm{eff}}(t)`` and output mode ``v_i^{\mathrm{eff}}(t)`` for the virtual cavity `i` are constructed via [`u_eff`](@ref) and [`v_eff`](@ref). +For multiple input/output modes the distortion of the pulse due to the subsequent/preceding virtual cavities needs to be taken into account. The effective input mode ``u_i^{\mathrm{eff}}(t)`` and output mode ``v_i^{\mathrm{eff}}(t)`` for the virtual cavity `i` are constructed via [`effective_input_mode`](@ref) and [`effective_output_mode`](@ref). ## Output modes and the correlation function @@ -94,13 +94,13 @@ The dominant output modes are extracted by computing the two-time correlation ma g^{(1)}(t_1, t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle ``` -and diagonalizing it. In this package, [`two_time_corr_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $L_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $L_s(t_1) \rho(t_1)$ and use this as the initial "state" for the propagation of $t_2$, with the same Hamiltonian and Lindblad terms. +and diagonalizing it. In this package, [`correlation_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $L_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $L_s(t_1) \rho(t_1)$ and use this as the initial "state" for the propagation of $t_2$, with the same Hamiltonian and Lindblad terms. The eigenvectors of the matrix $g^{(1)}(t_1, t_2)$ correspond to temporal modes and the eigenvalues to their mean photon-number weights. The full procedure is illustrated in the [Tutorial](@ref). ```julia Ls(t) = gu_t(t) * au_qo + √(1.0) * c_qo -g1 = two_time_corr_matrix(T, ρt, input_output_1, Ls) +g1 = correlation_matrix(T, ρt, input_output_1, Ls) F = eigen(g1) v_mode = F.vectors[:, end] / sqrt(T[2] - T[1]) ``` @@ -121,7 +121,7 @@ A(t) = \frac{1}{2}\begin{bmatrix} \end{bmatrix} ``` -The functions [`interaction_picture_A_2modes`](@ref), [`interaction_picture_A_3modes`](@ref), and [`interaction_picture_A_4modes`](@ref) build the coupling matrices for two, three and four modes. For two equal modes ($u(t) = v(t)$) the analytic solution of $M(t)$ is provided by [`interaction_picture_M_2modes_equal`](@ref). +The function [`coupling_matrix`](@ref) builds the coupling matrix for an arbitrary number of modes. The ODE for the coefficient matrix $M(t)$ is solved by [`solve_mode_evolution`](@ref). For two equal modes ($u(t) = v(t)$) the analytic solution of $M(t)$ is provided by [`solve_mode_evolution_symmetric`](@ref). Using the interaction picture for scattering with a Fock state $| n=20 \rangle$ on a two-level system, is demonstrated in the example [Interaction Picture Scattering with a Quantum Pulse](examples/06-1_interaction-picture__PRA2023_107-013706_fig2.md). @@ -129,8 +129,8 @@ Using the interaction picture for scattering with a Fock state $| n=20 \rangle$ Pulse propagation delays are modeled by a virtual delay cavity that absorbs a pulse `v(t)` while emitting a target pulse `u(t)` with a controlled delay. The functions -- [`uv_to_gin`](@ref) -- [`uv_to_gout`](@ref) +- [`coupling_delay_in`](@ref) +- [`coupling_delay_out`](@ref) compute the in-coupling and out-coupling strengths `g_in(t)` and `g_out(t)` for this delay cavity, according to @@ -140,6 +140,6 @@ compute the in-coupling and out-coupling strengths `g_in(t)` and `g_out(t)` for \tilde g_{\mathrm{in},v,u}(t) = -\frac{v^*(t)}{\sqrt{\int_0^t dt' \,|v(t')|^2 - \int_0^t dt' \,|u(t')|^2}}. ``` -The implementation mirrors [`u_to_gu`](@ref) and [`v_to_gv`](@ref): compute cumulative integrals on a time grid and return interpolated time-dependent couplings. +The implementation mirrors [`coupling_input`](@ref) and [`coupling_output`](@ref): compute cumulative integrals on a time grid and return interpolated time-dependent couplings. A simple single-pulse case is demonstrated in the example [Simple Pulse Delay with a Virtual Cavity](examples/08-1_pulse-delay__simple.md), where an input pulse is emitted, delayed by a virtual cavity, and captured into an output cavity. diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 77e84e5..f62a07b 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -48,11 +48,11 @@ nothing # hide ``` ```@example tutorial -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) ``` ```@example tutorial -L = get_lindblad(G_cas)[1] +L = lindblad(G_cas)[1] ``` ## 2. Numerical parameters and input pulse @@ -74,7 +74,7 @@ u1(t) = 1/(sqrt(σ)*π^(1/4)) * exp( -(t - 4σ)^2 / (2*σ^2) ) T = [0:0.002:1;]*T_end ΔT = T[2] - T[1] -gu_t = u_to_gu(u1, T) +gu_t = coupling_input(u1, T) dict_p_t = Dict(gu => gu_t) nothing # hide ``` @@ -119,7 +119,7 @@ c_qo = translate_qo(c, b) av_qo = translate_qo(av, b) Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo -g1_m = two_time_corr_matrix(T, ρt, input_output_1, Ls) +g1_m = correlation_matrix(T, ρt, input_output_1, Ls) nothing # hide ``` @@ -155,7 +155,7 @@ nothing # hide Finally, we treat the dominant output mode explicitly by providing `g_v(t)` as a time-dependent parameter, and propagate the system again. ```@example tutorial -gv_t = v_to_gv(v_mode, T) +gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]) dict_p_2 = Dict([γ, Δ] .=> [γ_, Δ_]) 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 1ec95a5..6d29045 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 @@ -40,11 +40,11 @@ nothing # hide # -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) # -L = get_lindblad(G_cas)[1] # only one Lindblad term in this example +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. @@ -64,7 +64,7 @@ u1(t) = 1/(sqrt(σ)*π^(1/4)) * exp(-(t - 4σ)^2 / (2*σ^2)) T = [0:0.002:1;]*T_end ΔT = T[2] - T[1] -gu_t = u_to_gu(u1, T) +gu_t = coupling_input(u1, T) dict_p_t = Dict(gu => gu_t) ## numeric bases @@ -110,7 +110,7 @@ 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$. Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo -g1_m = two_time_corr_matrix(T, ρt, input_output_1, Ls) +g1_m = correlation_matrix(T, ρt, input_output_1, Ls) nothing # hide # @@ -143,7 +143,7 @@ p_num_2 = [γ_, Δ_] dict_p_2 = Dict(p_sym_2 .=> p_num_2) ## time-dependent coupling for the output mode $v(t)$ -gv_t = v_to_gv(v_mode, T) +gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]) 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 index 6fff185..46ff752 100644 --- a/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl +++ b/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl @@ -37,11 +37,11 @@ nothing # hide # -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) # -L = get_lindblad(G_cas)[1] # only one Lindblad in this example +L = lindblad(G_cas)[1] # only one Lindblad in this example # @@ -61,7 +61,7 @@ u(t) = 1/(sqrt(τ)*π^(1/4)) * exp(-(t - tp)^2 / (2*τ^2)) T = [0:0.002:1;]*14τ ΔT = T[2] - T[1] -gu_ = u_to_gu(u, T) +gu_ = coupling_input(u, T) dict_p_t = Dict(gu => gu_) ## numeric bases @@ -93,7 +93,7 @@ nothing # hide # 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))'*au_qo + √(γ_)*c_qo -g1_m = two_time_corr_matrix(T, ρt, input_output, Ls); +g1_m = correlation_matrix(T, ρt, input_output, Ls); nothing # hide # diff --git a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl index ef8c0c1..f9f9814 100644 --- a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl +++ b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl @@ -47,8 +47,8 @@ G_v1 = SLH(1, gv1' * av1, 0) G_v2 = SLH(1, gv2' * av2, 0) G = cascade(G_s, G_v1, G_v2) -H = get_hamiltonian(G) -L = get_lindblad(G)[1] +H = hamiltonian(G) +L = lindblad(G)[1] nothing # hide # We use the parameters quoted in the paper and initialize the emitter in the excited state. @@ -93,7 +93,7 @@ t_1, ρt_1 = timeevolution.master_dynamic(T, ψ0, input_output_1) nothing # hide Ls = √(γ_) * a_qo -g1_m = two_time_corr_matrix(T, ρt_1, input_output_1, Ls) +g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) F = eigen(g1_m) n_avg = real.(F.values) * ΔT @@ -106,11 +106,11 @@ nothing # hide # After identifying the two modes, we add two cascaded virtual output cavities. # The coupling of the second cavity must be corrected for the reshaping caused by -# the first output cavity, which is done with [`v_eff`](@ref). +# the first output cavity, which is done with [`effective_output_mode`](@ref). -gv1_t = v_to_gv(v1_mode, T) -v2_eff = v_eff([v1_mode, v2_mode], T, 2) -gv2_t = v_to_gv(v2_eff, T) +gv1_t = coupling_output(v1_mode, T) +v2_eff = effective_output_mode([v1_mode, v2_mode], T, 2) +gv2_t = coupling_output(v2_eff, T) dict_p_2 = Dict([γ, g, ω12] .=> [γ_, g_, ω12_]) dict_p_t_2 = Dict(gv1 => gv1_t, gv2 => gv2_t) 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 index c94e7f4..7d2d021 100644 --- a/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl +++ b/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl @@ -45,11 +45,11 @@ nothing # hide # -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) # -L = get_lindblad(G_cas)[1] # only one Lindblad in this example +L = lindblad(G_cas)[1] # only one Lindblad in this example # Next, the numerical parameters and functions of the system are defined. @@ -67,7 +67,7 @@ u(t) = 1/(sqrt(τ)*π^(1/4)) * exp(-(t - tp)^2 / (2*τ^2)) T = [0:0.002:1;]*20 ΔT = T[2] - T[1] -gu_ = u_to_gu(u, T) +gu_ = coupling_input(u, T) dict_p_t = Dict(gu1 => gu_) nothing # hide @@ -107,7 +107,7 @@ au1_qo = translate_qo(au1, b) σ_qo(i, j) = translate_qo(σ(i, j), b) Ls(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) -g1_m = two_time_corr_matrix(T, ρt, input_output, Ls); +g1_m = correlation_matrix(T, ρt, input_output, Ls); p = heatmap( T, @@ -154,20 +154,20 @@ u2_new = conj.(reverse(v2_p)) nothing # hide # The pulse from input cavity $u_2$ is scattered on the cavity $u_1$. This distortion needs to be taken into account for the coupling of $u_2$, -# which is done with the function [`u_eff`](@ref). +# which is done with the function [`effective_input_mode`](@ref). # The coupling of the $u_1$ cavity needs no adaptation, since it directly couples to the two-level system. -gu1_ = u_to_gu(u1_new, T) +gu1_ = coupling_input(u1_new, T) u_new_data = [u1_new, u2_new] u_new_fct = [LinearInterpolation(u, T) for u in u_new_data] ## effective u2 mode and corresponding coupling -u2_for_gu2 = u_eff(u_new_fct, T, 2) -gu2_ = u_to_gu(u2_for_gu2, T) +u2_for_gu2 = effective_input_mode(u_new_fct, T, 2) +gu2_ = coupling_input(u2_for_gu2, T) ## coupling of the output mode -gv1_ = v_to_gv(v1_new, T) +gv1_ = coupling_output(v1_new, T) ## dictionary for the time-dependent functions g_sym = [gu1, gu2, gv1] 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 index 606bbbc..f348fce 100644 --- a/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl +++ b/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl @@ -41,15 +41,15 @@ nothing # hide # Note that one needs to be careful to not double-count the Hamiltonian terms with the concatenation rule. -H1 = get_hamiltonian(G_cav_L_R_drive) +H1 = hamiltonian(G_cav_L_R_drive) # -L1_L = get_lindblad(G_cav_L_R_drive)[1] +L1_L = lindblad(G_cav_L_R_drive)[1] # -L1_R = get_lindblad(G_cav_L_R_drive)[2] +L1_R = lindblad(G_cav_L_R_drive)[2] # Here, the usual classical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ appears as a combination of Hamiltonian and Lindblad term. # To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl). Since we do not want to include the basis of the atoms, we provide a dictionary of operators with the kwarg `operators` in the function [`translate_qo`](@ref). 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 index 85325da..acf2a93 100644 --- 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 @@ -40,15 +40,15 @@ nothing # hide # Note that one needs to be careful to not double-count the Hamiltonian terms with the concatenation rule. -H1 = get_hamiltonian(G_cav_L_R_drive) +H1 = hamiltonian(G_cav_L_R_drive) # -L1_L = get_lindblad(G_cav_L_R_drive)[1] +L1_L = lindblad(G_cav_L_R_drive)[1] # -L1_R = get_lindblad(G_cav_L_R_drive)[2] +L1_R = lindblad(G_cav_L_R_drive)[2] # The typical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ is a combination of Hamiltonian term and Lindblad. # We use the function `meanfield` to obtain the equation for the intra-cavity field, which leads to a closed set of equations in this particular case. 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 31609e5..56081fa 100644 --- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl +++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl @@ -49,11 +49,11 @@ G_L_t = G_L(2) ▷ G_ϕ(1, 2) ▷ G_L(1) G_t = G_R_t ⊞ G_L_t nothing # hide -H = get_hamiltonian(G_t) +H = hamiltonian(G_t) # -L = get_lindblad(G_t) +L = lindblad(G_t) L_R = L[1] # 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 index e59b26d..7275e3b 100644 --- a/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.jl +++ b/examples/05-2_N-QDs_bidirectional-waveguide_quantum-pulse_qo.jl @@ -1,6 +1,6 @@ # # Quantum Pulse Bi-Directional Waveguide # -# This example mirrors the example `Bi-Directional Waveguide` but uses the numeric SLH struct [`SLHqo`](@ref), to circumvent the symbolic derivation part. +# This example mirrors the example `Bi-Directional Waveguide` but passes numeric operators directly to [`SLH`](@ref), circumventing the symbolic derivation part. # Furthermore it drives the system with a *quantum* single-photon pulse via a virtual cavity. # As usual, we start by loading the packages and defining the operators and parameters of the system. @@ -43,15 +43,15 @@ T = [0:0.005:1;]*Tend ## u1(t) = sqrt(1 / (σt * √(2π)) * exp(-0.5 * (t - t0)^2 / σt^2)) # hide u1(t) = 1/(sqrt(σt)*π^(1/4)) * exp(-(t - t0)^2 / (2*σt^2)) -gu_t = u_to_gu(u1, T) +gu_t = coupling_input(u1, T) nothing # hide -# We use the [`SLHqo`](@ref) to directly use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions to the model the system. +# We use the [`SLH`](@ref) to directly use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions to the model the system. -G_u = SLHqo(1, t -> gu_t(t) * a_u, 0*one(b)) -G_ϕ(i, j) = SLHqo(exp(1im * ϕn[i]), 0*one(b), 0*one(b)) -G_R(i) = SLHqo(1, √(γRn[i]) * σ(i, 1, 2), -Δn[i] * σ(i, 2, 2)) -G_L(i) = SLHqo(1, √(γLn[i]) * σ(i, 1, 2), 0*one(b)) +G_u = SLH(1, t -> gu_t(t) * a_u, 0*one(b)) +G_ϕ(i, j) = SLH(exp(1im * ϕn[i]), 0*one(b), 0*one(b)) +G_R(i) = SLH(1, √(γRn[i]) * σ(i, 1, 2), -Δn[i] * σ(i, 2, 2)) +G_L(i) = SLH(1, √(γLn[i]) * σ(i, 1, 2), 0*one(b)) ## Cascade right-moving channel ## G_R_t = G_d ▷ cascade([G_R(i) ▷ G_ϕ(i, i + 1) for i=1:N-1]...) ▷ G_R(N) # for N > 1 # hide @@ -66,8 +66,8 @@ nothing # hide # The full Hamiltonian and Lindblad terms are extracted from the final SLH element. Note that as soon as one time-dependent function is involved in a cascade or concatenate, the returned $H$ and $L$ will also be time-dependent. -H = get_hamiltonian(G_t) -L = get_lindblad(G_t) +H = hamiltonian(G_t) +L = lindblad(G_t) L_R = L[1] L_L = L[2] diff --git a/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl b/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl index a7aecf2..0939f63 100644 --- a/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl +++ b/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl @@ -48,8 +48,8 @@ nothing # hide # first Lindblad operator corresponds to the reflected left-moving output and the # second one to the transmitted right-moving output. -H = get_hamiltonian(G_t) -L = get_lindblad(G_t) +H = hamiltonian(G_t) +L = lindblad(G_t) L_L = L[1] L_R = L[2] nothing # hide diff --git a/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl b/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl index 77bf436..3fae594 100644 --- a/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl +++ b/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl @@ -36,15 +36,15 @@ nothing # hide # -H = get_hamiltonian(G_cas) +H = hamiltonian(G_cas) # -L = get_lindblad(G_cas)[1] +L = lindblad(G_cas)[1] # Usually we deal with the above derived Hamiltonian and Lindblad. In this example, however, we transform the system into the interaction picture of the virtual cavity-cavity interaction Hamiltonian $H_{uv}$. -H_uv = get_hamiltonian(▷(G_u, G_v)) +H_uv = hamiltonian(▷(G_u, G_v)) # To do so, we first subtract $H_{uv}$ from $H$ and then replace the virtual cavity operators $a_v$ and $a_u$ as described in the [Theory](@ref) section. @@ -83,12 +83,12 @@ T_end = 12.0 T = [0:0.005:1;] * T_end ## virtual-cavity couplings for input/output modes -gu_t = u_to_gu(u, T) -gv_t = v_to_gv(u, T) # identical output mode v(t) = u(t) +gu_t = coupling_input(u, T) +gv_t = coupling_output(u, T) # identical output mode v(t) = u(t) ## interaction-picture coefficient matrix M(t) for u ↔ v -A_uv = interaction_picture_A_2modes(gu_t, gv_t) -M_t = interaction_picture_M(A_uv, T) +A_uv = coupling_matrix((gu_t, gv_t)) +M_t = solve_mode_evolution(A_uv, T) ## constant and time-dependent parameters dict_p = Dict(γ => γ_) diff --git a/examples/07-1_beamsplitter_loss__quantum-pulse.jl b/examples/07-1_beamsplitter_loss__quantum-pulse.jl index dcdd84e..29ea169 100644 --- a/examples/07-1_beamsplitter_loss__quantum-pulse.jl +++ b/examples/07-1_beamsplitter_loss__quantum-pulse.jl @@ -46,11 +46,11 @@ nothing # hide # -H = get_hamiltonian(G) +H = hamiltonian(G) # -L = get_lindblad(G) +L = lindblad(G) L[1] # @@ -69,8 +69,8 @@ T = [0:0.004:1;] * T_end ΔT = T[2] - T[1] ## time-dependent coupling for the virtual cavities -gu_t = u_to_gu(u, T) -gv_t = v_to_gv(u, T) # v(t) = u(t) +gu_t = coupling_input(u, T) +gv_t = coupling_output(u, T) # v(t) = u(t) dict_p_t = Dict(gu => gu_t, gv => gv_t) ## beam splitter parameters diff --git a/examples/07-2_hong-ou-mandel__quantum-pulse.jl b/examples/07-2_hong-ou-mandel__quantum-pulse.jl index 972a57c..c1c75a0 100644 --- a/examples/07-2_hong-ou-mandel__quantum-pulse.jl +++ b/examples/07-2_hong-ou-mandel__quantum-pulse.jl @@ -45,11 +45,11 @@ nothing # hide # -H = get_hamiltonian(G) +H = hamiltonian(G) # -L = get_lindblad(G) +L = lindblad(G) L[1] # @@ -71,13 +71,13 @@ T = [0:0.002:1;] * T_end ## time-dependent couplings for input and output modes u1 = u u2 = u -gu1_t = u_to_gu(u1, T) -gu2_t = u_to_gu(u2, T) +gu1_t = coupling_input(u1, T) +gu2_t = coupling_input(u2, T) v1(t) = u(t) v2(t) = u(t) -gv1_t = v_to_gv(v1, T) -gv2_t = v_to_gv(v2, T) +gv1_t = coupling_output(v1, T) +gv2_t = coupling_output(v2, T) dict_p_t = Dict(gu1 => gu1_t, gu2 => gu2_t, gv1 => gv1_t, gv2 => gv2_t) diff --git a/examples/08-1_pulse-delay__simple.jl b/examples/08-1_pulse-delay__simple.jl index 7b81ae5..7c532c9 100644 --- a/examples/08-1_pulse-delay__simple.jl +++ b/examples/08-1_pulse-delay__simple.jl @@ -46,8 +46,8 @@ 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) +H = hamiltonian(G_cas) +L = lindblad(G_cas) nothing # hide # @@ -63,10 +63,10 @@ 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) +gu_ = coupling_input(u, T) +gout_ = coupling_delay_out(u_del, u, T) +gin_ = coupling_delay_in(u_del, u, T) +gv_ = coupling_output(u_del, T) dict_p_t = Dict([gu, gout, gin, gv] .=> [gu_, gout_, gin_, gv_]) nothing # hide @@ -131,7 +131,7 @@ p # 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_ud = hamiltonian(cascade(G_u2, G_d_in)) H_int_sym_ = simplify(H - H_ud) M(i, j) = cnumber("M_{$(i)$(j)}") @@ -168,17 +168,17 @@ 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) +gu_ = coupling_input(u, T) +gout_ = coupling_delay_out(u_del, u, T) +gin_ = coupling_delay_in(u_del, u, T) +gv_ = coupling_output(u_del, T) nothing # hide # ## 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) +A_ud = coupling_matrix((gu_, gin_)) +M_t = solve_mode_evolution(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] diff --git a/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl b/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl index 85905ed..9fd2514 100644 --- a/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl +++ b/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl @@ -37,12 +37,12 @@ nothing # hide # -S_loop = get_scattering(G_loop) -L_loop = get_lindblad(G_loop)[1] +S_loop = scattering(G_loop) +L_loop = lindblad(G_loop)[1] # -H_loop = get_hamiltonian(G_loop) +H_loop = hamiltonian(G_loop) # The feedback loop leaves the OPO Hamiltonian unchanged but rescales the # coupling operator by $l = \eta / (1 + \sqrt{1-\eta^2})$. 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 0fc4c03..804c46b 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 @@ -37,8 +37,8 @@ nothing # hide H_s = p / 2 * (a'^2 + a^2) - K / 2 * (a'^2) * (a^2) + Δ * a' * a G_s = SLH(1, √(γ) * a, H_s) -H = get_hamiltonian(G_s) -L = get_lindblad(G_s)[1] +H = hamiltonian(G_s) +L = lindblad(G_s)[1] nothing # hide # Next, we define the numerical parameters from Sec. II.D. The classical pump @@ -103,7 +103,7 @@ nothing # hide a_qo = destroy(bc) Ls(t) = √(γ_) * a_qo -g1_m = two_time_corr_matrix(T, ρt_1, input_output_1, Ls) +g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) F = eigen(g1_m) n_avg = real.(F.values) * ΔT @@ -126,10 +126,10 @@ G_s2 = SLH(1, √(γ) * a2, H_s2) G_v = SLH(1, gv * av, 0) G = G_s2 ▷ G_v -H_2 = get_hamiltonian(G) -L_2 = get_lindblad(G)[1] +H_2 = hamiltonian(G) +L_2 = lindblad(G)[1] -gv_t = v_to_gv(v_mode, T) +gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict(p => p_t, gv => gv_t) bv = FockBasis(n_cut) 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 c20121f..5389b90 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 @@ -60,8 +60,8 @@ L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] G_atom = SLH(I4, L_atom, Δ*σ(2, 2)) G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom -H = get_hamiltonian(G_u_bs_d1_d2_atom) -L = get_lindblad(G_u_bs_d1_d2_atom) +H = hamiltonian(G_u_bs_d1_d2_atom) +L = lindblad(G_u_bs_d1_d2_atom) # @@ -113,27 +113,17 @@ t1 = 2tp + τ + 2tw rn = 1/√(2) tn = 1/√(2) -gu_ = u_to_gu(u, T) -gd1in_ = v_to_gv(u, T) -gd1out_ = u_to_gu(ud1, T) -gd2in_ = v_to_gv(u, T) -gd2out_ = u_to_gu(ud2, T) +gu_ = coupling_input(u, T) +gd1in_ = coupling_output(u, T) +gd1out_ = coupling_input(ud1, T) +gd2in_ = coupling_output(u, T) +gd2out_ = coupling_input(ud2, T) dict_p_t = Dict([gu, gd1in, gd1out, gd2in, gd2out] .=> [gu_, gd1in_, gd1out_, gd2in_, gd2out_]) nothing # hide -function interaction_picture_A_3modes(g1, g2, g3) - A(t) = - 0.5 * [ - 0 conj(g2(t)) * g1(t) g1(t) * conj(g3(t)); - -g2(t) * conj(g1(t)) 0 0; - -conj(g1(t)) * g3(t) 0 0 - ] - return A -end - -A_ud = interaction_picture_A_3modes(gu_, gd1in_, gd2in_) -M_t = interaction_picture_M(A_ud, T) +A_ud = coupling_matrix((gu_, gd1in_, gd2in_)) +M_t = solve_mode_evolution(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] diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl index c586b11..228e252 100644 --- a/src/QuantumInputOutput.jl +++ b/src/QuantumInputOutput.jl @@ -9,38 +9,43 @@ using Symbolics using SpecialFunctions: erf using DataInterpolations: LinearInterpolation, ExtrapolationType using NumericalIntegration: cumul_integrate -using LinearAlgebra: LinearAlgebra, I -using OrdinaryDiffEq # for u_eff +using LinearAlgebra: LinearAlgebra, I, mul! +using OrdinaryDiffEq +using StaticArrays +using FunctionWrappers import QuantumOpticsBase: expect import SecondQuantizedAlgebra: numeric_average const SQA = SecondQuantizedAlgebra -export SLH, # SLH.jl - SLHqo, - get_scattering, - get_lindblad, - get_hamiltonian, +export SLH, + Gaussian, + scattering, + lindblad, + hamiltonian, + # Composition ▷, cascade, ⊞, concatenate, feedback, - translate_qo, # translate.jl - u_to_gu, # utils.jl - v_to_gv, - u_to_gu_Gauss, - v_to_gv_Gauss, - u_eff, - v_eff, - uv_to_gout, - uv_to_gin, - two_time_corr_matrix, - interaction_picture_M, # interaction picture - interaction_picture_M_2modes_equal, - interaction_picture_A_2modes, - interaction_picture_A_3modes, - interaction_picture_A_4modes, + # Translation + translate_qo, + # Pulse coupling + coupling_input, + coupling_output, + coupling_delay_in, + coupling_delay_out, + # Effective modes + effective_input_mode, + effective_output_mode, + # Interaction picture + coupling_matrix, + solve_mode_evolution, + solve_mode_evolution_symmetric, + # Correlations + correlation_matrix, + # Operators substitute_operators include("SLH.jl") diff --git a/src/SLH.jl b/src/SLH.jl index 9baaf41..fd68d9d 100644 --- a/src/SLH.jl +++ b/src/SLH.jl @@ -1,429 +1,431 @@ -""" - SLH +using StaticArrays +using FunctionWrappers: FunctionWrapper -SLH triple with scattering matrix `S`, Lindblad term `L` and Hamiltonian `H`. -`S` and `L` can also be vectors of scattering matrices and Linblad terms. +# NOTE: FunctionWrapper is callable but NOT <: Function. +const _Callable = Union{Function,FunctionWrapper} -See also [`▷`](@ref) and [`⊞`](@ref) -""" -struct SLH{T,LT,H,S<:AbstractMatrix{T},L<:AbstractVector{LT}} - scattering::S - lindblad::L - hamiltonian::H - - function SLH( - scattering::S, - lindblad::L, - hamiltonian::H, - ) where {T,LT,H,S<:AbstractMatrix{T},L<:AbstractVector{LT}} - - @assert size(scattering, 1) == length(lindblad) - new{T,LT,H,S,L}(scattering, lindblad, hamiltonian) - end -end +# ────────────────────────────────────────────── +# Dispatch helpers: symbolic vs numeric +# ────────────────────────────────────────────── + +_adj(x::SQA.QNumber) = SQA._adjoint(x) +_adj(f::_Callable) = t -> adjoint(f(t)) +_adj(x) = adjoint(x) + +_post(x::BasicSymbolic) = simplify(x) +_post(x) = x + +# ────────────────────────────────────────────── +# Arithmetic helpers: static vs time-dependent +# ────────────────────────────────────────────── + +_is_time_dep(x) = x isa _Callable +_to_func(x) = _is_time_dep(x) ? x : (t -> x) + +_add(f::_Callable, g::_Callable) = t -> f(t) + g(t) +_add(f::_Callable, x) = iszero(x) ? f : (t -> f(t) + x) +_add(x, f::_Callable) = iszero(x) ? f : (t -> x + f(t)) +_add(x, y) = x + y -# ----------------------- -# QuantumOptics operators -# ----------------------- +_mul(f::_Callable, g::_Callable) = t -> f(t) * g(t) +_mul(s, f::_Callable) = isone(s) ? f : (t -> s * f(t)) +_mul(f::_Callable, s) = isone(s) ? f : (t -> f(t) * s) +_mul(x, y) = x * y + +# ────────────────────────────────────────────── +# FunctionWrapper with concrete return type +# ────────────────────────────────────────────── """ - SLHqo + _detect_operator_type(L, H) -SLH triple where `L` and `H` are QuantumOptics.jl operators or time-dependent functions -returning such operators. This is useful when you build models directly in a numeric basis, -without symbolic translation. +Determine the concrete return type for FunctionWrapper by inspecting +the elements of L and H. Checks FunctionWrapper type parameters first +(no evaluation needed), then static element types. +Errors if only plain closures are present (the caller should thread +the type through from input SLH objects instead). """ -struct SLHqo{S<:AbstractMatrix,L<:AbstractVector,H} - scattering::S - lindblad::L - hamiltonian::H - - function SLHqo( - scattering::S, - lindblad::L, - hamiltonian::H, - ) where {S<:AbstractMatrix,L<:AbstractVector,H} - @assert size(scattering, 1) == length(lindblad) - if !(hamiltonian isa Function || hamiltonian isa QuantumOpticsBase.AbstractOperator) - error( - "SLHqo expects H to be a QuantumOptics operator or a time-dependent function returning one; got $(typeof(hamiltonian))", - ) - end - for (i, l) in pairs(lindblad) - if !(l isa Function || l isa QuantumOpticsBase.AbstractOperator) - error( - "SLHqo expects L[$i] to be a QuantumOptics operator or a time-dependent function returning one; got $(typeof(l))", - ) - end - end - new{S,L,H}(scattering, lindblad, hamiltonian) +function _detect_operator_type(L, H) + # Check FunctionWrapper elements (carry explicit type info) + for l in L + l isa FunctionWrapper && return _fw_return_type(typeof(l)) + end + H isa FunctionWrapper && return _fw_return_type(typeof(H)) + # Check static (non-time-dep) elements + for l in L + _is_time_dep(l) || return typeof(l) end + _is_time_dep(H) || return typeof(H) + # No type information available + error( + "Cannot determine concrete operator type: all elements are untyped closures. " * + "Wrap at least one with FunctionWrapper{OpType, Tuple{Float64}} or include a static element.", + ) end -function SLHqo(S, L::AbstractVector, H) - S_ = Matrix(I, length(L), length(L)) * S - return SLHqo(S_, L, H) +_fw_return_type(::Type{FunctionWrapper{R,A}}) where {R,A} = R + +function _maybe_wrap_lindblad(L::SVector{N}, ::Type{OpType}) where {N,OpType} + if any(_is_time_dep, L) + fw_type = FunctionWrapper{OpType,Tuple{Float64}} + return SVector{N,fw_type}(ntuple(i -> fw_type(_to_func(L[i])), Val(N))) + end + return L end -function SLHqo(S, L, H) - S_ = ones(1, 1) * S - L_ = [L] - return SLHqo(S_, L_, H) + +function _maybe_wrap_hamiltonian(H, has_td::Bool, ::Type{OpType}) where {OpType} + if has_td + return FunctionWrapper{OpType,Tuple{Float64}}(_to_func(H)) + end + return H end -""" - get_scattering(G) +# ────────────────────────────────────────────── +# SLH type +# ────────────────────────────────────────────── -Return the scattering matrix `S` of an SLH or SLHqo object. """ -get_scattering(slh::SLHqo) = slh.scattering -""" - get_lindblad(G) + SLH{N, ST, LT, HT} + +SLH triple with scattering matrix `S`, Lindblad vector `L`, and Hamiltonian `H`. +`S` and `L` can also be vectors of scattering matrices and Lindblad terms. -Return the Lindblad vector `L` of an SLH or SLHqo object. +See also [`▷`](@ref), [`⊞`](@ref), [`feedback`](@ref) """ -get_lindblad(slh::SLHqo) = slh.lindblad +struct SLH{N,ST,LT,HT} + scattering::SMatrix{N,N,ST} + lindblad::SVector{N,LT} + hamiltonian::HT +end + """ - get_hamiltonian(G) + _op_type(G::SLH) -Return the Hamiltonian `H` of an SLH or SLHqo object. +Extract the operator return type from an SLH with FunctionWrapper elements. +Returns `nothing` if no FunctionWrapper type info is available. """ -get_hamiltonian(slh::SLHqo) = slh.hamiltonian +function _op_type(::SLH{N,ST,LT,HT}) where {N,ST,LT,HT} + LT <: FunctionWrapper && return _fw_return_type(LT) + HT <: FunctionWrapper && return _fw_return_type(HT) + return nothing +end + +# ────────────────────────────────────────────── +# Constructors +# ────────────────────────────────────────────── -function Base.isequal(slh1::SLHqo, slh2::SLHqo) - isequal(get_scattering(slh1), get_scattering(slh2)) && - isequal(get_lindblad(slh1), get_lindblad(slh2)) && - isequal(get_hamiltonian(slh1), get_hamiltonian(slh2)) +# Canonical: from SMatrix + SVector (handles FunctionWrapper wrapping) +function _build_slh(S::SMatrix{N,N}, L::SVector{N}, H) where {N} + has_td = any(_is_time_dep, L) || _is_time_dep(H) + if has_td + OpType = _detect_operator_type(L, H) + return _build_slh(S, L, H, OpType) + end + return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) end -_is_time_dep(x) = x isa Function -_to_func(x) = x isa Function ? x : (t -> x) +# With explicit operator type (skips detection — used by composition operations) +function _build_slh(S::SMatrix{N,N}, L::SVector{N}, H, ::Type{OpType}) where {N,OpType} + has_td = any(_is_time_dep, L) || _is_time_dep(H) + if has_td + L_w = _maybe_wrap_lindblad(L, OpType) + H_w = _maybe_wrap_hamiltonian(H, true, OpType) + return SLH{N,eltype(S),eltype(L_w),typeof(H_w)}(S, L_w, H_w) + end + return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) +end -# SLH(scattering::S, lindblad::L, hamiltonian::H) where {S,H,L} = SLH{S,H,L}(scattering,lindblad,hamiltonian) -function SLH(S, L::AbstractVector, H) - S_ = Matrix{typeof(S)}(I, length(L), length(L))*S - return SLH(S_, L, H) +# Nothing hint falls through to detection +_build_slh(S::SMatrix{N,N}, L::SVector{N}, H, ::Nothing) where {N} = _build_slh(S, L, H) + +# From AbstractMatrix + AbstractVector (includes SMatrix + SVector) +function SLH(S::AbstractMatrix, L::AbstractVector, H) + N = length(L) + @assert size(S, 1) == N && size(S, 2) == N + return _build_slh(SMatrix{N,N}(S), SVector{N}(L...), H) +end + +# Numeric scalar S + vector L → S * I_{NxN} +function SLH(S::Number, L::AbstractVector, H) + N = length(L) + S_mat = SMatrix{N,N}(S * LinearAlgebra.I) + return _build_slh(S_mat, SVector{N}(L...), H) end + +# Scalar S + scalar L → SLH{1} function SLH(S, L, H) - S_ = ones(1, 1)*S - L_ = [L] - return SLH(S_, L_, H) + S_mat = SMatrix{1,1}(S) + L_vec = SVector{1}(L) + return _build_slh(S_mat, L_vec, H) end -function Base.isequal(slh1::SLH, slh2::SLH) - isequal(get_scattering(slh1), get_scattering(slh2)) && - isequal(get_lindblad(slh1), get_lindblad(slh2)) && - isequal(get_hamiltonian(slh1), get_hamiltonian(slh2)) +# Symbolic/general scalar S + vector L → S * I +function SLH(S, L::AbstractVector, H) + N = length(L) + S_mat = SMatrix{N,N}([i == j ? S : 0 for i = 1:N, j = 1:N]) + return _build_slh(S_mat, SVector{N}(L...), H) end -get_scattering(slh::SLH) = slh.scattering -get_lindblad(slh::SLH) = slh.lindblad -get_hamiltonian(slh::SLH) = slh.hamiltonian +# ────────────────────────────────────────────── +# Accessors +# ────────────────────────────────────────────── """ - ▷(G::SLH...) + scattering(G::SLH) -Creates a new SLH triple by cascading the SLH triples from first to last according to -the rule +Return the scattering matrix `S` of an SLH object. +""" +scattering(G::SLH) = G.scattering -``SLH_1 \\triangleright SLH_2 = ( S_2 S_1, L_2 + S_2 L_1, H_1 + H_2 - \\frac{i}{2} L_2^\\dagger S_2 L_1 - L_1^\\dagger S_2^\\dagger L_2 ) `` +""" + lindblad(G::SLH) -Unicode `\\triangleright` alias of [`cascade`](@ref) +Return the Lindblad vector `L` of an SLH object. """ -function ▷(G1::SLH, G2::SLH) #\triangleright - # function cascade(G1::SLH,G2::SLH) # G1 ▷ G2 = G2 ◁ G1 - S1 = get_scattering(G1); - L1 = get_lindblad(G1); - H1 = get_hamiltonian(G1) - S2 = get_scattering(G2); - L2 = get_lindblad(G2); - H2 = get_hamiltonian(G2) - - lL1 = length(L1) - lL2 = length(L2) - @assert lL1==lL2 - - S_t = simplify.(S2*S1) - L_t = simplify.(L2 + S2*L1) - # H_t = (H1 + H2 - 1im/2*(L2'S2*L1 - L1'S2'L2)) # the adjoint also creates the transpose! - # H_t = (H1 + H2 - 1im/2*(SQA._adjoint(L2)*S2*L1 - SQA._adjoint(L1)*SQA._adjoint(S2)*L2)) - L1_ct = SQA._adjoint.(L1) - L2_ct = SQA._adjoint.(L2) - S2_ct = SQA._adjoint.(S2) - - H_t = simplify( - H1 + H2 - - 1im/2*( - sum(L2_ct[it]*(S2*L1)[it] for it = 1:lL1) - - sum(L1_ct[it]*(S2_ct*L2)[it] for it = 1:lL1) - ), - ) - return SLH(S_t, L_t, H_t) -end -▷(a::SLH, b::SLH, c::SLH...) = ▷(a▷b, c...) -# G's in the order as they are cascaded (from left to right) +lindblad(G::SLH) = G.lindblad """ - ▷(G::SLHqo...) + hamiltonian(G::SLH) -Cascades SLHqo triples from first to last, allowing time-dependent `L` or `H`. -If any input `L` or `H` is time-dependent, the result uses time-dependent functions. +Return the Hamiltonian `H` of an SLH object. """ -function ▷(G1::SLHqo, G2::SLHqo) - S1 = get_scattering(G1); - L1 = get_lindblad(G1); - H1 = get_hamiltonian(G1) - S2 = get_scattering(G2); - L2 = get_lindblad(G2); - H2 = get_hamiltonian(G2) - - lL1 = length(L1) - lL2 = length(L2) - @assert lL1 == lL2 - - S_t = S2 * S1 - - L1f = map(_to_func, L1) - L2f = map(_to_func, L2) - H1f = _to_func(H1) - H2f = _to_func(H2) - time_dep = - any(_is_time_dep, L1) || - any(_is_time_dep, L2) || - _is_time_dep(H1) || - _is_time_dep(H2) - - if time_dep - function L_t(i) - return t -> begin - L1_vec = [f(t) for f in L1f] - L2_vec = [f(t) for f in L2f] - (L2_vec+S2*L1_vec)[i] - end +hamiltonian(G::SLH) = G.hamiltonian + + +# ────────────────────────────────────────────── +# Equality +# ────────────────────────────────────────────── + +function Base.isequal(G1::SLH, G2::SLH) + isequal(scattering(G1), scattering(G2)) && + isequal(lindblad(G1), lindblad(G2)) && + isequal(hamiltonian(G1), hamiltonian(G2)) +end + +# ────────────────────────────────────────────── +# Matrix-vector helpers +# ────────────────────────────────────────────── + +@generated function _slh_matvec(S::SMatrix{N,N}, L::SVector{N}) where {N} + if N == 1 + return :(SVector{1}(_mul(S[1, 1], L[1]))) + end + exprs = [] + for i = 1:N + first_name = Symbol("tmp_$(i)_1") + terms = [:($first_name = _mul(S[$i, 1], L[1]))] + acc = first_name + for j = 2:N + tname = Symbol("tmp_$(i)_$(j)") + push!(terms, :($tname = _add($acc, _mul(S[$i, $j], L[$j])))) + acc = tname end - L_out = [L_t(i) for i = 1:lL1] - - H_out = - t -> begin - L1_vec = [f(t) for f in L1f] - L2_vec = [f(t) for f in L2f] - term1 = sum(adjoint(L2_vec[i]) * (S2*L1_vec)[i] for i = 1:lL1) - term2 = sum(adjoint(L1_vec[i]) * (adjoint(S2)*L2_vec)[i] for i = 1:lL1) - H1f(t) + H2f(t) - 1im / 2 * (term1 - term2) - end - return SLHqo(S_t, L_out, H_out) + push!(exprs, Expr(:block, terms..., acc)) end + return :(SVector($(exprs...))) +end - L_out = L2 + S2 * L1 - term1 = sum(adjoint(L2[i]) * (S2*L1)[i] for i = 1:lL1) - term2 = sum(adjoint(L1[i]) * (adjoint(S2)*L2)[i] for i = 1:lL1) - H_out = H1 + H2 - 1im / 2 * (term1 - term2) - return SLHqo(S_t, L_out, H_out) +@generated function _slh_dot(L1::SVector{N}, L2::SVector{N}) where {N} + if N == 1 + return :(_mul(L1[1], L2[1])) + end + expr = :(_mul(L1[1], L2[1])) + for i = 2:N + expr = :(_add($expr, _mul(L1[$i], L2[$i]))) + end + return expr end -▷(a::SLHqo, b::SLHqo, c::SLHqo...) = ▷(a▷b, c...) + +# ────────────────────────────────────────────── +# Cascade: ▷ +# Val(N*N) for ntuple in S2_adj construction +# ────────────────────────────────────────────── """ - cascade(G::SLH...) + ▷(G1::SLH{N}, G2::SLH{N}) where N -Creates a new SLH triple by cascading the SLH triples from first to last according to -the rule +Cascade two SLH triples: -```math -SLH_1 \\triangleright SLH_2 = ( S_2 S_1, L_2 + S_2 L_1, H_1 + H_2 - \\frac{i}{2} L_2^\\dagger S_2 L_1 - L_1^\\dagger S_2^\\dagger L_2 ) -``` +``G_1 \\triangleright G_2 = (S_2 S_1,\\; L_2 + S_2 L_1,\\; H_1 + H_2 - \\tfrac{i}{2}(L_2^\\dagger S_2 L_1 - L_1^\\dagger S_2^\\dagger L_2))`` -See also [`▷`](@ref). +Unicode `\\triangleright`. See also [`cascade`](@ref). """ -cascade(args...) = ▷(args...) -# ◁(G1::SLH,G2::SLH) = ▷(G2,G1) # unknown unicode character +function ▷(G1::SLH{N}, G2::SLH{N}) where {N} + S1, L1, H1 = scattering(G1), lindblad(G1), hamiltonian(G1) + S2, L2, H2 = scattering(G2), lindblad(G2), hamiltonian(G2) + + S_t = _post.(S2 * S1) + S2L1 = _slh_matvec(S2, L1) + L_t = SVector{N}(ntuple(i -> _post(_add(L2[i], S2L1[i])), Val(N))) + + # Cross terms for Hamiltonian + S2_adj = SMatrix{N,N}(ntuple(Val(N * N)) do k + i, j = divrem(k - 1, N) .+ (1, 1) + _adj(S2[i, j]) + end) + L1_adj = SVector{N}(ntuple(i -> _adj(L1[i]), Val(N))) + L2_adj = SVector{N}(ntuple(i -> _adj(L2[i]), Val(N))) + + cross1 = _slh_dot(L2_adj, S2L1) + S2adj_L2 = _slh_matvec(S2_adj, L2) + cross2 = _slh_dot(L1_adj, S2adj_L2) + + H_t = _post(_add(_add(H1, H2), _mul(-1im / 2, _add(cross1, _mul(-1, cross2))))) + + op_hint = _op_type(G1) + op_hint === nothing && (op_hint = _op_type(G2)) + return _build_slh(S_t, L_t, H_t, op_hint) +end + +▷(a::SLH, b::SLH, c::SLH...) = ▷(a ▷ b, c...) """ - ⊞(G::SLH...) + cascade(G::SLH...) -Creates a new SLH triple by concatenating the SLH triples according to -the rule +Cascade SLH triples from first to last. Alias for [`▷`](@ref). +""" +cascade(args...) = ▷(args...) -```math -SLH_1 \\boxplus SLH_2 = \\left( \\begin{pmatrix} S_1 & 0 \\\\ 0 & S_2 \\end{pmatrix}, \\begin{pmatrix} L_1 \\\\ L_2 \\end{pmatrix}, H_1 + H_2 \\right) -``` +# ────────────────────────────────────────────── +# Concatenate: ⊞ +# ────────────────────────────────────────────── -Unicode `\\boxplus` alias of [`concatenate`](@ref) """ -function ⊞(G1::SLH, G2::SLH) #\boxplus - S1 = get_scattering(G1); - L1 = get_lindblad(G1); - H1 = get_hamiltonian(G1) - S2 = get_scattering(G2); - L2 = get_lindblad(G2); - H2 = get_hamiltonian(G2) - - lS1 = size(S1, 1); - lS2 = size(S2, 1) - - S_t = Matrix{Any}(undef, lS1+lS2, lS1+lS2) # that is pretty ugly (but it works) - S_t .= zeros(lS1+lS2, lS1+lS2) - S_t[1:lS1, 1:lS1] .= S1 - S_t[(1+lS1):(lS1+lS2), (1+lS1):(lS1+lS2)] .= S2 - - L_t = [L1; L2] - H_t = H1 + H2 - - return SLH(S_t, L_t, H_t) -end -⊞(a::SLH, b::SLH, c::SLH...) = ⊞(a⊞b, c...) + ⊞(G1::SLH{N1}, G2::SLH{N2}) -""" - ⊞(G::SLHqo...) +Concatenate (parallel composition) of two SLH triples: + +``G_1 \\boxplus G_2 = \\left(\\begin{pmatrix} S_1 & 0 \\\\ 0 & S_2 \\end{pmatrix},\\; +\\begin{pmatrix} L_1 \\\\ L_2 \\end{pmatrix},\\; H_1 + H_2\\right)`` -Concatenates SLHqo triples, allowing time-dependent `L` or `H`. -If any input `L` or `H` is time-dependent, the result uses time-dependent functions. +Unicode `\\boxplus`. See also [`concatenate`](@ref). """ -function ⊞(G1::SLHqo, G2::SLHqo) - S1 = get_scattering(G1); - L1 = get_lindblad(G1); - H1 = get_hamiltonian(G1) - S2 = get_scattering(G2); - L2 = get_lindblad(G2); - H2 = get_hamiltonian(G2) - - lS1 = size(S1, 1); - lS2 = size(S2, 1) - T_S = promote_type(eltype(S1), eltype(S2)) - S_t = Matrix{T_S}(undef, lS1 + lS2, lS1 + lS2) - fill!(S_t, 0) - S_t[1:lS1, 1:lS1] .= S1 - S_t[(1+lS1):(lS1+lS2), (1+lS1):(lS1+lS2)] .= S2 - - L1f = map(_to_func, L1) - L2f = map(_to_func, L2) - time_dep = - any(_is_time_dep, L1) || - any(_is_time_dep, L2) || - _is_time_dep(H1) || - _is_time_dep(H2) - - if time_dep - L_out = vcat(L1f, L2f) - H_out = t -> _to_func(H1)(t) + _to_func(H2)(t) - return SLHqo(S_t, L_out, H_out) +@generated function ⊞(G1::SLH{N1}, G2::SLH{N2}) where {N1,N2} + N = N1 + N2 + s_exprs = [] + for j = 1:N, i = 1:N # column-major for SMatrix + if i <= N1 && j <= N1 + push!(s_exprs, :(S1[$i, $j])) + elseif i > N1 && j > N1 + push!(s_exprs, :(S2[$(i-N1), $(j-N1)])) + else + push!(s_exprs, 0) + end end - L_out = [L1; L2] - H_out = H1 + H2 - return SLHqo(S_t, L_out, H_out) + quote + S1, L1, H1 = scattering(G1), lindblad(G1), hamiltonian(G1) + S2, L2, H2 = scattering(G2), lindblad(G2), hamiltonian(G2) + S_t = SMatrix{$N,$N}($(s_exprs...)) + L_t = vcat(L1, L2) + H_t = _add(H1, H2) + op_hint = _op_type(G1) + op_hint === nothing && (op_hint = _op_type(G2)) + return _build_slh(S_t, L_t, H_t, op_hint) + end end -⊞(a::SLHqo, b::SLHqo, c::SLHqo...) = ⊞(a⊞b, c...) + +⊞(a::SLH, b::SLH, c::SLH...) = ⊞(a ⊞ b, c...) """ concatenate(G::SLH...) -Creates a new SLH triple by concatenating the SLH triples according to -the rule - -```math -SLH_1 \\boxplus SLH_2 = \\left( \\begin{pmatrix} S_1 & 0 \\\\ 0 & S_2 \\end{pmatrix}, \\begin{pmatrix} L_1 \\\\ L_2 \\end{pmatrix}, H_1 + H_2 \\right) -``` - -See also [`⊞`](@ref). +Concatenate (parallel composition) of SLH triples. Alias for [`⊞`](@ref). """ concatenate(args...) = ⊞(args...) -concatenation(args...) = concatenate(args...) # alias - -# feedback reduction rules -_dropindex(v::AbstractVector, i::Int) = [v[k] for k = 1:length(v) if k != i] - -function _dropindices(M::AbstractMatrix, i::Int, j::Int) - rows = [r for r = 1:size(M, 1) if r != i] - cols = [c for c = 1:size(M, 2) if c != j] - return M[rows, cols] +# ────────────────────────────────────────────── +# Feedback reduction +# dispatch through Val(N-1) so M is a type parameter +# all ntuple calls use Val +# ────────────────────────────────────────────── + +function _drop_row_col(S::SMatrix{N,N}, row::Int, col::Int, ::Val{M}) where {N,M} + SMatrix{M,M}(ntuple(Val(M * M)) do k + i, j = divrem(k - 1, M) .+ (1, 1) + ri = i >= row ? i + 1 : i + cj = j >= col ? j + 1 : j + S[ri, cj] + end) end -function _feedback_scalar(S::AbstractMatrix, L::AbstractVector, x::Int, y::Int) - n = size(S, 1) - @assert size(S, 2) == n - @assert length(L) == n - @assert 1 <= x <= n - @assert 1 <= y <= n - - S_xy = S[x, y] - S_xbar_ybar = _dropindices(S, x, y) - S_xbar_y = _dropindex(S[:, y], x) - S_x_ybar = permutedims(_dropindex(vec(S[x, :]), y)) - S_col_y = vec(S[:, y]) - L_xbar = _dropindex(L, x) - L_x = L[x] - - loop_gain = 1 / (1 - S_xy) - return S_xbar_ybar, S_xbar_y, S_x_ybar, S_col_y, L_xbar, L_x, loop_gain +function _drop_index(L::SVector{N}, idx::Int, ::Val{M}) where {N,M} + SVector{M}(ntuple(i -> L[i >= idx ? i + 1 : i], Val(M))) end -function _feedback_scattering_update(S_xbar_y::AbstractVector, loop_gain, S_x_ybar) - nrows = length(S_xbar_y) - ncols = length(S_x_ybar) - return [S_xbar_y[i] * loop_gain * S_x_ybar[j] for i = 1:nrows, j = 1:ncols] +function _get_col_dropped_row( + S::SMatrix{N,N}, + col::Int, + drop_row::Int, + ::Val{M}, +) where {N,M} + SVector{M}(ntuple(i -> S[i >= drop_row ? i + 1 : i, col], Val(M))) end -function _feedback_lindblad_update(S_xbar_y::AbstractVector, loop_gain, L_x) - return [S_xbar_y[i] * loop_gain * L_x for i = 1:length(S_xbar_y)] +function _get_row_dropped_col( + S::SMatrix{N,N}, + row::Int, + drop_col::Int, + ::Val{M}, +) where {N,M} + SVector{M}(ntuple(j -> S[row, j >= drop_col ? j + 1 : j], Val(M))) end -function _feedback_maps(n::Int, connections) - output_ports = collect(1:n) - input_ports = collect(1:n) - mapped = Tuple{Int,Int}[] +""" + feedback(G::SLH{N}, x::Int, y::Int) where N - for connection in connections - x = first(connection) - y = last(connection) - x_now = findfirst(==(x), output_ports) - y_now = findfirst(==(y), input_ports) - @assert x_now !== nothing "output port $x has already been eliminated" - @assert y_now !== nothing "input port $y has already been eliminated" - push!(mapped, (x_now, y_now)) - deleteat!(output_ports, x_now) - deleteat!(input_ports, y_now) - end +Apply the SLH feedback reduction rule: connect output port `x` to input port `y`. +Returns `SLH{N-1}`. - return mapped +See also [`SLH`](@ref), [`▷`](@ref), [`⊞`](@ref). +""" +function feedback(G::SLH{N}, x::Int, y::Int) where {N} + _feedback_impl(G, x, y, Val(N - 1)) end -""" - feedback(G::SLH, x::Int, y::Int) - feedback(G::SLH, connection::Pair{Int,Int}) - feedback(G::SLH, connections::Pair{Int,Int}...) +function _feedback_impl(G::SLH{N}, x::Int, y::Int, ::Val{M}) where {N,M} + S = scattering(G) + L = lindblad(G) + H = hamiltonian(G) -Apply the SLH feedback reduction rule to connect output port `x` to input port `y`. -Multiple connections can be passed as `x => y` pairs; in that form the pairs are -interpreted using the original port labels and reduced in sequence with the port -bookkeeping handled automatically. + @assert 1 <= x <= N && 1 <= y <= N -See also [`SLH`](@ref), [`▷`](@ref), and [`⊞`](@ref). -""" -function feedback(G::SLH, x::Int, y::Int) - S = get_scattering(G) - L = get_lindblad(G) - H = get_hamiltonian(G) + valM = Val(M) + S_xy = S[x, y] + loop_gain = 1 / (1 - S_xy) - S_xbar_ybar, S_xbar_y, S_x_ybar, S_col_y, L_xbar, L_x, loop_gain = - _feedback_scalar(S, L, x, y) + S_bar = _drop_row_col(S, x, y, valM) + S_col_y_no_x = _get_col_dropped_row(S, y, x, valM) + S_row_x_no_y = _get_row_dropped_col(S, x, y, valM) - S_red = simplify.( - S_xbar_ybar + _feedback_scattering_update(S_xbar_y, loop_gain, vec(S_x_ybar)), - ) - L_red = simplify.(L_xbar + _feedback_lindblad_update(S_xbar_y, loop_gain, L_x)) + S_update = SMatrix{M,M}(ntuple(Val(M * M)) do k + i, j = divrem(k - 1, M) .+ (1, 1) + _post(_mul(_mul(S_col_y_no_x[i], loop_gain), S_row_x_no_y[j])) + end) + S_red = SMatrix{M,M}(ntuple(Val(M * M)) do k + i, j = divrem(k - 1, M) .+ (1, 1) + _post(_add(S_bar[i, j], S_update[i, j])) + end) + + L_bar = _drop_index(L, x, valM) + L_x = L[x] + L_update = + SVector{M}(ntuple(i -> _post(_mul(_mul(S_col_y_no_x[i], loop_gain), L_x)), valM)) + L_red = SVector{M}(ntuple(i -> _post(_add(L_bar[i], L_update[i])), valM)) - L_ct = SQA._adjoint.(L) - term = sum(L_ct[i] * S_col_y[i] for i = 1:length(L)) * loop_gain * L_x - H_red = simplify(H + (term - SQA._adjoint(term)) / (2im)) + S_col_y_full = SVector{N}(ntuple(i -> S[i, y], Val(N))) + L_adj = SVector{N}(ntuple(i -> _adj(L[i]), Val(N))) + term = _mul(_slh_dot(L_adj, S_col_y_full), _mul(loop_gain, L_x)) + H_red = _post(_add(H, _mul(1 / (2im), _add(term, _mul(-1, _adj(term)))))) - return SLH(S_red, L_red, H_red) + return _build_slh(S_red, L_red, H_red, _op_type(G)) end feedback(G::SLH, connection::Pair{Int,Int}) = feedback(G, first(connection), last(connection)) function feedback(G::SLH, connections::Pair{Int,Int}...) - n = size(get_scattering(G), 1) + n = size(scattering(G), 1) G_red = G for (x, y) in _feedback_maps(n, connections) G_red = feedback(G_red, x, y) @@ -431,72 +433,23 @@ function feedback(G::SLH, connections::Pair{Int,Int}...) return G_red end -""" - feedback(G::SLHqo, x::Int, y::Int) - feedback(G::SLHqo, connection::Pair{Int,Int}) - feedback(G::SLHqo, connections::Pair{Int,Int}...) - -Apply the feedback reduction rule to an `SLHqo` triple. Time-dependent `L` or `H` -are preserved in the reduced model. -""" -function feedback(G::SLHqo, x::Int, y::Int) - S = get_scattering(G) - L = get_lindblad(G) - H = get_hamiltonian(G) - - S_xbar_ybar, S_xbar_y, S_x_ybar, S_col_y, L_xbar, _, loop_gain = - _feedback_scalar(S, L, x, y) - - S_red = S_xbar_ybar + _feedback_scattering_update(S_xbar_y, loop_gain, vec(S_x_ybar)) - - Lf = map(_to_func, L) - Hf = _to_func(H) - time_dep = any(_is_time_dep, L) || _is_time_dep(H) - - function reduced_L(i) - return t -> begin - L_vec = [f(t) for f in Lf] - _, S_xbar_y_t, _, _, L_xbar_t, L_x_t, loop_gain_t = - _feedback_scalar(S, L_vec, x, y) - (L_xbar_t+_feedback_lindblad_update(S_xbar_y_t, loop_gain_t, L_x_t))[i] - end - end - - function reduced_H(t) - L_vec = [f(t) for f in Lf] - _, _, _, S_col_y_t, _, L_x_t, loop_gain_t = _feedback_scalar(S, L_vec, x, y) - term = - sum(adjoint(L_vec[i]) * S_col_y_t[i] for i = 1:length(L_vec)) * - loop_gain_t * - L_x_t - Hf(t) + (term - adjoint(term)) / (2im) - end - - if time_dep - return SLHqo(S_red, [reduced_L(i) for i = 1:(length(L)-1)], reduced_H) - end - - L_red = begin - L_vec = copy(L) - _, S_xbar_y_t, _, _, L_xbar_t, L_x_t, loop_gain_t = - _feedback_scalar(S, L_vec, x, y) - L_xbar_t + _feedback_lindblad_update(S_xbar_y_t, loop_gain_t, L_x_t) - end - term = sum(adjoint(L[i]) * S_col_y[i] for i = 1:length(L)) * loop_gain * L[x] - H_red = H + (term - adjoint(term)) / (2im) - return SLHqo(S_red, L_red, H_red) -end - -feedback(G::SLHqo, connection::Pair{Int,Int}) = - feedback(G, first(connection), last(connection)) +function _feedback_maps(n::Int, connections) + output_ports = collect(1:n) + input_ports = collect(1:n) + mapped = Tuple{Int,Int}[] -function feedback(G::SLHqo, connections::Pair{Int,Int}...) - n = size(get_scattering(G), 1) - G_red = G - for (x, y) in _feedback_maps(n, connections) - G_red = feedback(G_red, x, y) + for connection in connections + x = first(connection) + y = last(connection) + x_now = findfirst(==(x), output_ports) + y_now = findfirst(==(y), input_ports) + @assert x_now !== nothing "output port $x has already been eliminated" + @assert y_now !== nothing "input port $y has already been eliminated" + push!(mapped, (x_now, y_now)) + deleteat!(output_ports, x_now) + deleteat!(input_ports, y_now) end - return G_red + return mapped end Base.length(h::SecondQuantizedAlgebra.ConcreteHilbertSpace) = 1 diff --git a/src/correlations.jl b/src/correlations.jl index 8f104c1..177dbc5 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -1,61 +1,54 @@ """ - two_time_corr_matrix(T, ρt, f, Ls; kwargs...) - two_time_corr_matrix(T, ρt, H, J, Ls; kwargs...) + correlation_matrix(T, ρt, f, Ls; kwargs...) + correlation_matrix(T, ρt, H, J, Ls; kwargs...) -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` for the operator `Ls`. -The first method supports time-dependent generators with the function `f`; the second is for time-independent `H` and `J`. +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. """ -function two_time_corr_matrix(T::Vector, ρt::Vector, f::Function, Ls::Function; kwargs...) - l_T = length(T) - @assert l_T == length(ρt) - Ls_ls = Ls.(T) - Ls_ls_dag = dagger.(Ls_ls) - ρ0_ = [Ls_ls[i] * ρt[i] for i = 1:l_T] - - g1_m = zeros(ComplexF64, l_T, l_T) - for it = 1:(l_T-1) - τ_, ρ_bar_τ = timeevolution.master_dynamic(T[it:end], ρ0_[it], f; kwargs...) - - g1 = [expect(Ls_ls_dag[it+i-1], ρ_bar_τ[i]) for i = 1:length(τ_)] - g1_m[it, it:end] = g1 - g1_m[it:end, it] = adjoint.(g1) +function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls::Function; kwargs...) + Ls_vec = Ls.(T) + Ls_dag_vec = dagger.(Ls_vec) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 + timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) end - return g1_m end -# time-dependent problem but constant Ls -function two_time_corr_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) - l_T = length(T) - @assert l_T == length(ρt) +function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) Ls_dag = dagger(Ls) - ρ0_ = [Ls * ρt[i] for i = 1:l_T] - - g1_m = zeros(ComplexF64, l_T, l_T) - for it = 1:(l_T-1) - τ_, ρ_bar_τ = timeevolution.master_dynamic(T[it:end], ρ0_[it], f; kwargs...) + l_T = length(T) + Ls_vec = fill(Ls, l_T) + Ls_dag_vec = fill(Ls_dag, l_T) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 + timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) + end +end - g1 = [expect(Ls_dag, ρ_bar_τ[i]) for i = 1:length(τ_)] - g1_m[it, it:end] = g1 - g1_m[it:end, it] = adjoint.(g1) +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) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 + timeevolution.master(T_slice, ρ0, H, J; kwargs...) end - return g1_m end -# two_time_corr_matrix for time-independent problems -function two_time_corr_matrix(T::Vector, ρt::Vector, H, J::Vector, Ls; kwargs...) +function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_dag_vec) l_T = length(T) @assert l_T == length(ρt) - Ls_dag = dagger(Ls) - ρ0_ = [Ls * ρt[i] for i = 1:l_T] g1_m = zeros(ComplexF64, l_T, l_T) - for it = 1:(l_T-1) - τ_, ρ_bar_τ = timeevolution.master(T[it:end], ρ0_[it], H, J; kwargs...) - - g1 = [expect(Ls_dag, ρ_bar_τ[i]) for i = 1:length(τ_)] - g1_m[it, it:end] = g1 - g1_m[it:end, it] = adjoint.(g1) + # Each iteration solves an independent master equation — parallelise + Threads.@threads for it = 1:(l_T-1) + ρ0_it = Ls_vec[it] * ρt[it] + τ_, ρ_bar_τ = solve_fn(@view(T[it:end]), ρ0_it) + + @inbounds for i in eachindex(ρ_bar_τ) + val = expect(Ls_dag_vec[it+i-1], ρ_bar_τ[i]) + g1_m[it, it+i-1] = val + g1_m[it+i-1, it] = conj(val) + end end return g1_m end diff --git a/src/interaction_picture.jl b/src/interaction_picture.jl index e22d806..2b20a7c 100644 --- a/src/interaction_picture.jl +++ b/src/interaction_picture.jl @@ -2,32 +2,74 @@ ### interaction picture ### ########################### +using StaticArrays: SMatrix + +_as_time_function(x::Number) = _ -> x +_as_time_function(x) = x # anything callable passes through + +""" + coupling_matrix(gs::Tuple) + +Build the antisymmetric coupling coefficient matrix `A(t)` from a tuple of +coupling functions/constants `gs = (g_1, ..., g_N)`. Returns a closure `t -> A(t)`. + +```math +A_{ij}(t) = \\frac{1}{2} \\begin{cases} +0 & i = j \\\\ +g_i(t)\\, g_j^*(t) & i < j \\\\ +-g_i^*(t)\\, g_j(t) & i > j +\\end{cases} +``` + +All couplings may be time-dependent or constant. +""" +function coupling_matrix(gs::NTuple{N}) where {N} + gfs = map(_as_time_function, gs) + function A(t) + SMatrix{N,N,ComplexF64}(ntuple(Val(N * N)) do k + i, j = divrem(k - 1, N) .+ (1, 1) + if i == j + zero(ComplexF64) + elseif j < i + val = 0.5 * gfs[j](t) * conj(gfs[i](t)) + val + else # j > i + val = 0.5 * gfs[i](t) * conj(gfs[j](t)) + -conj(val) + end + end) + end + return A +end + +coupling_matrix(g1, g2, gs...) = coupling_matrix((g1, g2, gs...)) + """ - interaction_picture_M(A, T; alg = Tsit5(), kwargs...) + solve_mode_evolution(A::Function, T; alg=Tsit5(), kwargs...) -Solve the Heisenberg coefficient-matrix equation for an interaction picture, -`dM/dt = A(t) * M(t)` with `M(0) = I`, and return the ODE solution. +Solve the interaction-picture coefficient-matrix ODE `dM/dt = A(t) M(t)` with `M(0) = I`. +Returns the ODE solution directly (callable as `sol(t)`). -The function `A(t)` must return an `N × N` complex matrix describing the linear -mode coupling in the interaction picture. The returned solution supports -interpolation `M(t)` for continuous `t`. +All kwargs are passed on to the ODE solver. """ -function interaction_picture_M(A::Function, T; alg = Tsit5(), kwargs...) +function solve_mode_evolution(A::Function, T; alg = Tsit5(), kwargs...) T0 = T[1] Tend = T[end] n = size(A(T0), 1) M0 = Matrix{ComplexF64}(I, n, n) - f_M(u, p, t) = A(t) * u - prob_M = ODEProblem(f_M, M0, (T0, Tend)) - sol = solve(prob_M, alg; saveat = T, kwargs...) - return t -> sol(t) + function f_M!(du, u, p, t) + mul!(du, A(t), u) + end + prob = ODEProblem(f_M!, M0, (T0, Tend)) + sol = solve(prob, alg; saveat = T, kwargs...) + return sol end """ - interaction_picture_M_2modes_equal(u, T) + solve_mode_evolution_symmetric(u, T) Analytic interaction-picture coefficient matrix for two modes when `u(t) = v(t)`. -The matrix is +Returns a callable `t -> M(t)` where ```math M(t) = \\begin{bmatrix} @@ -42,7 +84,7 @@ where \\sin^2 \\theta(t) = \\int_0^t |u(t')|^2\\,dt'. ``` """ -function interaction_picture_M_2modes_equal(u, T) +function solve_mode_evolution_symmetric(u, T) u_vals = u isa Function ? u.(T) : u sin2θ = cumul_integrate(T, abs2.(u_vals)) sin2θ = clamp.(real.(sin2θ), 0.0, 1.0) @@ -55,94 +97,5 @@ function interaction_picture_M_2modes_equal(u, T) M21 = LinearInterpolation(sθ, T; extrapolation = ExtrapolationType.Extension) M22 = LinearInterpolation(cθ, T; extrapolation = ExtrapolationType.Extension) - return t -> [M11(t) M12(t); M21(t) M22(t)] -end - -_as_time_function(x) = x isa Function ? x : (_ -> x) - -""" - interaction_picture_A_2modes(g1, g2) - -Coefficient matrix `A(t)` for two virtual modes `(1, 2)`, - -```math -A(t) = \\frac{1}{2} -\\begin{bmatrix} -0 & g_1(t) g_2^*(t) \\\\ --g_1^*(t) g_2(t) & 0 -\\end{bmatrix} -``` - -All couplings may be time-dependent or constant. -""" -function interaction_picture_A_2modes(g1, g2) - g1f = _as_time_function(g1) - g2f = _as_time_function(g2) - A(t) = 0.5 * [ - 0 g1f(t) * conj(g2f(t)); - -conj(g1f(t)) * g2f(t) 0 - ] - return A -end -interaction_picture_A_2modes(g_ls) = interaction_picture_A_2modes(g_ls...) - -""" - interaction_picture_A_3modes(g1, g2, g3) - -Coefficient matrix `A(t)` for three interacting modes ordered as `(1, 2, 3)` - -```math -A(t) = \\frac{1}{2} -\\begin{bmatrix} -0 & g_1(t) g_2^*(t) & g_1(t) g_3^*(t) \\\\ --g_1^*(t) g_2(t) & 0 & g_2(t) g_3^*(t) \\\\ --g_1^*(t) g_3(t) & -g_2^*(t) g_3(t) & 0 -\\end{bmatrix} -``` - -All couplings may be time-dependent or constant. -""" -function interaction_picture_A_3modes(g1, g2, g3) - g1f = _as_time_function(g1) - g2f = _as_time_function(g2) - g3f = _as_time_function(g3) - A(t) = - 0.5 * [ - 0 conj(g2f(t)) * g1f(t) g1f(t) * conj(g3f(t)); - -g2f(t) * conj(g1f(t)) 0 g2f(t) * conj(g3f(t)); - -conj(g1f(t)) * g3f(t) -conj(g2f(t)) * g3f(t) 0 - ] - return A -end - -""" - interaction_picture_A_4modes(g1, g2, g3, g4) - -Coefficient matrix `A(t)` for four interacting modes ordered as `(1, 2, 3, 4)` - -```math -A(t) = \\frac{1}{2} -\\begin{bmatrix} -0 & g_1(t) g_2^*(t) & g_1(t) g_3^*(t) & g_1(t) g_4^*(t) \\\\ --g_1^*(t) g_2(t) & 0 & g_2(t) g_3^*(t) & g_2(t) g_4^*(t) \\\\ --g_1^*(t) g_3(t) & -g_2^*(t) g_3(t) & 0 & g_3(t) g_4^*(t) \\\\ --g_1^*(t) g_4(t) & -g_2^*(t) g_4(t) & -g_3^*(t) g_4(t) & 0 -\\end{bmatrix} -``` - -All couplings may be time-dependent or constant. -""" -function interaction_picture_A_4modes(g1, g2, g3, g4) - g1f = _as_time_function(g1) - g2f = _as_time_function(g2) - g3f = _as_time_function(g3) - g4f = _as_time_function(g4) - A(t) = - 0.5 * [ - 0 g1f(t) * conj(g2f(t)) g1f(t) * conj(g3f(t)) g1f(t) * conj(g4f(t)); - -conj(g1f(t)) * g2f(t) 0 g2f(t) * conj(g3f(t)) g2f(t) * conj(g4f(t)); - -conj(g1f(t)) * g3f(t) -conj(g2f(t)) * g3f(t) 0 g3f(t) * conj(g4f(t)); - -conj(g1f(t)) * g4f(t) -conj(g2f(t)) * g4f(t) -conj(g3f(t)) * g4f(t) 0 - ] - return A + return t -> SMatrix{2,2}(M11(t), M21(t), M12(t), M22(t)) # column-major end diff --git a/src/pulses.jl b/src/pulses.jl index 6e3fbae..0b9b616 100644 --- a/src/pulses.jl +++ b/src/pulses.jl @@ -2,128 +2,165 @@ ### functions for virtual cavities ### ###################################### -_tol_div = 1e-10 #12 -_extrapolate = ExtrapolationType.Extension -_ϵu = 1e-10 -_ϵv = 1e-10 +const _tol_div = 1e-10 +const _extrapolate = ExtrapolationType.Extension +const _ϵ = 1e-10 _mode_interp(mode::AbstractVector, T::AbstractVector) = LinearInterpolation(mode, T; extrapolation = _extrapolate) _mode_interp(mode, T::AbstractVector) = mode +# ────────────────────────────────────────────── +# Gaussian pulse type +# ────────────────────────────────────────────── + """ - u_to_gu(u, T) + Gaussian(τ, σ; δ=0) -Compute the virtual-cavity coupling ``g_u(t)`` from an input mode `u(t)` sampled on `T` with a linear interpolation. -Returns a callable `t -> g_u(t)`. +Gaussian pulse shape descriptor with center time `τ`, width `σ`, and detuning `δ`. +Use with `coupling_input` and `coupling_output` for analytical coupling formulas. """ -function u_to_gu(u::Vector, T::Vector) - l_T = length(T) - ∫u2_t = cumul_integrate(T, abs2.(u)) .+ 0im - gu_t = zeros(ComplexF64, l_T) - for i = 2:l_T - if abs(sqrt(1 - ∫u2_t[i])) > _tol_div - gu_t[i] = u[i]' / sqrt(abs(1 - ∫u2_t[i]) + _ϵu) - end #else 0 +struct Gaussian{T} + τ::T + σ::T + δ::T +end +Gaussian(τ, σ; δ = zero(τ)) = Gaussian(promote(τ, σ, δ)...) + +function _gaussian_mode(g::Gaussian) + τ, σ, δ = g.τ, g.σ, g.δ + mode = if δ == 0 + t -> 1 / (√(σ) * π^(1 / 4)) * exp(-0.5 * (t - τ)^2 / σ^2) + else + t -> 1 / (√(σ) * π^(1 / 4)) * exp(-0.5 * (t - τ)^2 / σ^2) * exp(1im * δ * t) end - gu_int = LinearInterpolation(gu_t, T; extrapolation = _extrapolate) - return t -> gu_int(t) + ∫mode2 = t -> 0.5 * (erf((t - τ) / σ) + erf(τ / σ)) + return mode, ∫mode2 end -u_to_gu(u::Function, T::Vector) = u_to_gu(u.(T), T) -u_to_gu(u::LinearInterpolation, T::Vector) = u_to_gu(u.(T), T) -""" - v_to_gv(v, T) +# ────────────────────────────────────────────── +# Shared coupling core +# ────────────────────────────────────────────── -Compute the virtual-cavity coupling ``g_v(t)`` from an output mode `v(t)` sampled on `T` with a linear interpolation. -Returns a callable `t -> g_v(t)`. -""" -function v_to_gv(v::Vector, T::Vector) +function _compute_coupling(mode::Vector, T::Vector, denom_fn) l_T = length(T) - ∫v2_t = cumul_integrate(T, abs2.(v)) .+ 0im - gv_t = zeros(ComplexF64, l_T) - for i = 2:l_T - if abs(sqrt(∫v2_t[i])) > _tol_div - gv_t[i] = -v[i]' / sqrt(∫v2_t[i] + _ϵv) - end # else 0 + 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) + end end - gv_int = LinearInterpolation(gv_t, T; extrapolation = _extrapolate) - return t -> gv_int(t) + return LinearInterpolation(g, T; extrapolation = _extrapolate) end -v_to_gv(v::Function, T::Vector) = v_to_gv(v.(T), T) -v_to_gv(v::LinearInterpolation, T::Vector) = v_to_gv(v.(T), T) -""" - u_to_gu_Gauss(τ, σ; δ=0) +# ────────────────────────────────────────────── +# coupling_input / coupling_output +# ────────────────────────────────────────────── -Compute ``g_u(t)`` for a Gaussian input mode `u(t)` with delay `τ`, width `σ`, and detuning `δ`: +""" + coupling_input(u, T) -`` u(t) = 1/\\sqrt{ \\sigma \\sqrt{\\pi} } e^{-(t - \\tau)^2 / 2 \\sigma^2 )} e^{i \\delta 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)`). +""" +coupling_input(u::Vector, T::Vector) = _compute_coupling(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) +""" + coupling_input(g::Gaussian) -Returns a callable `t -> g_u(t)`. +Analytical input coupling ``g_u(t)`` for a Gaussian pulse. """ -function u_to_gu_Gauss(τ, σ; δ = 0) # slower than u_to_gu - if δ==0 # type stable - u = t -> 1/(√(σ)*π^(1/4)) * exp(-0.5*(t - τ)^2/σ^2) - else - u = t -> 1/(√(σ)*π^(1/4)) * exp(-0.5*(t - τ)^2/σ^2) * exp(1im*δ*t) - end - ∫u_2_t(t_) = 0.5 * (erf((t_ - τ) / σ) + erf(τ / σ)) - f = t_ -> u(t_)' / √(abs(1 - ∫u_2_t(t_)) + _ϵu) - return f +function coupling_input(g::Gaussian) + mode, ∫m2 = _gaussian_mode(g) + return t -> mode(t)' / √(abs(1 - ∫m2(t)) + _ϵ) end + """ - v_to_gv_Gauss(τ, σ; δ=0) + coupling_output(v, T) -Compute ``g_v(t)`` for a Gaussian output mode `v(t)` with delay `τ`, width `σ`, and detuning `δ`: +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)`). +""" +coupling_output(v::Vector, T::Vector) = _compute_coupling(-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) -`` v(t) = 1/\\sqrt{ \\sigma \\sqrt{\\pi} } e^{-(t - \\tau)^2 / 2 \\sigma^2 )} e^{i \\delta t} `` +""" + coupling_output(g::Gaussian) -Returns a callable `t -> g_v(t)`. +Analytical output coupling ``g_v(t)`` for a Gaussian pulse. """ -function v_to_gv_Gauss(τ, σ; δ = 0) # slower than v_to_gv - if δ==0 # type stable - v = t -> 1/(√(σ)*π^(1/4)) * exp(-0.5*(t - τ)^2/σ^2) - else - v = t -> 1/(√(σ)*π^(1/4)) * exp(-0.5*(t - τ)^2/σ^2) * exp(1im*δ*t) - end - ∫v_2_t(t_) = 0.5 * (erf((t_ - τ) / σ) + erf(τ / σ)) - f = t_ -> -v(t_)' / √(∫v_2_t(t_) + _ϵv) - return f +function coupling_output(g::Gaussian) + mode, ∫m2 = _gaussian_mode(g) + return t -> -mode(t)' / √(∫m2(t) + _ϵ) end +# ────────────────────────────────────────────── +# effective_output_mode (was v_eff) +# ────────────────────────────────────────────── + """ - v_eff(v_fcts, gv_fcts, T, i) - v_eff(v_fcts, T, i) + effective_output_mode(v_fcts, gv_fcts, T, i) + effective_output_mode(v_fcts, T, i) -Compute the effective output mode ``v_i^{\\mathrm{eff}}(t)`` for a system with multiple output modes, -due to the pulse distortion from the preceding output cavities. -The returned mode function is the relevant one for ``g_{v_i}``. -The output modes in the list `v_fcts` need to be sorted starting with the first -output cavity after the system. +Compute the effective output mode ``v_i^{\\mathrm{eff}}(t)`` for a system with multiple +output modes, due to the pulse distortion from the preceding output cavities. +The output modes in `v_fcts` must be sorted starting with the first output cavity after the system. -All kwargs are passed on to the ODE solver. +All kwargs are passed on to the ODE solver. """ -function v_eff(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...) +function effective_output_mode(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...) @assert i > 1 - function multiple_outputs_α(dα, α, p, t) # only for i>1 - for j = 1:(i-1) - dα[j] = - -gv_fcts[j](t) * - (v_fcts[i](t) + sum((gv_fcts[k](t))' * α[k] for k = 1:(j-1); init = 0.0)) - - 0.5 * abs2(gv_fcts[j](t)) * α[j] + n = i - 1 + # Capture as tuples for concrete closure types + gv_t = 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) + end + vi_t = v_i(t) + @inbounds for j = 1:n + coupling_sum = zero(ComplexF64) + for k = 1:(j-1) + coupling_sum += gv_buf[k]' * α[k] + end + dα[j] = -gv_buf[j] * (vi_t + coupling_sum) - 0.5 * abs2(gv_buf[j]) * α[j] end end - u0 = zeros(ComplexF64, i-1) + u0 = zeros(ComplexF64, n) + p = Vector{ComplexF64}(undef, n) tspan = (T[1], T[end]) - prob = ODEProblem(multiple_outputs_α, u0, tspan) + prob = ODEProblem(multiple_outputs_α!, u0, tspan, p) sol_α = solve(prob, alg; kwargs...) - v_i_im1(t) = v_fcts[i](t) + sum((gv_fcts[k](t))' * sol_α(t)[k] for k = 1:(i-1)) - return v_i_im1 + function v_i_eff(t) + α_t = sol_α(t) + result = v_i(t) + @inbounds for k = 1:n + result += gv_all[k](t)' * α_t[k] + end + return result + end + return v_i_eff end -v_eff(v_fcts, T, i) = v_eff(v_fcts, [v_to_gv(v_, T) for v_ in v_fcts], T, i) -function v_eff( + +effective_output_mode(v_fcts, T, i; kwargs...) = effective_output_mode( + v_fcts, + [coupling_output(v_, T) for v_ in v_fcts], + T, + i; + kwargs..., +) + +function effective_output_mode( v_data::AbstractVector{<:AbstractVector}, gv_data::AbstractVector{<:AbstractVector}, T::AbstractVector, @@ -133,48 +170,81 @@ function v_eff( ) v_fcts = [_mode_interp(v_, T) for v_ in v_data] gv_fcts = [_mode_interp(gv_, T) for gv_ in gv_data] - return v_eff(v_fcts, gv_fcts, T, i; alg, kwargs...) + return effective_output_mode(v_fcts, gv_fcts, T, i; alg, kwargs...) end -v_eff( + +effective_output_mode( v_data::AbstractVector{<:AbstractVector}, T::AbstractVector, i; alg = Tsit5(), kwargs..., -) = v_eff(v_data, [v_to_gv(v_, T).(T) for v_ in v_data], T, i; alg, kwargs...) +) = effective_output_mode( + v_data, + [coupling_output(v_, T).(T) for v_ in v_data], + T, + i; + alg, + kwargs..., +) + +# ────────────────────────────────────────────── +# effective_input_mode (was u_eff) +# ────────────────────────────────────────────── """ - u_eff(u_fcts, gu_fcts, T, i) - u_eff(u_fcts, T, i) + effective_input_mode(u_fcts, gu_fcts, T, i) + effective_input_mode(u_fcts, T, i) -Compute the effective input mode ``u_i^{\\mathrm{eff}}(t)`` for a system with multiple input modes, -due to the pulse distortion from the subsequent input cavities. -The returned mode function is the relevant one for ``g_{u_i}``. -The input modes in the list `u_fcts` need to be sorted starting with the first -input cavity before the system. +Compute the effective input mode ``u_i^{\\mathrm{eff}}(t)`` for a system with multiple +input modes, due to the pulse distortion from the subsequent input cavities. +The input modes in `u_fcts` must be sorted starting with the first input cavity before the system. -All kwargs are passed on to the ODE solver. +All kwargs are passed on to the ODE solver. """ -function u_eff(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...) +function effective_input_mode(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...) @assert i > 1 - function multiple_inputs_α(dα, α, p, t) # only for i>1 - for j = 1:(i-1) - dα[j] = - -gu_fcts[j](t) * - (u_fcts[i](t) - sum((gu_fcts[i](t))' * α[k] for k = 1:(j-1); init = 0.0)) + - 0.5 * abs2(gu_fcts[j](t)) * α[j] + n = i - 1 + # Capture as tuples for concrete closure types + gu_t = 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) + end + ui_t = u_i(t) + gui_t = gu_i(t) + @inbounds for j = 1:n + coupling_sum = zero(ComplexF64) + for k = 1:(j-1) + coupling_sum += gui_t' * α[k] + end + dα[j] = -gu_buf[j] * (ui_t - coupling_sum) + 0.5 * abs2(gu_buf[j]) * α[j] end - # 2 Typos in PRA2020-Kiilerich Eq.(A15): last term "-" → "+" and |g_ui|² → |g_uj|² end - u0 = zeros(ComplexF64, i-1) + u0 = zeros(ComplexF64, n) + p = Vector{ComplexF64}(undef, n) tspan = (T[1], T[end]) - prob = ODEProblem(multiple_inputs_α, u0, tspan) + prob = ODEProblem(multiple_inputs_α!, u0, tspan, p) sol_α = solve(prob, alg; kwargs...) - u_i_im1(t) = u_fcts[i](t) - sum((gu_fcts[k](t))' * sol_α(t)[k] for k = 1:(i-1)) - return u_i_im1 + function u_i_eff(t) + α_t = sol_α(t) + result = u_i(t) + @inbounds for k = 1:n + result -= gu_all[k](t)' * α_t[k] + end + return result + end + return u_i_eff end -u_eff(u_fcts, T, i) = u_eff(u_fcts, [u_to_gu(u_, T) for u_ in u_fcts], T, i) -function u_eff( + +effective_input_mode(u_fcts, T, i; kwargs...) = + effective_input_mode(u_fcts, [coupling_input(u_, T) for u_ in u_fcts], T, i; kwargs...) + +function effective_input_mode( u_data::AbstractVector{<:AbstractVector}, gu_data::AbstractVector{<:AbstractVector}, T::AbstractVector, @@ -184,62 +254,64 @@ function u_eff( ) u_fcts = [_mode_interp(u_, T) for u_ in u_data] gu_fcts = [_mode_interp(gu_, T) for gu_ in gu_data] - return u_eff(u_fcts, gu_fcts, T, i; alg, kwargs...) + return effective_input_mode(u_fcts, gu_fcts, T, i; alg, kwargs...) end -u_eff( + +effective_input_mode( u_data::AbstractVector{<:AbstractVector}, T::AbstractVector, i; alg = Tsit5(), kwargs..., -) = u_eff(u_data, [u_to_gu(u_, T).(T) for u_ in u_data], T, i; alg, kwargs...) +) = effective_input_mode( + u_data, + [coupling_input(u_, T).(T) for u_ in u_data], + T, + i; + alg, + kwargs..., +) -""" - uv_to_gout(u, v, T) +# ────────────────────────────────────────────── +# coupling_delay_out / coupling_delay_in +# ────────────────────────────────────────────── -Compute the out-coupling strength ``\\tilde g_{\\mathrm{out},u,v}(t)`` for a delay cavity -that absorbs an incoming pulse `v(t)` while simultaneously emitting the desired pulse `u(t)`. -Returns callable `t -> g_out(t)` based on samples on `T` with a linear interpolation. -""" -function uv_to_gout(u::Vector, v::Vector, T::Vector) +function _compute_coupling_delay(num_mode::Vector, u::Vector, v::Vector, T::Vector) l_T = length(T) - ∫u2_t = cumul_integrate(T, abs2.(u)) .+ 0im - ∫v2_t = cumul_integrate(T, abs2.(v)) .+ 0im - gout_t = zeros(ComplexF64, l_T) - for i = 2:l_T - denom = abs(∫v2_t[i] - ∫u2_t[i]) - if abs(sqrt(denom)) > _tol_div - gout_t[i] = u[i]' / sqrt(denom + _ϵu) + 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 + _ϵ) end end - gout_int = LinearInterpolation(gout_t, T; extrapolation = _extrapolate) - return t -> gout_int(t) + return LinearInterpolation(g, T; extrapolation = _extrapolate) end -uv_to_gout(u::Function, v::Function, T::Vector) = uv_to_gout(u.(T), v.(T), T) -uv_to_gout(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = - uv_to_gout(u.(T), v.(T), T) """ - uv_to_gin(u, v, T) + coupling_delay_out(u, v, T) -Compute the in-coupling strength ``\\tilde g_{\\mathrm{in},v,u}(t)`` for a delay cavity -that absorbs an incoming pulse `v(t)` while simultaneously emitting the desired pulse `u(t)`. -Returns callable `t -> g_in(t)` based on samples on `T` with a linear interpolation. +Compute the out-coupling strength for a delay cavity. +Returns `LinearInterpolation` directly. """ -function uv_to_gin(u::Vector, v::Vector, T::Vector) - l_T = length(T) - ∫u2_t = cumul_integrate(T, abs2.(u)) .+ 0im - ∫v2_t = cumul_integrate(T, abs2.(v)) .+ 0im - gin_t = zeros(ComplexF64, l_T) - for i = 2:l_T - denom = abs(∫v2_t[i] - ∫u2_t[i]) - if abs(sqrt(denom)) > _tol_div - gin_t[i] = -v[i]' / sqrt(denom + _ϵv) - end - end - gin_int = LinearInterpolation(gin_t, T; extrapolation = _extrapolate) - return t -> gin_int(t) -end -uv_to_gin(u::Function, v::Function, T::Vector) = uv_to_gin(u.(T), v.(T), T) -uv_to_gin(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = - uv_to_gin(u.(T), v.(T), T) +coupling_delay_out(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(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) = + coupling_delay_out(u.(T), v.(T), T) + +""" + coupling_delay_in(u, v, T) + +Compute the in-coupling strength for a delay cavity. +Returns `LinearInterpolation` directly. +""" +coupling_delay_in(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(-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/src/translate.jl b/src/translate.jl index 943fb9e..b9f0334 100644 --- a/src/translate.jl +++ b/src/translate.jl @@ -11,34 +11,27 @@ _translate_numeric_raw(op, b; level_map = nothing, operators = Dict(), op_type = op_type(to_numeric(op, b, operators; level_map = level_map)) _translate_one(b, operators, op_type) = - isempty(operators) ? op_type(one(b)) : - op_type(one(basis(collect(values(operators))[1]))) + isempty(operators) ? op_type(one(b)) : op_type(one(basis(first(values(operators))))) +# simplified, no applicable(), values are Number → wrap, anything else → pass through function _normalize_time_parameter(time_parameter) if isempty(time_parameter) return time_parameter end - out = Dict() + KT = keytype(time_parameter) + out = Dict{KT,Any}() for (k, v) in time_parameter - if isa(v, Function) - out[k] = v - elseif isa(v, Number) - out[k] = t -> v - else - error( - "time_parameter values must be callable or numeric, got $(typeof(v)) for key $(k)", - ) - end + out[k] = v isa Number ? (t -> v) : v end return out end +# use Tuple + map instead of generator splat to avoid per-call allocation function _translate_prefactor(arg_c, time_parameter) keys_ = collect(keys(time_parameter)) - values_ = collect(values(time_parameter)) - + values_tuple = Tuple(values(time_parameter)) pref_build = build_function(arg_c, keys_...; expression = Val(false)) - pref_f = t -> pref_build((v(t) for v in values_)...) + pref_f = t -> pref_build(map(v -> v(t), values_tuple)...) return true, pref_f end function _translate_prefactor(arg_c::Number, time_parameter) @@ -49,17 +42,17 @@ end translate_qo(op, b::QuantumOpticsBase.Basis; parameter=Dict(), time_parameter=Dict(), level_map=nothing, operators=Dict(), op_type=sparse) -Translate a symbolic operator `op` into a numeric QuantumOptics.jl operator with the corresponding basis `b`. -The dictionary `parameter` substitutes symbolic parameters with numbers. Time-dependent functions can be provide -with the dictionary `time_parameter`. -If `time_parameter` is non-empty, the result is a time-dependent function `t -> op(t)`. -The kwarg `level_map=nothing` is used to provide the names of levels for `transition` operators. -The operator type which should be returned can be set with the kwarg `op_type=sparse` and -a list of user-defined operators (e.g. on a different basis than `b`) can be provided with the -dictionary `operators=Dict()`. These operators will then be used to replace the symbolic expressions. +Translate a symbolic operator `op` into a numeric QuantumOptics.jl operator with the corresponding basis `b`. +The dictionary `parameter` substitutes symbolic parameters with numbers. Time-dependent functions can be provide +with the dictionary `time_parameter`. +If `time_parameter` is non-empty, the result is a time-dependent function `t -> op(t)`. +The kwarg `level_map=nothing` is used to provide the names of levels for `transition` operators. +The operator type which should be returned can be set with the kwarg `op_type=sparse` and +a list of user-defined operators (e.g. on a different basis than `b`) can be provided with the +dictionary `operators=Dict()`. These operators will then be used to replace the symbolic expressions. """ function translate_qo( - op::SQA.QMul, + op, b::QuantumOpticsBase.Basis; parameter = Dict(), time_parameter = Dict(), @@ -67,40 +60,45 @@ function translate_qo( operators = Dict(), op_type = sparse, ) + tp = _normalize_time_parameter(time_parameter) + return _translate_qo( + op, + b; + parameter, + time_parameter = tp, + level_map, + operators, + op_type, + ) +end + +# ── Internal dispatch: time_parameter already normalized ── - time_parameter = _normalize_time_parameter(time_parameter) +function _translate_qo( + op::SQA.QMul, + b::QuantumOpticsBase.Basis; + parameter = Dict(), + time_parameter = Dict(), + level_map = nothing, + operators = Dict(), + op_type = sparse, +) if isempty(time_parameter) - return _translate_numeric( - op, - b; - parameter = parameter, - level_map = level_map, - operators = operators, - op_type = op_type, - ) # TODO: test + return _translate_numeric(op, b; parameter, level_map, operators, op_type) end - op_ = substitute(op, parameter) # First substitute all symbolic parameters with the numerical values -> only time-dependent parameters left - - if isa(op_, QSym) # special case if parameter is 1.0: after substitute p*op -> 1.0*op -> op (has no fields arg_c, args_nc) - output_op_QMul_QO_ = _translate_numeric_raw( - op_, - b; - level_map = level_map, - operators = operators, - op_type = op_type, - ) # this is faster! - output_op_QMul_QO = t -> output_op_QMul_QO_ - return output_op_QMul_QO - elseif iszero(op_) # special case if parameter is 0.0 + op_ = substitute(op, parameter) + + if isa(op_, QSym) + output_op = _translate_numeric_raw(op_, b; level_map, operators, op_type) + return t -> output_op + elseif iszero(op_) return op_type(0 * one(b)) else arg_c = op_.arg_c args_nc = op_.args_nc - prod_args_nc = op_type( - prod((to_numeric(arg, b, operators; level_map = level_map)) for arg in args_nc), - ) - + prod_args_nc = + op_type(prod((to_numeric(arg, b, operators; level_map)) for arg in args_nc)) is_func, pref = _translate_prefactor(arg_c, time_parameter) if is_func return t -> pref(t) * prod_args_nc @@ -109,8 +107,8 @@ function translate_qo( return t -> prod_c_nc_num end end -# -function translate_qo( + +function _translate_qo( op::SQA.QAdd, b::QuantumOpticsBase.Basis; parameter = Dict(), @@ -119,40 +117,53 @@ function translate_qo( operators = Dict(), op_type = sparse, ) - - time_parameter = _normalize_time_parameter(time_parameter) - op = expand(op) # should ensure that only multiplications are in arguments - ### only static operators - returns a time-independent Hamiltonian and Lindblad + op = expand(op) if isempty(time_parameter) - return _translate_numeric( - op, - b; - parameter = parameter, - level_map = level_map, - operators = operators, - op_type = op_type, - ) # TODO: test + return _translate_numeric(op, b; parameter, level_map, operators, op_type) end - ### static and time-dependent operators args = arguments(substitute(op, parameter)) - args_translated = [ - translate_qo( - arg, - b; - parameter = parameter, - time_parameter = time_parameter, - level_map = level_map, - operators = operators, - op_type = sparse, - ) for arg in args - ] - - output_op_QAdd_QO = t -> sum(args_translated[i](t) for i = 1:length(args_translated)) - return output_op_QAdd_QO + # Translate first arg to determine concrete operator type + first_translated = _translate_qo( + args[1], + b; + parameter, + time_parameter, + level_map, + operators, + op_type = sparse, + ) + OpType = + typeof(first_translated isa Function ? first_translated(0.0) : first_translated) + FW = FunctionWrapper{OpType,Tuple{Float64}} + + args_wrapped = ntuple(length(args)) do k + a_k = if k == 1 + first_translated + else + _translate_qo( + args[k], + b; + parameter, + time_parameter, + level_map, + operators, + op_type = sparse, + ) + end + FW(a_k isa Function ? a_k : (_ -> a_k)) + end + + return t -> begin + result = args_wrapped[1](t) + for i = 2:length(args_wrapped) + result = result + args_wrapped[i](t) + end + result + end end -function translate_qo( +function _translate_qo( op::QSym, b::QuantumOpticsBase.Basis; parameter = Dict(), @@ -161,32 +172,14 @@ function translate_qo( operators = Dict(), op_type = sparse, ) - # QSym are only fundamental symbolic operators, e.g. a, σ_-, x, ... - - time_parameter = _normalize_time_parameter(time_parameter) if isempty(time_parameter) - return _translate_numeric( - op, - b; - parameter = parameter, - level_map = level_map, - operators = operators, - op_type = op_type, - ) # TODO: test + return _translate_numeric(op, b; parameter, level_map, operators, op_type) end - - output_op_QSym_QO_ = _translate_numeric_raw( - op, - b; - level_map = level_map, - operators = operators, - op_type = op_type, - ) # this is faster! - output_op_QSym_QO = t -> output_op_QSym_QO_ - return output_op_QSym_QO + output_op = _translate_numeric_raw(op, b; level_map, operators, op_type) + return t -> output_op end -function translate_qo( +function _translate_qo( arg_c_, b::QuantumOpticsBase.Basis; parameter = Dict(), @@ -195,12 +188,9 @@ function translate_qo( operators = Dict(), op_type = sparse, ) - # should only be needed for numbers and symbolic parameters - arg_c = substitute(arg_c_, parameter) one_b = _translate_one(b, operators, op_type) - time_parameter = _normalize_time_parameter(time_parameter) if isempty(time_parameter) return arg_c * one_b end @@ -216,9 +206,10 @@ end function translate_qo(ops::Vector, b::QuantumOpticsBase.Basis; kwargs...) return [translate_qo(op, b; kwargs...) for op in ops] end + function translate_qo(G::SLH, b::QuantumOpticsBase.Basis; kwargs...) - L = get_lindblad(G); - H = get_hamiltonian(G) + L = lindblad(G) + H = hamiltonian(G) H_QO = translate_qo(H, b; kwargs...) L_QO = [translate_qo(L_, b; kwargs...) for L_ in L] return H_QO, L_QO diff --git a/src/utils.jl b/src/utils.jl index a6c5dc5..ba06847 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -15,7 +15,7 @@ expect(op::SQA.QNumber, state; kwargs...) = numeric_average(op, state; kwargs... substitute_operators(op, dict::Dict; replace_adjoint=true) Like `substitute(op, dict::Dict)` but with special handling for `QMul` and `QAdd`. -This is needed if an operator is substitute by a `QMul` or `QAdd`, e.g. +This is needed if an operator is substitute by a `QMul` or `QAdd`, e.g. ``a_1 \\rightarrow g_2 a_2 + g_3 a_3`` If `replace_adjoint=true`, the dictionary is extended with adjoint substitutions @@ -33,19 +33,11 @@ function substitute_operators(op::SQA.QAdd, dict::Dict; replace_adjoint = true) end function substitute_operators(op::SQA.QMul, dict::Dict; replace_adjoint = true) dict_ = replace_adjoint ? _extend_with_adjoint(dict) : dict - return substitute( - op.arg_c, - dict_, - )*prod([substitute(arg, dict::Dict) for arg ∈ op.args_nc]) + return substitute(op.arg_c, dict_) * prod([substitute(arg, dict_) for arg ∈ op.args_nc]) end function _extend_with_adjoint(dict::Dict) - dict_ = Dict() - for (k, v) in dict - dict_[k] = v - ka = SQA.adjoint(k) - va = SQA.adjoint(v) - dict_[ka] = va - end - return dict_ + pairs_ = collect(dict) + adj_pairs = [SQA.adjoint(k) => SQA.adjoint(v) for (k, v) in pairs_] + return Dict(vcat(pairs_, adj_pairs)) end diff --git a/test/runtests.jl b/test/runtests.jl index a906722..dafa4ec 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ names = [ - # "test_code_quality.jl", "test_SLH.jl", - # "test_feedback.jl", # TODO + "test_feedback.jl", "test_translate.jl", "test_example_cavity_scattering.jl", "test_compare_example_05_1_05_2.jl", @@ -10,7 +9,7 @@ names = [ ] detected_tests = - filter(name->startswith(name, "test_") && endswith(name, ".jl"), readdir(".")) + filter(name -> startswith(name, "test_") && endswith(name, ".jl"), readdir(".")) unused_tests = setdiff(detected_tests, names) if length(unused_tests) != 0 diff --git a/test/test_SLH.jl b/test/test_SLH.jl index d4e7394..9e5ea20 100644 --- a/test/test_SLH.jl +++ b/test/test_SLH.jl @@ -1,6 +1,9 @@ using QuantumInputOutput using SecondQuantizedAlgebra +using QuantumOptics using SymbolicUtils +using FunctionWrappers: FunctionWrapper +using StaticArrays using Test @testset "SLH" begin @@ -16,34 +19,32 @@ using Test gu, Δ, γ = rnumbers("g_u Δ, γ") gv = cnumber("g_v"); - G_u = SLH(1, gu'*au, 0) # input cavity + G_u = SLH(1, gu'*au, 0) # input cavity G_c = SLH(1, √(γ)*c, Δ*c'c) # system cavity G_v = SLH(1, gv'*av, 0) # output cavity - G_c_S = get_scattering(G_c) - G_c_L = get_lindblad(G_c) - G_c_H = get_hamiltonian(G_c) + G_c_S = scattering(G_c) + G_c_L = lindblad(G_c) + G_c_H = hamiltonian(G_c) @test G_c_S == G_c.scattering @test G_c_L == G_c.lindblad @test G_c_H == G_c.hamiltonian + @test G_c_S isa SMatrix{1,1} + @test G_c_L isa SVector{1} + SLH(1, [√(γ)*c], Δ*c'c) @test isequal(G_c, SLH(1, [√(γ)*c], Δ*c'c)) - @test isequal(G_c, SLH(ones(1, 1), [√(γ)*c], Δ*c'c)) - - @test isequal(ones(1, 1), G_c_S) - @test isequal([√(γ)*c], G_c_L) - @test isequal(Δ*c'c, G_c_H) @testset "simple_cascade" begin G1 = G_u ▷ G_c - @test isequal(G1.scattering, ones(1, 1)) - @test isequal(G1.lindblad[1], simplify(gu'*au + √(γ)*c)) - @test isequal( - G1.hamiltonian, - simplify(G_c_H - 1im/2*((√(γ)*c)'*(1)*gu'*au - (gu'*au)'*(1)*(√(γ)*c))), + @test scattering(G1) isa SMatrix{1,1} + @test iszero(simplify(lindblad(G1)[1] - simplify(gu'*au + √(γ)*c))) + expected_H = simplify( + hamiltonian(G_c) - 1im/2*((√(γ)*c)'*(1)*gu'*au - (gu'*au)'*(1)*(√(γ)*c)), ) + @test iszero(simplify(hamiltonian(G1) - expected_H)) G2 = cascade(G_u, G_c, G_v) G3 = G_u ▷ G_c ▷ G_v @@ -52,7 +53,7 @@ using Test @test isequal(G2, ▷(G1, G_v)) @test isequal(G2, G3) - @test iszero(simplify(G2.lindblad[1] - (gu'*au + √(γ)*c + gv'*av))) + @test iszero(simplify(lindblad(G2)[1] - (gu'*au + √(γ)*c + gv'*av))) end @testset "simple_concatenate" begin @@ -61,15 +62,117 @@ using Test Gc = concatenate(G1, G2) - @test size(Gc.scattering) == (2, 2) - @test Gc.scattering == [1 0; 0 1] - @test length(Gc.lindblad) == 2 - @test isequal(Gc.lindblad[1], gu'*au) - @test isequal(Gc.lindblad[2], √(γ)*c) - @test isequal(Gc.hamiltonian, Δ*c'c) + @test scattering(Gc) isa SMatrix{2,2} + @test size(scattering(Gc)) == (2, 2) + @test length(lindblad(Gc)) == 2 + @test isequal(lindblad(Gc)[1], gu'*au) + @test isequal(lindblad(Gc)[2], √(γ)*c) + @test isequal(hamiltonian(Gc), Δ*c'c) Gc2 = G1 ⊞ G2 @test isequal(Gc, Gc2) end + @testset "two-port symbolic cascade" begin + hu2 = FockSpace(:u2) + hv2 = FockSpace(:v2) + h2 = hu1 ⊗ hu2 ⊗ hv1 ⊗ hv2 + + au1 = Destroy(h2, :a_u1, 1) + au2 = Destroy(h2, :a_u2, 2) + av1 = Destroy(h2, :a_v1, 3) + av2 = Destroy(h2, :a_v2, 4) + + gu1, gu2, gv1, gv2, t, r = rnumbers("g_u1 g_u2 g_v1 g_v2 t r") + + G_bs = SLH([r t; t -r], [0, 0], 0) + G_out = SLH(1, gv1' * av1, 0) ⊞ SLH(1, gv2' * av2, 0) + G_cas = G_bs ▷ G_out + + @test length(lindblad(G_cas)) == 2 + @test iszero(simplify(lindblad(G_cas)[1] - (gv1' * av1))) + @test iszero(simplify(lindblad(G_cas)[2] - (gv2' * av2))) + end + + @testset "numeric type stability" begin + bc = FockBasis(4) + a_op = destroy(bc) + H_s = sparse(0.5 * dagger(a_op) * a_op) + L_s = sparse(sqrt(1.0) * a_op) + gu_f(t) = exp(-t^2) * sparse(a_op) + gv_f(t) = exp(-(t - 2)^2) * sparse(a_op) + + @testset "static SLH concreteness" begin + G = SLH(1, L_s, H_s) + @test eltype(lindblad(G)) === typeof(L_s) + @test typeof(hamiltonian(G)) === typeof(H_s) + end + + @testset "time-dep SLH wraps into FunctionWrapper" begin + G_td = SLH(1, gu_f, H_s) + @test eltype(lindblad(G_td)) <: FunctionWrapper + @test typeof(hamiltonian(G_td)) <: FunctionWrapper + end + + @testset "cascade preserves FunctionWrapper" begin + G1 = SLH(1, gu_f, H_s) + G2 = SLH(1, gv_f, H_s) + G_cas = G1 ▷ G2 + @test eltype(lindblad(G_cas)) <: FunctionWrapper + @test typeof(hamiltonian(G_cas)) <: FunctionWrapper + @inferred lindblad(G_cas)[1](0.5) + end + + @testset "cascade mixed static/time-dep wraps uniformly" begin + G_cas = SLH(1, L_s, H_s) ▷ SLH(1, gu_f, H_s) + @test eltype(lindblad(G_cas)) <: FunctionWrapper + @test eltype(lindblad(G_cas)) !== Any + end + + @testset "concatenation mixed static/time-dep wraps uniformly" begin + G_cat = SLH(1, L_s, H_s) ⊞ SLH(1, gu_f, H_s) + LT = eltype(lindblad(G_cat)) + @test LT <: FunctionWrapper + @test LT !== Any + @inferred lindblad(G_cat)[1](0.5) + @inferred lindblad(G_cat)[2](0.5) + end + + @testset "concatenation static stays static" begin + G_cat = SLH(1, L_s, H_s) ⊞ SLH(1, L_s, H_s) + @test eltype(lindblad(G_cat)) === typeof(L_s) + @test !(eltype(lindblad(G_cat)) <: FunctionWrapper) + end + + @testset "FunctionWrapper call is inferred" begin + G_td = SLH(1, gu_f, H_s) + l = lindblad(G_td)[1] + @inferred l(0.5) + end + + @testset "_op_type extracts type from FunctionWrapper SLH" begin + G_td = SLH(1, gu_f, H_s) + @test QuantumInputOutput._op_type(G_td) === typeof(H_s) + end + + @testset "_op_type returns nothing for static SLH" begin + G_s = SLH(1, L_s, H_s) + @test QuantumInputOutput._op_type(G_s) === nothing + end + + @testset "SLH with only plain closures errors" begin + bare_f(t) = t * ones(ComplexF64, 5, 5) + bare_g(t) = (1 - t) * ones(ComplexF64, 5, 5) + @test_throws ErrorException SLH([1 0; 0 1], [bare_f, bare_g], bare_g) + end + + @testset "feedback preserves FunctionWrapper type" begin + G1 = SLH(1, gu_f, H_s) + G2 = SLH(1, gv_f, H_s) + G_cat = G1 ⊞ G2 + G_fb = feedback(G_cat, 2, 1) + @test eltype(lindblad(G_fb)) <: FunctionWrapper + @test typeof(hamiltonian(G_fb)) <: FunctionWrapper + end + end end diff --git a/test/test_compare_example_05_1_05_2.jl b/test/test_compare_example_05_1_05_2.jl index 76f1284..4df671e 100644 --- a/test/test_compare_example_05_1_05_2.jl +++ b/test/test_compare_example_05_1_05_2.jl @@ -1,6 +1,7 @@ using QuantumInputOutput using SecondQuantizedAlgebra using QuantumOptics +using FunctionWrappers: FunctionWrapper using Test @testset "compare_QDs_examples_05_1_05_2" begin @@ -12,9 +13,8 @@ using Test t0 = 4σt Tend = 3t0 T = collect(0.0:0.005:1.0) .* Tend - # u1(t) = sqrt(1 / (σt * √(2π)) * exp(-0.5 * (t - t0)^2 / σt^2)) u1(t) = 1/(sqrt(σt)*π^(1/4)) * exp(-(t - t0)^2 / (2*σt^2)) - gu_t = u_to_gu(u1, T) + gu_t = coupling_input(u1, T) # -------- Example 05-1 style (symbolic -> numeric) -------- ha(i) = NLevelSpace("a$(i)", 2) @@ -36,8 +36,8 @@ using Test G_L_t = G_L(2) ▷ G_ϕ(1, 2) ▷ G_L(1) G_t = G_R_t ⊞ G_L_t - H = get_hamiltonian(G_t) - L = get_lindblad(G_t) + H = hamiltonian(G_t) + L = lindblad(G_t) L_R = L[1] L_L = L[2] @@ -75,7 +75,7 @@ using Test I_R_1 = [real(expect(L_R_QO(ti)' * L_R_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] I_L_1 = [real(expect(L_L_QO(ti)' * L_L_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] - # -------- Example 05-2 style (SLHqo, quantum pulse) -------- + # -------- Example 05-2 style (numeric SLH, quantum pulse) -------- bu = FockBasis(4) bq = tensor([ba for _ = 1:N]...) b2 = bu ⊗ bq @@ -83,23 +83,26 @@ using Test σ_qds(α, i, j) = embed(b2, 1 + α, transition(ba, i, j)) a_u = destroy(bu) ⊗ one(bq) - G_u = SLHqo(1, t -> gu_t(t) * a_u, 0 * one(b2)) - G_ϕ_qo(i, j) = SLHqo(exp(1im * ϕn[i]), 0 * one(b2), 0 * one(b2)) - G_R_qo(i) = SLHqo(1, √(γRn[i]) * σ_qds(i, 1, 2), -Δn[i] * σ_qds(i, 2, 2)) - G_L_qo(i) = SLHqo(1, √(γLn[i]) * σ_qds(i, 1, 2), 0 * one(b2)) + # Now uses unified SLH instead of SLHqo + G_u = SLH(1, t -> gu_t(t) * a_u, 0 * one(b2)) + G_ϕ_qo(i, j) = SLH(exp(1im * ϕn[i]), 0 * one(b2), 0 * one(b2)) + G_R_qo(i) = SLH(1, √(γRn[i]) * σ_qds(i, 1, 2), -Δn[i] * σ_qds(i, 2, 2)) + G_L_qo(i) = SLH(1, √(γLn[i]) * σ_qds(i, 1, 2), 0 * one(b2)) G_R_t_qo = G_u ▷ G_R_qo(1) ▷ G_ϕ_qo(1, 2) ▷ G_R_qo(2) G_L_t_qo = G_L_qo(2) ▷ G_ϕ_qo(1, 2) ▷ G_L_qo(1) G_t_qo = G_R_t_qo ⊞ G_L_t_qo - H_qo = get_hamiltonian(G_t_qo) - L_qo = get_lindblad(G_t_qo) + H_qo = hamiltonian(G_t_qo) + L_qo = lindblad(G_t_qo) L_R_qo = L_qo[1] L_L_qo = L_qo[2] - Hf = H_qo isa Function ? H_qo : (t -> H_qo) - L_R_f = L_R_qo isa Function ? L_R_qo : (t -> L_R_qo) - L_L_f = L_L_qo isa Function ? L_L_qo : (t -> L_L_qo) + # FunctionWrapper is callable but not <: Function + _callable(x) = x isa Union{Function,FunctionWrapper} + Hf = _callable(H_qo) ? H_qo : (t -> H_qo) + L_R_f = _callable(L_R_qo) ? L_R_qo : (t -> L_R_qo) + L_L_f = _callable(L_L_qo) ? L_L_qo : (t -> L_L_qo) J_add_qo = [√(γ_add[i]) * σ_qds(i, 1, 2) for i = 1:N] diff --git a/test/test_example_cavity_scattering.jl b/test/test_example_cavity_scattering.jl index 84b2787..fbb9b9e 100644 --- a/test/test_example_cavity_scattering.jl +++ b/test/test_example_cavity_scattering.jl @@ -22,8 +22,8 @@ using Test G_v = SLH(1, gv' * av, 0) G_cas = ▷(G_u, G_c, G_v) - H = get_hamiltonian(G_cas) - L = get_lindblad(G_cas)[1] + H = hamiltonian(G_cas) + L = lindblad(G_cas)[1] γ_ = 1.0 Δ_ = 0.0 @@ -34,12 +34,11 @@ using Test T_p = 1/γ_ T_end = 12T_p σ = sqrt(0.5)*T_p - # u1(t) = sqrt(1/(σ*√(2π))*exp( -0.5*(t - 4T_p)^2/σ^2 )) u1(t) = 1/(sqrt(σ)*π^(1/4)) * exp(-(t - 4σ)^2 / (2*σ^2)) T = [0:0.004:1;]*T_end ΔT = T[2] - T[1] - gu_t = u_to_gu(u1, T) + gu_t = coupling_input(u1, T) dict_p_t = Dict(gu => gu_t) bu1 = FockBasis(2) @@ -71,7 +70,7 @@ using Test # correlation matrix Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo - g1_m = two_time_corr_matrix(T, ρt, input_output_1, Ls) + g1_m = correlation_matrix(T, ρt, input_output_1, Ls) F = eigen(g1_m) n_avg = round.(real.(F.values)*ΔT; digits = 3) @@ -86,11 +85,8 @@ using Test p_num_2 = [γ_, Δ_] dict_p_2 = Dict(p_sym_2 .=> p_num_2); - # time-dependent coupling for the output mode $v(t)$ - gv_t = v_to_gv(v_mode, T) + gv_t = coupling_output(v_mode, T) - # gvc_t = t -> conj(gv_t(t)) - # dict_p_t_2 = Dict([gu, gv, conj(gv)] .=> [gu_t, gv_t, gvc_t]); dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]); H_QO_2 = translate_qo(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) @@ -101,7 +97,6 @@ using Test return H, J, dagger.(J) end; - # time evolution for the system including the output cavity t_2, ρt_2 = timeevolution.master_dynamic(T, ψ0, input_output_2) n_v1_t = real.(expect(av_qo'*av_qo, ρt_2)) diff --git a/test/test_feedback.jl b/test/test_feedback.jl index a307724..57fa494 100644 --- a/test/test_feedback.jl +++ b/test/test_feedback.jl @@ -1,7 +1,10 @@ using QuantumInputOutput using SecondQuantizedAlgebra +using QuantumOptics using SymbolicUtils +using FunctionWrappers: FunctionWrapper using LinearAlgebra +using StaticArrays using Test @testset "feedback reduction" begin @@ -44,23 +47,24 @@ using Test G_red = feedback(G, 1, 1) loop_gain = (1 - s11)^(-1) - expected_S = reshape([simplify(s22 + s21 * loop_gain * s12)], 1, 1) - expected_L = [simplify(l2 + s21 * loop_gain * l1)] + expected_S = simplify(s22 + s21 * loop_gain * s12) + expected_L = simplify(l2 + s21 * loop_gain * l1) expected_term = simplify((l1' * s11 + l2' * s21) * loop_gain * l1) expected_H = simplify(h0 + (expected_term - expected_term') / (2im)) - @test isequal((G_red.scattering .- expected_S)[1, 1], 0) - @test isequal((G_red.lindblad .- expected_L)[1], 0) - @test isequal(G_red.hamiltonian - expected_H, 0) + @test scattering(G_red) isa SMatrix{1,1} + @test abs(scattering(G_red)[1, 1] - expected_S) < 1e-10 + @test abs(lindblad(G_red)[1] - expected_L) < 1e-10 + @test abs(hamiltonian(G_red) - expected_H) < 1e-10 @testset "coherent-feedback OPO loop" begin # TODO: problem with simplify of conj(conj(x)), conj((0+1im)*x) and fractions - κ = 0.7 - ϵ = 0.45 - η = 0.65 # κ = rnumber("κ") # ϵ = rnumber("ϵ") # η = rnumber("η") + κ = 0.7 + ϵ = 0.45 + η = 0.65 r = simplify(√(1 - η^2)) @@ -70,9 +74,9 @@ using Test G_loop = feedback(G_unconnected, 1 => 2, 2 => 1) l = simplify(η / (1 + r)) - @test isequal(G_loop.scattering, ones(1, 1)) - @test isequal(G_loop.lindblad, [simplify(l * √(κ) * a)]) - @test isequal(simplify(G_loop.hamiltonian - get_hamiltonian(G_opo)), 0) + @test scattering(G_loop) isa SMatrix{1,1} + @test iszero(simplify(lindblad(G_loop)[1] - simplify(l * √(κ) * a))) + @test isequal(simplify(hamiltonian(G_loop) - hamiltonian(G_opo)), 0) end @testset "bidirectional waveguide matches cascade model" begin @@ -103,9 +107,25 @@ using Test G_network = G_in ⊞ G_qd(1) ⊞ G_phase(1, 2) ⊞ G_qd(2) G_feedback = feedback(G_network, 1 => 3, 3 => 5, 5 => 7, 8 => 6, 6 => 4, 4 => 2) - # @test isequal([G_feedback.scattering - get_scattering(G_manual)], zeros(2,2)) # TODO - @test isequal(G_feedback.lindblad[1], get_lindblad(G_manual)[2]) - @test isequal(G_feedback.lindblad[2], get_lindblad(G_manual)[1]) - @test isequal(simplify(G_feedback.hamiltonian - get_hamiltonian(G_manual)), 0) + @test isequal(lindblad(G_feedback)[1], lindblad(G_manual)[2]) + @test isequal(lindblad(G_feedback)[2], lindblad(G_manual)[1]) + @test isequal(simplify(hamiltonian(G_feedback) - hamiltonian(G_manual)), 0) + end + + @testset "feedback preserves FunctionWrapper" begin + bc = FockBasis(4) + a_op = destroy(bc) + H_s = sparse(0.5 * dagger(a_op) * a_op) + L_s = sparse(sqrt(1.0) * a_op) + gu_f(t) = exp(-t^2) * sparse(a_op) + gv_f(t) = exp(-(t - 2)^2) * sparse(a_op) + + G_cat = SLH(1, gu_f, H_s) ⊞ SLH(1, gv_f, H_s) + G_fb = feedback(G_cat, 1, 1) + + LT = eltype(lindblad(G_fb)) + @test LT <: FunctionWrapper + @test LT !== Any + @inferred lindblad(G_fb)[1](0.5) end end diff --git a/test/test_interaction_picture.jl b/test/test_interaction_picture.jl index aa634de..0033832 100644 --- a/test/test_interaction_picture.jl +++ b/test/test_interaction_picture.jl @@ -35,9 +35,9 @@ using Test G_v = SLH(1, gv_sym' * av_sym, 0) G_cas = ▷(G_u, G_s, G_v) - H = get_hamiltonian(G_cas) - L = get_lindblad(G_cas)[1] - H_uv = get_hamiltonian(▷(G_u, G_v)) + H = hamiltonian(G_cas) + L = lindblad(G_cas)[1] + H_uv = hamiltonian(▷(G_u, G_v)) H_int_sym_ = simplify(H - H_uv) # Interaction-picture operator substitution @@ -45,25 +45,31 @@ using Test 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] - # int_dict = Dict([a0_ls; adjoint.(a0_ls)] .=> [a_int_ls; adjoint.(a_int_ls)]) int_dict = Dict(a0_ls .=> a_int_ls) H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict)) L_int_sym = simplify(substitute_operators(L, int_dict)) # Virtual-cavity couplings - gu_t = u_to_gu(u, T) - gv_t = v_to_gv(u, T) + gu_t = coupling_input(u, T) + gv_t = coupling_output(u, T) # Interaction-picture coefficient matrices - A_uv = interaction_picture_A_2modes(gu_t, gv_t) - M_num = interaction_picture_M(A_uv, T) - M_ana = interaction_picture_M_2modes_equal(u, T) + A_uv = coupling_matrix((gu_t, gv_t)) + sol_M = solve_mode_evolution(A_uv, T) + M_num = t -> sol_M(t) + M_ana = solve_mode_evolution_symmetric(u, T) @test abs(maximum([maximum(abs.(M_num(t))) for t in T]) - 1) < 1e-4 max_M_err = maximum([maximum(abs.(M_num(t) - M_ana(t))) for t in T]) @test max_M_err < 5e-4 + # coupling_matrix type stability + @test A_uv(0.5) isa SMatrix{2,2,ComplexF64} + @inferred A_uv(0.5) + A_uv(0.0) # warmup + @test (@allocated A_uv(0.5)) == 0 + # Numerical basis and operators bu = FockBasis(n_ph) ba = NLevelBasis(2) @@ -77,9 +83,6 @@ using Test dict_p = Dict(γ_sym => γ) M_ls = [M(i, j) for i = 1:la for j = 1:la] M_t_ls = [t -> M_num(t)[i, j] for i = 1:la for j = 1:la] - # M_t_c_ls = [t -> conj(M_num(t)[i, j]) for i = 1:la for j = 1:la] - # p_t_sym = [gu_sym, gv_sym, M_ls..., conj.(M_ls)...] - # p_t_num = [gu_t, gv_t, M_t_ls..., M_t_c_ls...] p_t_sym = [gu_sym, gv_sym, M_ls...] p_t_num = [gu_t, gv_t, M_t_ls...] dict_p_t = Dict(p_t_sym .=> p_t_num) diff --git a/test/test_translate.jl b/test/test_translate.jl index 42507d7..8d90dd2 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -13,8 +13,8 @@ using Test ha_ = NLevelSpace("a", 2) h = hc ⊗ ha_ - a = Destroy(h, :a, 1) # cavity - σ(i, j) = Transition(h, "σ", i, j, 2) # two-level atom + a = Destroy(h, :a, 1) # cavity + σ(i, j) = Transition(h, "σ", i, j, 2) # two-level atom bc1 = FockBasis(4) a_QO = destroy(bc1) @@ -31,7 +31,6 @@ using Test E_t(t) = 2*t + 1im E_t_c(t) = conj(E_t(t)) - # dict_p_t2 = Dict( [E, conj(E)] .=> [E_t, E_t_c] ) dict_p_t2 = Dict(E => E_t) @@ -99,6 +98,7 @@ using Test 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 F5_ = translate_qo(a*conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) @test sum(abs.((F5_(0.2) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8 F6 = translate_qo(conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) @@ -152,4 +152,21 @@ using Test @test isequal(simplify(y - (a_1*a_1*c1 + a_1*c3)), 0) @test isequal(simplify(y2 - (a_2*a_2*c1 + a_2*c3)), 0) end + + @testset "substitute_operators adjoint in args_nc" begin + h2 = FockSpace(:h2) + a = Destroy(h2, :a, 1) + b = Destroy(h2, :b, 1) + ad = SecondQuantizedAlgebra._adjoint(a) + bd = SecondQuantizedAlgebra._adjoint(b) + @cnumbers g + + # QMul with adjoint operator: g * a† * a + op = g * ad * a + dict_sub = Dict(a => b) + + result = substitute_operators(op, dict_sub) + expected = g * bd * b + @test isequal(simplify(result - expected), 0) + end end diff --git a/test/test_utils.jl b/test/test_utils.jl index 2ec9b69..20b5844 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -11,17 +11,17 @@ using Test T = [0.0:0.0001:1;]*6 - gu_num = u_to_gu(u, T) - gv_num = v_to_gv(v, T) + gu_num = coupling_input(u, T) + gv_num = coupling_output(v, T) - gu_num2 = u_to_gu(u.(T), T) - gv_num2 = v_to_gv(v.(T), T) + gu_num2 = coupling_input(u.(T), T) + gv_num2 = coupling_output(v.(T), T) @test sum(abs.(gu_num2.(T) - gu_num.(T))) < 1e-9 @test sum(abs.(gv_num2.(T) - gv_num.(T))) < 1e-9 - gu_ana = u_to_gu_Gauss(τ, σ; δ = δ) - gv_ana = v_to_gv_Gauss(τ, σ; δ = δ) + gu_ana = coupling_input(Gaussian(τ, σ; δ = δ)) + gv_ana = coupling_output(Gaussian(τ, σ; δ = δ)) 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]))) @@ -34,23 +34,23 @@ using Test u1(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - (τ - 2σ)) / σ)^2) u2(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - (τ + 2σ)) / σ)^2) - v_eff_f = v_eff([v1, v2], T, 2) - v_eff_data = v_eff([v1.(T), v2.(T)], T, 2) + v_eff_f = effective_output_mode([v1, v2], T, 2) + v_eff_data = effective_output_mode([v1.(T), v2.(T)], T, 2) @test maximum(abs.(v_eff_f.(T) .- v_eff_data.(T))) < 1e-7 - gv1 = v_to_gv(v1, T) - gv2 = v_to_gv(v2, T) - v_eff_data2 = v_eff([v1.(T), v2.(T)], [gv1.(T), gv2.(T)], T, 2) + gv1 = coupling_output(v1, T) + gv2 = coupling_output(v2, T) + v_eff_data2 = effective_output_mode([v1.(T), v2.(T)], [gv1.(T), gv2.(T)], T, 2) @test maximum(abs.(v_eff_f.(T) .- v_eff_data2.(T))) < 1e-7 - u_eff_f = u_eff([u1, u2], T, 2) - u_eff_data = u_eff([u1.(T), u2.(T)], T, 2) + u_eff_f = effective_input_mode([u1, u2], T, 2) + u_eff_data = effective_input_mode([u1.(T), u2.(T)], T, 2) @test maximum(abs.(u_eff_f.(T) .- u_eff_data.(T))) / maximum(abs.(u_eff_f.(T) .+ u_eff_data.(T))) < 1e-5 - gu1 = u_to_gu(u1, T) - gu2 = u_to_gu(u2, T) - u_eff_data2 = u_eff([u1.(T), u2.(T)], [gu1.(T), gu2.(T)], T, 2) + gu1 = coupling_input(u1, T) + gu2 = coupling_input(u2, T) + 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-5 end