diff --git a/TODO.md b/TODO.md index 6ba13c92..4d503b83 100644 --- a/TODO.md +++ b/TODO.md @@ -86,7 +86,7 @@ - [x] feedback reduction - [x] interaction picture - [x] general expression for A(t) - - [x] two_time_corr_matrix for Ls const (see example 02-2) + - [x] two_time_corr_matrix for Js const (see example 02-2) - [x] Formatter/SpellCheck - [x] JET diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl index 17a44b6b..9f9cfa71 100644 --- a/benchmarks/correlations.jl +++ b/benchmarks/correlations.jl @@ -23,7 +23,7 @@ function benchmark_correlations!(SUITE) G_cas = ▷(G_u, G_c, G_v) H_sym = hamiltonian(G_cas) - L_sym = jump_operator(G_cas)[1] + J_sym = jump_operator(G_cas)[1] γ_ = 1.0 σ_pulse = 1 / γ_ @@ -40,17 +40,17 @@ function benchmark_correlations!(SUITE) dict_p_t = Dict(g_u => gu_t) H_QO = to_numeric(H_sym, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = to_numeric(L_sym, b; parameter = dict_p, time_parameter = dict_p_t) + J_QO = to_numeric(J_sym, b; parameter = dict_p, time_parameter = dict_p_t) ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0) - _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [L_QO]) + _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [J_QO]) au_qo = to_numeric(au, b) c_qo = to_numeric(c, b) - Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo + Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo ## --- Two-time correlation --- - # Pass the time-dependent operators (`H_QO`, `[L_QO]`) directly to the solver rather + # Pass the time-dependent operators (`H_QO`, `[J_QO]`) directly to the solver rather # than wrapping them in a `(t, ρ)` closure: since v0.10 `to_numeric` returns a lazy # `TimeDependentSum`, the direct path lets the solver build its integrator once and is # markedly faster than re-reading the operator from a function at every step. @@ -58,6 +58,6 @@ function benchmark_correlations!(SUITE) SUITE["Correlations"]["two-time"] = BenchmarkGroup() SUITE["Correlations"]["two-time"]["single photon cavity"] = - @benchmarkable correlation_matrix($T, $ρt, $H_QO, [$L_QO], $Ls) + @benchmarkable correlation_matrix($T, $ρt, $H_QO, [$J_QO], $Js) return nothing end diff --git a/benchmarks/interaction_picture.jl b/benchmarks/interaction_picture.jl index c4833c9a..ae5596b3 100644 --- a/benchmarks/interaction_picture.jl +++ b/benchmarks/interaction_picture.jl @@ -67,7 +67,7 @@ function benchmark_interaction_picture!(SUITE) G_cas = ▷(G_u, G_s, G_v) H = hamiltonian(G_cas) - L = jump_operator(G_cas)[1] + J = jump_operator(G_cas)[1] H_uv = hamiltonian(▷(G_u, G_v)) H_int_ = simplify(H - H_uv) @@ -82,7 +82,7 @@ function benchmark_interaction_picture!(SUITE) SUITE["Interaction Picture"]["operator substitution"]["TLS cascade"] = @benchmarkable begin simplify(substitute($H_int_, $int_dict)) - simplify(substitute($L, $int_dict)) + simplify(substitute($J, $int_dict)) end return nothing diff --git a/benchmarks/slh_algebra.jl b/benchmarks/slh_algebra.jl index 384427d5..2e35a91f 100644 --- a/benchmarks/slh_algebra.jl +++ b/benchmarks/slh_algebra.jl @@ -99,20 +99,20 @@ function benchmark_slh_algebra!(SUITE) G_t = G_R_t ⊞ G_L_t H_f = hamiltonian(G_t) - L_f = jump_operator(G_t) + J_f = jump_operator(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] + J_callables = [_callable(Ji) ? Ji : (t -> Ji) for Ji in J_f] SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide H(t)"] = @benchmarkable $Hf($t_mid) - SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide L(t)"] = @benchmarkable begin - for Li in $L_callables - Li($t_mid) + SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide J(t)"] = @benchmarkable begin + for Ji in $J_callables + Ji($t_mid) end end diff --git a/benchmarks/translation.jl b/benchmarks/translation.jl index ffae277f..150cb71d 100644 --- a/benchmarks/translation.jl +++ b/benchmarks/translation.jl @@ -15,7 +15,7 @@ function benchmark_translation!(SUITE) a = Destroy(h, :a, 1) σ(i, j) = Transition(h, Symbol("σ"), i, j, 2) - ## Derive H and L from cascade + ## Derive H and J from cascade @variables gu_sym::Real γ_sym::Real gv_sym::Real hu_ = FockSpace(:u) hc_ = FockSpace(:c) @@ -30,7 +30,7 @@ function benchmark_translation!(SUITE) G_v = SLH(1, gv_sym * av_, 0) G_cas = ▷(G_u, G_c, G_v) H_sym = hamiltonian(G_cas) - L_sym = jump_operator(G_cas)[1] + J_sym = jump_operator(G_cas)[1] ## --- Static translation (no time dependence) --- @@ -63,7 +63,7 @@ function benchmark_translation!(SUITE) time_parameter = $dict_p_t, ) - # Full cascade H and L translation + # Full cascade H and J translation γ_ = 1.0 σ_pulse = 1 / γ_ T = [0:0.002:1;] * 12σ_pulse @@ -79,9 +79,9 @@ function benchmark_translation!(SUITE) dict_p_cav = Dict([γ_sym, Δ_sym, gv_sym] .=> [γ_, 0.0, 0]) dict_p_t_cav = Dict([gu_sym, gv_sym] .=> [gu_t, gv_t]) - SUITE["Translation"]["time-dependent"]["3-cavity H+L"] = @benchmarkable begin + SUITE["Translation"]["time-dependent"]["3-cavity H+J"] = @benchmarkable begin to_numeric($H_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) - to_numeric($L_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) + to_numeric($J_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) end ## --- Closure evaluation (the ODE hot loop) --- @@ -89,15 +89,15 @@ function benchmark_translation!(SUITE) SUITE["Translation"]["closure evaluation"] = BenchmarkGroup() H_QO = to_numeric(H_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) - L_QO = to_numeric(L_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) + J_QO = to_numeric(J_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) t_mid = T[length(T)÷2] SUITE["Translation"]["closure evaluation"]["3-cavity H(t)"] = @benchmarkable $H_QO($t_mid) - SUITE["Translation"]["closure evaluation"]["3-cavity L(t)"] = - @benchmarkable $L_QO($t_mid) + SUITE["Translation"]["closure evaluation"]["3-cavity J(t)"] = + @benchmarkable $J_QO($t_mid) return nothing end diff --git a/docs/src/implementation.md b/docs/src/implementation.md index 51809215..99427231 100644 --- a/docs/src/implementation.md +++ b/docs/src/implementation.md @@ -19,7 +19,7 @@ av = Destroy(h, :a_v, 3) @variables γ::Real g_u::Complex g_v::Complex ``` -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. +An SLH component is represented as `(S, J, 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 [`hamiltonian`](@ref) and [`jump_operator`](@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. @@ -27,12 +27,12 @@ The resulting effective operators are accessed by [`hamiltonian`](@ref) and [`ju G_cas = ▷(G_u, G_s, G_v) H = hamiltonian(G_cas) -L = jump_operator(G_cas) +J = jump_operator(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, 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. +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 `J` or `H` while still using the same cascade, concatenate, and feedback rules. This can be much faster. ## Translate to numerics @@ -55,7 +55,7 @@ gu_t = coupling_input(t -> exp(-t^2), 0:0.01:5) dict_p_t = Dict(g_u => gu_t) H_QO = to_numeric(H, b; parameter=dict_p, time_parameter=dict_p_t) -L_QO = to_numeric(L, b; parameter=dict_p, time_parameter=dict_p_t) +J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t) ``` In some cases it can be useful to define your own set of numeric operators which should replace the symbolic expressions, e.g. to reduce the Hilbert space if the output cavities are not analyzed but they are already included in the symbolic derivation. Such a list of operators can be provide with the dictionary `operators`. @@ -105,16 +105,16 @@ gv2_eff = coupling_output(v2_eff, T) The dominant output modes are extracted by computing the two-time correlation matrix ```math -g^{(1)}(t_1, t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle +g^{(1)}(t_1, t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle ``` -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. +and diagonalizing it. In this package, [`correlation_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $J_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $J_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 = correlation_matrix(T, ρt, input_output_1, Ls) +Js(t) = gu_t(t) * au_qo + √(1.0) * c_qo +g1 = correlation_matrix(T, ρt, input_output_1, Js) F = eigen(g1) v_mode = F.vectors[:, end] / sqrt(T[2] - T[1]) ``` diff --git a/docs/src/theory.md b/docs/src/theory.md index e1ba6c5f..acbae397 100644 --- a/docs/src/theory.md +++ b/docs/src/theory.md @@ -86,10 +86,10 @@ This is the multi-mode extension used by [`effective_input_mode`](@ref) and ## SLH Networks for Cascaded Pulses -In the SLH formalism, each component is specified by a triple `(S, L, H)`: +In the SLH formalism, each component is specified by a triple `(S, J, H)`: - `S`cattering matrix -- `L`indblad operators +- `J`ump (Lindblad) operators - `H`amiltonian Networks are built using composition rules: @@ -99,8 +99,8 @@ Networks are built using composition rules: ```math G_1 \triangleright G_2 = \left( S_2 S_1,\; -L_2 + S_2 L_1,\; -H_1 + H_2 + \frac{1}{2i}\left(L_2^\dagger S_2 L_1 - L_1^\dagger S_2^\dagger L_2\right) +J_2 + S_2 J_1,\; +H_1 + H_2 + \frac{1}{2i}\left(J_2^\dagger S_2 J_1 - J_1^\dagger S_2^\dagger J_2\right) \right). ``` @@ -113,27 +113,27 @@ S_1 & 0 \\ 0 & S_2 \end{bmatrix},\; \begin{bmatrix} -L_1 \\ -L_2 +J_1 \\ +J_2 \end{bmatrix},\; H_1 + H_2 \right). ``` 3. **Feedback reduction**: connecting output port $x$ back into input port $y$ of an $n$-port component -$G = (S, L, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde L, \tilde H)$ with +$G = (S, J, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde J, \tilde H)$ with ```math \tilde S = S_{\bar x,\bar y} + S_{\bar x,y}(1 - S_{x,y})^{-1}S_{x,\bar y}, ``` ```math -\tilde L = L_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}L_x, +\tilde J = J_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}J_x, ``` ```math \tilde H = H + \frac{1}{2i}\left[ -\left(\sum_{j=1}^n L_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}L_x - \mathrm{h.c.} +\left(\sum_{j=1}^n J_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}J_x - \mathrm{h.c.} \right]. ``` @@ -144,7 +144,7 @@ obtain an effective SLH triple for the full problem. This describes a master equ ## Common SLH Elements -Appendix 1 of [Combes et al. (2017)](https://doi.org/10.1080/23746149.2017.1343097) lists SLH triples for commonly used components. The most frequently used elements in this package are summarized below using the standard ``(S, L, H)`` ordering. +Appendix 1 of [Combes et al. (2017)](https://doi.org/10.1080/23746149.2017.1343097) lists SLH triples for commonly used components. The most frequently used elements in this package are summarized below using the standard ``(S, J, H)`` ordering. Phase shifter (single input/output): @@ -210,10 +210,10 @@ The output field is generally multimode. To determine a proper **basis of tempor occupations, we compute the first-order correlation function ```math -g^{(1)}(t_1, t_2) = \langle \hat L_s^\dagger(t_1)\, \hat L_s(t_2) \rangle, +g^{(1)}(t_1, t_2) = \langle \hat J_s^\dagger(t_1)\, \hat J_s(t_2) \rangle, ``` -where $\hat L_s$ is the output operator (e.g., $\hat L_s = g_u(t) \hat{a}_u + \sqrt{\kappa} \hat c$ for a single sided cavity interacting with a pulse$u(t)$). +where $\hat J_s$ is the output operator (e.g., $\hat J_s = g_u(t) \hat{a}_u + \sqrt{\kappa} \hat c$ for a single sided cavity interacting with a pulse$u(t)$). The eigen-decomposition ```math diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 9e6b9eb5..26556eed 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -52,7 +52,7 @@ H = hamiltonian(G_cas) ``` ```@example tutorial -L = jump_operator(G_cas)[1] +J = jump_operator(G_cas)[1] ``` ## 2. Numerical parameters and input pulse @@ -88,7 +88,7 @@ bv1 = FockBasis(2) b = bu1 ⊗ bc1 ⊗ bv1 H_QO = to_numeric(H, b; parameter=dict_p, time_parameter=dict_p_t) -L_QO = to_numeric(L, b; parameter=dict_p, time_parameter=dict_p_t) +J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t) nothing # hide ``` @@ -99,7 +99,7 @@ We now solve the master equation. The required callback for `timeevolution.maste ```@example tutorial function input_output_1(t, ρ) Ht = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return Ht, J, dagger.(J) end @@ -111,15 +111,15 @@ nothing # hide ## 4. Two-time correlation function To extract the dominant output mode, we compute the two-time correlation matrix -``g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators. +``g^{(1)}(t_1,t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators. ```@example tutorial au_qo = to_numeric(au, b) c_qo = to_numeric(c, b) av_qo = to_numeric(av, b) -Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo -g1_m = correlation_matrix(T, ρt, input_output_1, Ls) +Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo +g1_m = correlation_matrix(T, ρt, input_output_1, Js) nothing # hide ``` @@ -161,11 +161,11 @@ dict_p_t_2 = Dict([g_u, g_v] .=> [gu_t, gv_t]) dict_p_2 = Dict([γ, Δ] .=> [γ_, Δ_]) H_QO_2 = to_numeric(H, b; parameter=dict_p_2, time_parameter=dict_p_t_2) -L_QO_2 = to_numeric(L, b; parameter=dict_p_2, time_parameter=dict_p_t_2) +J_QO_2 = to_numeric(J, b; parameter=dict_p_2, time_parameter=dict_p_t_2) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return Ht, J, dagger.(J) end nothing # hide diff --git a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl index 26c558b2..3141fc82 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 @@ -44,7 +44,7 @@ H = hamiltonian(G_cas) # -L = jump_operator(G_cas)[1] # only one Lindblad term in this example +J = jump_operator(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. @@ -77,7 +77,7 @@ nothing # hide # We now use the function [`to_numeric`](@ref) to create the numeric operators. If the kwarg `time_parameter` is provided the created operator is a time-dependent function. H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) nothing # hide @@ -86,7 +86,7 @@ nothing # hide ## time-dependent function for timeevolution.master_dynamic that returns H(t), J(t) and Jd(t) function input_output_1(t, ρ) H = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return H, J, dagger.(J) end; @@ -107,10 +107,10 @@ n_c_t = real.(expect(c_qo'*c_qo, ρt)) n_u1_t = real.(expect(au_qo'*au_qo, ρt)) nothing # hide -# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$. +# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$. -Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo -g1_m = correlation_matrix(T, ρt, H_QO, [L_QO], Ls) +Js(t) = gu_t(t)*au_qo + √(γ_)*c_qo +g1_m = correlation_matrix(T, ρt, H_QO, [J_QO], Js) nothing # hide # @@ -148,10 +148,10 @@ gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([g_u, g_v] .=> [gu_t, gv_t]) H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +J_QO_2 = to_numeric(J, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) H = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return H, J, dagger.(J) end; @@ -214,10 +214,10 @@ cdc_3 = c_3'c_3 ## we use the same Hamiltonian as before but add a depasing term to the dissipation H_QO_3 = to_numeric(H, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_3 = to_numeric(L, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) +J_QO_3 = to_numeric(J, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_3(t, ρ) H = H_QO_3(t) - J = [L_QO_3(t), √(γ_)*cdc_3] + J = [J_QO_3(t), √(γ_)*cdc_3] return H, J, dagger.(J) end; @@ -231,8 +231,8 @@ nothing # hide # -L0(t) = √(γ_)*c_3 + gu_t(t)*au_3 + gv_t(t)*av_3 -I_out = [real(expect(dagger(L0(t_3[i]))*L0(t_3[i]), ρt_3[i])) for i = 1:length(t_3)] +J0(t) = √(γ_)*c_3 + gu_t(t)*au_3 + gv_t(t)*av_3 +I_out = [real(expect(dagger(J0(t_3[i]))*J0(t_3[i]), ρt_3[i])) for i = 1:length(t_3)] n_u1_t_3 = real.(expect(au'*au, ρt_3)) n_v1_t_3 = real.(expect(av'*av, ρt_3)) diff --git a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl index d9af02d1..94966b3d 100644 --- a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl +++ b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl @@ -43,7 +43,7 @@ H = G_cas.hamiltonian # -L = G_cas.jump_operator[1] # only one Lindblad in this example +J = G_cas.jump_operator[1] # only one Lindblad in this example # @@ -71,11 +71,11 @@ bv1 = FockBasis(2) b = bu1 ⊗ ba1 ⊗ bv1 H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) H = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return H, J, dagger.(J) end; @@ -99,8 +99,8 @@ popu_e = real.(expect(σ_qo(2, 2), ρt)) popu_v_n1 = real.(expect(proj_v(1), ρt)) popu_v_n2 = real.(expect(proj_v(2), ρt)) -L0(t) = √(γ_)*σ_qo(1, 2) + √(Γ_)*au_qo + gv_(t)*av_qo -I_out = [real(expect(dagger(L0(t_[i]))*L0(t_[i]), ρt[i])) for i = 1:length(t_)] +J0(t) = √(γ_)*σ_qo(1, 2) + √(Γ_)*au_qo + gv_(t)*av_qo +I_out = [real(expect(dagger(J0(t_[i]))*J0(t_[i]), ρt[i])) for i = 1:length(t_)] I_out_int = [0.0] for i = 2:length(I_out) push!(I_out_int, I_out_int[end]+I_out[i]*Δ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 4f9daf75..ae82e75a 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 @@ -41,7 +41,7 @@ H = hamiltonian(G_cas) # -L = jump_operator(G_cas)[1] # only one jump operator in this example +J = jump_operator(G_cas)[1] # only one jump operator in this example # @@ -75,13 +75,13 @@ cdc_qo = c_qo'c_qo ## translate to numeric Hamiltonian and Lindblad H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) nothing # hide # We additionally include a cavity dephasing term and solve the dynamics. function input_output(t, ρ) H = H_QO(t) - J = [L_QO(t), √(γ_p_)*cdc_qo] + J = [J_QO(t), √(γ_p_)*cdc_qo] return H, J, dagger.(J) end; @@ -90,10 +90,10 @@ end; t_, ρt = timeevolution.master_dynamic(T, ψ0, input_output) 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$. +# We calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle J_s^\dagger(t_1) J_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 = correlation_matrix(T, ρt, input_output, Ls); +Js(t) = (gu_(t))'*au_qo + √(γ_)*c_qo +g1_m = correlation_matrix(T, ρt, input_output, Js); 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 92c5188f..a45a8083 100644 --- a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl +++ b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl @@ -48,7 +48,7 @@ G_v2 = SLH(1, g_v2' * av2, 0) G = cascade(G_s, G_v1, G_v2) H = hamiltonian(G) -L = jump_operator(G)[1] +J = jump_operator(G)[1] nothing # hide # We use the parameters quoted in the paper and initialize the emitter in the excited state. @@ -81,10 +81,10 @@ nothing # hide # correlation matrix g⁽¹⁾(t₁, t₂). H_QO_1 = to_numeric(H, b; parameter = dict_p_1) -L_QO_1 = to_numeric(L, b; parameter = dict_p_1) +J_QO_1 = to_numeric(J, b; parameter = dict_p_1) function input_output_1(t, ρ) - J = [L_QO_1] + J = [J_QO_1] return H_QO_1, J, dagger.(J) end @@ -92,8 +92,8 @@ end 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) +Js = √(γ_) * a_qo +g1_m = correlation_matrix(T, ρt_1, input_output_1, Js) F = eigen(g1_m) n_avg = real.(F.values) * ΔT @@ -116,10 +116,10 @@ dict_p_2 = Dict([γ, g, ω12] .=> [γ_, g_, ω12_]) dict_p_t_2 = Dict(g_v1 => gv1_t, g_v2 => gv2_t) H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +J_QO_2 = to_numeric(J, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return H_QO_2(t), J, dagger.(J) end 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 88ce932b..97dfd53b 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 @@ -50,7 +50,7 @@ H = hamiltonian(G_cas) # -L = jump_operator(G_cas)[1] # only one Lindblad in this example +J = jump_operator(G_cas)[1] # only one Lindblad in this example # Next, the numerical parameters and functions of the system are defined. @@ -87,10 +87,10 @@ bv1 = FockBasis(2) b = bu2 ⊗ bu1 ⊗ bs1 ⊗ bv1; H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) H = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return H, J, dagger.(J) end nothing # hide @@ -102,13 +102,13 @@ nothing # hide t_, ρt = timeevolution.master_dynamic(T, ψ0, input_output) nothing # hide -# Now we analyze the output modes with the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$. +# Now we analyze the output modes with the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle$. au1_qo = to_numeric(au1, b) σ_qo(i, j) = to_numeric(σ(i, j), b) -Ls(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) -g1_m = correlation_matrix(T, ρt, input_output, Ls); +Js(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) +g1_m = correlation_matrix(T, ρt, input_output, Js); p = heatmap( T, @@ -184,10 +184,10 @@ nothing # hide # The time-dependent couplings are used to define the numeric Hamiltonian and Lindblad term, and then solve the dynamics of the system. H_QO_2 = to_numeric(H, b; parameter = dict_p_out, time_parameter = dict_p_t_out) -L_QO_2 = to_numeric(L, b; parameter = dict_p_out, time_parameter = dict_p_t_out) +J_QO_2 = to_numeric(J, b; parameter = dict_p_out, time_parameter = dict_p_t_out) function input_output_2(t, ρ) H = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return H, J, dagger.(J) end 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 673a9df6..4ea18209 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 @@ -46,11 +46,11 @@ H1 = hamiltonian(G_cav_L_R_drive) # -L1_L = jump_operator(G_cav_L_R_drive)[1] +J1_L = jump_operator(G_cav_L_R_drive)[1] # -L1_R = jump_operator(G_cav_L_R_drive)[2] +J1_R = jump_operator(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 [`to_numeric`](@ref). @@ -73,9 +73,9 @@ a_QO = destroy(bc1) ops_dict = Dict([a, a'] .=> [a_QO, dagger(a_QO)]) H1_QO = to_numeric(H1, bc1; parameter = dict_p1, operators = ops_dict) -L1_L_QO = to_numeric(L1_L, bc1; parameter = dict_p1, operators = ops_dict) -L1_R_QO = to_numeric(L1_R, bc1; parameter = dict_p1, operators = ops_dict) -J1_QO = [L1_L_QO, L1_R_QO] +J1_L_QO = to_numeric(J1_L, bc1; parameter = dict_p1, operators = ops_dict) +J1_R_QO = to_numeric(J1_R, bc1; parameter = dict_p1, operators = ops_dict) +J1_QO = [J1_L_QO, J1_R_QO] nothing # hide # @@ -140,11 +140,11 @@ H2 = G_ac_drive.hamiltonian # -L2_L = G_ac_drive.jump_operator[1] +J2_L = G_ac_drive.jump_operator[1] # -L2_R = G_ac_drive.jump_operator[2] +J2_R = G_ac_drive.jump_operator[2] # @@ -186,15 +186,15 @@ a_QO2 = to_numeric(a, b) ## translate to numeric Hamiltonian and Lindblad H_QO = to_numeric(H2, b; parameter = dict_p2, time_parameter = dict_p_t2) -L2_L_QO = to_numeric(L2_L, b; parameter = dict_p2, time_parameter = dict_p_t2) -L2_R_QO = to_numeric(L2_R, b; parameter = dict_p2) +J2_L_QO = to_numeric(J2_L, b; parameter = dict_p2, time_parameter = dict_p_t2) +J2_R_QO = to_numeric(J2_R, b; parameter = dict_p2) ## additional atomic decay into free space J_add = [√(γn)*σ_QO(α, 1, 2) for α = 1:Natoms] function input_output(t, ρ) H = H_QO(t) - J = [L2_L_QO(t), L2_R_QO, J_add...] + J = [J2_L_QO(t), J2_R_QO, J_add...] return H, J, dagger.(J) end nothing # hide @@ -208,13 +208,13 @@ nothing # hide # -L2_L_QO_dag(t) = dagger(L2_L_QO(t)) +J2_L_QO_dag(t) = dagger(J2_L_QO(t)) l_t = length(t2_) n_trans2 = zeros(l_t) n_ref2 = zeros(l_t) for it = 1:l_t - n_trans2[it] = abs(expect(dagger(L2_R_QO)*L2_R_QO, ρt2[it])) - n_ref2[it] = abs(expect(L2_L_QO_dag(t2_[it])*L2_L_QO(t2_[it]), ρt2[it])) + n_trans2[it] = abs(expect(dagger(J2_R_QO)*J2_R_QO, ρt2[it])) + n_ref2[it] = abs(expect(J2_L_QO_dag(t2_[it])*J2_L_QO(t2_[it]), ρt2[it])) end nothing # hide 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 117fd513..544b50f1 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 @@ -45,16 +45,16 @@ H1 = hamiltonian(G_cav_L_R_drive) # -L1_L = jump_operator(G_cav_L_R_drive)[1] +J1_L = jump_operator(G_cav_L_R_drive)[1] # -L1_R = jump_operator(G_cav_L_R_drive)[2] +J1_R = jump_operator(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. -eqs_a = meanfield([a], H1, [L1_L, L1_R]) +eqs_a = meanfield([a], H1, [J1_L, J1_R]) complete!(eqs_a) nothing # hide @@ -153,16 +153,16 @@ H2 = G_ac_drive.hamiltonian # -L2_L = G_ac_drive.jump_operator[1] +J2_L = G_ac_drive.jump_operator[1] # -L2_R = G_ac_drive.jump_operator[2] +J2_R = G_ac_drive.jump_operator[2] # We derive the equations of motion for system with a second-order mean-field approximation. J_add = [√(γ)*σ(α, 1, 2) for α = 1:Natoms] -eqs2 = meanfield([a'a, σ(1, 2, 2)], H2, [L2_L, L2_R, J_add...]; order = 2) +eqs2 = meanfield([a'a, σ(1, 2, 2)], H2, [J2_L, J2_R, J_add...]; order = 2) nothing # hide # ```math 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 09b72af3..4bffd223 100644 --- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl +++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl @@ -55,12 +55,12 @@ H = hamiltonian(G_t) # -L = jump_operator(G_t) -L_R = L[1] +J = jump_operator(G_t) +J_R = J[1] # -L_L = L[2] +J_L = J[2] # Note that this Hamiltonian and Lindblad terms (without the drive) describe the collective decay of the quantum dots. @@ -100,15 +100,15 @@ ba = NLevelBasis(2) b = tensor([ba for i = 1:N]...) H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) -L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) +J_R_QO = to_numeric(J_R, b; parameter = dict_p, time_parameter = dict_p_t) +J_L_QO = to_numeric(J_L, b; parameter = dict_p, time_parameter = dict_p_t) σ_qo(α, i, j) = to_numeric(σ(α, 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...] + J = [J_R_QO(t), J_L_QO(t), J_add...] return Ht, J, dagger.(J) end nothing # hide @@ -128,10 +128,10 @@ 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])) + JR = J_R_QO(ti) + JL = J_L_QO(ti) + I_R[i] = real(expect(JR'JR, ρt[i])) + I_L[i] = real(expect(JL'JL, ρt[i])) end nothing # hide @@ -162,34 +162,34 @@ nothing # hide # Materialize the lazy `TimeDependentSum` to a concrete operator at each time, so the # quantum-regression products below give a plain operator usable as the solver's initial state. -L0(t) = dense(static_operator(L_R_QO(t))) -L0_dag(t) = dagger(L0(t)) -L0_ref(t) = dense(static_operator(L_L_QO(t))) -L0_ref_dag(t) = dagger(L0_ref(t)) +J0(t) = dense(static_operator(J_R_QO(t))) +J0_dag(t) = dagger(J0(t)) +J0_ref(t) = dense(static_operator(J_L_QO(t))) +J0_ref_dag(t) = dagger(J0_ref(t)) for it1 = 1:(lT-1) ρ_t1 = ρt[it1] t_2, ρ_2 = timeevolution.master_dynamic( T[it1:end], - L0(T[it1]) * ρ_t1 * L0_dag(T[it1]), + J0(T[it1]) * ρ_t1 * J0_dag(T[it1]), input_output, ) ## transmission - G2_ls = real.([expect(L0_dag(t_2[j]) * L0(t_2[j]), ρ_2[j]) for j = 1:length(t_2)]) + G2_ls = real.([expect(J0_dag(t_2[j]) * J0(t_2[j]), ρ_2[j]) for j = 1:length(t_2)]) G2[it1, it1:end] = G2_ls G2[it1:end, it1] = G2_ls t_2_r, ρ_2_r = timeevolution.master_dynamic( T[it1:end], - L0_ref(T[it1]) * ρ_t1 * L0_ref_dag(T[it1]), + J0_ref(T[it1]) * ρ_t1 * J0_ref_dag(T[it1]), input_output, ) ## reflection G2_ls_r = real.([ - expect(L0_ref_dag(t_2_r[j]) * L0_ref(t_2_r[j]), ρ_2_r[j]) for j = 1:length(t_2_r) + expect(J0_ref_dag(t_2_r[j]) * J0_ref(t_2_r[j]), ρ_2_r[j]) for j = 1:length(t_2_r) ]) G2_ref[it1, it1:end] = G2_ls_r G2_ref[it1:end, it1] = G2_ls_r 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 4ed14413..914da658 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 @@ -65,18 +65,18 @@ G_L_t = G_L(2) ▷ G_ϕ(1, 2) ▷ G_L(1) G_t = G_R_t ⊞ G_L_t 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. +# 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 $J$ will also be time-dependent. H = hamiltonian(G_t) -L = jump_operator(G_t) -L_R = L[1] -L_L = L[2] +J = jump_operator(G_t) +J_R = J[1] +J_L = J[2] J_add = [√(γ_add[i]) * σ(i, 1, 2) for i = 1:N] function input_output(t, ρ) Ht = H(t) - J = [L_R(t), L_L(t), J_add...] + J = [J_R(t), J_L(t), J_add...] return Ht, J, dagger.(J) end nothing # hide @@ -95,10 +95,10 @@ nothing # hide I_R = zeros(length(t)) I_L = zeros(length(t)) for (i, ti) in enumerate(t) - LR = L_R(ti) - LL = L_L(ti) - I_R[i] = real(expect(LR' * LR, ρt[i])) - I_L[i] = real(expect(LL' * LL, ρt[i])) + JR = J_R(ti) + JL = J_L(ti) + I_R[i] = real(expect(JR' * JR, ρt[i])) + I_L[i] = real(expect(JL' * JL, ρt[i])) end nothing # hide 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 e2db6850..804283a9 100644 --- a/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl +++ b/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl @@ -51,9 +51,9 @@ nothing # hide # second one to the transmitted right-moving output. H = hamiltonian(G_t) -L = jump_operator(G_t) -L_L = L[1] -L_R = L[2] +J = jump_operator(G_t) +J_L = J[1] +J_R = J[2] nothing # hide # @@ -90,15 +90,15 @@ ba = NLevelBasis(2) b = tensor([ba for _ = 1:N]...) H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) -L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) +J_R_QO = to_numeric(J_R, b; parameter = dict_p, time_parameter = dict_p_t) +J_L_QO = to_numeric(J_L, b; parameter = dict_p, time_parameter = dict_p_t) σ_qo(α, i, j) = to_numeric(σ(α, 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...] + J = [J_R_QO(t), J_L_QO(t), J_add...] return Ht, J, dagger.(J) end nothing # hide @@ -116,10 +116,10 @@ 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])) + JR = J_R_QO(ti) + JL = J_L_QO(ti) + I_R[i] = real(expect(JR' * JR, ρt[i])) + I_L[i] = real(expect(JL' * JL, ρt[i])) end 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 da518eb4..a1ee55ff 100644 --- a/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl +++ b/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl @@ -42,7 +42,7 @@ H = hamiltonian(G_cas) # -L = jump_operator(G_cas)[1] +J = jump_operator(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}$. @@ -68,7 +68,7 @@ H_int_sym = simplify(substitute(H_int_sym_, int_dict)) # -L_int_sym = simplify(substitute(L, int_dict)) +J_int_sym = simplify(substitute(J, int_dict)) # The above Hamiltonian and Lindblad operator are the ones in the interaction picture of the virtual cavity interaction. # We define the numerical parameters of the system, calculate the solution for the coefficient matrix $M(t)$ and solve the time evolution of the system. @@ -110,11 +110,11 @@ bv = FockBasis(5) b = bu ⊗ ba ⊗ bv H_int_QO = to_numeric(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) Ht = H_int_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return Ht, J, dagger.(J) end nothing # hide diff --git a/examples/07-1_beamsplitter_loss__quantum-pulse.jl b/examples/07-1_beamsplitter_loss__quantum-pulse.jl index 85eeaf9f..e053ae8c 100644 --- a/examples/07-1_beamsplitter_loss__quantum-pulse.jl +++ b/examples/07-1_beamsplitter_loss__quantum-pulse.jl @@ -51,12 +51,12 @@ H = hamiltonian(G) # -L = jump_operator(G) -L[1] +J = jump_operator(G) +J[1] # -L[2] +J[2] # @@ -93,11 +93,11 @@ av_qo = one(bu) ⊗ destroy(bv) ## translate to numeric operators H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] +J_QO = [to_numeric(Ji, b; parameter = dict_p, time_parameter = dict_p_t) for Ji in J] function input_output(t, ρ) Ht = H_QO(t) - J = [L_QO[1](t), L_QO[2](t)] + J = [J_QO[1](t), J_QO[2](t)] return Ht, J, dagger.(J) end nothing # hide diff --git a/examples/07-2_hong-ou-mandel__quantum-pulse.jl b/examples/07-2_hong-ou-mandel__quantum-pulse.jl index 03dbae19..a0903740 100644 --- a/examples/07-2_hong-ou-mandel__quantum-pulse.jl +++ b/examples/07-2_hong-ou-mandel__quantum-pulse.jl @@ -50,12 +50,12 @@ H = hamiltonian(G) # -L = jump_operator(G) -L[1] +J = jump_operator(G) +J[1] # -L[2] +J[2] # @@ -106,11 +106,11 @@ av2_qo = one(bu1) ⊗ one(bu2) ⊗ one(bv1) ⊗ destroy(bv2) ## translate to numeric operators H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] +J_QO = [to_numeric(Ji, b; parameter = dict_p, time_parameter = dict_p_t) for Ji in J] function input_output(t, ρ) Ht = H_QO(t) - J = [L_QO[1](t), L_QO[2](t)] + J = [J_QO[1](t), J_QO[2](t)] return Ht, J, dagger.(J) end nothing # hide @@ -147,7 +147,7 @@ n_v2_coll(t) = (t > 0.9*T[end])*R*av2_qo'av2_qo function input_output_mc(t, ρ) Ht = H_QO(t) - J = [L_QO[1](t), L_QO[2](t), n_v1_coll(t), n_v2_coll(t)] + J = [J_QO[1](t), J_QO[2](t), n_v1_coll(t), n_v2_coll(t)] return Ht, J, dagger.(J) end diff --git a/examples/08-1_pulse-delay__simple.jl b/examples/08-1_pulse-delay__simple.jl index ac01fe24..fbf3a893 100644 --- a/examples/08-1_pulse-delay__simple.jl +++ b/examples/08-1_pulse-delay__simple.jl @@ -49,7 +49,7 @@ G_v2 = concatenate(SLH(1, 0, 0), G_v) G_cas = cascade(G_u2, G_d, G_v2) H = hamiltonian(G_cas) -L = jump_operator(G_cas) +J = jump_operator(G_cas) nothing # hide # @@ -83,11 +83,11 @@ bv = FockBasis(n) b = bu ⊗ bd ⊗ bv H_QO = to_numeric(H, b; time_parameter = dict_p_t) -L_QO = [to_numeric(L[i], b; time_parameter = dict_p_t) for i = 1:length(L)] +J_QO = [to_numeric(J[i], b; time_parameter = dict_p_t) for i = 1:length(J)] function input_output(t, ρ) Ht = H_QO(t) - Jt = [L_QO[i](t) for i = 1:length(L_QO)] + Jt = [J_QO[i](t) for i = 1:length(J_QO)] return Ht, Jt, dagger.(Jt) end nothing # hide @@ -150,12 +150,12 @@ H_int_sym = simplify(substitute(H_int_sym_, int_dict)) # -L_int_sym = simplify.(substitute.(L, Ref(int_dict))) -L_int_sym[1] +J_int_sym = simplify.(substitute.(J, Ref(int_dict))) +J_int_sym[1] # -L_int_sym[2] +J_int_sym[2] # @@ -203,19 +203,19 @@ operators = Dict( ) H_int_QO = to_numeric(H_int_sym, b_int; time_parameter = dict_p_t_int, operators) -L_int_QO = [ - to_numeric(L_int_sym[i], b_int; time_parameter = dict_p_t_int, operators) for - i = 1:length(L_int_sym) +J_int_QO = [ + to_numeric(J_int_sym[i], b_int; time_parameter = dict_p_t_int, operators) for + i = 1:length(J_int_sym) ] ## H_int_QO = to_numeric(H_int_sym, b_int; time_parameter=dict_p_t_int) # hide -## L_int_QO = [to_numeric(L_int_sym[i], b_int; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide +## J_int_QO = [to_numeric(J_int_sym[i], b_int; time_parameter=dict_p_t_int) for i=1:length(J_int_sym)] # hide nothing # hide # function input_output_int(t, ρ) Ht = H_int_QO(t) - Jt = [L_int_QO[i](t) for i = 1:length(L_int_QO)] + Jt = [J_int_QO[i](t) for i = 1:length(J_int_QO)] return Ht, Jt, dagger.(Jt) end 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 fa308860..84b6b5f7 100644 --- a/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl +++ b/examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.jl @@ -36,7 +36,7 @@ nothing # hide # S_loop = scattering(G_loop) -L_loop = jump_operator(G_loop)[1] +J_loop = jump_operator(G_loop)[1] # diff --git a/examples/10-1_SUPER_excitation.jl b/examples/10-1_SUPER_excitation.jl index 4c1b9f9f..8b52b997 100644 --- a/examples/10-1_SUPER_excitation.jl +++ b/examples/10-1_SUPER_excitation.jl @@ -1,403 +1,403 @@ -# # Input-Output Analysis of Quantum Dot SUPER Excitation -# -# This example analyzes the SUPER excitation scheme for quantum dots with the input-output formalism [J. Kerber et al., TODO](https://arxiv.org/). Two red-detuned pulses allow for a close to 100% excitation of a two-level quantum emitter. At the microscopic level, the SUPER mechanism exhibits its nonlinear three-photon Raman-type character, leading to a net photon-number change of −2 in one mode and +1 in the other. - -# In the first part we describe the dynamics within a cumulant expansion approach for coherent light. We then transform into the interaction-picture of the input and output cavities, which allows us to describe the interaction with large Fock states. - -# We start by loading the needed packages and specifying the model. - -using QuantumInputOutput -using SecondQuantizedAlgebra -using QuantumCumulants -using ModelingToolkitBase -using OrdinaryDiffEq -using QuantumOptics -using Plots -using LaTeXStrings -sol_values(sol, op, eqs) = get_solution(sol, op, eqs).(sol.t) - -## Hilbert spaces -hu2 = FockSpace(:u2) # virtual input cavity (u2) -hu1 = FockSpace(:u1) # virtual input cavity (u1) -hs1 = NLevelSpace(:atom, 2) # TLS -hv1 = FockSpace(:v1) # virtual output cavity (v1) -hv2 = FockSpace(:v2) # virtual output cavity (v2) -h = tensor(hu2, hu1, hs1, hv1, hv2) - -## Operators -au2 = Destroy(h, :a_u2, 1) -au1 = Destroy(h, :a_u1, 2) -s(i, j) = Transition(h, :s, i, j) -av1 = Destroy(h, :a_v1, 4) -av2 = Destroy(h, :a_v2, 5) - -## Symbolic parameters: decay rate and virtual cavity couplings -@variables γ::Real gu1::Number gu2::Number gv1::Number gv2::Number -@independent_variables t # Symbolic time variable - -## SLH triplets -G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 -G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 -G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system -G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 -G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 - -## cascade SLH triplets -G_cas = ▷(G_u2, G_u1, G_2lvl, G_v1, G_v2) -nothing # hide - -# - -## Hamiltonian and Lindbladian -Hcas = hamiltonian(G_cas) -Lcas = jump_operator(G_cas)[1] -Lcasd = adjoint(Lcas) -nothing # hide - -# To deal with time-dependent functions in QuantumCumulants, we need to register them. Furthermore, due to a problem for the conjugate of registered functions (conj is ignored), we first need to create the adjoint of the jump operators and then substitute the time-dependent functions. - -## Time-dependent couplings -@register_symbolic gu1_t(t) -@register_symbolic gu2_t(t) -@register_symbolic gv1_t(t) -@register_symbolic gv2_t(t) - -g_ls = [gu2, gu1, gv1, gv2] -gt_ls = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] -dict_gt = Dict(g_ls .=> gt_ls) - -## Insert time-dependence -Hcas_t = substitute(Hcas, dict_gt) -Lcas_t = substitute(Lcas, dict_gt) -Lcasd_t = substitute(Lcasd, dict_gt) -nothing # hide - -# We calculate the coupling for the input and output cavities. The modified couplings of the second input and output modes are obtained with the function [`effective_input_mode`](@ref) and [`effective_output_mode`](@ref), respectively. Due to the fast oscillations, the tolerance of the numeric solver needs to be improved. The parameters are taken from [T. K. Bracht et al., PRX Quantum 2, 040354 (2021)](https://doi.org/10.1103/PRXQuantum.2.040354). - -## Time grid -dt = 1e-4 -Tend = 12*2 -T = [dt:dt:Tend;] - -## Numerical pulse parameters in THz and ps -γ_ = 1e-2 # TLS decay rate -Δ1_ = -2π*1.934 # detuning pulse 1 -Δ2_ = -2π*4.634 # detuning pulse 2 -Α1_ = 22.65π # pulse area 1 -Α2_ = 19.29π # pulse area 2 -σ1_ = 2.4 # temporal FWHM pulse 1 -σ2_ = 3.04 # temporal FWHM pulse 2 -τ1_ = 0.0 + 12 # time shift pulse 1 -τ2_ = -0.73 + 12 # time shift pulse 2 - -## Normalized input modes -u1(t_) = 1/(√(σ1_)*π^(1/4)) * exp(-(t_ - τ1_)^2 / (2*σ1_^2)) * exp(-1im*Δ1_*t_) -u2(t_) = 1/(√(σ2_)*π^(1/4)) * exp(-(t_ - τ2_)^2 / (2*σ2_^2)) * exp(-1im*Δ2_*t_) - -## Coupling functions -gu1_t_ = coupling_input(u1, T) -gu1_t(t) = gu1_t_(t) -gv1_t_ = coupling_output(u1, T) -gv1_t(t) = gv1_t_(t) - -## Cascade-modified couplings (effective modes) -abstol = 1e-10 -reltol = 1e-10 -u_fcts = [u1, u2] -u2_eff = effective_input_mode(u_fcts, T, 2; abstol, reltol) -gu2_t_ = coupling_input(u2_eff, T) -gu2_t(t) = gu2_t_(t) - -v_fcts = [u1, u2] # output modes = input modes -v2_eff = effective_output_mode(v_fcts, T, 2; abstol, reltol) -gv2_t_ = coupling_output(v2_eff, T) -gv2_t(t) = gv2_t_(t) -nothing # hide - -# After deriving the mean-field equations we define the initial state, create the ODE problem and solve the dynamics. - -## First-order cumulant expansion -order = 1 -ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] -eqs = meanfield(ops, Hcas_t, [Lcas_t]; Jdagger = [Lcasd_t], order = order, iv = t) - -## Coherent-state amplitudes -α1 = Α1_ / (2*√(2)*π^(1/4)*√(σ1_*γ_)) # field 1 -α2 = Α2_ / (2*√(2)*π^(1/4)*√(σ2_*γ_)) # field 2 -u0 = [α1, α2, 0, 0, 0, 0.0im] - -## Solve ODE system -sys = mtkcompile(System(eqs; name = :sys)) -u0_p_map_cas = Dict([unknowns(sys); γ] .=> [u0; γ_]) -prob_cas = ODEProblem(sys, u0_p_map_cas, (dt, Tend)) -sol = solve(prob_cas, Tsit5(); abstol, reltol) -nothing # hide - -## Expectation values -t_cas = sol.t # time vector -s22_cas = sol_values(sol, s(2, 2), eqs) -nu1_cas = abs2.(sol_values(sol, au1, eqs)) -nu2_cas = abs2.(sol_values(sol, au2, eqs)) -nv1_cas = abs2.(sol_values(sol, av1, eqs)) -nv2_cas = abs2.(sol_values(sol, av2, eqs)) -nothing # hide - -# - -common = (; - xlims = (t_cas[1]-0.01, t_cas[end]), - tickfontsize = 18, - guidefontsize = 18, - legendfontsize = 18, -) -p1 = plot(t_cas, real.(s22_cas); color = :red, label = L"\mathrm{cascade}") -p2 = plot( - t_cas, - nu1_cas; - color = :blue, - label = L"\langle \hat{n}_{u_1} \rangle", - ylabel = L"\langle\hat{n}_i\rangle", - xticks = ([0, 10, 20], ["", "", ""]), - yticks = ([0, 5e3, 10e3, 15e3], [L"0", L"5\cdot10^3", L"10\cdot10^3", L"15\cdot10^3"]), - ylims = (nv1_cas[1]-0.5e3, nu1_cas[1]+0.5e3), - legend = :left, - common..., -) -plot!(p2, t_cas, nv1_cas; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") -plot!( - p2, - t_cas, - nu2_cas; - color = :blue, - ls = :dash, - label = L"\langle \hat{n}_{u_2} \rangle", -) -plot!( - p2, - t_cas, - nv2_cas; - color = :red, - ls = :dash, - label = L"\langle \hat{n}_{v_2} \rangle", -) -p3 = plot( - t_cas, - nu1_cas .+ nv1_cas .- nu1_cas[1]; - color = :blue, - label = L"\mathrm{mode~1}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - xticks = ([0, 10, 20], latexstring.([0, 10, 20])), - yticks = ([-4, -2, 0, 2, 4], latexstring.([-4, -2, 0, 2, 4])), - legend = :topleft, - common..., -) -plot!(p3, t_cas, nu2_cas .+ nv2_cas .- nu2_cas[1]; color = :red, label = L"\mathrm{mode~2}") -pl1 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) -display(pl1) - -# ## Interaction picture - -## In the following, we will transform into the interaction picture of the virtual cavities. - -## Interaction picture: cavity dynamics -H_uv = hamiltonian(▷(G_u2, G_u1, G_v1, G_v2)) -H_int_ = simplify(Hcas - H_uv) - -M(i, j) = Symbolics.variable(Symbol("M_{$(i)$(j)}"); T = Number) -a0_ls = [au2, au1, av1, av2] -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 .=> a_int_ls) - -H_int = substitute(H_int_, int_dict) -L_int = simplify(substitute(Lcas, int_dict)) -Ld_int = simplify(substitute(Lcasd, int_dict)) -nothing # hide - -## Coefficient matrix M -Mat = Matrix{Any}(undef, la, la) -mod = @__MODULE__ # hide -for i = 1:la, j = 1:la - name = Symbol("Ma_$(i)$(j)") - @eval @register_symbolic $name(t) - Mat[i, j] = getfield(mod, name)(t) -end -Mat_ls = [Mat[i, j] for i = 1:la for j = 1:la] -M_ls = [M(i, j) for i = 1:la for j = 1:la] - -## Time-evolution of the matrix M(t) -A_uv = coupling_matrix((gu2_t_, gu1_t_, gv1_t_, gv2_t_)) -M_t = solve_mode_evolution(A_uv, T) - -for i = 1:la, j = 1:la - fname = Symbol("Ma_$(i)$(j)") - @eval begin - $fname(t) = M_t(t)[$i, $j] - end -end -dict_Mt = Dict(M_ls .=> Mat_ls) -dict_gt_Mt = merge(dict_gt, dict_Mt) - -H_int_t = substitute(H_int, dict_gt_Mt) -L_int_t = substitute(L_int, dict_gt_Mt) -Ld_int_t = substitute(Ld_int, dict_gt_Mt) - -eqs_int = meanfield(ops, H_int_t, [L_int_t]; Jdagger = [Ld_int_t], order = order, iv = t); - -## Solve ODE system in interaction picture -sys_int = mtkcompile(System(eqs_int; name = :sysI)) -u0_p_map_int = Dict([unknowns(sys_int); γ] .=> [u0; γ_]) -prob_int = ODEProblem(sys_int, u0_p_map_int, (dt, Tend)) -sol_int = solve(prob_int, Tsit5(); abstol, reltol) -nothing # hide - -## Expectation values -t_int = sol_int.t -s22_int = real.(sol_values(sol_int, s(2, 2), eqs_int)) -nu1_int = abs2.(sol_values(sol_int, au1, eqs_int)) -nu2_int = abs2.(sol_values(sol_int, au2, eqs_int)) -nv1_int = abs2.(sol_values(sol_int, av1, eqs_int)) -nv2_int = abs2.(sol_values(sol_int, av2, eqs_int)) - -# - -pl4 = plot( - t_int, - nu1_int .- nu1_int[1]; - color = :blue, - label = L"\mathrm{mode~1~(int.)}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - xlims = (t_int[1]-0.01, t_int[end]), - yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), - legend = :right, - tickfontsize = 18, - guidefontsize = 18, - legendfontsize = 18, - size = (800, 400), -) -plot!(pl4, t_int, nu2_int .- nu2_int[1]; color = :red, label = L"\mathrm{mode~2~(int.)}") -display(pl4) - -# ## Fock state input - -# Let us now compare the dynamics for coherent input pulses with the case of incident non-classical photon number eigenstates (Fock states), where we choose states with the same mean photon numbers as for the coherent pulses. Since in the atom-field interaction only a few photons are exchanged, the quantum states of the excitation pulses are only changed by a couple of photons, we only need to keep a couple of nearby Fock states in the computational basis. - -# We define the basis of the system and create the dictionary for the time-dependent variables to translate the Hamiltonian and Lindblad operator to a QuantumOptics.jl operator. - -n1_fock = round(Int, abs2(α1)) -n2_fock = round(Int, abs2(α2)) - -bu1 = FockBasis(n1_fock+2, n1_fock-6) -bu2 = FockBasis(n2_fock+3, n2_fock-3) -bs1 = NLevelBasis(2) -bv1 = FockBasis(1) -bv2 = FockBasis(1) -b = tensor([bu2, bu1, bs1, bv1, bv2]...) - -g_t_ls = [gu2_t_, gu1_t_, gv1_t_, gv2_t_] -M_t_ls = [t -> M_t(t)[i, j] for i = 1:la for j = 1:la] -dict_fock = Dict([g_ls; M_ls] .=> [g_t_ls; M_t_ls]) - -H_int_fock = to_numeric(H_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) -L_int_fock = to_numeric(L_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) - -## To solve the dynamics, we create the time-dependent function for the open quantum system and define the initial state. - -function input_output(t, ρ) - Ht = H_int_fock(t) - J = [L_int_fock(t)] - return Ht, J, QuantumOptics.dagger.(J) -end - -# initial state -ψu2 = fockstate(bu2, n2_fock) -ψu1 = fockstate(bu1, n1_fock) -ψs1 = nlevelstate(bs1, 1) -ψv1 = fockstate(bv1, 0) -ψv2 = fockstate(bv2, 0) -ψ0 = tensor(ψu2, ψu1, ψs1, ψv1, ψv2) - -T_fock = [0:0.001:1;]*T[end] -# t_fock, ρt_fock = timeevolution.master_dynamic(T_fock, ψ0, input_output; abstol, reltol) -using Random -Random.seed!(1) # hide -t_fock, ρt_fock = timeevolution.mcwf_dynamic(T_fock, ψ0, input_output; abstol, reltol) - -# Due to the relatively long computation time of timeevolution.master_dynamic, we simulate a single trajectory with timeevolution.mcwf_dynamic. - -## Expectation values -s22_fock = real.(expect(s(2, 2), ρt_fock)) -nu1_fock = real.(expect(au1'au1, ρt_fock)) -nu2_fock = real.(expect(au2'au2, ρt_fock)) -nv1_fock = real.(expect(av1'av1, ρt_fock)) -nv2_fock = real.(expect(av2'av2, ρt_fock)) - -common = (; - xlims = (t_fock[1]-0.01, t_fock[end]), - tickfontsize = 18, - guidefontsize = 18, - legendfontsize = 18, -) -p3_1 = plot(t_int, s22_int; color = :blue, label = L"\mathrm{coherent state}") -plot!(p3_1, t_fock, s22_fock; color = :red, label = L"\mathrm{Fock state}") -p3_2 = plot( - t_fock, - nu1_fock .- nu1_fock[1]; - color = :green, - label = L"\mathrm{Fock: mode~1}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - xlims = (t_fock[1]-0.01, t_fock[end]), - yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), - legend = :right, - common..., -) -plot!( - p3_2, - t_fock, - nu2_fock .- nu2_fock[1]; - color = :yellow, - label = L"\mathrm{Fock: mode~2}", -) -plot!( - p3_2, - t_int, - nu1_int .- nu1_int[1]; - color = :blue, - label = L"\mathrm{Coherent: mode~1}", -) -plot!( - p3_2, - t_int, - nu2_int .- nu2_int[1]; - color = :red, - label = L"\mathrm{Coherent: mode~2}", -) -pl3 = plot(p3_1, p3_2; layout = (2, 1), size = (800, 600)) -display(pl3) - -# Due to the vanishing relative phase of the Fock states, the oscillations disappear. - -# ## Package versions - -# These results were obtained using the following versions: - -using InteractiveUtils -versioninfo() - -using Pkg -Pkg.status( - [ - "QuantumInputOutput", - "SecondQuantizedAlgebra", - "QuantumCumulants", - "ModelingToolkitBase", - "OrdinaryDiffEq", - "QuantumOptics", - "Plots", - "LaTeXStrings", - ], - mode = PKGMODE_MANIFEST, -) +# # Input-Output Analysis of Quantum Dot SUPER Excitation +# +# This example analyzes the SUPER excitation scheme for quantum dots with the input-output formalism [J. Kerber et al., TODO](https://arxiv.org/). Two red-detuned pulses allow for a close to 100% excitation of a two-level quantum emitter. At the microscopic level, the SUPER mechanism exhibits its nonlinear three-photon Raman-type character, leading to a net photon-number change of −2 in one mode and +1 in the other. + +# In the first part we describe the dynamics within a cumulant expansion approach for coherent light. We then transform into the interaction-picture of the input and output cavities, which allows us to describe the interaction with large Fock states. + +# We start by loading the needed packages and specifying the model. + +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumCumulants +using ModelingToolkitBase +using OrdinaryDiffEq +using QuantumOptics +using Plots +using LaTeXStrings +sol_values(sol, op, eqs) = get_solution(sol, op, eqs).(sol.t) + +## Hilbert spaces +hu2 = FockSpace(:u2) # virtual input cavity (u2) +hu1 = FockSpace(:u1) # virtual input cavity (u1) +hs1 = NLevelSpace(:atom, 2) # TLS +hv1 = FockSpace(:v1) # virtual output cavity (v1) +hv2 = FockSpace(:v2) # virtual output cavity (v2) +h = tensor(hu2, hu1, hs1, hv1, hv2) + +## Operators +au2 = Destroy(h, :a_u2, 1) +au1 = Destroy(h, :a_u1, 2) +s(i, j) = Transition(h, :s, i, j) +av1 = Destroy(h, :a_v1, 4) +av2 = Destroy(h, :a_v2, 5) + +## Symbolic parameters: decay rate and virtual cavity couplings +@variables γ::Real gu1::Number gu2::Number gv1::Number gv2::Number +@independent_variables t # Symbolic time variable + +## SLH triplets +G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 +G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 +G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system +G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 +G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 + +## cascade SLH triplets +G_cas = ▷(G_u2, G_u1, G_2lvl, G_v1, G_v2) +nothing # hide + +# + +## Hamiltonian and Lindbladian +Hcas = hamiltonian(G_cas) +Jcas = jump_operator(G_cas)[1] +Jcasd = adjoint(Jcas) +nothing # hide + +# To deal with time-dependent functions in QuantumCumulants, we need to register them. Furthermore, due to a problem for the conjugate of registered functions (conj is ignored), we first need to create the adjoint of the jump operators and then substitute the time-dependent functions. + +## Time-dependent couplings +@register_symbolic gu1_t(t) +@register_symbolic gu2_t(t) +@register_symbolic gv1_t(t) +@register_symbolic gv2_t(t) + +g_ls = [gu2, gu1, gv1, gv2] +gt_ls = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] +dict_gt = Dict(g_ls .=> gt_ls) + +## Insert time-dependence +Hcas_t = substitute(Hcas, dict_gt) +Jcas_t = substitute(Jcas, dict_gt) +Jcasd_t = substitute(Jcasd, dict_gt) +nothing # hide + +# We calculate the coupling for the input and output cavities. The modified couplings of the second input and output modes are obtained with the function [`effective_input_mode`](@ref) and [`effective_output_mode`](@ref), respectively. Due to the fast oscillations, the tolerance of the numeric solver needs to be improved. The parameters are taken from [T. K. Bracht et al., PRX Quantum 2, 040354 (2021)](https://doi.org/10.1103/PRXQuantum.2.040354). + +## Time grid +dt = 1e-4 +Tend = 12*2 +T = [dt:dt:Tend;] + +## Numerical pulse parameters in THz and ps +γ_ = 1e-2 # TLS decay rate +Δ1_ = -2π*1.934 # detuning pulse 1 +Δ2_ = -2π*4.634 # detuning pulse 2 +Α1_ = 22.65π # pulse area 1 +Α2_ = 19.29π # pulse area 2 +σ1_ = 2.4 # temporal FWHM pulse 1 +σ2_ = 3.04 # temporal FWHM pulse 2 +τ1_ = 0.0 + 12 # time shift pulse 1 +τ2_ = -0.73 + 12 # time shift pulse 2 + +## Normalized input modes +u1(t_) = 1/(√(σ1_)*π^(1/4)) * exp(-(t_ - τ1_)^2 / (2*σ1_^2)) * exp(-1im*Δ1_*t_) +u2(t_) = 1/(√(σ2_)*π^(1/4)) * exp(-(t_ - τ2_)^2 / (2*σ2_^2)) * exp(-1im*Δ2_*t_) + +## Coupling functions +gu1_t_ = coupling_input(u1, T) +gu1_t(t) = gu1_t_(t) +gv1_t_ = coupling_output(u1, T) +gv1_t(t) = gv1_t_(t) + +## Cascade-modified couplings (effective modes) +abstol = 1e-10 +reltol = 1e-10 +u_fcts = [u1, u2] +u2_eff = effective_input_mode(u_fcts, T, 2; abstol, reltol) +gu2_t_ = coupling_input(u2_eff, T) +gu2_t(t) = gu2_t_(t) + +v_fcts = [u1, u2] # output modes = input modes +v2_eff = effective_output_mode(v_fcts, T, 2; abstol, reltol) +gv2_t_ = coupling_output(v2_eff, T) +gv2_t(t) = gv2_t_(t) +nothing # hide + +# After deriving the mean-field equations we define the initial state, create the ODE problem and solve the dynamics. + +## First-order cumulant expansion +order = 1 +ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] +eqs = meanfield(ops, Hcas_t, [Jcas_t]; Jdagger = [Jcasd_t], order = order, iv = t) + +## Coherent-state amplitudes +α1 = Α1_ / (2*√(2)*π^(1/4)*√(σ1_*γ_)) # field 1 +α2 = Α2_ / (2*√(2)*π^(1/4)*√(σ2_*γ_)) # field 2 +u0 = [α1, α2, 0, 0, 0, 0.0im] + +## Solve ODE system +sys = mtkcompile(System(eqs; name = :sys)) +u0_p_map_cas = Dict([unknowns(sys); γ] .=> [u0; γ_]) +prob_cas = ODEProblem(sys, u0_p_map_cas, (dt, Tend)) +sol = solve(prob_cas, Tsit5(); abstol, reltol) +nothing # hide + +## Expectation values +t_cas = sol.t # time vector +s22_cas = sol_values(sol, s(2, 2), eqs) +nu1_cas = abs2.(sol_values(sol, au1, eqs)) +nu2_cas = abs2.(sol_values(sol, au2, eqs)) +nv1_cas = abs2.(sol_values(sol, av1, eqs)) +nv2_cas = abs2.(sol_values(sol, av2, eqs)) +nothing # hide + +# + +common = (; + xlims = (t_cas[1]-0.01, t_cas[end]), + tickfontsize = 18, + guidefontsize = 18, + legendfontsize = 18, +) +p1 = plot(t_cas, real.(s22_cas); color = :red, label = L"\mathrm{cascade}") +p2 = plot( + t_cas, + nu1_cas; + color = :blue, + label = L"\langle \hat{n}_{u_1} \rangle", + ylabel = L"\langle\hat{n}_i\rangle", + xticks = ([0, 10, 20], ["", "", ""]), + yticks = ([0, 5e3, 10e3, 15e3], [L"0", L"5\cdot10^3", L"10\cdot10^3", L"15\cdot10^3"]), + ylims = (nv1_cas[1]-0.5e3, nu1_cas[1]+0.5e3), + legend = :left, + common..., +) +plot!(p2, t_cas, nv1_cas; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") +plot!( + p2, + t_cas, + nu2_cas; + color = :blue, + ls = :dash, + label = L"\langle \hat{n}_{u_2} \rangle", +) +plot!( + p2, + t_cas, + nv2_cas; + color = :red, + ls = :dash, + label = L"\langle \hat{n}_{v_2} \rangle", +) +p3 = plot( + t_cas, + nu1_cas .+ nv1_cas .- nu1_cas[1]; + color = :blue, + label = L"\mathrm{mode~1}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + xticks = ([0, 10, 20], latexstring.([0, 10, 20])), + yticks = ([-4, -2, 0, 2, 4], latexstring.([-4, -2, 0, 2, 4])), + legend = :topleft, + common..., +) +plot!(p3, t_cas, nu2_cas .+ nv2_cas .- nu2_cas[1]; color = :red, label = L"\mathrm{mode~2}") +pl1 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) +display(pl1) + +# ## Interaction picture + +## In the following, we will transform into the interaction picture of the virtual cavities. + +## Interaction picture: cavity dynamics +H_uv = hamiltonian(▷(G_u2, G_u1, G_v1, G_v2)) +H_int_ = simplify(Hcas - H_uv) + +M(i, j) = Symbolics.variable(Symbol("M_{$(i)$(j)}"); T = Number) +a0_ls = [au2, au1, av1, av2] +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 .=> a_int_ls) + +H_int = substitute(H_int_, int_dict) +J_int = simplify(substitute(Jcas, int_dict)) +Jd_int = simplify(substitute(Jcasd, int_dict)) +nothing # hide + +## Coefficient matrix M +Mat = Matrix{Any}(undef, la, la) +mod = @__MODULE__ # hide +for i = 1:la, j = 1:la + name = Symbol("Ma_$(i)$(j)") + @eval @register_symbolic $name(t) + Mat[i, j] = getfield(mod, name)(t) +end +Mat_ls = [Mat[i, j] for i = 1:la for j = 1:la] +M_ls = [M(i, j) for i = 1:la for j = 1:la] + +## Time-evolution of the matrix M(t) +A_uv = coupling_matrix((gu2_t_, gu1_t_, gv1_t_, gv2_t_)) +M_t = solve_mode_evolution(A_uv, T) + +for i = 1:la, j = 1:la + fname = Symbol("Ma_$(i)$(j)") + @eval begin + $fname(t) = M_t(t)[$i, $j] + end +end +dict_Mt = Dict(M_ls .=> Mat_ls) +dict_gt_Mt = merge(dict_gt, dict_Mt) + +H_int_t = substitute(H_int, dict_gt_Mt) +J_int_t = substitute(J_int, dict_gt_Mt) +Jd_int_t = substitute(Jd_int, dict_gt_Mt) + +eqs_int = meanfield(ops, H_int_t, [J_int_t]; Jdagger = [Jd_int_t], order = order, iv = t); + +## Solve ODE system in interaction picture +sys_int = mtkcompile(System(eqs_int; name = :sysI)) +u0_p_map_int = Dict([unknowns(sys_int); γ] .=> [u0; γ_]) +prob_int = ODEProblem(sys_int, u0_p_map_int, (dt, Tend)) +sol_int = solve(prob_int, Tsit5(); abstol, reltol) +nothing # hide + +## Expectation values +t_int = sol_int.t +s22_int = real.(sol_values(sol_int, s(2, 2), eqs_int)) +nu1_int = abs2.(sol_values(sol_int, au1, eqs_int)) +nu2_int = abs2.(sol_values(sol_int, au2, eqs_int)) +nv1_int = abs2.(sol_values(sol_int, av1, eqs_int)) +nv2_int = abs2.(sol_values(sol_int, av2, eqs_int)) + +# + +pl4 = plot( + t_int, + nu1_int .- nu1_int[1]; + color = :blue, + label = L"\mathrm{mode~1~(int.)}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + xlims = (t_int[1]-0.01, t_int[end]), + yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), + legend = :right, + tickfontsize = 18, + guidefontsize = 18, + legendfontsize = 18, + size = (800, 400), +) +plot!(pl4, t_int, nu2_int .- nu2_int[1]; color = :red, label = L"\mathrm{mode~2~(int.)}") +display(pl4) + +# ## Fock state input + +# Let us now compare the dynamics for coherent input pulses with the case of incident non-classical photon number eigenstates (Fock states), where we choose states with the same mean photon numbers as for the coherent pulses. Since in the atom-field interaction only a few photons are exchanged, the quantum states of the excitation pulses are only changed by a couple of photons, we only need to keep a couple of nearby Fock states in the computational basis. + +# We define the basis of the system and create the dictionary for the time-dependent variables to translate the Hamiltonian and Lindblad operator to a QuantumOptics.jl operator. + +n1_fock = round(Int, abs2(α1)) +n2_fock = round(Int, abs2(α2)) + +bu1 = FockBasis(n1_fock+2, n1_fock-6) +bu2 = FockBasis(n2_fock+3, n2_fock-3) +bs1 = NLevelBasis(2) +bv1 = FockBasis(1) +bv2 = FockBasis(1) +b = tensor([bu2, bu1, bs1, bv1, bv2]...) + +g_t_ls = [gu2_t_, gu1_t_, gv1_t_, gv2_t_] +M_t_ls = [t -> M_t(t)[i, j] for i = 1:la for j = 1:la] +dict_fock = Dict([g_ls; M_ls] .=> [g_t_ls; M_t_ls]) + +H_int_fock = to_numeric(H_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) +J_int_fock = to_numeric(J_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) + +## To solve the dynamics, we create the time-dependent function for the open quantum system and define the initial state. + +function input_output(t, ρ) + Ht = H_int_fock(t) + J = [J_int_fock(t)] + return Ht, J, QuantumOptics.dagger.(J) +end + +# initial state +ψu2 = fockstate(bu2, n2_fock) +ψu1 = fockstate(bu1, n1_fock) +ψs1 = nlevelstate(bs1, 1) +ψv1 = fockstate(bv1, 0) +ψv2 = fockstate(bv2, 0) +ψ0 = tensor(ψu2, ψu1, ψs1, ψv1, ψv2) + +T_fock = [0:0.001:1;]*T[end] +# t_fock, ρt_fock = timeevolution.master_dynamic(T_fock, ψ0, input_output; abstol, reltol) +using Random +Random.seed!(1) # hide +t_fock, ρt_fock = timeevolution.mcwf_dynamic(T_fock, ψ0, input_output; abstol, reltol) + +# Due to the relatively long computation time of timeevolution.master_dynamic, we simulate a single trajectory with timeevolution.mcwf_dynamic. + +## Expectation values +s22_fock = real.(expect(s(2, 2), ρt_fock)) +nu1_fock = real.(expect(au1'au1, ρt_fock)) +nu2_fock = real.(expect(au2'au2, ρt_fock)) +nv1_fock = real.(expect(av1'av1, ρt_fock)) +nv2_fock = real.(expect(av2'av2, ρt_fock)) + +common = (; + xlims = (t_fock[1]-0.01, t_fock[end]), + tickfontsize = 18, + guidefontsize = 18, + legendfontsize = 18, +) +p3_1 = plot(t_int, s22_int; color = :blue, label = L"\mathrm{coherent state}") +plot!(p3_1, t_fock, s22_fock; color = :red, label = L"\mathrm{Fock state}") +p3_2 = plot( + t_fock, + nu1_fock .- nu1_fock[1]; + color = :green, + label = L"\mathrm{Fock: mode~1}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + xlims = (t_fock[1]-0.01, t_fock[end]), + yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), + legend = :right, + common..., +) +plot!( + p3_2, + t_fock, + nu2_fock .- nu2_fock[1]; + color = :yellow, + label = L"\mathrm{Fock: mode~2}", +) +plot!( + p3_2, + t_int, + nu1_int .- nu1_int[1]; + color = :blue, + label = L"\mathrm{Coherent: mode~1}", +) +plot!( + p3_2, + t_int, + nu2_int .- nu2_int[1]; + color = :red, + label = L"\mathrm{Coherent: mode~2}", +) +pl3 = plot(p3_1, p3_2; layout = (2, 1), size = (800, 600)) +display(pl3) + +# Due to the vanishing relative phase of the Fock states, the oscillations disappear. + +# ## Package versions + +# These results were obtained using the following versions: + +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status( + [ + "QuantumInputOutput", + "SecondQuantizedAlgebra", + "QuantumCumulants", + "ModelingToolkitBase", + "OrdinaryDiffEq", + "QuantumOptics", + "Plots", + "LaTeXStrings", + ], + mode = PKGMODE_MANIFEST, +) 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 a11db69b..b8ba2af4 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 @@ -38,7 +38,7 @@ H_s = p / 2 * (a'^2 + a^2) - K / 2 * (a'^2) * (a^2) + Δ * a' * a G_s = SLH(1, √(γ) * a, H_s) H = hamiltonian(G_s) -L = jump_operator(G_s)[1] +J = jump_operator(G_s)[1] nothing # hide # Next, we define the numerical parameters from Sec. II.D. The classical pump @@ -89,11 +89,11 @@ dict_p = Dict([γ, K, Δ] .=> [γ_, K_, Δ_]) dict_p_t = Dict(p => p_t) H_QO = to_numeric(H, b1; parameter = dict_p, time_parameter = dict_p_t) -L_QO = to_numeric(L, b1; parameter = dict_p, time_parameter = dict_p_t) +J_QO = to_numeric(J, b1; parameter = dict_p, time_parameter = dict_p_t) function input_output_1(t, ρ) Ht = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return Ht, J, dagger.(J) end @@ -102,8 +102,8 @@ t_1, ρt_1 = timeevolution.master_dynamic(T, ψ0, input_output_1) nothing # hide a_qo = destroy(bc) -Ls(t) = √(γ_) * a_qo -g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) +Js(t) = √(γ_) * a_qo +g1_m = correlation_matrix(T, ρt_1, input_output_1, Js) F = eigen(g1_m) n_avg = real.(F.values) * ΔT @@ -127,7 +127,7 @@ G_v = SLH(1, g_v * av, 0) G = G_s2 ▷ G_v H_2 = hamiltonian(G) -L_2 = jump_operator(G)[1] +J_2 = jump_operator(G)[1] gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict(p => p_t, g_v => gv_t) @@ -136,11 +136,11 @@ bv = FockBasis(n_cut) b2 = bc ⊗ bv H_QO_2 = to_numeric(H_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) -L_QO_2 = to_numeric(L_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) +J_QO_2 = to_numeric(J_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return Ht, J, dagger.(J) end 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 3aa4d771..3b8079c8 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 @@ -61,12 +61,12 @@ nothing #hide # -L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] -G_atom = SLH(I4, L_atom, Δ*σ(2, 2)) +J_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] +G_atom = SLH(I4, J_atom, Δ*σ(2, 2)) G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom H = hamiltonian(G_u_bs_d1_d2_atom) -L = jump_operator(G_u_bs_d1_d2_atom) +J = jump_operator(G_u_bs_d1_d2_atom) # @@ -97,7 +97,7 @@ int_dict = Dict(a_ls .=> a_int_ls) ## substitute interaction picture operators H_int = simplify(substitute(H_int_, int_dict)) -L_int = [simplify(substitute(L_, int_dict)) for L_ in L] +J_int = [simplify(substitute(J_, int_dict)) for J_ in J] ## Pulse parameters n = 5#9 @@ -150,12 +150,12 @@ b = bu ⊗ bd1 ⊗ bd2 ⊗ bs dict_p_Δ(Δn) = Dict([γ, Δ, r, t] .=> [γ_, Δn, rn, tn]) H_QO_Δ(Δn) = to_numeric(H_int, b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) -L_QO_Δ(Δn) = [ - to_numeric(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for - i = 1:length(L) +J_QO_Δ(Δn) = [ + to_numeric(J_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for + i = 1:length(J) ] # H_QO_Δ(Δn) = to_numeric(H, b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) -# L_QO_Δ(Δn) = [to_numeric(L[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(L)] +# J_QO_Δ(Δn) = [to_numeric(J[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(J)] nothing # hide # @@ -170,13 +170,13 @@ ad2_QO = to_numeric(ad2, b) Δn = 0 H_QO = H_QO_Δ(Δn) -L_QO = L_QO_Δ(Δn) +J_QO = J_QO_Δ(Δn) input_output = (t, ρ) -> ( H_QO(t), - [L_QO[i](t) for i = 1:length(L_QO)], - [dagger(L_QO[i](t)) for i = 1:length(L_QO)], + [J_QO[i](t) for i = 1:length(J_QO)], + [dagger(J_QO[i](t)) for i = 1:length(J_QO)], ) t_, ρt = timeevolution.master_dynamic(T, ψ0_fock, input_output) @@ -203,13 +203,13 @@ for (it, Δn) in enumerate(Δ_ls) Δn = Δ_ls[it] H_QO = H_QO_Δ(Δn) - L_QO = L_QO_Δ(Δn) + J_QO = J_QO_Δ(Δn) input_output = (t, ρ) -> ( H_QO(t), - [L_QO[i](t) for i = 1:length(L_QO)], - [dagger(L_QO[i](t)) for i = 1:length(L_QO)], + [J_QO[i](t) for i = 1:length(J_QO)], + [dagger(J_QO[i](t)) for i = 1:length(J_QO)], ) t_, ρt = timeevolution.master_dynamic([T[1], T[end]], ψ0_fock, input_output) diff --git a/examples/drafts/10-1_SUPER_excitation.jl b/examples/drafts/10-1_SUPER_excitation.jl index 4f764184..be8ac5c5 100644 --- a/examples/drafts/10-1_SUPER_excitation.jl +++ b/examples/drafts/10-1_SUPER_excitation.jl @@ -1,1435 +1,1435 @@ -### #TODO: - -# # TODO: same in other file: delete conj fcts -# # TODO: delete nv term in other file! not correct in Δn! - - -# # Input-Output Analysis of Quantum Dot SUPER Excitation -# -# This script collects and demonstrates several approaches used in [PAPER REF]. -# -# We provide a quantum pulse description of the SUPER excitation scheme [T. K. Bracht et al., PRX Quantum 2, 040354 (2021)] (https://doi.org/10.1103/PRXQuantum.2.040354) -# and show that SUPER is a three-photon process [Q. W. Richter et al., Physical Review Research 7, 013079 (2025)] (https://doi.org/10.1103/PhysRevResearch.7.013079), -# [L. Vannucci and N. Gregersen, Optics Express Vol. 32, Issue 20, pp. 35381-35394 (2024)] (https://doi.org/10.1364/OE.533998). -# Due to large photon numbers, we use QuantumCumulants.jl [D. Plankensteiner et al., Quantum 6, 617 (2022)] (https://doi.org/10.22331/q-2022-01-04-617). -# -# We also propose full master-equation treatments for large-α coherent states in [PAPER REF]; here we show the “Displacement Frame” approach. -# -# Additionally we demonstrate how suitable output modes can be determined utilizing a correlation-based approach -# -# 1) Two-Input-Two-Output (2I-2O) SUPER Cascade -# -# Two Gaussian coherent-state pulses are emitted by two virtual input cavities, drive a TLS, and are collected by matched output cavities (v_i(t) = u_i(t)). -# The difference between input and output photon numbers per pulse reveals the three-photon nature. -# -# 2) Two-Input-Two-Output (2I-2O) SUPER Concatenation -# -# The same physics can be captured by concatenating channels: for each pulse i, input cavity i -> TLS -> output cavity i. -# The SLH description has two inputs and two outputs; photon and system dynamics recover the multi-photon character. -# -# 3) Interaction Picture and Displacement Frame of 2I-2O Cascade -# -# We transform the 2I-2O cascade into the cavity–cavity interaction picture [V. R. Christiansen et al., Physical Review A 107, 013706 (2023)] (https://doi.org/10.1103/PhysRevA.107.013706). -# We then enter a displacement frame that splits coherent backgrounds from quantum fluctuations, reducing the Fock spaces and enabling full master-equation treatment. -# -# 4) Input-Output SUPER Cascade -# -# We find a single coherent mode representing the SUPER superposition. This compressed description hides the explicit multi-photon nature, but greatly speeds up the QuantumCumulants.jl-based evaluation of the two-time correlation matrix g^(1)(t_1, t_2). -# Diagonalizing g^(1) yields the dominant output modes of the Input–TLS SLH triplet and of the full I–O SUPER cascade. -# # TODO: re-order - -using QuantumInputOutput -using SecondQuantizedAlgebra -using QuantumCumulants -using ModelingToolkitBase -using OrdinaryDiffEq -using LinearAlgebra -using QuantumOptics -using NumericalIntegration -using JLD2 -using Plots -using LaTeXStrings - -abstol = 1e-10 -reltol = 1e-10 -cd(@__DIR__) -solution_values(sol, op, eqs) = get_solution(sol, op, eqs).(sol.t) -# parameter_map(sys, ps, vals) = -# Dict(p => v for (p, v) in zip(ps, vals) if any(isequal(p), parameters(sys))) - -################################################### -############# 1) 2I-2O SUPER Cascade ############## -################################################### - -# Hilbert spaces -hu2 = FockSpace(:u2) # virtual input cavity (u2) -hu1 = FockSpace(:u1) # virtual input cavity (u1) -hs1 = NLevelSpace(:atom, 2) # TLS -hv1 = FockSpace(:v1) # virtual output cavity (v1) -hv2 = FockSpace(:v2) # virtual output cavity (v2) -h = tensor(hu2, hu1, hs1, hv1, hv2) - -# Operators -au2 = Destroy(h, :a_u2, 1) -au1 = Destroy(h, :a_u1, 2) -s(i, j) = Transition(h, :s, i, j) -av1 = Destroy(h, :a_v1, 4) -av2 = Destroy(h, :a_v2, 5) - -# Symbolic parameters: decay rate and virtual cavity couplings -@variables γ::Real gu1::Number gu2::Number gv1::Number gv2::Number -@independent_variables t # Symbolic time variable - -################################################### -############# SUPER Pulse Parameters ############## -################################################### - -# Time grid -dt = 1e-4; -Tend = 12*2 -T = [dt:dt:Tend;]; -ΔT = T[2] - T[1] - -# Numerical pulse parameters -γ_ = 1e-2 # TLS decay rate -Δ1_ = -2π*1.934 # detuning pulse 1 -Δ2_ = -2π*4.634 # detuning pulse 2 -α1_ = 22.65π # pulse area 1 -α2_ = 19.29π # pulse area 2 -σ1_ = 2.4 # temporal FWHM pulse 1 -σ2_ = 3.04 # temporal FWHM pulse 2 -τ1_ = 0.0 + 12 # time shift pulse 1 -τ2_ = -0.73 + 12 # time shift pulse 2 - -# Classical Gaussian pulses -Ω1(t) = 0.5*1/(√(2π))*α1_/σ1_*exp(-(t-τ1_)^2/(2*σ1_^2))*exp(-1im*Δ1_*t) -Ω2(t) = 0.5*1/(√(2π))*α2_/σ2_*exp(-(t-τ2_)^2/(2*σ2_^2))*exp(-1im*Δ2_*t) - -################################################### -##### Classical SUPER Scheme (reference) ########## -################################################### - -@register_symbolic Ω(t) -@register_symbolic Ω_c(t) - -Ω(t) = Ω1(t) + Ω2(t) -Ω_c(t) = conj(Ω(t)) - -H_cl = -(Ω_c(t)*s(1, 2) + Ω(t)*s(2, 1)) -J_cl = [s(1, 2)]; -rates_cl = [γ] - -ops_cl = [s(2, 2), s(1, 2)] -eqs_cl = meanfield(ops_cl, H_cl, J_cl; rates = rates_cl, iv = t) -sys_cl = mtkcompile(System(eqs_cl; name = :sys_cl)) - -u0_cl = zeros(ComplexF64, length(eqs_cl)) -u0_p_map_cl = merge(Dict(unknowns(sys_cl) .=> u0_cl), Dict(γ => γ_)) -prob_cl = ODEProblem(sys_cl, u0_p_map_cl, (dt, Tend)) -sol_cl = solve(prob_cl, Tsit5(); abstol, reltol) -t_cl = sol_cl.t -s22_cl = real.(solution_values(sol_cl, s(2, 2), eqs_cl)); - -using PyPlot; -pygui(true) - -close(123) -figure(123) -subplot(211) -PyPlot.plot(t_cl, abs.(Ω.(t_cl))) -PyPlot.grid(true) -subplot(212) -PyPlot.plot(t_cl, s22_cl) -PyPlot.grid(true) - -### effective: 3photon process### - -Ω1_max = Α1_/σ1_/√(2π) -Ω2_max = Α2_/σ2_/√(2π) - -Ω_3ph = Ω1_max^2*Ω2_max/(4*Δ1_*Δ2_) -T_3ph = π/Ω_3ph -0.01/3.99 - -Ω_max = maximum(abs.(Ω1.(T)) + abs.(Ω2.(T))) -Ω_eff(t) = 0.5 * Ω_3ph * abs.(abs.(Ω1.(t)) + abs.(Ω2.(t)))/Ω_max - -b = NLevelBasis(2) -s(i, j) = transition(b, i, j) - -H_(t, psi) = Ω_eff(t)*(s(1, 2) + s(2, 1)) - -t_, psi_ = timeevolution.schroedinger_dynamic(T, nlevelstate(b, 1), H_) - -s22_eff = real(expect(s(2, 2), psi_)) -subplot(212) -PyPlot.plot(t_, s22_eff) - - -111 - -################################################### -####### IOT and QuantumCumulants.jl section ####### -################################################### - -# SLH triplets -G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 -G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 -G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system -G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 -G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 - -# 2I-2O cascade SLH triplet -G_cas = ▷(G_u2, G_u1, G_2lvl, G_v1, G_v2) # cascade - -# Hamiltonian and Lindbladian -Hcas = hamiltonian(G_cas) -Lcas = jump_operator(G_cas)[1] -Lcasd = adjoint(Lcas) - -# Time-dependent couplings -@register_symbolic gu1_t(t) -@register_symbolic gu2_t(t) -@register_symbolic gv1_t(t) -@register_symbolic gv2_t(t) -# Complex conjugates -@register_symbolic gu1_c_t(t) -@register_symbolic gu2_c_t(t) -@register_symbolic gv1_c_t(t) -@register_symbolic gv2_c_t(t) - -g_ls = [gu2, gu1, gv1, gv2] -gt_ls = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] -gct_ls = [gu2_c_t(t), gu1_c_t(t), gv1_c_t(t), gv2_c_t(t)] -dict_gt = Dict([g_ls; conj.(g_ls)] .=> [gt_ls; gct_ls]); - -# Insert time-dependence -Hcas_t = substitute(Hcas, dict_gt); -Lcas_t = substitute(Lcas, dict_gt); -Lcasd_t = substitute(Lcasd, dict_gt); - -# Normalized input modes -u1(t_) = 1/(√(σ1_)*π^(1/4)) * exp(-(t_ - τ1_)^2 / (2*σ1_^2)) * exp(-1im*Δ1_*t_) -u2(t_) = 1/(√(σ2_)*π^(1/4)) * exp(-(t_ - τ2_)^2 / (2*σ2_^2)) * exp(-1im*Δ2_*t_) - -# Coupling functions -gu1_t_ = coupling_input(u1, T) -gu1_t(t) = gu1_t_(t) -gv1_t_ = coupling_output(u1, T) -gv1_t(t) = gv1_t_(t) -gu1_c_t(t) = conj(gu1_t(t)) -gv1_c_t(t) = conj(gv1_t(t)) - -# Cascade-modified couplings (effective modes) -u_fcts = [u1, u2] -u2_eff = effective_input_mode(u_fcts, T, 2; abstol, reltol) -gu2_t_ = coupling_input(u2_eff, T) -gu2_t(t) = gu2_t_(t) -gu2_c_t(t) = conj(gu2_t(t)); - -v_fcts = [u1, u2] # output modes = input modes -v2_eff = effective_output_mode(v_fcts, T, 2; abstol, reltol) -gv2_t_ = coupling_output(v2_eff, T) -gv2_t(t) = gv2_t_(t) -gv2_c_t(t) = conj(gv2_t(t)); - -# First-order cumulant expansion -order = 1 -ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] -eqs = meanfield(ops, Hcas_t, [Lcas_t]; Jdagger = [Lcasd_t], order = order, iv = t) - -# Classical-to-coherent amplitude relation -α1_io = α1_ / (2*√(2)*π^(1/4)*√(σ1_*γ_)) # field 1 -α2_io = α2_ / (2*√(2)*π^(1/4)*√(σ2_*γ_)) # field 2 -u0 = [α1_io, α2_io, 0, 0, 0, 0.0im] - -# Solve ODE system -sys = mtkcompile(System(eqs; name = :sys)) -u0_p_map_cas = Dict([unknowns(sys); γ] .=> [u0; γ_]) -prob_cas = ODEProblem(sys, u0_p_map_cas, (dt, Tend)) -sol = solve(prob_cas, Tsit5(); abstol, reltol); - -# Expectation values -t_cas = sol.t # time vector -s22_cas = solution_values(sol, s(2, 2), eqs); -nu1_cas = abs2.(solution_values(sol, au1, eqs)) -nu2_cas = abs2.(solution_values(sol, au2, eqs)) -nv1_cas = abs2.(solution_values(sol, av1, eqs)) -nv2_cas = abs2.(solution_values(sol, av2, eqs)) - -# Plots -plot_font = 18 -common = (; - xlims = (t_cas[1]-0.01, t_cas[end]), - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, -) -p1 = plot( - t_cl, - s22_cl; - color = :blue, - label = L"\mathrm{classical}", - ylabel = L"P", - xticks = ([0, 10, 20], ["", "", ""]), - yticks = ([0, 0.5, 1], [L"0.0", L"0.5", L"1.0"]), - legend = :topleft, - common..., -) -plot!(p1, t_cas, real.(s22_cas); color = :red, ls = :dash, label = L"\mathrm{cascade}") -p2 = plot( - t_cas, - nu1_cas; - color = :blue, - label = L"\langle \hat{n}_{u_1} \rangle", - ylabel = L"\langle\hat{n}_i\rangle", - xticks = ([0, 10, 20], ["", "", ""]), - yticks = ([0, 5e3, 10e3, 15e3], [L"0", L"5\cdot10^3", L"10\cdot10^3", L"15\cdot10^3"]), - ylims = (nv1_cas[1]-0.5e3, nu1_cas[1]+0.5e3), - legend = :left, - common..., -) -plot!(p2, t_cas, nv1_cas; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") -plot!( - p2, - t_cas, - nu2_cas; - color = :blue, - ls = :dash, - label = L"\langle \hat{n}_{u_2} \rangle", -) -plot!( - p2, - t_cas, - nv2_cas; - color = :red, - ls = :dash, - label = L"\langle \hat{n}_{v_2} \rangle", -) -p3 = plot( - t_cas, - nu1_cas .+ nv1_cas .- nu1_cas[1]; - color = :blue, - label = L"\mathrm{mode~1}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - xticks = ([0, 10, 20], latexstring.([0, 10, 20])), - yticks = ([-4, -2, 0, 2, 4], latexstring.([-4, -2, 0, 2, 4])), - legend = :topleft, - common..., -) -plot!(p3, t_cas, nu2_cas .+ nv2_cas .- nu2_cas[1]; color = :red, label = L"\mathrm{mode~2}") -pl1 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) -display(pl1) -# savefig(pl1, "plot_2I2O_cascade.pdf") - - -#################################################################### -##### 3) 2I-2O SUPER Cascade interaction picture and Displacement ## -#################################################################### - -# Interaction picture: cavity dynamics -H_uv = hamiltonian(▷(G_u2, G_u1, G_v1, G_v2)) -H_int_sym_ = simplify(Hcas - H_uv) - -M(i, j) = Symbolics.variable(Symbol("M_{$(i)$(j)}"); T = Number) -Mc(i, j) = Symbolics.variable(Symbol("Mc_{$(i)$(j)}"); T = Number) - -a0_ls = [au2, au1, av1, av2] -la = length(a0_ls) -a_int_ls = [sum(M(i, j)*a0_ls[j] for j = 1:la) for i = 1:la] -# a_c_int_ls = [sum(Mc(i,j)*a0_ls[j]' for j = 1:la) for i = 1:la] # # TODO -a_c_int_ls = [sum(conj(M(i, j))*a0_ls[j]' for j = 1:la) for i = 1:la] # # TODO - -int_dict = Dict([a0_ls; adjoint.(a0_ls)] .=> [a_int_ls; a_c_int_ls]) - -H_int_sym = substitute(H_int_sym_, int_dict) -L_int_sym = simplify(substitute(Lcas, int_dict)) -Ld_int_sym = simplify(substitute(Lcasd, int_dict)) - -# M-matrix -Mat = Matrix{Any}(undef, la, la) -Matc = Matrix{Any}(undef, la, la) -mod = @__MODULE__ -for i = 1:la, j = 1:la - name = Symbol("Ma_$(i)$(j)") - namec = Symbol("Mac_$(i)$(j)") - - @eval @register_symbolic $name(t) - @eval @register_symbolic $namec(t) - - Mat[i, j] = getfield(mod, name)(t) - Matc[i, j] = getfield(mod, namec)(t) -end -Mat_ls = [Mat[i, j] for i = 1:la for j = 1:la] -Mat_conls = [Matc[i, j] for i = 1:la for j = 1:la] - -M_ls = [M(i, j) for i = 1:la for j = 1:la] -Mc_ls = [Mc(i, j) for i = 1:la for j = 1:la] - -# Time-evolution matrix M(t) -A_uv = coupling_matrix((gu2_t_, gu1_t_, gv1_t_, gv2_t_)) -M_t = solve_mode_evolution(A_uv, T) -Mc_t(t) = adjoint(M_t(t)) - -for i = 1:la, j = 1:la - fname = Symbol("Ma_$(i)$(j)") - fnamec = Symbol("Mac_$(i)$(j)") - - @eval begin - $fname(t) = M_t(t)[$i, $j] - $fnamec(t) = Mc_t(t)[$i, $j] - end -end - -dict_Mt = Dict([M_ls..., Mc_ls...] .=> [Mat_ls..., Mat_conls...]) -dict_gt_Mt = merge(dict_gt, dict_Mt) - -H_int_sym_t = substitute(H_int_sym, dict_gt_Mt) -L_int_sym_t = substitute(L_int_sym, dict_gt_Mt) -Ld_int_sym_t = substitute(Ld_int_sym, dict_gt_Mt) - -eqs_int = meanfield( - ops, - H_int_sym_t, - [L_int_sym_t]; - Jdagger = [Ld_int_sym_t], - order = order, - iv = t, -); - -# Solve ODE system in interaction picture -sys_int = mtkcompile(System(eqs_int; name = :sysI)) -u0_p_map_int = Dict([unknowns(sys_int); γ] .=> [u0; γ_]) -prob_int = ODEProblem(sys_int, u0_p_map_int, (dt, Tend)) -sol_int = solve(prob_int, Tsit5(); abstol, reltol); - -# Expectation values -t_int = sol_int.t -s22_int = solution_values(sol_int, s(2, 2), eqs_int); -nu1_int = abs2.(solution_values(sol_int, au1, eqs_int)) -nu2_int = abs2.(solution_values(sol_int, au2, eqs_int)) -nv1_int = abs2.(solution_values(sol_int, av1, eqs_int)) -nv2_int = abs2.(solution_values(sol_int, av2, eqs_int)) - -### Displacement frame ### - -# QuantumOptics.jl Hilbert space for fluctuations # TODO: Why is a cutoff of 1 enough? -bu1 = FockBasis(1) # vacuum fluctuation mode 1 -bu2 = FockBasis(1) # vacuum fluctuation mode 2 -bs1 = NLevelBasis(2) # TLS -bv1 = FockBasis(1) # vacuum mode 1 -bv2 = FockBasis(1) # vacuum mode 2 -b = tensor([bu2, bu1, bs1, bv1, bv2]...) - -ψu1 = fockstate(bu1, 0) -ψu2 = fockstate(bu2, 0) -ψs1 = nlevelstate(bs1, 1) -ψv1 = fockstate(bv1, 0) -ψv2 = fockstate(bv2, 0) -ψ0 = tensor(ψu2, ψu1, ψs1, ψv1, ψv2) - -# Operators -anu2 = destroy(bu2) -anu1 = destroy(bu1) -anv1 = destroy(bv1) -anv2 = destroy(bv2) -s12 = embed(b, 3, transition(bs1, 1, 2)) -s22 = embed(b, 3, transition(bs1, 2, 2)) - -# Mode operator vector -avn = [embed(b, 1, anu2), embed(b, 2, anu1), embed(b, 4, anv1), embed(b, 5, anv2)] - -# Coherent amplitude shifts (input modes only) -αvec = [α2_io, α1_io, 0.0, 0.0] - -# Displacement-frame Hamiltonian and jumps -function Ham_displaced(t, ρ) - Mall = M_t(t) - gall = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] - - # Modified coupling vectors for H and L terms - gn = [gall[1], gall[2], -gall[3], -gall[4]] - gnc = conj.(gn) - gl = [gall[1], gall[2], gall[3], gall[4]] - glc = conj.(gl) - - # Subspaces and local annihilation operators - B = [bu2, bu1, bv1, bv2] - Av = [anu2, anu1, anv1, anv2] - - # Displaced operators: a -> b + α (b = fluctuation mode) - avn_disp = [] - for k = 1:4 - if k < 3 - push!(avn_disp, embed(b, k, (Av[k] + αvec[k]*one(B[k])))) - else - push!(avn_disp, embed(b, k+1, (Av[k] + αvec[k]*one(B[k])))) - end - end - - # Interaction pieces (IP + displacement) - A_gn = transpose(gn) * Mall * avn_disp - A_gnc = avn_disp' * Mall' * gnc - - A_gl = transpose(gl) * Mall * avn_disp - A_glc = avn_disp' * Mall' * glc - - # Displacement Hamiltonian - Ht = sqrt(γ_) * 1im/2 * (A_gn * s12' - s12 * A_gnc) - - # Jumps (displacement Lindbladian) - J = [sqrt(γ_) * s12 + A_gl] - Jd = [sqrt(γ_) * s12' + A_glc] - - return Ht, J, Jd -end - -# Master equation in displacement frame -t_df, ρt_df = timeevolution.master_dynamic(T, ψ0, Ham_displaced; abstol, reltol); - -# Expectation values -s22_df = real.(expect(s22, ρt_df)) -nu2_df = real.( - expect( - embed( - b, - 1, - (anu2 + α2_io*identityoperator(bu2))'*(anu2 + α2_io*identityoperator(bu2)), - ), - ρt_df, - ), -) -nu1_df = real.( - expect( - embed( - b, - 2, - (anu1 + α1_io*identityoperator(bu1))'*(anu1 + α1_io*identityoperator(bu1)), - ), - ρt_df, - ), -) -nv1_df = real.(expect(embed(b, 4, anv1'*anv1), ρt_df)) -nv2_df = real.(expect(embed(b, 5, anv2'*anv2), ρt_df)) - - -# Plot net loss/gain in interaction picture and displacement frames -pl4 = plot( - t_int, - nu1_int .+ nv1_int .- nu1_int[1]; - color = :blue, - label = L"\mathrm{mode~1~(int.)}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - xlims = (t_df[1]-0.01, t_df[end]), - yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), - legend = :right, - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, - size = (800, 400), -) -plot!( - pl4, - t_int, - nu2_int .+ nv2_int .- nu2_int[1]; - color = :red, - label = L"\mathrm{mode~2~(int.)}", -) -plot!( - pl4, - t_df, - nu1_df .+ nv1_df .- nu1_df[1]; - color = :green, - ls = :dash, - label = L"\mathrm{mode~1~(dis.)}", -) -plot!( - pl4, - t_df, - nu2_df .+ nv2_df .- nu2_df[1]; - color = :orange, - ls = :dash, - label = L"\mathrm{mode~2~(dis.)}", -) -display(pl4) -# savefig(pl4, "plot_int_dis.pdf") - - -################################################### -######### 2) 2I-2O SUPER Concatenation ############ -################################################### - -# TODO: No substitution with γ/2! - interpretation: decay rate twice - -# SLH triplets -G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 -G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 -G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system # TODO -G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 -G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 - -# Channel-wise cascades -G_ch1 = ▷(G_u2, G_2lvl, G_v2) # pulse 1 -G_ch2 = ▷(G_u1, G_2lvl, G_v1) # pulse 2 -Gcon = ⊞(G_ch1, G_ch2) # concatenation - -# Hamiltonian and Lindblad operators (two channels) -Hcon = hamiltonian(Gcon) -Lch1g = jump_operator(Gcon)[1] -Lch2g = jump_operator(Gcon)[2] -Ldch1g = adjoint(Lch1g) -Ldch2g = adjoint(Lch2g) - -# Split decay γ -> γ/2 (symmetric into two channels) -# dict_gam = Dict(γ => γ/2); # TODO -dict_gam = Dict(γ => γ); - -# Apply decay adjustment -Lch1 = substitute(Lch1g, dict_gam) # TODO -Lch2 = substitute(Lch2g, dict_gam) -Ldch1 = substitute(Ldch1g, dict_gam) -Ldch2 = substitute(Ldch2g, dict_gam) - -# Insert time-dependence -Hcon_t = substitute(Hcon, dict_gt); -Lch1_t = substitute(Lch1, dict_gt); -Lch2_t = substitute(Lch2, dict_gt); -Ldch1_t = substitute(Ldch1, dict_gt); -Ldch2_t = substitute(Ldch2, dict_gt); - -gu2_t_ = coupling_input(u2, T) -gu2_t(t) = gu2_t_(t) -gv2_t_ = coupling_output(u2, T) -gv2_t(t) = gv2_t_(t) - -gu2_c_t(t) = conj(gu2_t(t)) -gv2_c_t(t) = conj(gv2_t(t)) - -# First-order cumulant expansion -order = 1 -ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] -eqs_con = meanfield( - ops, - Hcon_t, - [Lch1_t, Lch2_t]; - Jdagger = [Ldch1_t, Ldch2_t], - order = order, - iv = t, -) - -# Solve ODE system -sys_con = mtkcompile(System(eqs_con; name = :sys_con)) - -u0_p_map_con = Dict([unknowns(sys_con); γ] .=> [u0; γ_]) # TODO -prob_con = ODEProblem(sys_con, u0_p_map_con, (dt, Tend)) -sol_con = solve(prob_con, Tsit5(); abstol, reltol); - -# Expectation values -t_con = sol_con.t -s22_con = solution_values(sol_con, s(2, 2), eqs_con) -nu1_con = abs2.(solution_values(sol_con, au1, eqs_con)) -nu2_con = abs2.(solution_values(sol_con, au2, eqs_con)) -nv1_con = abs2.(solution_values(sol_con, av1, eqs_con)) -nv2_con = abs2.(solution_values(sol_con, av2, eqs_con)) - - -common = (; - xlims = (t_con[1]-0.01, t_con[end]), - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, -) -p1 = plot( - t_cl, - real.(s22_cl); - color = :blue, - label = L"\mathrm{cascade}", - ylabel = L"P", - xticks = ([0, 10, 20], ["", "", ""]), - legend = :topleft, - common..., -) -plot!(p1, t_con, real.(s22_con); color = :red, ls = :dash, label = L"\mathrm{concat.}") -p2 = plot( - t_con, - nu1_con; - color = :blue, - label = L"\langle \hat{n}_{u_1} \rangle", - ylabel = L"\langle\hat{n}_i\rangle", - xticks = ([0, 10, 20], ["", "", ""]), - ylims = (nv1_con[1]-0.5e3, nu1_con[1]+0.5e3), - legend = :left, - common..., -) -plot!(p2, t_con, nv1_con; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") -plot!( - p2, - t_con, - nu2_con; - color = :blue, - ls = :dash, - label = L"\langle \hat{n}_{u_2} \rangle", -) -plot!( - p2, - t_con, - nv2_con; - color = :red, - ls = :dash, - label = L"\langle \hat{n}_{v_2} \rangle", -) -p3 = plot( - t_con, - nu1_con .+ nv1_con .- nu1_con[1]; - color = :blue, - label = L"\mathrm{mode~1}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - yticks = ([-2, -1, 0, 1, 2], latexstring.([-2, -1, 0, 1, 2])), - legend = :bottomleft, - common..., -) -plot!(p3, t_con, nu2_con .+ nv2_con .- nu2_con[1]; color = :red, label = L"\mathrm{mode~2}") -pl2 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) -display(pl2) -# savefig(pl2, "plot_concat.pdf") - - -################################################### -########## 4) I-O SUPER Cascade ################### -################################################### - -# Single-mode SUPER description -dt = 1e-5 -Tend = 12*2 -T = [dt:dt:Tend;] -ΔT = T[2] - T[1] - -# Hilbert space -hu = FockSpace(:u) # virtual input cavity -hv = FockSpace(:v) # virtual output cavity -hs = tensor(hu, hs1, hv) - -# Operators -aun = Destroy(hs, :a_u, 1) -σ(i, j) = Transition(hs, :σ, i, j) -avn = Destroy(hs, :a_v, 3) - -# Couplings -@variables gu::Number gv::Number - -################################################### -# QuantumCumulants.jl section -################################################### - -# SLH triplets -G_u = SLH(1, gu*aun, 0) # input cavity -G_2lvls = SLH(1, √(γ)*σ(1, 2), 0) # 2-level system -G_v = SLH(1, gv*avn, 0) # output cavity - -# I-O cascade -G_cas_1m = ▷(G_u, G_2lvls, G_v) - -# Hamiltonian and Lindbladian -Hcas_1m = hamiltonian(G_cas_1m) -Lcas_1m = jump_operator(G_cas_1m)[1] -Lcas_1md = adjoint(Lcas_1m) - -# Time-dependent couplings -@register_symbolic gu_t(t) -@register_symbolic gv_t(t) -# Complex conjugates -@register_symbolic gu_c_t(t) -@register_symbolic gv_c_t(t) - -gs_ls = [gu, gv] -gst_ls = [gu_t(t), gv_t(t)] -gsct_ls = [gu_c_t(t), gv_c_t(t)] - -dict_gts = Dict([gs_ls; conj.(gs_ls)] .=> [gst_ls; gsct_ls]); - -# Insert time-dependence -Hcas_1m_t = substitute(Hcas_1m, dict_gts); -Lcas_1m_t = substitute(Lcas_1m, dict_gts); -Lcas_1md_t = substitute(Lcas_1md, dict_gts); - -0.5*1/(√(2π))*α1_/σ1_*exp(-(t-τ1_)^2/(2*σ1_^2))*exp(-1im*Δ1_*t) - -# Pulses and composite mode -Ω1(t_) = α1_/(sqrt(2π*σ1_^2))*exp(-(t_-τ1_)^2/(2σ1_^2))*exp(-1im*Δ1_*t_) -Ω2(t_) = α2_/(sqrt(2π*σ2_^2))*exp(-(t_-τ2_)^2/(2σ2_^2))*exp(-1im*Δ2_*t_) -Ωall(t_) = Ω1(t_) + Ω2(t_) -αeff = sqrt(cumul_integrate(T, abs2.(Ωall.(T)))[end]) -u(t_) = Ωall(t_) / αeff - -gu_t_ = coupling_input(u, T) -gu_t(t) = gu_t_(t) -gv_t_ = coupling_output(u, T) -gv_t(t) = gv_t_(t) - -gu_c_t(t) = conj(gu_t(t)) -gv_c_t(t) = conj(gv_t(t)) - -# Second-order expansion for correlation-based mode extraction -order = 2 -# Minimal operator set to capture second-order dynamics -ops_1m = [aun, aun'aun, aun*aun, σ(1, 2), σ(2, 2), avn, avn'avn] -# ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); # TODO: needed for corr? -# ops_all = union(ops1,ops_all1[1:1],ops_all1[3:14], ops_all1[16:16], ops_all1[18:26], ops_all1[28:29]) - -eqs_1m = - meanfield(ops_1m, Hcas_1m_t, [Lcas_1m_t]; Jdagger = [Lcas_1md_t], order = order, iv = t); -complete!(eqs_1m); -length(eqs_1m) - -# Classical-to-coherent relation -α_1m = αeff/(2*sqrt(γ_)) # may change if cross-terms are significant - -u0_1m = zeros(ComplexF64, length(eqs_1m)) -u0_1m[1] = α_1m -u0_1m[2] = abs2(α_1m) -u0_1m[3] = α_1m*α_1m - -# Solve ODE system -sys_1m = mtkcompile(System(eqs_1m; name = :sys_1m)) -ps = [γ] -p0 = [γ_] -u0_p_maps_1m = Dict([unknowns(sys_1m); γ] .=> [u0_1m; γ_]) -prob_1m = ODEProblem(sys_1m, u0_p_maps_1m, (dt, Tend)) -sol_1m = solve(prob_1m, Tsit5(), saveat = dt; abstol, reltol); - -t_1m = sol_1m.t; -s22_1m = solution_values(sol_1m, σ(2, 2), eqs_1m); -nu_1m = real.(solution_values(sol_1m, aun'aun, eqs_1m)) -nv_1m = real.(solution_values(sol_1m, avn'avn, eqs_1m)) - -common = (; - xlims = (t_1m[1]-0.01, t_1m[end]), - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, -) -p1 = plot( - t_cl, - s22_cl; - color = :blue, - label = L"\mathrm{classical}", - ylabel = L"P", - xticks = ([0, 10, 20], ["", "", ""]), - legend = :topleft, - common..., -) -plot!(p1, t_1m, real.(s22_1m); color = :red, ls = :dash, label = L"\mathrm{single~SUPER}") -p2 = plot( - t_1m, - nu_1m; - color = :blue, - label = L"\langle \hat{n}_{u} \rangle", - ylabel = L"\langle\hat{n}_i\rangle", - xticks = ([0, 10, 20], ["", "", ""]), - ylims = (nv_1m[1]-0.5e3, nu_1m[1]+0.5e3), - legend = :left, - common..., -) -plot!(p2, t_1m, nv_1m; color = :red, label = L"\langle \hat{n}_{v} \rangle") -p3 = plot( - t_1m, - nu_1m .+ nv_1m .- nu_1m[1]; - color = :blue, - label = L"\mathrm{mode}", - ylabel = L"\langle\Delta\hat{n}_i\rangle", - xlabel = L"\gamma t", - yticks = ([-1, -0.5, 0], latexstring.([-1, -0.5, 0])), - legend = :bottomleft, - common..., -) -pl4 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) -display(pl4) -# savefig(pl4, "plot_IO_cascade.pdf") - - -################################################### -########### I-O SUPER Correlation ################# -################################################### - -# Same as above but with sparser saveat (va = 0.05) and only Input–System cascade - -dt = 1e-5 -Tend = 12*2 -T = [dt:dt:Tend;] -ΔT = T[2] - T[1] - -@register_symbolic gu_t(t) -@register_symbolic gv_t(t) -@register_symbolic gu_c_t(t) -@register_symbolic gv_c_t(t) - -G_cas_1m0 = ▷(G_u, G_2lvls, G_v) - -Hcas_1m0 = hamiltonian(G_cas_1m0) -Lcas_1m0 = jump_operator(G_cas_1m0)[1] -Lcas_1m0d = adjoint(Lcas_1m0) - -# Set v-coupling to zero (collect only input-system correlations) -gs0_ls = [gu, gv] -gs0t_ls = [gu_t(t), 0.0*gv_t(t)] -gsc0t_ls = [gu_c_t(t), 0.0*gv_c_t(t)] - -dict_gts0 = Dict([gs0_ls; conj.(gs0_ls)] .=> [gs0t_ls; gsc0t_ls]); - -Hcas_1m0_t = substitute(Hcas_1m0, dict_gts0) -Lcas_1m0_t = substitute(Lcas_1m0, dict_gts0) -Lcas_1m0d_t = substitute(Lcas_1m0d, dict_gts0) - -gu_t_ = coupling_input(u, T) -gu_t(t) = gu_t_(t) -gv_t_ = coupling_output(u, T) -gv_t(t) = 0.0*gv_t_(t) - -gu_c_t(t) = conj(gu_t(t)) -gv_c_t(t) = conj(gv_t(t)) - -# order = 2 -# ops1 = [aun, aun', σ(1,2), σ(2,2),σ(1,2)', avn, avn'] -# ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); -# ops_all = union(ops1,ops_all1[1:1],ops_all1[3:14], ops_all1[16:16], ops_all1[18:26], ops_all1[28:29]) -# eqs_1m0 = meanfield(ops_all, Hcas_1m0_t, [Lcas_1m0_t]; Jdagger=[Lcas_1m0d_t], order=order, iv=t) - -# complete!(eqs_1m0) -# length(eqs_1m0) - -ops_1m0 = [aun, aun'aun, aun*aun, σ(1, 2), σ(2, 2)] -eqs_1m0 = meanfield( - ops_1m0, - Hcas_1m0_t, - [Lcas_1m0_t]; - Jdagger = [Lcas_1m0d_t], - order = order, - iv = t, -); -complete!(eqs_1m0); -length(eqs_1m0) - -u0_1m0 = zeros(ComplexF64, length(eqs_1m0)) -u0_1m0[1] = α_1m -u0_1m0[2] = abs2(α_1m) -u0_1m0[3] = α_1m*α_1m - -va = 0.05 - -# Solve ODE system -sys_1m0 = mtkcompile(System(eqs_1m0; name = :sys_1m0)) -u0_p_maps_1m0 = Dict([unknowns(sys_1m0); γ] .=> [u0_1m0; γ_]) -prob_1m0 = ODEProblem(sys_1m0, u0_p_maps_1m0, (dt, Tend)) -sol_1m0 = solve(prob_1m0, Tsit5(), saveat = va; abstol, reltol); - -################################################### -########### Correlation Part ###################### -################################################### - -# Lindblad vector for correlation functions -Lvec=[aun, σ(1, 2)] - -# Discrete time grid -N = length(sol_1m0.t) - -# g^(1)(t1,t2) container for each operator pair -Gmat_bar = zeros(ComplexF64, N, length(Lvec), length(Lvec), N) -t1 = sol_1m0.t; -t2 = sol_1m0.t; -τ0 = 0.0 - -# All operator correlation combinations -Corrmat = [ - [ - CorrelationFunction(Lvec[k]', Lvec[l], eqs_1m0; steady_state = false, iv0 = t) - for l = 1:length(Lvec) - ] for k = 1:length(Lvec) -] -# Minimal set containing full dynamics -Corrvec = [Corrmat[end][i] for i = 1:length(Lvec)] - -xx = CorrelationFunction(σ(1, 2), σ(1, 2), eqs_1m0; steady_state = false, iv0 = t) -xx.eqs.states - -################################################### -###### QC.jl extended correlation utilities ####### -################################################### - -# Structural equality of two lhs expressions -function same_lhs(a, b) - operation(a) === operation(b) && begin - aa = Tuple(arguments(a)) - bb = Tuple(arguments(b)) - length(aa) == length(bb) && all(isequal.(aa, bb)) - end -end - -# Index mapping between Corrvec and Corrmat -function idx_tuple(Lvec, Corrvec, Corrmat) - idxvec = zeros(Int8, length(Lvec), length(Lvec)) - for l = 1:length(Lvec) - eqs1 = equations(Corrvec[l].de) - for k = 1:length(Lvec) - target_lhs = equations(Corrmat[k][l].de)[1].lhs - idx = findfirst(eq -> same_lhs(eq.lhs, target_lhs), eqs1) - idxvec[k, l] = idx - end - end - idx_tuple = Tuple(ntuple(i -> c[i], length(Lvec)) for c in eachcol(idxvec)) - return idx_tuple -end - -# Operators included -oper_tup = (Lvec[1], Lvec[2]) -# Index tuple for selecting correlation components -idx_tup = idx_tuple(Lvec, Corrvec, Corrmat) - -# Compute full Gmat_bar on the discrete time grid -function compute_Gmat( - sol_1m0, - eqs_1m0, - Lvec, - Corrvec; - p0, - ps, - saveat, - τ0 = 0.0, - alg = OrdinaryDiffEq.Tsit5(), - state_syms = oper_tup, - extract_idxs = idx_tup, -) - - N = length(sol_1m0.t) - Gmat_bar = zeros(ComplexF64, N, length(Lvec), length(Lvec), N) - - padleftN(v::AbstractVector{<:Number}, N::Int) = ( - length(v) < N ? vcat(zeros(ComplexF64, N - length(v)), ComplexF64.(v)) : - ComplexF64.(v) - ) - - for i1 = 1:(N-1) - τend = sol_1m0.t[end] - sol_1m0.t[i1] - z = i1 - - for l = 1:length(Lvec) - corr = Corrvec[l] - complete(corr.de) - - iv = ModelingToolkit.get_iv(corr.de) - de_shift = substitute(corr.de, Dict(t => iv + sol_1m0.t[z])) - - bc_var = complete(corr.de)[end].lhs - rhs_val = sol_1m0[state_syms[l]][z] - corrn_fin = substitute(de_shift, Dict(bc_var => rhs_val)) - - @named csys = System(corrn_fin) - - u0_c = correlation_u0(corr, sol_1m0.u[z]) - p0_c = correlation_p0(corr, sol_1m0.u[z], ps .=> p0) - dict = merge(Dict(u0_c), Dict(p0_c)) - - prob_c = ODEProblem(csys, dict, (τ0, τend)) - sol_c = solve(prob_c, alg; saveat = saveat) - - idx_v = extract_idxs[l] - u_flat = [ComplexF64.(getindex.(sol_c.u, idx_v[i])) for i = 1:length(Lvec)] - - if i1 == 1 - for k = 1:length(Lvec) - Gmat_bar[i1, k, l, :] .= u_flat[k] - end - else - for k = 1:length(Lvec) - Gmat_bar[i1, k, l, :] .= padleftN(u_flat[k], N) - end - end - end - end - - return Gmat_bar -end - - -# Optionally compute and save Gmat_bar -#Gmat_bar = compute_Gmat( -# sol_1m0, eqs_1m0, Lvec, Corrvec; -# p0 = p0, -# ps = ps, -# saveat = va, -# τ0 = 0.0, -# alg = Tsit5(), -# state_syms = oper_tup, -# extract_idxs = idx_tup -#) - -# Load precomputed correlation data -#JLD2.save_object("Corr_Single.jld2",[Gmat_bar,t1]) -Gmat_bar = JLD2.load_object("Corr_Single.jld2")[1]; -t1 = JLD2.load_object("Corr_Single.jld2")[2]; -N = length(t1) -# Build g^(1)(t1,t2) -gcorr1 = zeros(ComplexF64, N, N) -gcorr = zeros(ComplexF64, N, N) - -for i1 = 1:N - for j1 = 1:N - gcorr1[j1, i1] = - gu_c_t(t1[j1])*gu_t(t1[i1])*Gmat_bar[i1, 1, 1, j1] + - gu_c_t(t1[j1])*sqrt(γ_)*Gmat_bar[i1, 1, 2, j1] + - gu_t(t1[i1])*sqrt(γ_)*Gmat_bar[i1, 2, 1, j1] + - γ_*Gmat_bar[i1, 2, 2, j1] - end -end -gcorr1 -Λ = zeros(ComplexF64, N, N) -for i = 1:N - Λ[i, i] = diag(real.(gcorr1))[i] -end - -# Symmetrized final g^(1) -gcorr = gcorr1 + tril(gcorr1, 1)' - Λ - - -################################################### -########### Plots: I-O SUPER Corr. ################ -################################################### - -################################################### -# Heat plot -################################################### - -# Normalization for visualization -g_abs = abs.(gcorr) -g_max = maximum(g_abs) -g_vis = g_abs ./ (g_max > 0 ? g_max : 1) - -tick_vals = collect(0.0:0.2:1.0) -tick_labs = latexstring.(string.(round.(tick_vals; digits = 1))) -pl6 = heatmap( - t1, - t1, - g_vis; - color = :thermal, - clims = (0, 1), - xlabel = L"\gamma t_2", - ylabel = L"\gamma t_1", - colorbar_title = L"|g^{(1)}|\ \mathrm{(norm)}", - colorbar_ticks = (tick_vals, tick_labs), - aspect_ratio = :equal, - tickfontsize = 16, - guidefontsize = 18, - size = (600, 500), -) -display(pl6) -# savefig(pl6, "heat_IO.pdf") - - -##################################################################### -# Eigenvalues and Eigenmodes -##################################################################### - -F = eigen(gcorr); - -ΔT = t1[2]-t1[1] -ns_avg = F.values -modess = F.vectors -vs_mode = (modess[:, end]) / sqrt(ΔT) -ns_avgt = ns_avg*(t1[2] - t1[1]) -ns = ns_avgt[end] -ns_avgt - -################################################### -########### Plot I-O SUPER Cascade ################ -################################################### - -# Dominant mode magnitude vs composite input -pl7 = plot( - t1, - abs.(u.(t1)); - color = :blue, - label = L"v(t) = u(t)", - ylabel = L"|v(t)|", - xlims = (t1[1]-0.01, t1[end]), - xticks = ([0, 10, 20], ["", "", ""]), - legend = :topleft, - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, - size = (800, 200), -) -plot!(pl7, t1, abs.(vs_mode); color = :red, ls = :dash, label = L"\mathrm{corr.}") -display(pl7) -# savefig(pl7, "mode_IO.pdf") - - -################################################### -## Loss mode: I-O SUPER with corr-based approach ## -################################################### -dt = 1e-5 -Tend = 12*2 -T = [dt:dt:Tend;] -ΔT = T[2] - T[1] - -# Hilbert space -hu = FockSpace(:u) # virtual input cavity -hv = FockSpace(:v) # virtual output cavity -hs = tensor(hu, hs1, hv) - -# Operators -aun = Destroy(hs, :a_u, 1) -σ(i, j) = Transition(hs, :ss, i, j) -avn = Destroy(hs, :a_v, 3) - -# Couplings -gu, gv = cnumbers("g_{u} g_{v}") - - -################################################### -# QuantumCumulants.jl section -################################################### - -# SLH triplets -G_u = SLH(1, gu*aun, 0) # input cavity -G_2lvls = SLH(1, √(γ)*σ(1, 2), 0) # 2-level system -G_v = SLH(1, gv*avn, 0) # output cavity - -# I-O cascade -G_cas_1m = ▷(G_u, G_2lvls, G_v) - -# Hamiltonian and Lindbladian -Hcas_1m = hamiltonian(G_cas_1m) -Lcas_1m = jump_operator(G_cas_1m)[1] -Lcas_1md = adjoint(Lcas_1m) - -# Time-dependent couplings -@register_symbolic gu_t(t) -@register_symbolic gv_t(t) -# Complex conjugates -@register_symbolic gu_c_t(t) -@register_symbolic gv_c_t(t) - -gs_ls = [gu, gv] -gst_ls = [gu_t(t), gv_t(t)] -gsct_ls = [gu_c_t(t), gv_c_t(t)] - -dict_gts = Dict([gs_ls; conj.(gs_ls)] .=> [gst_ls; gsct_ls]); - -# Insert time-dependence -Hcas_1m_t = substitute(Hcas_1m, dict_gts); -Lcas_1m_t = substitute(Lcas_1m, dict_gts); -Lcas_1md_t = substitute(Lcas_1md, dict_gts); - -# Pulses and composite mode -αeff = 0.0 - -Ω1(t_) = α1_/(sqrt(2π*σ1_^2)) * exp(-(t_-τ1_)^2/(2σ1_^2)) * exp(-1im*Δ1_*t_) - -Ω2(t_) = α2_/(sqrt(2π*σ2_^2)) * exp(-(t_-τ2_)^2/(2σ2_^2)) * exp(-1im*Δ2_*t_) - -Ωall(t_) = Ω1(t_) + Ω2(t_) - -αeff = sqrt(cumul_integrate(T, abs2.(Ωall.(T)))[end]) - -u(t_) = Ωall(t_) / αeff - -gu_t_ = coupling_input(u, T) -gu_t(t) = gu_t_(t) -gv_t_ = coupling_output(u, T) -gv_t(t) = gv_t_(t) - - -gu_c_t(t) = conj(gu_t(t)) -gv_c_t(t) = conj(gv_t(t)) - - -order = 2 -ops1 = [aun, aun', σ(1, 2), σ(2, 2), σ(1, 2)', avn, avn'] -ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); -ops_all = union( - ops1, - ops_all1[1:1], - ops_all1[3:14], - ops_all1[16:16], - ops_all1[18:26], - ops_all1[28:29], -) -eqs_1m = meanfield( - ops_all, - Hcas_1m_t, - [Lcas_1m_t]; - Jdagger = [Lcas_1md_t], - order = order, - iv = t, -) - -complete!(eqs_1m) -length(eqs_1m) - -# Classical-to-coherent relation -α_1m = αeff/(2*sqrt(γ_)) # may change if cross-terms are significant - -function ψ0_us1v(α, eqs) - bu = FockBasis(rd(abs2(α) + 20abs(α)), rd(abs2(α) - 20abs(α))) # include ±20α photon fluctuations - bs1 = NLevelBasis(2) - bv = FockBasis(1) - b = tensor([bu, bs1, bv]...) - ψu = coherentstate(bu, α) - ψs1 = nlevelstate(bs1, 1) - ψv = fockstate(bv, 0) - ψ0 = LazyKet(b, (ψu, ψs1, ψv)) - return initial_values(eqs, ψ0) -end - -# Initial values -u0_1m = ψ0_us1v(α_1m, eqs_1m) - -va = 0.05 - -# Solve ODE system -sys_1m = mtkcompile(System(eqs_1m; name = :sys_1m)) -ps = [γ] -p0 = [γ_] -u0maps = Dict(unknowns(sys_1m) .=> u0_1m) -pmaps = parameter_map(sys_1m, ps, p0) -probs = ODEProblem(sys_1m, merge(u0maps, pmaps), (dt, Tend)) -sol_1m0v = solve(probs, Tsit5(), saveat = va; abstol, reltol); - - -################################################### -########### Correlation Part ###################### -################################################### - -# Lindblad vector (includes output) -Lvec=[aun, σ(1, 2), avn] - -N = length(sol_1m0v.t) -Gmat_bar = zeros(ComplexF64, N, length(Lvec), length(Lvec), N) -t1 = sol_1m0v.t; -t2 = sol_1m0v.t; -τ0 = 0.0 - -Corrmat = [ - [ - CorrelationFunction(Lvec[k]', Lvec[l], eqs_1m; steady_state = false) for - l = 1:length(Lvec) - ] for k = 1:length(Lvec) -] -Corrvec = [Corrmat[end][i] for i = 1:length(Lvec)] - -oper_tup = (Lvec[1], Lvec[2], Lvec[3]) -idx_tup = idx_tuple(Lvec, Corrvec, Corrmat) - -# Optionally compute and save -#Gmat_bar = compute_Gmat( -# sol_1m0v, eqs_1m, Lvec, Corrvec; -# p0 = p0, -# ps = ps, -# saveat = va, -# τ0 = 0.0, -# alg = Tsit5(), -# state_syms = oper_tup, -# extract_idxs = idx_tup -#) - -# Load precomputed data -#JLD2.save_object("Corr_v.jld2",[Gmat_bar,t1]) -Gmat_bar = JLD2.load_object("Corr_v.jld2")[1]; -t1 = JLD2.load_object("Corr_v.jld2")[2]; -N = length(t1) -gcorr1 = zeros(ComplexF64, N, N) -gcorr = zeros(ComplexF64, N, N) - -gvec(t) = [gu_t(t), sqrt(γ_), gv_t(t)] -gvecc(t) = [gu_c_t(t), sqrt(γ_), gv_c_t(t)] - -for i1 = 1:N - for j1 = 1:N - gcorr1[i1, j1] = sum( - gvecc(t1[i1])[k]*gvec(t1[j1])[l]*(Gmat_bar[j1, k, l, i1]) for k = 1:3 for - l = 1:3 - ) - end -end -gcorr1 -Λ = zeros(ComplexF64, N, N) -for i = 1:N - Λ[i, i] = diag(real.(gcorr1))[i] -end - -# Final g^(1)(t1,t2) -gcorr = gcorr1 + tril(gcorr1, 1)' - Λ - - - -################################################### -########### Plot I-O SUPER Corr. loss ############# -################################################### - -#################################################################### -# Heat plot -#################################################################### - -g_abs = abs.(gcorr) -g_max = maximum(g_abs) -g_vis = g_abs ./ (g_max > 0 ? g_max : 1) - -tick_vals = collect(0.0:0.2:1.0) -tick_labs = latexstring.(string.(round.(tick_vals; digits = 1))) -pl8 = heatmap( - t1, - t1, - g_vis; - color = :thermal, - clims = (0, 1), - xlabel = L"\gamma t_2", - ylabel = L"\gamma t_1", - colorbar_title = L"|g^{(1)}|\ \mathrm{(norm)}", - colorbar_ticks = (tick_vals, tick_labs), - aspect_ratio = :equal, - tickfontsize = 16, - guidefontsize = 18, - size = (600, 500), -) -display(pl8) -# savefig(pl8, "heat_IO_loss.pdf") - - -##################################################################### -# Eigenvalues and eigenmodes -##################################################################### - -F = eigen(gcorr); - -ΔT = t1[2]-t1[1] -ns_avg = F.values -modess = F.vectors -vs_mode = (modess[:, end]) / sqrt(ΔT) -ns_avgt = ns_avg*(t1[2] - t1[1]) -ns = ns_avgt[end] -ns_avgt - -################################################### -########### Plot I-O SUPER Cascade mode ########### -################################################### - -# Dominant mode magnitude (loss-inclusive) -pl9 = plot( - t1, - abs.(u.(t1)); - color = :blue, - label = L"v(t) = u(t)", - ylabel = L"|v(t)|", - xlims = (t1[1]-0.01, t1[end]), - xticks = ([0, 10, 20], ["", "", ""]), - legend = :topleft, - tickfontsize = plot_font, - guidefontsize = plot_font, - legendfontsize = plot_font, - size = (800, 200), -) -plot!(pl9, t1, abs.(vs_mode); label = L"\mathrm{loss}") -display(pl9) -# savefig(pl9, "mode_IO_loss.pdf") +### #TODO: + +# # TODO: same in other file: delete conj fcts +# # TODO: delete nv term in other file! not correct in Δn! + + +# # Input-Output Analysis of Quantum Dot SUPER Excitation +# +# This script collects and demonstrates several approaches used in [PAPER REF]. +# +# We provide a quantum pulse description of the SUPER excitation scheme [T. K. Bracht et al., PRX Quantum 2, 040354 (2021)] (https://doi.org/10.1103/PRXQuantum.2.040354) +# and show that SUPER is a three-photon process [Q. W. Richter et al., Physical Review Research 7, 013079 (2025)] (https://doi.org/10.1103/PhysRevResearch.7.013079), +# [L. Vannucci and N. Gregersen, Optics Express Vol. 32, Issue 20, pp. 35381-35394 (2024)] (https://doi.org/10.1364/OE.533998). +# Due to large photon numbers, we use QuantumCumulants.jl [D. Plankensteiner et al., Quantum 6, 617 (2022)] (https://doi.org/10.22331/q-2022-01-04-617). +# +# We also propose full master-equation treatments for large-α coherent states in [PAPER REF]; here we show the “Displacement Frame” approach. +# +# Additionally we demonstrate how suitable output modes can be determined utilizing a correlation-based approach +# +# 1) Two-Input-Two-Output (2I-2O) SUPER Cascade +# +# Two Gaussian coherent-state pulses are emitted by two virtual input cavities, drive a TLS, and are collected by matched output cavities (v_i(t) = u_i(t)). +# The difference between input and output photon numbers per pulse reveals the three-photon nature. +# +# 2) Two-Input-Two-Output (2I-2O) SUPER Concatenation +# +# The same physics can be captured by concatenating channels: for each pulse i, input cavity i -> TLS -> output cavity i. +# The SLH description has two inputs and two outputs; photon and system dynamics recover the multi-photon character. +# +# 3) Interaction Picture and Displacement Frame of 2I-2O Cascade +# +# We transform the 2I-2O cascade into the cavity–cavity interaction picture [V. R. Christiansen et al., Physical Review A 107, 013706 (2023)] (https://doi.org/10.1103/PhysRevA.107.013706). +# We then enter a displacement frame that splits coherent backgrounds from quantum fluctuations, reducing the Fock spaces and enabling full master-equation treatment. +# +# 4) Input-Output SUPER Cascade +# +# We find a single coherent mode representing the SUPER superposition. This compressed description hides the explicit multi-photon nature, but greatly speeds up the QuantumCumulants.jl-based evaluation of the two-time correlation matrix g^(1)(t_1, t_2). +# Diagonalizing g^(1) yields the dominant output modes of the Input–TLS SLH triplet and of the full I–O SUPER cascade. +# # TODO: re-order + +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumCumulants +using ModelingToolkitBase +using OrdinaryDiffEq +using LinearAlgebra +using QuantumOptics +using NumericalIntegration +using JLD2 +using Plots +using LaTeXStrings + +abstol = 1e-10 +reltol = 1e-10 +cd(@__DIR__) +solution_values(sol, op, eqs) = get_solution(sol, op, eqs).(sol.t) +# parameter_map(sys, ps, vals) = +# Dict(p => v for (p, v) in zip(ps, vals) if any(isequal(p), parameters(sys))) + +################################################### +############# 1) 2I-2O SUPER Cascade ############## +################################################### + +# Hilbert spaces +hu2 = FockSpace(:u2) # virtual input cavity (u2) +hu1 = FockSpace(:u1) # virtual input cavity (u1) +hs1 = NLevelSpace(:atom, 2) # TLS +hv1 = FockSpace(:v1) # virtual output cavity (v1) +hv2 = FockSpace(:v2) # virtual output cavity (v2) +h = tensor(hu2, hu1, hs1, hv1, hv2) + +# Operators +au2 = Destroy(h, :a_u2, 1) +au1 = Destroy(h, :a_u1, 2) +s(i, j) = Transition(h, :s, i, j) +av1 = Destroy(h, :a_v1, 4) +av2 = Destroy(h, :a_v2, 5) + +# Symbolic parameters: decay rate and virtual cavity couplings +@variables γ::Real gu1::Number gu2::Number gv1::Number gv2::Number +@independent_variables t # Symbolic time variable + +################################################### +############# SUPER Pulse Parameters ############## +################################################### + +# Time grid +dt = 1e-4; +Tend = 12*2 +T = [dt:dt:Tend;]; +ΔT = T[2] - T[1] + +# Numerical pulse parameters +γ_ = 1e-2 # TLS decay rate +Δ1_ = -2π*1.934 # detuning pulse 1 +Δ2_ = -2π*4.634 # detuning pulse 2 +α1_ = 22.65π # pulse area 1 +α2_ = 19.29π # pulse area 2 +σ1_ = 2.4 # temporal FWHM pulse 1 +σ2_ = 3.04 # temporal FWHM pulse 2 +τ1_ = 0.0 + 12 # time shift pulse 1 +τ2_ = -0.73 + 12 # time shift pulse 2 + +# Classical Gaussian pulses +Ω1(t) = 0.5*1/(√(2π))*α1_/σ1_*exp(-(t-τ1_)^2/(2*σ1_^2))*exp(-1im*Δ1_*t) +Ω2(t) = 0.5*1/(√(2π))*α2_/σ2_*exp(-(t-τ2_)^2/(2*σ2_^2))*exp(-1im*Δ2_*t) + +################################################### +##### Classical SUPER Scheme (reference) ########## +################################################### + +@register_symbolic Ω(t) +@register_symbolic Ω_c(t) + +Ω(t) = Ω1(t) + Ω2(t) +Ω_c(t) = conj(Ω(t)) + +H_cl = -(Ω_c(t)*s(1, 2) + Ω(t)*s(2, 1)) +J_cl = [s(1, 2)]; +rates_cl = [γ] + +ops_cl = [s(2, 2), s(1, 2)] +eqs_cl = meanfield(ops_cl, H_cl, J_cl; rates = rates_cl, iv = t) +sys_cl = mtkcompile(System(eqs_cl; name = :sys_cl)) + +u0_cl = zeros(ComplexF64, length(eqs_cl)) +u0_p_map_cl = merge(Dict(unknowns(sys_cl) .=> u0_cl), Dict(γ => γ_)) +prob_cl = ODEProblem(sys_cl, u0_p_map_cl, (dt, Tend)) +sol_cl = solve(prob_cl, Tsit5(); abstol, reltol) +t_cl = sol_cl.t +s22_cl = real.(solution_values(sol_cl, s(2, 2), eqs_cl)); + +using PyPlot; +pygui(true) + +close(123) +figure(123) +subplot(211) +PyPlot.plot(t_cl, abs.(Ω.(t_cl))) +PyPlot.grid(true) +subplot(212) +PyPlot.plot(t_cl, s22_cl) +PyPlot.grid(true) + +### effective: 3photon process### + +Ω1_max = Α1_/σ1_/√(2π) +Ω2_max = Α2_/σ2_/√(2π) + +Ω_3ph = Ω1_max^2*Ω2_max/(4*Δ1_*Δ2_) +T_3ph = π/Ω_3ph +0.01/3.99 + +Ω_max = maximum(abs.(Ω1.(T)) + abs.(Ω2.(T))) +Ω_eff(t) = 0.5 * Ω_3ph * abs.(abs.(Ω1.(t)) + abs.(Ω2.(t)))/Ω_max + +b = NLevelBasis(2) +s(i, j) = transition(b, i, j) + +H_(t, psi) = Ω_eff(t)*(s(1, 2) + s(2, 1)) + +t_, psi_ = timeevolution.schroedinger_dynamic(T, nlevelstate(b, 1), H_) + +s22_eff = real(expect(s(2, 2), psi_)) +subplot(212) +PyPlot.plot(t_, s22_eff) + + +111 + +################################################### +####### IOT and QuantumCumulants.jl section ####### +################################################### + +# SLH triplets +G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 +G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 +G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system +G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 +G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 + +# 2I-2O cascade SLH triplet +G_cas = ▷(G_u2, G_u1, G_2lvl, G_v1, G_v2) # cascade + +# Hamiltonian and Lindbladian +Hcas = hamiltonian(G_cas) +Jcas = jump_operator(G_cas)[1] +Jcasd = adjoint(Jcas) + +# Time-dependent couplings +@register_symbolic gu1_t(t) +@register_symbolic gu2_t(t) +@register_symbolic gv1_t(t) +@register_symbolic gv2_t(t) +# Complex conjugates +@register_symbolic gu1_c_t(t) +@register_symbolic gu2_c_t(t) +@register_symbolic gv1_c_t(t) +@register_symbolic gv2_c_t(t) + +g_ls = [gu2, gu1, gv1, gv2] +gt_ls = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] +gct_ls = [gu2_c_t(t), gu1_c_t(t), gv1_c_t(t), gv2_c_t(t)] +dict_gt = Dict([g_ls; conj.(g_ls)] .=> [gt_ls; gct_ls]); + +# Insert time-dependence +Hcas_t = substitute(Hcas, dict_gt); +Jcas_t = substitute(Jcas, dict_gt); +Jcasd_t = substitute(Jcasd, dict_gt); + +# Normalized input modes +u1(t_) = 1/(√(σ1_)*π^(1/4)) * exp(-(t_ - τ1_)^2 / (2*σ1_^2)) * exp(-1im*Δ1_*t_) +u2(t_) = 1/(√(σ2_)*π^(1/4)) * exp(-(t_ - τ2_)^2 / (2*σ2_^2)) * exp(-1im*Δ2_*t_) + +# Coupling functions +gu1_t_ = coupling_input(u1, T) +gu1_t(t) = gu1_t_(t) +gv1_t_ = coupling_output(u1, T) +gv1_t(t) = gv1_t_(t) +gu1_c_t(t) = conj(gu1_t(t)) +gv1_c_t(t) = conj(gv1_t(t)) + +# Cascade-modified couplings (effective modes) +u_fcts = [u1, u2] +u2_eff = effective_input_mode(u_fcts, T, 2; abstol, reltol) +gu2_t_ = coupling_input(u2_eff, T) +gu2_t(t) = gu2_t_(t) +gu2_c_t(t) = conj(gu2_t(t)); + +v_fcts = [u1, u2] # output modes = input modes +v2_eff = effective_output_mode(v_fcts, T, 2; abstol, reltol) +gv2_t_ = coupling_output(v2_eff, T) +gv2_t(t) = gv2_t_(t) +gv2_c_t(t) = conj(gv2_t(t)); + +# First-order cumulant expansion +order = 1 +ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] +eqs = meanfield(ops, Hcas_t, [Jcas_t]; Jdagger = [Jcasd_t], order = order, iv = t) + +# Classical-to-coherent amplitude relation +α1_io = α1_ / (2*√(2)*π^(1/4)*√(σ1_*γ_)) # field 1 +α2_io = α2_ / (2*√(2)*π^(1/4)*√(σ2_*γ_)) # field 2 +u0 = [α1_io, α2_io, 0, 0, 0, 0.0im] + +# Solve ODE system +sys = mtkcompile(System(eqs; name = :sys)) +u0_p_map_cas = Dict([unknowns(sys); γ] .=> [u0; γ_]) +prob_cas = ODEProblem(sys, u0_p_map_cas, (dt, Tend)) +sol = solve(prob_cas, Tsit5(); abstol, reltol); + +# Expectation values +t_cas = sol.t # time vector +s22_cas = solution_values(sol, s(2, 2), eqs); +nu1_cas = abs2.(solution_values(sol, au1, eqs)) +nu2_cas = abs2.(solution_values(sol, au2, eqs)) +nv1_cas = abs2.(solution_values(sol, av1, eqs)) +nv2_cas = abs2.(solution_values(sol, av2, eqs)) + +# Plots +plot_font = 18 +common = (; + xlims = (t_cas[1]-0.01, t_cas[end]), + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, +) +p1 = plot( + t_cl, + s22_cl; + color = :blue, + label = L"\mathrm{classical}", + ylabel = L"P", + xticks = ([0, 10, 20], ["", "", ""]), + yticks = ([0, 0.5, 1], [L"0.0", L"0.5", L"1.0"]), + legend = :topleft, + common..., +) +plot!(p1, t_cas, real.(s22_cas); color = :red, ls = :dash, label = L"\mathrm{cascade}") +p2 = plot( + t_cas, + nu1_cas; + color = :blue, + label = L"\langle \hat{n}_{u_1} \rangle", + ylabel = L"\langle\hat{n}_i\rangle", + xticks = ([0, 10, 20], ["", "", ""]), + yticks = ([0, 5e3, 10e3, 15e3], [L"0", L"5\cdot10^3", L"10\cdot10^3", L"15\cdot10^3"]), + ylims = (nv1_cas[1]-0.5e3, nu1_cas[1]+0.5e3), + legend = :left, + common..., +) +plot!(p2, t_cas, nv1_cas; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") +plot!( + p2, + t_cas, + nu2_cas; + color = :blue, + ls = :dash, + label = L"\langle \hat{n}_{u_2} \rangle", +) +plot!( + p2, + t_cas, + nv2_cas; + color = :red, + ls = :dash, + label = L"\langle \hat{n}_{v_2} \rangle", +) +p3 = plot( + t_cas, + nu1_cas .+ nv1_cas .- nu1_cas[1]; + color = :blue, + label = L"\mathrm{mode~1}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + xticks = ([0, 10, 20], latexstring.([0, 10, 20])), + yticks = ([-4, -2, 0, 2, 4], latexstring.([-4, -2, 0, 2, 4])), + legend = :topleft, + common..., +) +plot!(p3, t_cas, nu2_cas .+ nv2_cas .- nu2_cas[1]; color = :red, label = L"\mathrm{mode~2}") +pl1 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) +display(pl1) +# savefig(pl1, "plot_2I2O_cascade.pdf") + + +#################################################################### +##### 3) 2I-2O SUPER Cascade interaction picture and Displacement ## +#################################################################### + +# Interaction picture: cavity dynamics +H_uv = hamiltonian(▷(G_u2, G_u1, G_v1, G_v2)) +H_int_sym_ = simplify(Hcas - H_uv) + +M(i, j) = Symbolics.variable(Symbol("M_{$(i)$(j)}"); T = Number) +Mc(i, j) = Symbolics.variable(Symbol("Mc_{$(i)$(j)}"); T = Number) + +a0_ls = [au2, au1, av1, av2] +la = length(a0_ls) +a_int_ls = [sum(M(i, j)*a0_ls[j] for j = 1:la) for i = 1:la] +# a_c_int_ls = [sum(Mc(i,j)*a0_ls[j]' for j = 1:la) for i = 1:la] # # TODO +a_c_int_ls = [sum(conj(M(i, j))*a0_ls[j]' for j = 1:la) for i = 1:la] # # TODO + +int_dict = Dict([a0_ls; adjoint.(a0_ls)] .=> [a_int_ls; a_c_int_ls]) + +H_int_sym = substitute(H_int_sym_, int_dict) +J_int_sym = simplify(substitute(Jcas, int_dict)) +Jd_int_sym = simplify(substitute(Jcasd, int_dict)) + +# M-matrix +Mat = Matrix{Any}(undef, la, la) +Matc = Matrix{Any}(undef, la, la) +mod = @__MODULE__ +for i = 1:la, j = 1:la + name = Symbol("Ma_$(i)$(j)") + namec = Symbol("Mac_$(i)$(j)") + + @eval @register_symbolic $name(t) + @eval @register_symbolic $namec(t) + + Mat[i, j] = getfield(mod, name)(t) + Matc[i, j] = getfield(mod, namec)(t) +end +Mat_ls = [Mat[i, j] for i = 1:la for j = 1:la] +Mat_conls = [Matc[i, j] for i = 1:la for j = 1:la] + +M_ls = [M(i, j) for i = 1:la for j = 1:la] +Mc_ls = [Mc(i, j) for i = 1:la for j = 1:la] + +# Time-evolution matrix M(t) +A_uv = coupling_matrix((gu2_t_, gu1_t_, gv1_t_, gv2_t_)) +M_t = solve_mode_evolution(A_uv, T) +Mc_t(t) = adjoint(M_t(t)) + +for i = 1:la, j = 1:la + fname = Symbol("Ma_$(i)$(j)") + fnamec = Symbol("Mac_$(i)$(j)") + + @eval begin + $fname(t) = M_t(t)[$i, $j] + $fnamec(t) = Mc_t(t)[$i, $j] + end +end + +dict_Mt = Dict([M_ls..., Mc_ls...] .=> [Mat_ls..., Mat_conls...]) +dict_gt_Mt = merge(dict_gt, dict_Mt) + +H_int_sym_t = substitute(H_int_sym, dict_gt_Mt) +J_int_sym_t = substitute(J_int_sym, dict_gt_Mt) +Jd_int_sym_t = substitute(Jd_int_sym, dict_gt_Mt) + +eqs_int = meanfield( + ops, + H_int_sym_t, + [J_int_sym_t]; + Jdagger = [Jd_int_sym_t], + order = order, + iv = t, +); + +# Solve ODE system in interaction picture +sys_int = mtkcompile(System(eqs_int; name = :sysI)) +u0_p_map_int = Dict([unknowns(sys_int); γ] .=> [u0; γ_]) +prob_int = ODEProblem(sys_int, u0_p_map_int, (dt, Tend)) +sol_int = solve(prob_int, Tsit5(); abstol, reltol); + +# Expectation values +t_int = sol_int.t +s22_int = solution_values(sol_int, s(2, 2), eqs_int); +nu1_int = abs2.(solution_values(sol_int, au1, eqs_int)) +nu2_int = abs2.(solution_values(sol_int, au2, eqs_int)) +nv1_int = abs2.(solution_values(sol_int, av1, eqs_int)) +nv2_int = abs2.(solution_values(sol_int, av2, eqs_int)) + +### Displacement frame ### + +# QuantumOptics.jl Hilbert space for fluctuations # TODO: Why is a cutoff of 1 enough? +bu1 = FockBasis(1) # vacuum fluctuation mode 1 +bu2 = FockBasis(1) # vacuum fluctuation mode 2 +bs1 = NLevelBasis(2) # TLS +bv1 = FockBasis(1) # vacuum mode 1 +bv2 = FockBasis(1) # vacuum mode 2 +b = tensor([bu2, bu1, bs1, bv1, bv2]...) + +ψu1 = fockstate(bu1, 0) +ψu2 = fockstate(bu2, 0) +ψs1 = nlevelstate(bs1, 1) +ψv1 = fockstate(bv1, 0) +ψv2 = fockstate(bv2, 0) +ψ0 = tensor(ψu2, ψu1, ψs1, ψv1, ψv2) + +# Operators +anu2 = destroy(bu2) +anu1 = destroy(bu1) +anv1 = destroy(bv1) +anv2 = destroy(bv2) +s12 = embed(b, 3, transition(bs1, 1, 2)) +s22 = embed(b, 3, transition(bs1, 2, 2)) + +# Mode operator vector +avn = [embed(b, 1, anu2), embed(b, 2, anu1), embed(b, 4, anv1), embed(b, 5, anv2)] + +# Coherent amplitude shifts (input modes only) +αvec = [α2_io, α1_io, 0.0, 0.0] + +# Displacement-frame Hamiltonian and jumps +function Ham_displaced(t, ρ) + Mall = M_t(t) + gall = [gu2_t(t), gu1_t(t), gv1_t(t), gv2_t(t)] + + # Modified coupling vectors for H and J terms + gn = [gall[1], gall[2], -gall[3], -gall[4]] + gnc = conj.(gn) + gl = [gall[1], gall[2], gall[3], gall[4]] + glc = conj.(gl) + + # Subspaces and local annihilation operators + B = [bu2, bu1, bv1, bv2] + Av = [anu2, anu1, anv1, anv2] + + # Displaced operators: a -> b + α (b = fluctuation mode) + avn_disp = [] + for k = 1:4 + if k < 3 + push!(avn_disp, embed(b, k, (Av[k] + αvec[k]*one(B[k])))) + else + push!(avn_disp, embed(b, k+1, (Av[k] + αvec[k]*one(B[k])))) + end + end + + # Interaction pieces (IP + displacement) + A_gn = transpose(gn) * Mall * avn_disp + A_gnc = avn_disp' * Mall' * gnc + + A_gl = transpose(gl) * Mall * avn_disp + A_glc = avn_disp' * Mall' * glc + + # Displacement Hamiltonian + Ht = sqrt(γ_) * 1im/2 * (A_gn * s12' - s12 * A_gnc) + + # Jumps (displacement Lindbladian) + J = [sqrt(γ_) * s12 + A_gl] + Jd = [sqrt(γ_) * s12' + A_glc] + + return Ht, J, Jd +end + +# Master equation in displacement frame +t_df, ρt_df = timeevolution.master_dynamic(T, ψ0, Ham_displaced; abstol, reltol); + +# Expectation values +s22_df = real.(expect(s22, ρt_df)) +nu2_df = real.( + expect( + embed( + b, + 1, + (anu2 + α2_io*identityoperator(bu2))'*(anu2 + α2_io*identityoperator(bu2)), + ), + ρt_df, + ), +) +nu1_df = real.( + expect( + embed( + b, + 2, + (anu1 + α1_io*identityoperator(bu1))'*(anu1 + α1_io*identityoperator(bu1)), + ), + ρt_df, + ), +) +nv1_df = real.(expect(embed(b, 4, anv1'*anv1), ρt_df)) +nv2_df = real.(expect(embed(b, 5, anv2'*anv2), ρt_df)) + + +# Plot net loss/gain in interaction picture and displacement frames +pl4 = plot( + t_int, + nu1_int .+ nv1_int .- nu1_int[1]; + color = :blue, + label = L"\mathrm{mode~1~(int.)}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + xlims = (t_df[1]-0.01, t_df[end]), + yticks = ([-2, -1, 0, 1], latexstring.([-2, -1, 0, 1])), + legend = :right, + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, + size = (800, 400), +) +plot!( + pl4, + t_int, + nu2_int .+ nv2_int .- nu2_int[1]; + color = :red, + label = L"\mathrm{mode~2~(int.)}", +) +plot!( + pl4, + t_df, + nu1_df .+ nv1_df .- nu1_df[1]; + color = :green, + ls = :dash, + label = L"\mathrm{mode~1~(dis.)}", +) +plot!( + pl4, + t_df, + nu2_df .+ nv2_df .- nu2_df[1]; + color = :orange, + ls = :dash, + label = L"\mathrm{mode~2~(dis.)}", +) +display(pl4) +# savefig(pl4, "plot_int_dis.pdf") + + +################################################### +######### 2) 2I-2O SUPER Concatenation ############ +################################################### + +# TODO: No substitution with γ/2! - interpretation: decay rate twice + +# SLH triplets +G_u2 = SLH(1, gu2*au2, 0) # input cavity 2 +G_u1 = SLH(1, gu1*au1, 0) # input cavity 1 +G_2lvl = SLH(1, √(γ)*s(1, 2), 0) # 2-level system # TODO +G_v1 = SLH(1, gv1*av1, 0) # output cavity 1 +G_v2 = SLH(1, gv2*av2, 0) # output cavity 2 + +# Channel-wise cascades +G_ch1 = ▷(G_u2, G_2lvl, G_v2) # pulse 1 +G_ch2 = ▷(G_u1, G_2lvl, G_v1) # pulse 2 +Gcon = ⊞(G_ch1, G_ch2) # concatenation + +# Hamiltonian and Lindblad operators (two channels) +Hcon = hamiltonian(Gcon) +Jch1g = jump_operator(Gcon)[1] +Jch2g = jump_operator(Gcon)[2] +Jdch1g = adjoint(Jch1g) +Jdch2g = adjoint(Jch2g) + +# Split decay γ -> γ/2 (symmetric into two channels) +# dict_gam = Dict(γ => γ/2); # TODO +dict_gam = Dict(γ => γ); + +# Apply decay adjustment +Jch1 = substitute(Jch1g, dict_gam) # TODO +Jch2 = substitute(Jch2g, dict_gam) +Jdch1 = substitute(Jdch1g, dict_gam) +Jdch2 = substitute(Jdch2g, dict_gam) + +# Insert time-dependence +Hcon_t = substitute(Hcon, dict_gt); +Jch1_t = substitute(Jch1, dict_gt); +Jch2_t = substitute(Jch2, dict_gt); +Jdch1_t = substitute(Jdch1, dict_gt); +Jdch2_t = substitute(Jdch2, dict_gt); + +gu2_t_ = coupling_input(u2, T) +gu2_t(t) = gu2_t_(t) +gv2_t_ = coupling_output(u2, T) +gv2_t(t) = gv2_t_(t) + +gu2_c_t(t) = conj(gu2_t(t)) +gv2_c_t(t) = conj(gv2_t(t)) + +# First-order cumulant expansion +order = 1 +ops = [au1, au2, s(2, 2), s(2, 1), av1, av2] +eqs_con = meanfield( + ops, + Hcon_t, + [Jch1_t, Jch2_t]; + Jdagger = [Jdch1_t, Jdch2_t], + order = order, + iv = t, +) + +# Solve ODE system +sys_con = mtkcompile(System(eqs_con; name = :sys_con)) + +u0_p_map_con = Dict([unknowns(sys_con); γ] .=> [u0; γ_]) # TODO +prob_con = ODEProblem(sys_con, u0_p_map_con, (dt, Tend)) +sol_con = solve(prob_con, Tsit5(); abstol, reltol); + +# Expectation values +t_con = sol_con.t +s22_con = solution_values(sol_con, s(2, 2), eqs_con) +nu1_con = abs2.(solution_values(sol_con, au1, eqs_con)) +nu2_con = abs2.(solution_values(sol_con, au2, eqs_con)) +nv1_con = abs2.(solution_values(sol_con, av1, eqs_con)) +nv2_con = abs2.(solution_values(sol_con, av2, eqs_con)) + + +common = (; + xlims = (t_con[1]-0.01, t_con[end]), + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, +) +p1 = plot( + t_cl, + real.(s22_cl); + color = :blue, + label = L"\mathrm{cascade}", + ylabel = L"P", + xticks = ([0, 10, 20], ["", "", ""]), + legend = :topleft, + common..., +) +plot!(p1, t_con, real.(s22_con); color = :red, ls = :dash, label = L"\mathrm{concat.}") +p2 = plot( + t_con, + nu1_con; + color = :blue, + label = L"\langle \hat{n}_{u_1} \rangle", + ylabel = L"\langle\hat{n}_i\rangle", + xticks = ([0, 10, 20], ["", "", ""]), + ylims = (nv1_con[1]-0.5e3, nu1_con[1]+0.5e3), + legend = :left, + common..., +) +plot!(p2, t_con, nv1_con; color = :red, label = L"\langle \hat{n}_{v_1} \rangle") +plot!( + p2, + t_con, + nu2_con; + color = :blue, + ls = :dash, + label = L"\langle \hat{n}_{u_2} \rangle", +) +plot!( + p2, + t_con, + nv2_con; + color = :red, + ls = :dash, + label = L"\langle \hat{n}_{v_2} \rangle", +) +p3 = plot( + t_con, + nu1_con .+ nv1_con .- nu1_con[1]; + color = :blue, + label = L"\mathrm{mode~1}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + yticks = ([-2, -1, 0, 1, 2], latexstring.([-2, -1, 0, 1, 2])), + legend = :bottomleft, + common..., +) +plot!(p3, t_con, nu2_con .+ nv2_con .- nu2_con[1]; color = :red, label = L"\mathrm{mode~2}") +pl2 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) +display(pl2) +# savefig(pl2, "plot_concat.pdf") + + +################################################### +########## 4) I-O SUPER Cascade ################### +################################################### + +# Single-mode SUPER description +dt = 1e-5 +Tend = 12*2 +T = [dt:dt:Tend;] +ΔT = T[2] - T[1] + +# Hilbert space +hu = FockSpace(:u) # virtual input cavity +hv = FockSpace(:v) # virtual output cavity +hs = tensor(hu, hs1, hv) + +# Operators +aun = Destroy(hs, :a_u, 1) +σ(i, j) = Transition(hs, :σ, i, j) +avn = Destroy(hs, :a_v, 3) + +# Couplings +@variables gu::Number gv::Number + +################################################### +# QuantumCumulants.jl section +################################################### + +# SLH triplets +G_u = SLH(1, gu*aun, 0) # input cavity +G_2lvls = SLH(1, √(γ)*σ(1, 2), 0) # 2-level system +G_v = SLH(1, gv*avn, 0) # output cavity + +# I-O cascade +G_cas_1m = ▷(G_u, G_2lvls, G_v) + +# Hamiltonian and Lindbladian +Hcas_1m = hamiltonian(G_cas_1m) +Jcas_1m = jump_operator(G_cas_1m)[1] +Jcas_1md = adjoint(Jcas_1m) + +# Time-dependent couplings +@register_symbolic gu_t(t) +@register_symbolic gv_t(t) +# Complex conjugates +@register_symbolic gu_c_t(t) +@register_symbolic gv_c_t(t) + +gs_ls = [gu, gv] +gst_ls = [gu_t(t), gv_t(t)] +gsct_ls = [gu_c_t(t), gv_c_t(t)] + +dict_gts = Dict([gs_ls; conj.(gs_ls)] .=> [gst_ls; gsct_ls]); + +# Insert time-dependence +Hcas_1m_t = substitute(Hcas_1m, dict_gts); +Jcas_1m_t = substitute(Jcas_1m, dict_gts); +Jcas_1md_t = substitute(Jcas_1md, dict_gts); + +0.5*1/(√(2π))*α1_/σ1_*exp(-(t-τ1_)^2/(2*σ1_^2))*exp(-1im*Δ1_*t) + +# Pulses and composite mode +Ω1(t_) = α1_/(sqrt(2π*σ1_^2))*exp(-(t_-τ1_)^2/(2σ1_^2))*exp(-1im*Δ1_*t_) +Ω2(t_) = α2_/(sqrt(2π*σ2_^2))*exp(-(t_-τ2_)^2/(2σ2_^2))*exp(-1im*Δ2_*t_) +Ωall(t_) = Ω1(t_) + Ω2(t_) +αeff = sqrt(cumul_integrate(T, abs2.(Ωall.(T)))[end]) +u(t_) = Ωall(t_) / αeff + +gu_t_ = coupling_input(u, T) +gu_t(t) = gu_t_(t) +gv_t_ = coupling_output(u, T) +gv_t(t) = gv_t_(t) + +gu_c_t(t) = conj(gu_t(t)) +gv_c_t(t) = conj(gv_t(t)) + +# Second-order expansion for correlation-based mode extraction +order = 2 +# Minimal operator set to capture second-order dynamics +ops_1m = [aun, aun'aun, aun*aun, σ(1, 2), σ(2, 2), avn, avn'avn] +# ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); # TODO: needed for corr? +# ops_all = union(ops1,ops_all1[1:1],ops_all1[3:14], ops_all1[16:16], ops_all1[18:26], ops_all1[28:29]) + +eqs_1m = + meanfield(ops_1m, Hcas_1m_t, [Jcas_1m_t]; Jdagger = [Jcas_1md_t], order = order, iv = t); +complete!(eqs_1m); +length(eqs_1m) + +# Classical-to-coherent relation +α_1m = αeff/(2*sqrt(γ_)) # may change if cross-terms are significant + +u0_1m = zeros(ComplexF64, length(eqs_1m)) +u0_1m[1] = α_1m +u0_1m[2] = abs2(α_1m) +u0_1m[3] = α_1m*α_1m + +# Solve ODE system +sys_1m = mtkcompile(System(eqs_1m; name = :sys_1m)) +ps = [γ] +p0 = [γ_] +u0_p_maps_1m = Dict([unknowns(sys_1m); γ] .=> [u0_1m; γ_]) +prob_1m = ODEProblem(sys_1m, u0_p_maps_1m, (dt, Tend)) +sol_1m = solve(prob_1m, Tsit5(), saveat = dt; abstol, reltol); + +t_1m = sol_1m.t; +s22_1m = solution_values(sol_1m, σ(2, 2), eqs_1m); +nu_1m = real.(solution_values(sol_1m, aun'aun, eqs_1m)) +nv_1m = real.(solution_values(sol_1m, avn'avn, eqs_1m)) + +common = (; + xlims = (t_1m[1]-0.01, t_1m[end]), + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, +) +p1 = plot( + t_cl, + s22_cl; + color = :blue, + label = L"\mathrm{classical}", + ylabel = L"P", + xticks = ([0, 10, 20], ["", "", ""]), + legend = :topleft, + common..., +) +plot!(p1, t_1m, real.(s22_1m); color = :red, ls = :dash, label = L"\mathrm{single~SUPER}") +p2 = plot( + t_1m, + nu_1m; + color = :blue, + label = L"\langle \hat{n}_{u} \rangle", + ylabel = L"\langle\hat{n}_i\rangle", + xticks = ([0, 10, 20], ["", "", ""]), + ylims = (nv_1m[1]-0.5e3, nu_1m[1]+0.5e3), + legend = :left, + common..., +) +plot!(p2, t_1m, nv_1m; color = :red, label = L"\langle \hat{n}_{v} \rangle") +p3 = plot( + t_1m, + nu_1m .+ nv_1m .- nu_1m[1]; + color = :blue, + label = L"\mathrm{mode}", + ylabel = L"\langle\Delta\hat{n}_i\rangle", + xlabel = L"\gamma t", + yticks = ([-1, -0.5, 0], latexstring.([-1, -0.5, 0])), + legend = :bottomleft, + common..., +) +pl4 = plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) +display(pl4) +# savefig(pl4, "plot_IO_cascade.pdf") + + +################################################### +########### I-O SUPER Correlation ################# +################################################### + +# Same as above but with sparser saveat (va = 0.05) and only Input–System cascade + +dt = 1e-5 +Tend = 12*2 +T = [dt:dt:Tend;] +ΔT = T[2] - T[1] + +@register_symbolic gu_t(t) +@register_symbolic gv_t(t) +@register_symbolic gu_c_t(t) +@register_symbolic gv_c_t(t) + +G_cas_1m0 = ▷(G_u, G_2lvls, G_v) + +Hcas_1m0 = hamiltonian(G_cas_1m0) +Jcas_1m0 = jump_operator(G_cas_1m0)[1] +Jcas_1m0d = adjoint(Jcas_1m0) + +# Set v-coupling to zero (collect only input-system correlations) +gs0_ls = [gu, gv] +gs0t_ls = [gu_t(t), 0.0*gv_t(t)] +gsc0t_ls = [gu_c_t(t), 0.0*gv_c_t(t)] + +dict_gts0 = Dict([gs0_ls; conj.(gs0_ls)] .=> [gs0t_ls; gsc0t_ls]); + +Hcas_1m0_t = substitute(Hcas_1m0, dict_gts0) +Jcas_1m0_t = substitute(Jcas_1m0, dict_gts0) +Jcas_1m0d_t = substitute(Jcas_1m0d, dict_gts0) + +gu_t_ = coupling_input(u, T) +gu_t(t) = gu_t_(t) +gv_t_ = coupling_output(u, T) +gv_t(t) = 0.0*gv_t_(t) + +gu_c_t(t) = conj(gu_t(t)) +gv_c_t(t) = conj(gv_t(t)) + +# order = 2 +# ops1 = [aun, aun', σ(1,2), σ(2,2),σ(1,2)', avn, avn'] +# ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); +# ops_all = union(ops1,ops_all1[1:1],ops_all1[3:14], ops_all1[16:16], ops_all1[18:26], ops_all1[28:29]) +# eqs_1m0 = meanfield(ops_all, Hcas_1m0_t, [Jcas_1m0_t]; Jdagger=[Jcas_1m0d_t], order=order, iv=t) + +# complete!(eqs_1m0) +# length(eqs_1m0) + +ops_1m0 = [aun, aun'aun, aun*aun, σ(1, 2), σ(2, 2)] +eqs_1m0 = meanfield( + ops_1m0, + Hcas_1m0_t, + [Jcas_1m0_t]; + Jdagger = [Jcas_1m0d_t], + order = order, + iv = t, +); +complete!(eqs_1m0); +length(eqs_1m0) + +u0_1m0 = zeros(ComplexF64, length(eqs_1m0)) +u0_1m0[1] = α_1m +u0_1m0[2] = abs2(α_1m) +u0_1m0[3] = α_1m*α_1m + +va = 0.05 + +# Solve ODE system +sys_1m0 = mtkcompile(System(eqs_1m0; name = :sys_1m0)) +u0_p_maps_1m0 = Dict([unknowns(sys_1m0); γ] .=> [u0_1m0; γ_]) +prob_1m0 = ODEProblem(sys_1m0, u0_p_maps_1m0, (dt, Tend)) +sol_1m0 = solve(prob_1m0, Tsit5(), saveat = va; abstol, reltol); + +################################################### +########### Correlation Part ###################### +################################################### + +# Lindblad vector for correlation functions +Jvec=[aun, σ(1, 2)] + +# Discrete time grid +N = length(sol_1m0.t) + +# g^(1)(t1,t2) container for each operator pair +Gmat_bar = zeros(ComplexF64, N, length(Jvec), length(Jvec), N) +t1 = sol_1m0.t; +t2 = sol_1m0.t; +τ0 = 0.0 + +# All operator correlation combinations +Corrmat = [ + [ + CorrelationFunction(Jvec[k]', Jvec[l], eqs_1m0; steady_state = false, iv0 = t) + for l = 1:length(Jvec) + ] for k = 1:length(Jvec) +] +# Minimal set containing full dynamics +Corrvec = [Corrmat[end][i] for i = 1:length(Jvec)] + +xx = CorrelationFunction(σ(1, 2), σ(1, 2), eqs_1m0; steady_state = false, iv0 = t) +xx.eqs.states + +################################################### +###### QC.jl extended correlation utilities ####### +################################################### + +# Structural equality of two lhs expressions +function same_lhs(a, b) + operation(a) === operation(b) && begin + aa = Tuple(arguments(a)) + bb = Tuple(arguments(b)) + length(aa) == length(bb) && all(isequal.(aa, bb)) + end +end + +# Index mapping between Corrvec and Corrmat +function idx_tuple(Jvec, Corrvec, Corrmat) + idxvec = zeros(Int8, length(Jvec), length(Jvec)) + for l = 1:length(Jvec) + eqs1 = equations(Corrvec[l].de) + for k = 1:length(Jvec) + target_lhs = equations(Corrmat[k][l].de)[1].lhs + idx = findfirst(eq -> same_lhs(eq.lhs, target_lhs), eqs1) + idxvec[k, l] = idx + end + end + idx_tuple = Tuple(ntuple(i -> c[i], length(Jvec)) for c in eachcol(idxvec)) + return idx_tuple +end + +# Operators included +oper_tup = (Jvec[1], Jvec[2]) +# Index tuple for selecting correlation components +idx_tup = idx_tuple(Jvec, Corrvec, Corrmat) + +# Compute full Gmat_bar on the discrete time grid +function compute_Gmat( + sol_1m0, + eqs_1m0, + Jvec, + Corrvec; + p0, + ps, + saveat, + τ0 = 0.0, + alg = OrdinaryDiffEq.Tsit5(), + state_syms = oper_tup, + extract_idxs = idx_tup, +) + + N = length(sol_1m0.t) + Gmat_bar = zeros(ComplexF64, N, length(Jvec), length(Jvec), N) + + padleftN(v::AbstractVector{<:Number}, N::Int) = ( + length(v) < N ? vcat(zeros(ComplexF64, N - length(v)), ComplexF64.(v)) : + ComplexF64.(v) + ) + + for i1 = 1:(N-1) + τend = sol_1m0.t[end] - sol_1m0.t[i1] + z = i1 + + for l = 1:length(Jvec) + corr = Corrvec[l] + complete(corr.de) + + iv = ModelingToolkit.get_iv(corr.de) + de_shift = substitute(corr.de, Dict(t => iv + sol_1m0.t[z])) + + bc_var = complete(corr.de)[end].lhs + rhs_val = sol_1m0[state_syms[l]][z] + corrn_fin = substitute(de_shift, Dict(bc_var => rhs_val)) + + @named csys = System(corrn_fin) + + u0_c = correlation_u0(corr, sol_1m0.u[z]) + p0_c = correlation_p0(corr, sol_1m0.u[z], ps .=> p0) + dict = merge(Dict(u0_c), Dict(p0_c)) + + prob_c = ODEProblem(csys, dict, (τ0, τend)) + sol_c = solve(prob_c, alg; saveat = saveat) + + idx_v = extract_idxs[l] + u_flat = [ComplexF64.(getindex.(sol_c.u, idx_v[i])) for i = 1:length(Jvec)] + + if i1 == 1 + for k = 1:length(Jvec) + Gmat_bar[i1, k, l, :] .= u_flat[k] + end + else + for k = 1:length(Jvec) + Gmat_bar[i1, k, l, :] .= padleftN(u_flat[k], N) + end + end + end + end + + return Gmat_bar +end + + +# Optionally compute and save Gmat_bar +#Gmat_bar = compute_Gmat( +# sol_1m0, eqs_1m0, Jvec, Corrvec; +# p0 = p0, +# ps = ps, +# saveat = va, +# τ0 = 0.0, +# alg = Tsit5(), +# state_syms = oper_tup, +# extract_idxs = idx_tup +#) + +# Load precomputed correlation data +#JLD2.save_object("Corr_Single.jld2",[Gmat_bar,t1]) +Gmat_bar = JLD2.load_object("Corr_Single.jld2")[1]; +t1 = JLD2.load_object("Corr_Single.jld2")[2]; +N = length(t1) +# Build g^(1)(t1,t2) +gcorr1 = zeros(ComplexF64, N, N) +gcorr = zeros(ComplexF64, N, N) + +for i1 = 1:N + for j1 = 1:N + gcorr1[j1, i1] = + gu_c_t(t1[j1])*gu_t(t1[i1])*Gmat_bar[i1, 1, 1, j1] + + gu_c_t(t1[j1])*sqrt(γ_)*Gmat_bar[i1, 1, 2, j1] + + gu_t(t1[i1])*sqrt(γ_)*Gmat_bar[i1, 2, 1, j1] + + γ_*Gmat_bar[i1, 2, 2, j1] + end +end +gcorr1 +Λ = zeros(ComplexF64, N, N) +for i = 1:N + Λ[i, i] = diag(real.(gcorr1))[i] +end + +# Symmetrized final g^(1) +gcorr = gcorr1 + tril(gcorr1, 1)' - Λ + + +################################################### +########### Plots: I-O SUPER Corr. ################ +################################################### + +################################################### +# Heat plot +################################################### + +# Normalization for visualization +g_abs = abs.(gcorr) +g_max = maximum(g_abs) +g_vis = g_abs ./ (g_max > 0 ? g_max : 1) + +tick_vals = collect(0.0:0.2:1.0) +tick_labs = latexstring.(string.(round.(tick_vals; digits = 1))) +pl6 = heatmap( + t1, + t1, + g_vis; + color = :thermal, + clims = (0, 1), + xlabel = L"\gamma t_2", + ylabel = L"\gamma t_1", + colorbar_title = L"|g^{(1)}|\ \mathrm{(norm)}", + colorbar_ticks = (tick_vals, tick_labs), + aspect_ratio = :equal, + tickfontsize = 16, + guidefontsize = 18, + size = (600, 500), +) +display(pl6) +# savefig(pl6, "heat_IO.pdf") + + +##################################################################### +# Eigenvalues and Eigenmodes +##################################################################### + +F = eigen(gcorr); + +ΔT = t1[2]-t1[1] +ns_avg = F.values +modess = F.vectors +vs_mode = (modess[:, end]) / sqrt(ΔT) +ns_avgt = ns_avg*(t1[2] - t1[1]) +ns = ns_avgt[end] +ns_avgt + +################################################### +########### Plot I-O SUPER Cascade ################ +################################################### + +# Dominant mode magnitude vs composite input +pl7 = plot( + t1, + abs.(u.(t1)); + color = :blue, + label = L"v(t) = u(t)", + ylabel = L"|v(t)|", + xlims = (t1[1]-0.01, t1[end]), + xticks = ([0, 10, 20], ["", "", ""]), + legend = :topleft, + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, + size = (800, 200), +) +plot!(pl7, t1, abs.(vs_mode); color = :red, ls = :dash, label = L"\mathrm{corr.}") +display(pl7) +# savefig(pl7, "mode_IO.pdf") + + +################################################### +## Loss mode: I-O SUPER with corr-based approach ## +################################################### +dt = 1e-5 +Tend = 12*2 +T = [dt:dt:Tend;] +ΔT = T[2] - T[1] + +# Hilbert space +hu = FockSpace(:u) # virtual input cavity +hv = FockSpace(:v) # virtual output cavity +hs = tensor(hu, hs1, hv) + +# Operators +aun = Destroy(hs, :a_u, 1) +σ(i, j) = Transition(hs, :ss, i, j) +avn = Destroy(hs, :a_v, 3) + +# Couplings +gu, gv = cnumbers("g_{u} g_{v}") + + +################################################### +# QuantumCumulants.jl section +################################################### + +# SLH triplets +G_u = SLH(1, gu*aun, 0) # input cavity +G_2lvls = SLH(1, √(γ)*σ(1, 2), 0) # 2-level system +G_v = SLH(1, gv*avn, 0) # output cavity + +# I-O cascade +G_cas_1m = ▷(G_u, G_2lvls, G_v) + +# Hamiltonian and Lindbladian +Hcas_1m = hamiltonian(G_cas_1m) +Jcas_1m = jump_operator(G_cas_1m)[1] +Jcas_1md = adjoint(Jcas_1m) + +# Time-dependent couplings +@register_symbolic gu_t(t) +@register_symbolic gv_t(t) +# Complex conjugates +@register_symbolic gu_c_t(t) +@register_symbolic gv_c_t(t) + +gs_ls = [gu, gv] +gst_ls = [gu_t(t), gv_t(t)] +gsct_ls = [gu_c_t(t), gv_c_t(t)] + +dict_gts = Dict([gs_ls; conj.(gs_ls)] .=> [gst_ls; gsct_ls]); + +# Insert time-dependence +Hcas_1m_t = substitute(Hcas_1m, dict_gts); +Jcas_1m_t = substitute(Jcas_1m, dict_gts); +Jcas_1md_t = substitute(Jcas_1md, dict_gts); + +# Pulses and composite mode +αeff = 0.0 + +Ω1(t_) = α1_/(sqrt(2π*σ1_^2)) * exp(-(t_-τ1_)^2/(2σ1_^2)) * exp(-1im*Δ1_*t_) + +Ω2(t_) = α2_/(sqrt(2π*σ2_^2)) * exp(-(t_-τ2_)^2/(2σ2_^2)) * exp(-1im*Δ2_*t_) + +Ωall(t_) = Ω1(t_) + Ω2(t_) + +αeff = sqrt(cumul_integrate(T, abs2.(Ωall.(T)))[end]) + +u(t_) = Ωall(t_) / αeff + +gu_t_ = coupling_input(u, T) +gu_t(t) = gu_t_(t) +gv_t_ = coupling_output(u, T) +gv_t(t) = gv_t_(t) + + +gu_c_t(t) = conj(gu_t(t)) +gv_c_t(t) = conj(gv_t(t)) + + +order = 2 +ops1 = [aun, aun', σ(1, 2), σ(2, 2), σ(1, 2)', avn, avn'] +ops_all1 = unique([ops1[i]*ops1[j] for i = 1:length(ops1) for j = 1:length(ops1)]); +ops_all = union( + ops1, + ops_all1[1:1], + ops_all1[3:14], + ops_all1[16:16], + ops_all1[18:26], + ops_all1[28:29], +) +eqs_1m = meanfield( + ops_all, + Hcas_1m_t, + [Jcas_1m_t]; + Jdagger = [Jcas_1md_t], + order = order, + iv = t, +) + +complete!(eqs_1m) +length(eqs_1m) + +# Classical-to-coherent relation +α_1m = αeff/(2*sqrt(γ_)) # may change if cross-terms are significant + +function ψ0_us1v(α, eqs) + bu = FockBasis(rd(abs2(α) + 20abs(α)), rd(abs2(α) - 20abs(α))) # include ±20α photon fluctuations + bs1 = NLevelBasis(2) + bv = FockBasis(1) + b = tensor([bu, bs1, bv]...) + ψu = coherentstate(bu, α) + ψs1 = nlevelstate(bs1, 1) + ψv = fockstate(bv, 0) + ψ0 = LazyKet(b, (ψu, ψs1, ψv)) + return initial_values(eqs, ψ0) +end + +# Initial values +u0_1m = ψ0_us1v(α_1m, eqs_1m) + +va = 0.05 + +# Solve ODE system +sys_1m = mtkcompile(System(eqs_1m; name = :sys_1m)) +ps = [γ] +p0 = [γ_] +u0maps = Dict(unknowns(sys_1m) .=> u0_1m) +pmaps = parameter_map(sys_1m, ps, p0) +probs = ODEProblem(sys_1m, merge(u0maps, pmaps), (dt, Tend)) +sol_1m0v = solve(probs, Tsit5(), saveat = va; abstol, reltol); + + +################################################### +########### Correlation Part ###################### +################################################### + +# Lindblad vector (includes output) +Jvec=[aun, σ(1, 2), avn] + +N = length(sol_1m0v.t) +Gmat_bar = zeros(ComplexF64, N, length(Jvec), length(Jvec), N) +t1 = sol_1m0v.t; +t2 = sol_1m0v.t; +τ0 = 0.0 + +Corrmat = [ + [ + CorrelationFunction(Jvec[k]', Jvec[l], eqs_1m; steady_state = false) for + l = 1:length(Jvec) + ] for k = 1:length(Jvec) +] +Corrvec = [Corrmat[end][i] for i = 1:length(Jvec)] + +oper_tup = (Jvec[1], Jvec[2], Jvec[3]) +idx_tup = idx_tuple(Jvec, Corrvec, Corrmat) + +# Optionally compute and save +#Gmat_bar = compute_Gmat( +# sol_1m0v, eqs_1m, Jvec, Corrvec; +# p0 = p0, +# ps = ps, +# saveat = va, +# τ0 = 0.0, +# alg = Tsit5(), +# state_syms = oper_tup, +# extract_idxs = idx_tup +#) + +# Load precomputed data +#JLD2.save_object("Corr_v.jld2",[Gmat_bar,t1]) +Gmat_bar = JLD2.load_object("Corr_v.jld2")[1]; +t1 = JLD2.load_object("Corr_v.jld2")[2]; +N = length(t1) +gcorr1 = zeros(ComplexF64, N, N) +gcorr = zeros(ComplexF64, N, N) + +gvec(t) = [gu_t(t), sqrt(γ_), gv_t(t)] +gvecc(t) = [gu_c_t(t), sqrt(γ_), gv_c_t(t)] + +for i1 = 1:N + for j1 = 1:N + gcorr1[i1, j1] = sum( + gvecc(t1[i1])[k]*gvec(t1[j1])[l]*(Gmat_bar[j1, k, l, i1]) for k = 1:3 for + l = 1:3 + ) + end +end +gcorr1 +Λ = zeros(ComplexF64, N, N) +for i = 1:N + Λ[i, i] = diag(real.(gcorr1))[i] +end + +# Final g^(1)(t1,t2) +gcorr = gcorr1 + tril(gcorr1, 1)' - Λ + + + +################################################### +########### Plot I-O SUPER Corr. loss ############# +################################################### + +#################################################################### +# Heat plot +#################################################################### + +g_abs = abs.(gcorr) +g_max = maximum(g_abs) +g_vis = g_abs ./ (g_max > 0 ? g_max : 1) + +tick_vals = collect(0.0:0.2:1.0) +tick_labs = latexstring.(string.(round.(tick_vals; digits = 1))) +pl8 = heatmap( + t1, + t1, + g_vis; + color = :thermal, + clims = (0, 1), + xlabel = L"\gamma t_2", + ylabel = L"\gamma t_1", + colorbar_title = L"|g^{(1)}|\ \mathrm{(norm)}", + colorbar_ticks = (tick_vals, tick_labs), + aspect_ratio = :equal, + tickfontsize = 16, + guidefontsize = 18, + size = (600, 500), +) +display(pl8) +# savefig(pl8, "heat_IO_loss.pdf") + + +##################################################################### +# Eigenvalues and eigenmodes +##################################################################### + +F = eigen(gcorr); + +ΔT = t1[2]-t1[1] +ns_avg = F.values +modess = F.vectors +vs_mode = (modess[:, end]) / sqrt(ΔT) +ns_avgt = ns_avg*(t1[2] - t1[1]) +ns = ns_avgt[end] +ns_avgt + +################################################### +########### Plot I-O SUPER Cascade mode ########### +################################################### + +# Dominant mode magnitude (loss-inclusive) +pl9 = plot( + t1, + abs.(u.(t1)); + color = :blue, + label = L"v(t) = u(t)", + ylabel = L"|v(t)|", + xlims = (t1[1]-0.01, t1[end]), + xticks = ([0, 10, 20], ["", "", ""]), + legend = :topleft, + tickfontsize = plot_font, + guidefontsize = plot_font, + legendfontsize = plot_font, + size = (800, 200), +) +plot!(pl9, t1, abs.(vs_mode); label = L"\mathrm{loss}") +display(pl9) +# savefig(pl9, "mode_IO_loss.pdf") diff --git a/examples/drafts/11-1_two-level-system_spontaneous-emission.jl b/examples/drafts/11-1_two-level-system_spontaneous-emission.jl index af631c51..0bbab00e 100644 --- a/examples/drafts/11-1_two-level-system_spontaneous-emission.jl +++ b/examples/drafts/11-1_two-level-system_spontaneous-emission.jl @@ -36,7 +36,7 @@ G_v = SLH(1, g_v' * av, 0) G = cascade(G_s, G_v) H = hamiltonian(G) -L = jump_operator(G)[1] +J = jump_operator(G)[1] nothing # hide # We use the spontaneous-emission rate as the frequency unit. At @@ -67,18 +67,18 @@ nothing # hide # # ```math # g^{(1)}(t_1,t_2) = -# \langle L_s^\dagger(t_1)L_s(t_2)\rangle, -# \qquad L_s = \sqrt{\gamma}\,\sigma_-. +# \langle J_s^\dagger(t_1)J_s(t_2)\rangle, +# \qquad J_s = \sqrt{\gamma}\,\sigma_-. # ``` dict_p_1 = Dict(γ => γ_, g_v => 0.0) H_QO_1 = to_numeric(H, b; parameter = dict_p_1) -L_QO_1 = to_numeric(L, b; parameter = dict_p_1) +J_QO_1 = to_numeric(J, b; parameter = dict_p_1) -t_1, ρt_1 = timeevolution.master(T, ψ0, H_QO_1, [L_QO_1]) +t_1, ρt_1 = timeevolution.master(T, ψ0, H_QO_1, [J_QO_1]) -Ls = √(γ_) * σm_qo -g1 = correlation_matrix(T, ρt_1, H_QO_1, [L_QO_1], Ls) +Js = √(γ_) * σm_qo +g1 = correlation_matrix(T, ρt_1, H_QO_1, [J_QO_1], Js) # We diagonalize the quadrature-weighted integral kernel. Its eigenvalues are # the mean photon numbers in the corresponding temporal modes. `Hermitian` @@ -116,11 +116,11 @@ gv_t = coupling_output(v_mode_interpolation, T_2) dict_p_2 = Dict(γ => γ_) dict_p_t_2 = Dict(g_v => gv_t) H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +J_QO_2 = to_numeric(J, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return Ht, J, dagger.(J) end diff --git a/src/SLH.jl b/src/SLH.jl index 65cea3b3..bc4b3429 100644 --- a/src/SLH.jl +++ b/src/SLH.jl @@ -40,23 +40,23 @@ _mul(x, y) = _isunit(x) ? y : _isunit(y) ? x : x * y # ────────────────────────────────────────────── """ - _detect_operator_type(L, H) + _detect_operator_type(J, H) Determine the concrete return type for FunctionWrapper by inspecting -the elements of L and H. Checks FunctionWrapper type parameters first +the elements of J 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). """ -function _detect_operator_type(L, H) +function _detect_operator_type(J, H) # Check FunctionWrapper elements (carry explicit type info) - for l in L - l isa FunctionWrapper && return _fw_return_type(typeof(l)) + for j in J + j isa FunctionWrapper && return _fw_return_type(typeof(j)) 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) + for j in J + _is_time_dep(j) || return typeof(j) end _is_time_dep(H) || return typeof(H) # No type information available @@ -68,12 +68,12 @@ end _fw_return_type(::Type{FunctionWrapper{R,A}}) where {R,A} = R -function _maybe_wrap_jump_operators(L::SVector{N}, ::Type{OpType}) where {N,OpType} - if any(_is_time_dep, L) +function _maybe_wrap_jump_operators(J::SVector{N}, ::Type{OpType}) where {N,OpType} + if any(_is_time_dep, J) fw_type = FunctionWrapper{OpType,Tuple{Float64}} - return SVector{N,fw_type}(ntuple(i -> fw_type(_to_func(L[i])), Val(N))) + return SVector{N,fw_type}(ntuple(i -> fw_type(_to_func(J[i])), Val(N))) end - return L + return J end function _maybe_wrap_hamiltonian(H, has_td::Bool, ::Type{OpType}) where {OpType} @@ -88,23 +88,23 @@ end # ────────────────────────────────────────────── """ - SLH{N, ST, LT, HT} + SLH{N, ST, JT, HT} -SLH triple with scattering matrix `S`, jump-operator vector `L`, and Hamiltonian `H`. -`S` and `L` can also be vectors of scattering matrices and jump operators. +SLH triple with scattering matrix `S`, jump-operator vector `J`, and Hamiltonian `H`. +`S` and `J` can also be vectors of scattering matrices and jump operators. See also [`▷`](@ref), [`⊞`](@ref), [`feedback`](@ref) """ -struct SLH{N,ST,LT,HT,L} +struct SLH{N,ST,JT,HT,L} scattering::SMatrix{N,N,ST,L} - jump_operator::SVector{N,LT} + jump_operator::SVector{N,JT} hamiltonian::HT - function SLH{N,ST,LT,HT}( + function SLH{N,ST,JT,HT}( S::SMatrix{N,N,ST,L}, - jump_operator::SVector{N,LT}, + jump_operator::SVector{N,JT}, H::HT, - ) where {N,ST,LT,HT,L} - return new{N,ST,LT,HT,L}(S, jump_operator, H) + ) where {N,ST,JT,HT,L} + return new{N,ST,JT,HT,L}(S, jump_operator, H) end end @@ -114,8 +114,8 @@ end Extract the operator return type from an SLH with FunctionWrapper elements. Returns `nothing` if no FunctionWrapper type info is available. """ -function _op_type(::SLH{N,ST,LT,HT}) where {N,ST,LT,HT} - LT <: FunctionWrapper && return _fw_return_type(LT) +function _op_type(::SLH{N,ST,JT,HT}) where {N,ST,JT,HT} + JT <: FunctionWrapper && return _fw_return_type(JT) HT <: FunctionWrapper && return _fw_return_type(HT) return nothing end @@ -125,55 +125,55 @@ end # ────────────────────────────────────────────── # 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) +function _build_slh(S::SMatrix{N,N}, J::SVector{N}, H) where {N} + has_td = any(_is_time_dep, J) || _is_time_dep(H) if has_td - OpType = _detect_operator_type(L, H) - return _build_slh(S, L, H, OpType) + OpType = _detect_operator_type(J, H) + return _build_slh(S, J, H, OpType) end - return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) + return SLH{N,eltype(S),eltype(J),typeof(H)}(S, J, H) end # 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) +function _build_slh(S::SMatrix{N,N}, J::SVector{N}, H, ::Type{OpType}) where {N,OpType} + has_td = any(_is_time_dep, J) || _is_time_dep(H) if has_td - L_w = _maybe_wrap_jump_operators(L, OpType) + J_w = _maybe_wrap_jump_operators(J, 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) + return SLH{N,eltype(S),eltype(J_w),typeof(H_w)}(S, J_w, H_w) end - return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) + return SLH{N,eltype(S),eltype(J),typeof(H)}(S, J, H) end # Nothing hint falls through to detection -_build_slh(S::SMatrix{N,N}, L::SVector{N}, H, ::Nothing) where {N} = _build_slh(S, L, H) +_build_slh(S::SMatrix{N,N}, J::SVector{N}, H, ::Nothing) where {N} = _build_slh(S, J, H) # From AbstractMatrix + AbstractVector (includes SMatrix + SVector) -function SLH(S::AbstractMatrix, L::AbstractVector, H) - N = length(L) +function SLH(S::AbstractMatrix, J::AbstractVector, H) + N = length(J) @assert size(S, 1) == N && size(S, 2) == N - return _build_slh(SMatrix{N,N}(S), SVector{N}(L...), H) + return _build_slh(SMatrix{N,N}(S), SVector{N}(J...), H) end -# Numeric scalar S + vector L → S * I_{NxN} -function SLH(S::Number, L::AbstractVector, H) - N = length(L) +# Numeric scalar S + vector J → S * I_{NxN} +function SLH(S::Number, J::AbstractVector, H) + N = length(J) S_mat = SMatrix{N,N}(S * LinearAlgebra.I) - return _build_slh(S_mat, SVector{N}(L...), H) + return _build_slh(S_mat, SVector{N}(J...), H) end -# Scalar S + scalar L → SLH{1} -function SLH(S, L, H) +# Scalar S + scalar J → SLH{1} +function SLH(S, J, H) S_mat = SMatrix{1,1}(S) - L_vec = SVector{1}(L) - return _build_slh(S_mat, L_vec, H) + J_vec = SVector{1}(J) + return _build_slh(S_mat, J_vec, H) end -# Symbolic/general scalar S + vector L → S * I -function SLH(S, L::AbstractVector, H) - N = length(L) +# Symbolic/general scalar S + vector J → S * I +function SLH(S, J::AbstractVector, H) + N = length(J) 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) + return _build_slh(S_mat, SVector{N}(J...), H) end # ────────────────────────────────────────────── @@ -190,7 +190,7 @@ scattering(G::SLH) = G.scattering """ jump_operator(G::SLH) -Return the jump-operator vector `L` of an SLH object. +Return the jump-operator vector `J` of an SLH object. """ jump_operator(G::SLH) = G.jump_operator @@ -226,18 +226,18 @@ end # Matrix-vector helpers # ────────────────────────────────────────────── -@generated function _slh_matvec(S::SMatrix{N,N}, L::SVector{N}) where {N} +@generated function _slh_matvec(S::SMatrix{N,N}, J::SVector{N}) where {N} if N == 1 - return :(SVector{1}(_mul(S[1, 1], L[1]))) + return :(SVector{1}(_mul(S[1, 1], J[1]))) end exprs = [] for i = 1:N first_name = Symbol("tmp_$(i)_1") - terms = [:($first_name = _mul(S[$i, 1], L[1]))] + terms = [:($first_name = _mul(S[$i, 1], J[1]))] acc = first_name for j = 2:N tname = Symbol("tmp_$(i)_$(j)") - push!(terms, :($tname = _add($acc, _mul(S[$i, $j], L[$j])))) + push!(terms, :($tname = _add($acc, _mul(S[$i, $j], J[$j])))) acc = tname end push!(exprs, Expr(:block, terms..., acc)) @@ -245,13 +245,13 @@ end return :(SVector($(exprs...))) end -@generated function _slh_dot(L1::SVector{N}, L2::SVector{N}) where {N} +@generated function _slh_dot(J1::SVector{N}, J2::SVector{N}) where {N} if N == 1 - return :(_mul(L1[1], L2[1])) + return :(_mul(J1[1], J2[1])) end - expr = :(_mul(L1[1], L2[1])) + expr = :(_mul(J1[1], J2[1])) for i = 2:N - expr = :(_add($expr, _mul(L1[$i], L2[$i]))) + expr = :(_add($expr, _mul(J1[$i], J2[$i]))) end return expr end @@ -265,27 +265,27 @@ end Cascade two SLH triples: -``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))`` +``G_1 \\triangleright G_2 = (S_2 S_1,\\; J_2 + S_2 J_1,\\; H_1 + H_2 - \\tfrac{i}{2}(J_2^\\dagger S_2 J_1 - J_1^\\dagger S_2^\\dagger J_2))`` Unicode `\\triangleright`. See also [`cascade`](@ref). """ function ▷(G1::SLH{N}, G2::SLH{N}) where {N} - S1, L1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) - S2, L2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) + S1, J1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) + S2, J2, H2 = scattering(G2), jump_operator(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))) + S2J1 = _slh_matvec(S2, J1) + J_t = SVector{N}(ntuple(i -> _post(_add(J2[i], S2J1[i])), Val(N))) - L2_adj = SVector{N}(ntuple(i -> _adj(L2[i]), Val(N))) - cross1 = _slh_dot(L2_adj, S2L1) + J2_adj = SVector{N}(ntuple(i -> _adj(J2[i]), Val(N))) + cross1 = _slh_dot(J2_adj, S2J1) X = _mul(-1im / 2, cross1) H_t = _post(_add(_add(H1, H2), _add(X, _adj(X)))) op_hint = _op_type(G1) op_hint === nothing && (op_hint = _op_type(G2)) - return _build_slh(S_t, L_t, H_t, op_hint) + return _build_slh(S_t, J_t, H_t, op_hint) end function ▷(::SLH{N1}, ::SLH{N2}) where {N1,N2} @@ -315,7 +315,7 @@ cascade(args...) = ▷(args...) 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)`` +\\begin{pmatrix} J_1 \\\\ J_2 \\end{pmatrix},\\; H_1 + H_2\\right)`` Unicode `\\boxplus`. See also [`concatenate`](@ref). """ @@ -333,14 +333,14 @@ Unicode `\\boxplus`. See also [`concatenate`](@ref). end quote - S1, L1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) - S2, L2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) + S1, J1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) + S2, J2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) S_t = SMatrix{$N,$N}($(s_exprs...)) - L_t = vcat(L1, L2) + J_t = vcat(J1, J2) 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) + return _build_slh(S_t, J_t, H_t, op_hint) end end @@ -368,8 +368,8 @@ function _drop_row_col(S::SMatrix{N,N}, row::Int, col::Int, ::Val{M}) where {N,M end) end -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))) +function _drop_index(J::SVector{N}, idx::Int, ::Val{M}) where {N,M} + SVector{M}(ntuple(i -> J[i >= idx ? i + 1 : i], Val(M))) end function _get_col_dropped_row( @@ -404,7 +404,7 @@ end function _feedback_impl(G::SLH{N}, x::Int, y::Int, ::Val{M}) where {N,M} S = scattering(G) - L = jump_operator(G) + J = jump_operator(G) H = hamiltonian(G) @assert 1 <= x <= N && 1 <= y <= N @@ -426,18 +426,18 @@ function _feedback_impl(G::SLH{N}, x::Int, y::Int, ::Val{M}) where {N,M} _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)) + J_bar = _drop_index(J, x, valM) + J_x = J[x] + J_update = + SVector{M}(ntuple(i -> _post(_mul(_mul(S_col_y_no_x[i], loop_gain), J_x)), valM)) + J_red = SVector{M}(ntuple(i -> _post(_add(J_bar[i], J_update[i])), valM)) 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)) + J_adj = SVector{N}(ntuple(i -> _adj(J[i]), Val(N))) + term = _mul(_slh_dot(J_adj, S_col_y_full), _mul(loop_gain, J_x)) H_red = _post(_add(H, _mul(1 / (2im), _add(term, _mul(-1, _adj(term)))))) - return _build_slh(S_red, L_red, H_red, _op_type(G)) + return _build_slh(S_red, J_red, H_red, _op_type(G)) end feedback(G::SLH, connection::Pair{Int,Int}) = @@ -480,13 +480,13 @@ end Translate the Hamiltonian and Lindblad operators of an SLH object `G` into numeric [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators on the basis `b`. -Returns the tuple `(H_QO, L_QO)`, where `L_QO` is a vector holding one translated operator +Returns the tuple `(H_QO, J_QO)`, where `J_QO` is a vector holding one translated operator per jump operator in `jump_operator(G)`. All keyword arguments (`parameter`, `time_parameter`, `operators`, `adjoint_ops`, `op_type`) are forwarded to [`SecondQuantizedAlgebra.to_numeric`](@ref). """ function SQA.to_numeric(G::SLH, b::QuantumOpticsBase.Basis; kwargs...) H_QO = SQA.to_numeric(hamiltonian(G), b; kwargs...) - L_QO = [SQA.to_numeric(L_, b; kwargs...) for L_ in jump_operator(G)] - return H_QO, L_QO + J_QO = [SQA.to_numeric(J_, b; kwargs...) for J_ in jump_operator(G)] + return H_QO, J_QO end diff --git a/src/correlations.jl b/src/correlations.jl index 3ea313a7..80f1a05c 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -1,20 +1,20 @@ """ - correlation_matrix(T, ρt, f, Ls; kwargs...) - correlation_matrix(T, ρt, H, J, Ls; kwargs...) + correlation_matrix(T, ρt, f, Js; kwargs...) + correlation_matrix(T, ρt, H, J, Js; kwargs...) Compute the two-time correlation matrix -``g^{(1)}(t_1, t_2) = \\langle L_s^\\dagger(t_1) L_s(t_2) \\rangle`` +``g^{(1)}(t_1, t_2) = \\langle J_s^\\dagger(t_1) J_s(t_2) \\rangle`` on the time grid `T`. Writes directly into output matrix. Supply the dynamics either as a `master_dynamic`-style function `f(t, ρ)`, or as operators passed straight to the solver: a time-dependent `H` (e.g. the `TimeDependentSum` from [`to_numeric`](@ref)) with jump operators `J`, or a constant `H` with constant `J`. The operator form is much faster for time-dependent problems (the integrator is built once). -`Ls` is either a constant operator or a function `Ls(t)` returning the operator at `t`. +`Js` is either a constant operator or a function `Js(t)` returning the operator at `t`. """ -function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) - Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) - _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 +function correlation_matrix(T::Vector, ρt::Vector, f::Function, Js; kwargs...) + Js_vec, Js_dag_vec = _sample_operator_and_adjoint(T, Js) + _correlation_loop(T, ρt, Js_vec, Js_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) end end @@ -24,18 +24,18 @@ function correlation_matrix( ρt::Vector, H::QuantumOpticsBase.AbstractTimeDependentOperator, J::AbstractVector, - Ls; + Js; kwargs..., ) - Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) - _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 + Js_vec, Js_dag_vec = _sample_operator_and_adjoint(T, Js) + _correlation_loop(T, ρt, Js_vec, Js_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, copy(H), [copy(j) for j in J]; kwargs...) end end -function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Ls; kwargs...) - Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) - _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 +function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Js; kwargs...) + Js_vec, Js_dag_vec = _sample_operator_and_adjoint(T, Js) + _correlation_loop(T, ρt, Js_vec, Js_dag_vec) do T_slice, ρ0 timeevolution.master(T_slice, ρ0, H, J; kwargs...) end end @@ -50,8 +50,8 @@ function _sample_operator_and_adjoint( ) throw( ArgumentError( - "`Ls` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " * - "function `Ls(t)` returning the concrete operator at time `t`.", + "`Js` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " * + "function `Js(t)` returning the concrete operator at time `t`.", ), ) end @@ -60,18 +60,18 @@ function _sample_operator_and_adjoint(T::Vector, op) return fill(op, nt), fill(dagger(op), nt) end -function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_dag_vec) +function _correlation_loop(solve_fn, T, ρt, Js_vec, Js_dag_vec) l_T = length(T) @assert l_T == length(ρt) g1_m = zeros(ComplexF64, l_T, l_T) # Each iteration solves an independent master equation — parallelise Threads.@threads for it = 1:(l_T-1) - ρ0_it = Ls_vec[it] * ρt[it] + ρ0_it = Js_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]) + val = expect(Js_dag_vec[it+i-1], ρ_bar_τ[i]) g1_m[it, it+i-1] = val g1_m[it+i-1, it] = conj(val) end diff --git a/test/test_SLH.jl b/test/test_SLH.jl index 1ac400e7..f8218047 100644 --- a/test/test_SLH.jl +++ b/test/test_SLH.jl @@ -25,16 +25,16 @@ using Test G_v = SLH(1, gv'*av, 0) # output cavity G_c_S = scattering(G_c) - G_c_L = jump_operator(G_c) + G_c_J = jump_operator(G_c) G_c_H = hamiltonian(G_c) @test G_c_S == G_c.scattering - @test G_c_L == G_c.jump_operator + @test G_c_J == G_c.jump_operator @test G_c_H == G_c.hamiltonian @test_deprecated lindblad(G_c) @test G_c_S isa SMatrix{1,1} - @test G_c_L isa SVector{1} + @test G_c_J isa SVector{1} SLH(1, [√(γ)*c], Δ*c'c) @test isequal(G_c, SLH(1, [√(γ)*c], Δ*c'c)) @@ -116,13 +116,13 @@ using Test 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) + J_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(jump_operator(G)) === typeof(L_s) + G = SLH(1, J_s, H_s) + @test eltype(jump_operator(G)) === typeof(J_s) @test typeof(hamiltonian(G)) === typeof(H_s) end @@ -142,30 +142,30 @@ using Test end @testset "cascade mixed static/time-dep wraps uniformly" begin - G_cas = SLH(1, L_s, H_s) ▷ SLH(1, gu_f, H_s) + G_cas = SLH(1, J_s, H_s) ▷ SLH(1, gu_f, H_s) @test eltype(jump_operator(G_cas)) <: FunctionWrapper @test eltype(jump_operator(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(jump_operator(G_cat)) - @test LT <: FunctionWrapper - @test LT !== Any + G_cat = SLH(1, J_s, H_s) ⊞ SLH(1, gu_f, H_s) + JT = eltype(jump_operator(G_cat)) + @test JT <: FunctionWrapper + @test JT !== Any @inferred jump_operator(G_cat)[1](0.5) @inferred jump_operator(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(jump_operator(G_cat)) === typeof(L_s) + G_cat = SLH(1, J_s, H_s) ⊞ SLH(1, J_s, H_s) + @test eltype(jump_operator(G_cat)) === typeof(J_s) @test !(eltype(jump_operator(G_cat)) <: FunctionWrapper) end @testset "FunctionWrapper call is inferred" begin G_td = SLH(1, gu_f, H_s) - l = jump_operator(G_td)[1] - @inferred l(0.5) + j = jump_operator(G_td)[1] + @inferred j(0.5) end @testset "_op_type extracts type from FunctionWrapper SLH" begin @@ -174,7 +174,7 @@ using Test end @testset "_op_type returns nothing for static SLH" begin - G_s = SLH(1, L_s, H_s) + G_s = SLH(1, J_s, H_s) @test QuantumInputOutput._op_type(G_s) === nothing end diff --git a/test/test_compare_example_05_1_05_2.jl b/test/test_compare_example_05_1_05_2.jl index d36d5c6a..09246e90 100644 --- a/test/test_compare_example_05_1_05_2.jl +++ b/test/test_compare_example_05_1_05_2.jl @@ -39,9 +39,9 @@ using Test G_t = G_R_t ⊞ G_L_t H = hamiltonian(G_t) - L = jump_operator(G_t) - L_R = L[1] - L_L = L[2] + Jop = jump_operator(G_t) + J_R = Jop[1] + J_L = Jop[2] γ_ = 1.0 β = 0.9 @@ -59,23 +59,23 @@ using Test ba = NLevelBasis(2) b = tensor([ba for _ = 1:N]...) H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) - L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) - L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) + J_R_QO = to_numeric(J_R, b; parameter = dict_p, time_parameter = dict_p_t) + J_L_QO = to_numeric(J_L, b; parameter = dict_p, time_parameter = dict_p_t) σ_qo(α, i, j) = to_numeric(σ(α, i, j), b) J_add = [√(γ_add[i]) * σ_qo(i, 1, 2) for i = 1:N] function input_output_1(t, ρ) Ht = H_QO(t) - J = [L_R_QO(t), L_L_QO(t), J_add...] + J = [J_R_QO(t), J_L_QO(t), J_add...] return Ht, J, dagger.(J) end ψ0 = tensor([nlevelstate(ba, 1) for _ = 1:N]...) t1, ρt1 = timeevolution.master_dynamic(T, ψ0, input_output_1) - 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)] + I_R_1 = [real(expect(J_R_QO(ti)' * J_R_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] + I_L_1 = [real(expect(J_L_QO(ti)' * J_L_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] # -------- Example 05-2 style (numeric SLH, quantum pulse) -------- bu = FockBasis(4) @@ -96,29 +96,29 @@ using Test G_t_qo = G_R_t_qo ⊞ G_L_t_qo H_qo = hamiltonian(G_t_qo) - L_qo = jump_operator(G_t_qo) - L_R_qo = L_qo[1] - L_L_qo = L_qo[2] + J_qo = jump_operator(G_t_qo) + J_R_qo = J_qo[1] + J_L_qo = J_qo[2] # 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_R_f = _callable(J_R_qo) ? J_R_qo : (t -> J_R_qo) + J_L_f = _callable(J_L_qo) ? J_L_qo : (t -> J_L_qo) J_add_qo = [√(γ_add[i]) * σ_qds(i, 1, 2) for i = 1:N] function input_output_2(t, ρ) Ht = Hf(t) - J = [L_R_f(t), L_L_f(t), J_add_qo...] + J = [J_R_f(t), J_L_f(t), J_add_qo...] return Ht, J, dagger.(J) end ψ0_qo = coherentstate(bu, α0) ⊗ tensor([nlevelstate(ba, 1) for _ = 1:N]...) t2, ρt2 = timeevolution.master_dynamic(T, ψ0_qo, input_output_2) - I_R_2 = [real(expect(L_R_f(ti)' * L_R_f(ti), ρt2[i])) for (i, ti) in enumerate(t2)] - I_L_2 = [real(expect(L_L_f(ti)' * L_L_f(ti), ρt2[i])) for (i, ti) in enumerate(t2)] + I_R_2 = [real(expect(J_R_f(ti)' * J_R_f(ti), ρt2[i])) for (i, ti) in enumerate(t2)] + I_L_2 = [real(expect(J_L_f(ti)' * J_L_f(ti), ρt2[i])) for (i, ti) in enumerate(t2)] @test maximum(abs.(I_R_1 .- I_R_2)) < 5e-5 @test maximum(abs.(I_L_1 .- I_L_2)) < 5e-5 diff --git a/test/test_correlations.jl b/test/test_correlations.jl index eb4ca207..81c30934 100644 --- a/test/test_correlations.jl +++ b/test/test_correlations.jl @@ -12,7 +12,7 @@ using Test H = 0.0 * one(ba) J = [sqrt(γ) * σm] - Ls = J[1] + Js = J[1] T_end = 10 / γ T = collect(0.0:0.01:T_end) @@ -23,10 +23,10 @@ using Test end _, ρt_dyn = timeevolution.master_dynamic(T, ψ0, f_const) - g1_dyn = correlation_matrix(T, ρt_dyn, f_const, Ls) + g1_dyn = correlation_matrix(T, ρt_dyn, f_const, Js) _, ρt_static = timeevolution.master(T, ψ0, H, J) - g1_static = correlation_matrix(T, ρt_static, H, J, Ls) + g1_static = correlation_matrix(T, ρt_static, H, J, Js) @test maximum(abs.(g1_dyn .- g1_static)) < 1e-8 @@ -46,16 +46,16 @@ using Test overlap_aligned = abs(sum(conj.(mode_exp) .* mode_num_aligned) * ΔT) @test overlap_aligned > 0.99 - @testset "time-dependent Ls rejected" begin - td_Ls = TimeDependentSum([t -> 1.0 + 0im], J) - @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Ls) + @testset "time-dependent Js rejected" begin + td_Js = TimeDependentSum([t -> 1.0 + 0im], J) + @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Js) end @testset "time-dependent operator path matches static" begin H_td = TimeDependentSum([t -> 0.0 + 0im], [one(ba)]) J_td = [TimeDependentSum([t -> sqrt(γ) + 0im], [σm])] _, ρt_td = timeevolution.master_dynamic(T, ψ0, H_td, J_td) - g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Ls) + g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Js) @test maximum(abs.(g1_td .- g1_static)) < 1e-8 end end diff --git a/test/test_example_cavity_scattering.jl b/test/test_example_cavity_scattering.jl index d74c61ce..138b5918 100644 --- a/test/test_example_cavity_scattering.jl +++ b/test/test_example_cavity_scattering.jl @@ -24,7 +24,7 @@ using Test G_cas = ▷(G_u, G_c, G_v) H = hamiltonian(G_cas) - L = jump_operator(G_cas)[1] + Jop = jump_operator(G_cas)[1] γ_ = 1.0 Δ_ = 0.0 @@ -48,11 +48,11 @@ using Test b = bu1 ⊗ bc1 ⊗ bv1 H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) + J_QO = to_numeric(Jop, b; parameter = dict_p, time_parameter = dict_p_t) function input_output_1(t, ρ) Ht = H_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return Ht, J, dagger.(J) end @@ -70,8 +70,8 @@ using Test @test abs(sum(n_u_t .* abs2.(gu_t.(T)))*ΔT - 1) < 1e-2 # correlation matrix - Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo - g1_m = correlation_matrix(T, ρt, input_output_1, Ls) + Js(t) = gu_t(t)*au_qo + √(γ_)*c_qo + g1_m = correlation_matrix(T, ρt, input_output_1, Js) F = eigen(g1_m) n_avg = round.(real.(F.values)*ΔT; digits = 3) @@ -91,10 +91,10 @@ using Test dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]) H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) - L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) + J_QO_2 = to_numeric(Jop, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) H = H_QO_2(t) - J = [L_QO_2(t)] + J = [J_QO_2(t)] return H, J, dagger.(J) end diff --git a/test/test_feedback.jl b/test/test_feedback.jl index 024d73dd..5dac0e72 100644 --- a/test/test_feedback.jl +++ b/test/test_feedback.jl @@ -18,22 +18,22 @@ using Test # s12_r = rnumber("s12_r") # s21_r = rnumber("s21_r") # s22_r = rnumber("s22_r") - # l1_r = rnumber("l1_r") - # l2_r = rnumber("l2_r") + # j1_r = rnumber("j1_r") + # j2_r = rnumber("j2_r") # # # s11_i = rnumber("s11_i") # s12_i = rnumber("s12_i") # s21_i = rnumber("s21_i") # s22_i = rnumber("s22_i") - # l1_i = rnumber("l1_i") - # l2_i = rnumber("l2_i") + # j1_i = rnumber("j1_i") + # j2_i = rnumber("j2_i") # # # s11 = s11_r + 1im*s11_i # s12 = s12_r + 1im*s12_i # s21 = s21_r + 1im*s21_i # s22 = s22_r + 1im*s22_i - # l1 = l1_r + 1im*l1_i - # l2 = l2_r + 1im*l2_i + # j1 = j1_r + 1im*j1_i + # j2 = j2_r + 1im*j2_i # h0 = rnumber("h0") @@ -41,22 +41,22 @@ using Test s12 = 0.3 + 1im*0.24 s21 = 0.5 + 1im*0.23 s22 = 0.12 + 1im*0.29 - l1 = 0.18 + 1im*0.21 - l2 = 0.15 + 1im*0.25 + j1 = 0.18 + 1im*0.21 + j2 = 0.15 + 1im*0.25 h0 = 1.0 - G = SLH([s11 s12; s21 s22], [l1, l2], h0) + G = SLH([s11 s12; s21 s22], [j1, j2], h0) G_red = feedback(G, 1, 1) loop_gain = (1 - s11)^(-1) 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_J = simplify(j2 + s21 * loop_gain * j1) + expected_term = simplify((j1' * s11 + j2' * s21) * loop_gain * j1) expected_H = simplify(h0 + (expected_term - expected_term') / (2im)) @test scattering(G_red) isa SMatrix{1,1} @test abs(scattering(G_red)[1, 1] - expected_S) < 1e-10 - @test abs(jump_operator(G_red)[1] - expected_L) < 1e-10 + @test abs(jump_operator(G_red)[1] - expected_J) < 1e-10 @test abs(hamiltonian(G_red) - expected_H) < 1e-10 @testset "coherent-feedback OPO loop" begin @@ -118,16 +118,16 @@ using Test 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) + J_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(jump_operator(G_fb)) - @test LT <: FunctionWrapper - @test LT !== Any + JT = eltype(jump_operator(G_fb)) + @test JT <: FunctionWrapper + @test JT !== Any @inferred jump_operator(G_fb)[1](0.5) end end diff --git a/test/test_interaction_picture.jl b/test/test_interaction_picture.jl index 9db62b54..ba76b13c 100644 --- a/test/test_interaction_picture.jl +++ b/test/test_interaction_picture.jl @@ -39,7 +39,7 @@ using Test G_cas = ▷(G_u, G_s, G_v) H = hamiltonian(G_cas) - L = jump_operator(G_cas)[1] + Jop = jump_operator(G_cas)[1] H_uv = hamiltonian(▷(G_u, G_v)) H_int_sym_ = simplify(H - H_uv) @@ -51,7 +51,7 @@ using Test int_dict = Dict(a0_ls .=> a_int_ls) H_int_sym = simplify(substitute(H_int_sym_, int_dict)) - L_int_sym = simplify(substitute(L, int_dict)) + J_int_sym = simplify(substitute(Jop, int_dict)) # Virtual-cavity couplings gu_t = coupling_input(u, T) @@ -116,11 +116,11 @@ using Test dict_p_t = Dict(p_t_sym .=> p_t_num) H_int_QO = to_numeric(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = to_numeric(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) + J_QO = to_numeric(J_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) function input_output_I(t, ρ) Ht = H_int_QO(t) - J = [L_QO(t)] + J = [J_QO(t)] return Ht, J, dagger.(J) end @@ -140,11 +140,11 @@ using Test dict_p_t_s = Dict(gu_sym => gu_t, gv_sym => gv_t) H_QO = to_numeric(H, b; parameter = dict_p_s, time_parameter = dict_p_t_s) - L_QO_S = to_numeric(L, b; parameter = dict_p_s, time_parameter = dict_p_t_s) + J_QO_S = to_numeric(Jop, b; parameter = dict_p_s, time_parameter = dict_p_t_s) function input_output_S(t, ρ) Ht = H_QO(t) - J = [L_QO_S(t)] + J = [J_QO_S(t)] return Ht, J, dagger.(J) end @@ -158,12 +158,12 @@ using Test dict_p_small = Dict(γ_sym => γ_small) H_int_QO_small = to_numeric(H_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) - L_QO_small = - to_numeric(L_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) + J_QO_small = + to_numeric(J_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) function input_output_I_small(t, ρ) Ht = H_int_QO_small(t) - J = [L_QO_small(t)] + J = [J_QO_small(t)] return Ht, J, dagger.(J) end diff --git a/test/test_translate.jl b/test/test_translate.jl index 5fe1168c..fe33fed7 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -169,7 +169,7 @@ mat(F, t) = dense(static_operator(F(t))) @testset "to_numeric SLH overload" begin G_sym = SLH(1, [sqrt(κ_R) * a * E], Δ * a' * a) - H_QO, L_QO = to_numeric( + H_QO, J_QO = to_numeric( G_sym, bc1; parameter = dict_p1, @@ -178,10 +178,10 @@ mat(F, t) = dense(static_operator(F(t))) ) @test H_QO isa TimeDependentSum - @test length(L_QO) == 1 - @test L_QO[1] isa TimeDependentSum + @test length(J_QO) == 1 + @test J_QO[1] isa TimeDependentSum @test sum(abs.((mat(H_QO, 0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8 - @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < + @test sum(abs.((mat(J_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 end