diff --git a/TODO.md b/TODO.md index 4d503b83..6ba13c92 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 Js const (see example 02-2) + - [x] two_time_corr_matrix for Ls const (see example 02-2) - [x] Formatter/SpellCheck - [x] JET diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl index 9f9cfa71..17a44b6b 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) - J_sym = jump_operator(G_cas)[1] + L_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) - J_QO = to_numeric(J_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) ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0) - _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [J_QO]) + _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [L_QO]) au_qo = to_numeric(au, b) c_qo = to_numeric(c, b) - Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo + Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo ## --- Two-time correlation --- - # Pass the time-dependent operators (`H_QO`, `[J_QO]`) directly to the solver rather + # Pass the time-dependent operators (`H_QO`, `[L_QO]`) directly to the solver rather # than wrapping them in a `(t, ρ)` closure: since v0.10 `to_numeric` returns a lazy # `TimeDependentSum`, the direct path lets the solver build its integrator once and is # markedly faster than re-reading the operator from a function at every step. @@ -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, [$J_QO], $Js) + @benchmarkable correlation_matrix($T, $ρt, $H_QO, [$L_QO], $Ls) return nothing end diff --git a/benchmarks/interaction_picture.jl b/benchmarks/interaction_picture.jl index ae5596b3..c4833c9a 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) - J = jump_operator(G_cas)[1] + L = 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($J, $int_dict)) + simplify(substitute($L, $int_dict)) end return nothing diff --git a/benchmarks/slh_algebra.jl b/benchmarks/slh_algebra.jl index 2e35a91f..384427d5 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) - J_f = jump_operator(G_t) + L_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) - J_callables = [_callable(Ji) ? Ji : (t -> Ji) for Ji in J_f] + L_callables = [_callable(Li) ? Li : (t -> Li) for Li in L_f] SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide H(t)"] = @benchmarkable $Hf($t_mid) - SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide J(t)"] = @benchmarkable begin - for Ji in $J_callables - Ji($t_mid) + SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide L(t)"] = @benchmarkable begin + for Li in $L_callables + Li($t_mid) end end diff --git a/benchmarks/translation.jl b/benchmarks/translation.jl index 150cb71d..ffae277f 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 J from cascade + ## Derive H and L 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) - J_sym = jump_operator(G_cas)[1] + L_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 J translation + # Full cascade H and L 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+J"] = @benchmarkable begin + SUITE["Translation"]["time-dependent"]["3-cavity H+L"] = @benchmarkable begin to_numeric($H_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) + to_numeric($L_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) - J_QO = to_numeric(J_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) 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 J(t)"] = - @benchmarkable $J_QO($t_mid) + SUITE["Translation"]["closure evaluation"]["3-cavity L(t)"] = + @benchmarkable $L_QO($t_mid) return nothing end diff --git a/docs/src/implementation.md b/docs/src/implementation.md index 99427231..51809215 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, 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. +An SLH component is represented as `(S, L, H)` by the [`SLH`](@ref) type. The cascade [`▷`](@ref), concatenation [`⊞`](@ref), and feedback reduction [`feedback`](@ref) rules implement the standard network composition from the SLH framework. The resulting effective operators are accessed by [`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) -J = jump_operator(G_cas) +L = 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 `J` 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 `L` 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) -J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t) +L_QO = to_numeric(L, 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 J_s^\dagger(t_1) J_s(t_2) \rangle +g^{(1)}(t_1, t_2) = \langle L_s^\dagger(t_1) L_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 $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. +and diagonalizing it. In this package, [`correlation_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $L_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $L_s(t_1) \rho(t_1)$ and use this as the initial "state" for the propagation of $t_2$, with the same Hamiltonian and Lindblad terms. The eigenvectors of the matrix $g^{(1)}(t_1, t_2)$ correspond to temporal modes and the eigenvalues to their mean photon-number weights. The full procedure is illustrated in the [Tutorial](@ref). ```julia -Js(t) = gu_t(t) * au_qo + √(1.0) * c_qo -g1 = correlation_matrix(T, ρt, input_output_1, Js) +Ls(t) = gu_t(t) * au_qo + √(1.0) * c_qo +g1 = correlation_matrix(T, ρt, input_output_1, Ls) 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 acbae397..e1ba6c5f 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, J, H)`: +In the SLH formalism, each component is specified by a triple `(S, L, H)`: - `S`cattering matrix -- `J`ump (Lindblad) operators +- `L`indblad 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,\; -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) +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) \right). ``` @@ -113,27 +113,27 @@ S_1 & 0 \\ 0 & S_2 \end{bmatrix},\; \begin{bmatrix} -J_1 \\ -J_2 +L_1 \\ +L_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, J, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde J, \tilde H)$ with +$G = (S, L, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde L, \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 J = J_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}J_x, +\tilde L = L_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}L_x, ``` ```math \tilde H = H + \frac{1}{2i}\left[ -\left(\sum_{j=1}^n J_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}J_x - \mathrm{h.c.} +\left(\sum_{j=1}^n L_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}L_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, J, 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, L, 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 J_s^\dagger(t_1)\, \hat J_s(t_2) \rangle, +g^{(1)}(t_1, t_2) = \langle \hat L_s^\dagger(t_1)\, \hat L_s(t_2) \rangle, ``` -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)$). +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)$). The eigen-decomposition ```math diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 26556eed..9e6b9eb5 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -52,7 +52,7 @@ H = hamiltonian(G_cas) ``` ```@example tutorial -J = jump_operator(G_cas)[1] +L = 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) -J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t) +L_QO = to_numeric(L, b; parameter=dict_p, time_parameter=dict_p_t) nothing # hide ``` @@ -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 = [J_QO(t)] + J = [L_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 J_s^\dagger(t_1) J_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators. +``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. ```@example tutorial au_qo = to_numeric(au, b) c_qo = to_numeric(c, b) av_qo = to_numeric(av, b) -Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo -g1_m = correlation_matrix(T, ρt, input_output_1, Js) +Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo +g1_m = correlation_matrix(T, ρt, input_output_1, Ls) 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) -J_QO_2 = to_numeric(J, 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) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_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 85fddf91..3f7ed085 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) # -J = jump_operator(G_cas)[1] # only one Lindblad term in this example +L = 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) -J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) nothing # hide @@ -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 = [J_QO(t)] + J = [L_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 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$. +# 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$. -Js(t) = gu_t(t)*au_qo + √(γ_)*c_qo -g1_m = correlation_matrix(T, ρt, H_QO, [J_QO], Js) +Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo +g1_m = correlation_matrix(T, ρt, H_QO, [L_QO], Ls) nothing # hide # @@ -149,10 +149,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) -J_QO_2 = to_numeric(J, 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) function input_output_2(t, ρ) H = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_QO_2(t)] return H, J, dagger.(J) end; @@ -215,10 +215,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) -J_QO_3 = to_numeric(J, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) +L_QO_3 = to_numeric(L, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_3(t, ρ) H = H_QO_3(t) - J = [J_QO_3(t), √(γ_)*cdc_3] + J = [L_QO_3(t), √(γ_)*cdc_3] return H, J, dagger.(J) end; @@ -232,8 +232,8 @@ nothing # hide # -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)] +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)] 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 94966b3d..d9af02d1 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 # -J = G_cas.jump_operator[1] # only one Lindblad in this example +L = 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) -J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) H = H_QO(t) - J = [J_QO(t)] + J = [L_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)) -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_)] +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_)] 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 ae82e75a..4f9daf75 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) # -J = jump_operator(G_cas)[1] # only one jump operator in this example +L = 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) -J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) nothing # hide # We additionally include a cavity dephasing term and solve the dynamics. function input_output(t, ρ) H = H_QO(t) - J = [J_QO(t), √(γ_p_)*cdc_qo] + J = [L_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 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$. +# 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$. -Js(t) = (gu_(t))'*au_qo + √(γ_)*c_qo -g1_m = correlation_matrix(T, ρt, input_output, Js); +Ls(t) = (gu_(t))'*au_qo + √(γ_)*c_qo +g1_m = correlation_matrix(T, ρt, input_output, Ls); nothing # hide # diff --git a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl index a45a8083..92c5188f 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) -J = jump_operator(G)[1] +L = 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) -J_QO_1 = to_numeric(J, b; parameter = dict_p_1) +L_QO_1 = to_numeric(L, b; parameter = dict_p_1) function input_output_1(t, ρ) - J = [J_QO_1] + J = [L_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 -Js = √(γ_) * a_qo -g1_m = correlation_matrix(T, ρt_1, input_output_1, Js) +Ls = √(γ_) * a_qo +g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) 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) -J_QO_2 = to_numeric(J, 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) function input_output_2(t, ρ) - J = [J_QO_2(t)] + J = [L_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 97dfd53b..88ce932b 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) # -J = jump_operator(G_cas)[1] # only one Lindblad in this example +L = 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) -J_QO = to_numeric(J, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) H = H_QO(t) - J = [J_QO(t)] + J = [L_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 J_s^\dagger(t_1) J_s(t_2) \rangle$. +# 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$. au1_qo = to_numeric(au1, b) σ_qo(i, j) = to_numeric(σ(i, j), b) -Js(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) -g1_m = correlation_matrix(T, ρt, input_output, Js); +Ls(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) +g1_m = correlation_matrix(T, ρt, input_output, Ls); 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) -J_QO_2 = to_numeric(J, 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) function input_output_2(t, ρ) H = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_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 4ea18209..673a9df6 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) # -J1_L = jump_operator(G_cav_L_R_drive)[1] +L1_L = jump_operator(G_cav_L_R_drive)[1] # -J1_R = jump_operator(G_cav_L_R_drive)[2] +L1_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) -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] +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] nothing # hide # @@ -140,11 +140,11 @@ H2 = G_ac_drive.hamiltonian # -J2_L = G_ac_drive.jump_operator[1] +L2_L = G_ac_drive.jump_operator[1] # -J2_R = G_ac_drive.jump_operator[2] +L2_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) -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) +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) ## additional atomic decay into free space J_add = [√(γn)*σ_QO(α, 1, 2) for α = 1:Natoms] function input_output(t, ρ) H = H_QO(t) - J = [J2_L_QO(t), J2_R_QO, J_add...] + J = [L2_L_QO(t), L2_R_QO, J_add...] return H, J, dagger.(J) end nothing # hide @@ -208,13 +208,13 @@ nothing # hide # -J2_L_QO_dag(t) = dagger(J2_L_QO(t)) +L2_L_QO_dag(t) = dagger(L2_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(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])) + 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])) 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 544b50f1..117fd513 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) # -J1_L = jump_operator(G_cav_L_R_drive)[1] +L1_L = jump_operator(G_cav_L_R_drive)[1] # -J1_R = jump_operator(G_cav_L_R_drive)[2] +L1_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, [J1_L, J1_R]) +eqs_a = meanfield([a], H1, [L1_L, L1_R]) complete!(eqs_a) nothing # hide @@ -153,16 +153,16 @@ H2 = G_ac_drive.hamiltonian # -J2_L = G_ac_drive.jump_operator[1] +L2_L = G_ac_drive.jump_operator[1] # -J2_R = G_ac_drive.jump_operator[2] +L2_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, [J2_L, J2_R, J_add...]; order = 2) +eqs2 = meanfield([a'a, σ(1, 2, 2)], H2, [L2_L, L2_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 4bffd223..09b72af3 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) # -J = jump_operator(G_t) -J_R = J[1] +L = jump_operator(G_t) +L_R = L[1] # -J_L = J[2] +L_L = L[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) -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) +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) σ_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 = [J_R_QO(t), J_L_QO(t), J_add...] + J = [L_R_QO(t), L_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) - 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])) + LR = L_R_QO(ti) + LL = L_L_QO(ti) + I_R[i] = real(expect(LR'LR, ρt[i])) + I_L[i] = real(expect(LL'LL, ρt[i])) end nothing # hide @@ -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. -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)) +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)) for it1 = 1:(lT-1) ρ_t1 = ρt[it1] t_2, ρ_2 = timeevolution.master_dynamic( T[it1:end], - J0(T[it1]) * ρ_t1 * J0_dag(T[it1]), + L0(T[it1]) * ρ_t1 * L0_dag(T[it1]), input_output, ) ## transmission - G2_ls = real.([expect(J0_dag(t_2[j]) * J0(t_2[j]), ρ_2[j]) for j = 1:length(t_2)]) + G2_ls = real.([expect(L0_dag(t_2[j]) * L0(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], - J0_ref(T[it1]) * ρ_t1 * J0_ref_dag(T[it1]), + L0_ref(T[it1]) * ρ_t1 * L0_ref_dag(T[it1]), input_output, ) ## reflection G2_ls_r = real.([ - expect(J0_ref_dag(t_2_r[j]) * J0_ref(t_2_r[j]), ρ_2_r[j]) for j = 1:length(t_2_r) + expect(L0_ref_dag(t_2_r[j]) * L0_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 914da658..4ed14413 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 $J$ 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 $L$ will also be time-dependent. H = hamiltonian(G_t) -J = jump_operator(G_t) -J_R = J[1] -J_L = J[2] +L = jump_operator(G_t) +L_R = L[1] +L_L = L[2] J_add = [√(γ_add[i]) * σ(i, 1, 2) for i = 1:N] function input_output(t, ρ) Ht = H(t) - J = [J_R(t), J_L(t), J_add...] + J = [L_R(t), L_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) - 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])) + 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])) 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 804283a9..e2db6850 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) -J = jump_operator(G_t) -J_L = J[1] -J_R = J[2] +L = jump_operator(G_t) +L_L = L[1] +L_R = L[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) -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) +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) σ_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 = [J_R_QO(t), J_L_QO(t), J_add...] + J = [L_R_QO(t), L_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) - 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])) + LR = L_R_QO(ti) + LL = L_L_QO(ti) + I_R[i] = real(expect(LR' * LR, ρt[i])) + I_L[i] = real(expect(LL' * LL, ρt[i])) end nothing # hide 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 a1ee55ff..da518eb4 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) # -J = jump_operator(G_cas)[1] +L = 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)) # -J_int_sym = simplify(substitute(J, int_dict)) +L_int_sym = simplify(substitute(L, 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) -J_QO = to_numeric(J_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) function input_output(t, ρ) Ht = H_int_QO(t) - J = [J_QO(t)] + J = [L_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 e053ae8c..85eeaf9f 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) # -J = jump_operator(G) -J[1] +L = jump_operator(G) +L[1] # -J[2] +L[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) -J_QO = [to_numeric(Ji, b; parameter = dict_p, time_parameter = dict_p_t) for Ji in J] +L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] function input_output(t, ρ) Ht = H_QO(t) - J = [J_QO[1](t), J_QO[2](t)] + J = [L_QO[1](t), L_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 a0903740..03dbae19 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) # -J = jump_operator(G) -J[1] +L = jump_operator(G) +L[1] # -J[2] +L[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) -J_QO = [to_numeric(Ji, b; parameter = dict_p, time_parameter = dict_p_t) for Ji in J] +L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] function input_output(t, ρ) Ht = H_QO(t) - J = [J_QO[1](t), J_QO[2](t)] + J = [L_QO[1](t), L_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 = [J_QO[1](t), J_QO[2](t), n_v1_coll(t), n_v2_coll(t)] + J = [L_QO[1](t), L_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 fbf3a893..ac01fe24 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) -J = jump_operator(G_cas) +L = 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) -J_QO = [to_numeric(J[i], b; time_parameter = dict_p_t) for i = 1:length(J)] +L_QO = [to_numeric(L[i], b; time_parameter = dict_p_t) for i = 1:length(L)] function input_output(t, ρ) Ht = H_QO(t) - Jt = [J_QO[i](t) for i = 1:length(J_QO)] + Jt = [L_QO[i](t) for i = 1:length(L_QO)] return Ht, Jt, dagger.(Jt) end nothing # hide @@ -150,12 +150,12 @@ H_int_sym = simplify(substitute(H_int_sym_, int_dict)) # -J_int_sym = simplify.(substitute.(J, Ref(int_dict))) -J_int_sym[1] +L_int_sym = simplify.(substitute.(L, Ref(int_dict))) +L_int_sym[1] # -J_int_sym[2] +L_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) -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) +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) ] ## H_int_QO = to_numeric(H_int_sym, b_int; time_parameter=dict_p_t_int) # 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 +## L_int_QO = [to_numeric(L_int_sym[i], b_int; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide nothing # hide # function input_output_int(t, ρ) Ht = H_int_QO(t) - Jt = [J_int_QO[i](t) for i = 1:length(J_int_QO)] + Jt = [L_int_QO[i](t) for i = 1:length(L_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 84b6b5f7..fa308860 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) -J_loop = jump_operator(G_loop)[1] +L_loop = jump_operator(G_loop)[1] # diff --git a/examples/10-1_SUPER_excitation.jl b/examples/10-1_SUPER_excitation.jl index fd6cfb0e..2af9d054 100644 --- a/examples/10-1_SUPER_excitation.jl +++ b/examples/10-1_SUPER_excitation.jl @@ -50,8 +50,8 @@ nothing # hide ## Hamiltonian and Lindbladian Hcas = hamiltonian(G_cas) -Jcas = jump_operator(G_cas)[1] -Jcasd = adjoint(Jcas) +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. @@ -68,8 +68,8 @@ 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) +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). @@ -119,7 +119,7 @@ nothing # hide ## 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) +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 @@ -144,13 +144,15 @@ nothing # hide # -common = (; - xlims = (t_cas[1]-0.01, t_cas[end]), - tickfontsize = 18, - guidefontsize = 18, - legendfontsize = 18, +common = (; xlims = (t_cas[1]-0.01, t_cas[end])) +p1 = plot( + t_cas, + real.(s22_cas); + color = :red, + label = false, + ylabel = L"\langle\hat\sigma^{ee}\rangle", + common..., ) -p1 = plot(t_cas, real.(s22_cas); color = :red, label = L"\mathrm{cascade}") p2 = plot( t_cas, nu1_cas; @@ -193,7 +195,7 @@ p3 = plot( common..., ) plot!(p3, t_cas, nu2_cas .+ nv2_cas .- nu2_cas[1]; color = :red, label = L"\mathrm{mode~2}") -plot(p1, p2, p3; layout = (3, 1), size = (800, 800)) +plot(p1, p2, p3; layout = (3, 1), size = (600, 700)) # ## Interaction picture @@ -210,8 +212,8 @@ 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)) +L_int = simplify(substitute(Lcas, int_dict)) +Ld_int = simplify(substitute(Lcasd, int_dict)) nothing # hide ## Coefficient matrix M @@ -239,10 +241,10 @@ 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) +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, [J_int_t]; Jdagger = [Jd_int_t], order = order, iv = t); +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)) @@ -272,10 +274,7 @@ pl4 = plot( 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), + size = (500, 350), ) plot!(pl4, t_int, nu2_int .- nu2_int[1]; color = :red, label = L"\mathrm{mode~2~(int.)}") @@ -300,13 +299,13 @@ 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) +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 = [J_int_fock(t)] + J = [L_int_fock(t)] return Ht, J, QuantumOptics.dagger.(J) end @@ -334,13 +333,15 @@ 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, +common = (; xlims = (t_fock[1]-0.01, t_fock[end])) +p3_1 = plot( + t_int, + s22_int; + color = :blue, + label = L"\mathrm{coherent state}", + ylabel = L"\langle\hat\sigma^{ee}\rangle", + common..., ) -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, @@ -375,7 +376,7 @@ plot!( color = :red, label = L"\mathrm{Coherent: mode~2}", ) -pl3 = plot(p3_1, p3_2; layout = (2, 1), size = (800, 600)) +pl3 = plot(p3_1, p3_2; layout = (2, 1), size = (600, 500)) # Due to the vanishing relative phase of the Fock states, the oscillations disappear. 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 b8ba2af4..a11db69b 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) -J = jump_operator(G_s)[1] +L = 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) -J_QO = to_numeric(J, b1; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b1; parameter = dict_p, time_parameter = dict_p_t) function input_output_1(t, ρ) Ht = H_QO(t) - J = [J_QO(t)] + J = [L_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) -Js(t) = √(γ_) * a_qo -g1_m = correlation_matrix(T, ρt_1, input_output_1, Js) +Ls(t) = √(γ_) * a_qo +g1_m = correlation_matrix(T, ρt_1, input_output_1, Ls) 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) -J_2 = jump_operator(G)[1] +L_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) -J_QO_2 = to_numeric(J_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) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_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 3b8079c8..3aa4d771 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 # -J_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] -G_atom = SLH(I4, J_atom, Δ*σ(2, 2)) +L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)] +G_atom = SLH(I4, L_atom, Δ*σ(2, 2)) G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom H = hamiltonian(G_u_bs_d1_d2_atom) -J = jump_operator(G_u_bs_d1_d2_atom) +L = 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)) -J_int = [simplify(substitute(J_, int_dict)) for J_ in J] +L_int = [simplify(substitute(L_, int_dict)) for L_ in L] ## 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) -J_QO_Δ(Δn) = [ - to_numeric(J_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for - i = 1:length(J) +L_QO_Δ(Δn) = [ + to_numeric(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for + i = 1:length(L) ] # H_QO_Δ(Δn) = to_numeric(H, b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) -# J_QO_Δ(Δn) = [to_numeric(J[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(J)] +# L_QO_Δ(Δn) = [to_numeric(L[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(L)] nothing # hide # @@ -170,13 +170,13 @@ ad2_QO = to_numeric(ad2, b) Δn = 0 H_QO = H_QO_Δ(Δn) -J_QO = J_QO_Δ(Δn) +L_QO = L_QO_Δ(Δn) input_output = (t, ρ) -> ( H_QO(t), - [J_QO[i](t) for i = 1:length(J_QO)], - [dagger(J_QO[i](t)) for i = 1:length(J_QO)], + [L_QO[i](t) for i = 1:length(L_QO)], + [dagger(L_QO[i](t)) for i = 1:length(L_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) - J_QO = J_QO_Δ(Δn) + L_QO = L_QO_Δ(Δn) input_output = (t, ρ) -> ( H_QO(t), - [J_QO[i](t) for i = 1:length(J_QO)], - [dagger(J_QO[i](t)) for i = 1:length(J_QO)], + [L_QO[i](t) for i = 1:length(L_QO)], + [dagger(L_QO[i](t)) for i = 1:length(L_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 44c4f569..4a493b0e 100644 --- a/examples/drafts/10-1_SUPER_excitation.jl +++ b/examples/drafts/10-1_SUPER_excitation.jl @@ -184,8 +184,8 @@ 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) +Lcas = jump_operator(G_cas)[1] +Lcasd = adjoint(Lcas) # Time-dependent couplings @register_symbolic gu1_t(t) @@ -205,8 +205,8 @@ 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); +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_) @@ -236,7 +236,7 @@ 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) +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 @@ -344,8 +344,8 @@ 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)) +L_int_sym = simplify(substitute(Lcas, int_dict)) +Ld_int_sym = simplify(substitute(Lcasd, int_dict)) # M-matrix Mat = Matrix{Any}(undef, la, la) @@ -386,14 +386,14 @@ 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) +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, - [J_int_sym_t]; - Jdagger = [Jd_int_sym_t], + [L_int_sym_t]; + Jdagger = [Ld_int_sym_t], order = order, iv = t, ); @@ -448,7 +448,7 @@ 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 + # 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]] @@ -577,27 +577,27 @@ 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) +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 -Jch1 = substitute(Jch1g, dict_gam) # TODO -Jch2 = substitute(Jch2g, dict_gam) -Jdch1 = substitute(Jdch1g, dict_gam) -Jdch2 = substitute(Jdch2g, dict_gam) +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); -Jch1_t = substitute(Jch1, dict_gt); -Jch2_t = substitute(Jch2, dict_gt); -Jdch1_t = substitute(Jdch1, dict_gt); -Jdch2_t = substitute(Jdch2, 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) @@ -613,8 +613,8 @@ 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], + [Lch1_t, Lch2_t]; + Jdagger = [Ldch1_t, Ldch2_t], order = order, iv = t, ) @@ -734,8 +734,8 @@ 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) +Lcas_1m = jump_operator(G_cas_1m)[1] +Lcas_1md = adjoint(Lcas_1m) # Time-dependent couplings @register_symbolic gu_t(t) @@ -752,8 +752,8 @@ 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); +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) @@ -780,7 +780,7 @@ ops_1m = [aun, aun'aun, aun*aun, σ(1, 2), σ(2, 2), avn, avn'avn] # 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); + meanfield(ops_1m, Hcas_1m_t, [Lcas_1m_t]; Jdagger = [Lcas_1md_t], order = order, iv = t); complete!(eqs_1m); length(eqs_1m) @@ -869,8 +869,8 @@ T = [dt:dt:Tend;] 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) +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] @@ -880,8 +880,8 @@ 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) +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) @@ -895,7 +895,7 @@ gv_c_t(t) = conj(gv_t(t)) # 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) +# 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) @@ -904,8 +904,8 @@ 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], + [Lcas_1m0_t]; + Jdagger = [Lcas_1m0d_t], order = order, iv = t, ); @@ -930,13 +930,13 @@ sol_1m0 = solve(prob_1m0, Tsit5(), saveat = va; abstol, reltol); ################################################### # Lindblad vector for correlation functions -Jvec=[aun, σ(1, 2)] +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(Jvec), length(Jvec), N) +Gmat_bar = zeros(ComplexF64, N, length(Lvec), length(Lvec), N) t1 = sol_1m0.t; t2 = sol_1m0.t; τ0 = 0.0 @@ -944,12 +944,12 @@ t2 = sol_1m0.t; # 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) + 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(Jvec)] +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 @@ -968,30 +968,30 @@ function same_lhs(a, b) 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) +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(Jvec) + 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(Jvec)) for c in eachcol(idxvec)) + idx_tuple = Tuple(ntuple(i -> c[i], length(Lvec)) for c in eachcol(idxvec)) return idx_tuple end # Operators included -oper_tup = (Jvec[1], Jvec[2]) +oper_tup = (Lvec[1], Lvec[2]) # Index tuple for selecting correlation components -idx_tup = idx_tuple(Jvec, Corrvec, Corrmat) +idx_tup = idx_tuple(Lvec, Corrvec, Corrmat) # Compute full Gmat_bar on the discrete time grid function compute_Gmat( sol_1m0, eqs_1m0, - Jvec, + Lvec, Corrvec; p0, ps, @@ -1003,7 +1003,7 @@ function compute_Gmat( ) N = length(sol_1m0.t) - Gmat_bar = zeros(ComplexF64, N, length(Jvec), length(Jvec), N) + 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)) : @@ -1014,7 +1014,7 @@ function compute_Gmat( τend = sol_1m0.t[end] - sol_1m0.t[i1] z = i1 - for l = 1:length(Jvec) + for l = 1:length(Lvec) corr = Corrvec[l] complete(corr.de) @@ -1035,14 +1035,14 @@ function compute_Gmat( 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)] + u_flat = [ComplexF64.(getindex.(sol_c.u, idx_v[i])) for i = 1:length(Lvec)] if i1 == 1 - for k = 1:length(Jvec) + for k = 1:length(Lvec) Gmat_bar[i1, k, l, :] .= u_flat[k] end else - for k = 1:length(Jvec) + for k = 1:length(Lvec) Gmat_bar[i1, k, l, :] .= padleftN(u_flat[k], N) end end @@ -1055,7 +1055,7 @@ end # Optionally compute and save Gmat_bar #Gmat_bar = compute_Gmat( -# sol_1m0, eqs_1m0, Jvec, Corrvec; +# sol_1m0, eqs_1m0, Lvec, Corrvec; # p0 = p0, # ps = ps, # saveat = va, @@ -1201,8 +1201,8 @@ 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) +Lcas_1m = jump_operator(G_cas_1m)[1] +Lcas_1md = adjoint(Lcas_1m) # Time-dependent couplings @register_symbolic gu_t(t) @@ -1219,8 +1219,8 @@ 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); +Lcas_1m_t = substitute(Lcas_1m, dict_gts); +Lcas_1md_t = substitute(Lcas_1md, dict_gts); # Pulses and composite mode αeff = 0.0 @@ -1259,8 +1259,8 @@ ops_all = union( eqs_1m = meanfield( ops_all, Hcas_1m_t, - [Jcas_1m_t]; - Jdagger = [Jcas_1md_t], + [Lcas_1m_t]; + Jdagger = [Lcas_1md_t], order = order, iv = t, ) @@ -1303,28 +1303,28 @@ sol_1m0v = solve(probs, Tsit5(), saveat = va; abstol, reltol); ################################################### # Lindblad vector (includes output) -Jvec=[aun, σ(1, 2), avn] +Lvec=[aun, σ(1, 2), avn] N = length(sol_1m0v.t) -Gmat_bar = zeros(ComplexF64, N, length(Jvec), length(Jvec), N) +Gmat_bar = zeros(ComplexF64, N, length(Lvec), length(Lvec), 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) + 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(Jvec)] +Corrvec = [Corrmat[end][i] for i = 1:length(Lvec)] -oper_tup = (Jvec[1], Jvec[2], Jvec[3]) -idx_tup = idx_tuple(Jvec, Corrvec, Corrmat) +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, Jvec, Corrvec; +# sol_1m0v, eqs_1m, Lvec, Corrvec; # p0 = p0, # ps = ps, # saveat = va, 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 0bbab00e..af631c51 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) -J = jump_operator(G)[1] +L = 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 J_s^\dagger(t_1)J_s(t_2)\rangle, -# \qquad J_s = \sqrt{\gamma}\,\sigma_-. +# \langle L_s^\dagger(t_1)L_s(t_2)\rangle, +# \qquad L_s = \sqrt{\gamma}\,\sigma_-. # ``` dict_p_1 = Dict(γ => γ_, g_v => 0.0) H_QO_1 = to_numeric(H, b; parameter = dict_p_1) -J_QO_1 = to_numeric(J, b; parameter = dict_p_1) +L_QO_1 = to_numeric(L, b; parameter = dict_p_1) -t_1, ρt_1 = timeevolution.master(T, ψ0, H_QO_1, [J_QO_1]) +t_1, ρt_1 = timeevolution.master(T, ψ0, H_QO_1, [L_QO_1]) -Js = √(γ_) * σm_qo -g1 = correlation_matrix(T, ρt_1, H_QO_1, [J_QO_1], Js) +Ls = √(γ_) * σm_qo +g1 = correlation_matrix(T, ρt_1, H_QO_1, [L_QO_1], Ls) # 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) -J_QO_2 = to_numeric(J, 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) function input_output_2(t, ρ) Ht = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_QO_2(t)] return Ht, J, dagger.(J) end diff --git a/src/SLH.jl b/src/SLH.jl index bc4b3429..65cea3b3 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(J, H) + _detect_operator_type(L, H) Determine the concrete return type for FunctionWrapper by inspecting -the elements of J and H. Checks FunctionWrapper type parameters first +the elements of L and H. Checks FunctionWrapper type parameters first (no evaluation needed), then static element types. Errors if only plain closures are present (the caller should thread the type through from input SLH objects instead). """ -function _detect_operator_type(J, H) +function _detect_operator_type(L, H) # Check FunctionWrapper elements (carry explicit type info) - for j in J - j isa FunctionWrapper && return _fw_return_type(typeof(j)) + for l in L + l isa FunctionWrapper && return _fw_return_type(typeof(l)) end H isa FunctionWrapper && return _fw_return_type(typeof(H)) # Check static (non-time-dep) elements - for j in J - _is_time_dep(j) || return typeof(j) + for l in L + _is_time_dep(l) || return typeof(l) 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(J::SVector{N}, ::Type{OpType}) where {N,OpType} - if any(_is_time_dep, J) +function _maybe_wrap_jump_operators(L::SVector{N}, ::Type{OpType}) where {N,OpType} + if any(_is_time_dep, L) fw_type = FunctionWrapper{OpType,Tuple{Float64}} - return SVector{N,fw_type}(ntuple(i -> fw_type(_to_func(J[i])), Val(N))) + return SVector{N,fw_type}(ntuple(i -> fw_type(_to_func(L[i])), Val(N))) end - return J + return L end function _maybe_wrap_hamiltonian(H, has_td::Bool, ::Type{OpType}) where {OpType} @@ -88,23 +88,23 @@ end # ────────────────────────────────────────────── """ - SLH{N, ST, JT, HT} + SLH{N, ST, LT, HT} -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. +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. See also [`▷`](@ref), [`⊞`](@ref), [`feedback`](@ref) """ -struct SLH{N,ST,JT,HT,L} +struct SLH{N,ST,LT,HT,L} scattering::SMatrix{N,N,ST,L} - jump_operator::SVector{N,JT} + jump_operator::SVector{N,LT} hamiltonian::HT - function SLH{N,ST,JT,HT}( + function SLH{N,ST,LT,HT}( S::SMatrix{N,N,ST,L}, - jump_operator::SVector{N,JT}, + jump_operator::SVector{N,LT}, H::HT, - ) where {N,ST,JT,HT,L} - return new{N,ST,JT,HT,L}(S, jump_operator, H) + ) where {N,ST,LT,HT,L} + return new{N,ST,LT,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,JT,HT}) where {N,ST,JT,HT} - JT <: FunctionWrapper && return _fw_return_type(JT) +function _op_type(::SLH{N,ST,LT,HT}) where {N,ST,LT,HT} + LT <: FunctionWrapper && return _fw_return_type(LT) HT <: FunctionWrapper && return _fw_return_type(HT) return nothing end @@ -125,55 +125,55 @@ end # ────────────────────────────────────────────── # Canonical: from SMatrix + SVector (handles FunctionWrapper wrapping) -function _build_slh(S::SMatrix{N,N}, J::SVector{N}, H) where {N} - has_td = any(_is_time_dep, J) || _is_time_dep(H) +function _build_slh(S::SMatrix{N,N}, L::SVector{N}, H) where {N} + has_td = any(_is_time_dep, L) || _is_time_dep(H) if has_td - OpType = _detect_operator_type(J, H) - return _build_slh(S, J, H, OpType) + OpType = _detect_operator_type(L, H) + return _build_slh(S, L, H, OpType) end - return SLH{N,eltype(S),eltype(J),typeof(H)}(S, J, H) + return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) end # With explicit operator type (skips detection — used by composition operations) -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) +function _build_slh(S::SMatrix{N,N}, L::SVector{N}, H, ::Type{OpType}) where {N,OpType} + has_td = any(_is_time_dep, L) || _is_time_dep(H) if has_td - J_w = _maybe_wrap_jump_operators(J, OpType) + L_w = _maybe_wrap_jump_operators(L, OpType) H_w = _maybe_wrap_hamiltonian(H, true, OpType) - return SLH{N,eltype(S),eltype(J_w),typeof(H_w)}(S, J_w, H_w) + return SLH{N,eltype(S),eltype(L_w),typeof(H_w)}(S, L_w, H_w) end - return SLH{N,eltype(S),eltype(J),typeof(H)}(S, J, H) + return SLH{N,eltype(S),eltype(L),typeof(H)}(S, L, H) end # Nothing hint falls through to detection -_build_slh(S::SMatrix{N,N}, J::SVector{N}, H, ::Nothing) where {N} = _build_slh(S, J, H) +_build_slh(S::SMatrix{N,N}, L::SVector{N}, H, ::Nothing) where {N} = _build_slh(S, L, H) # From AbstractMatrix + AbstractVector (includes SMatrix + SVector) -function SLH(S::AbstractMatrix, J::AbstractVector, H) - N = length(J) +function SLH(S::AbstractMatrix, L::AbstractVector, H) + N = length(L) @assert size(S, 1) == N && size(S, 2) == N - return _build_slh(SMatrix{N,N}(S), SVector{N}(J...), H) + return _build_slh(SMatrix{N,N}(S), SVector{N}(L...), H) end -# Numeric scalar S + vector J → S * I_{NxN} -function SLH(S::Number, J::AbstractVector, H) - N = length(J) +# Numeric scalar S + vector L → S * I_{NxN} +function SLH(S::Number, L::AbstractVector, H) + N = length(L) S_mat = SMatrix{N,N}(S * LinearAlgebra.I) - return _build_slh(S_mat, SVector{N}(J...), H) + return _build_slh(S_mat, SVector{N}(L...), H) end -# Scalar S + scalar J → SLH{1} -function SLH(S, J, H) +# Scalar S + scalar L → SLH{1} +function SLH(S, L, H) S_mat = SMatrix{1,1}(S) - J_vec = SVector{1}(J) - return _build_slh(S_mat, J_vec, H) + L_vec = SVector{1}(L) + return _build_slh(S_mat, L_vec, H) end -# Symbolic/general scalar S + vector J → S * I -function SLH(S, J::AbstractVector, H) - N = length(J) +# Symbolic/general scalar S + vector L → S * I +function SLH(S, L::AbstractVector, H) + N = length(L) S_mat = SMatrix{N,N}([i == j ? S : 0 for i = 1:N, j = 1:N]) - return _build_slh(S_mat, SVector{N}(J...), H) + return _build_slh(S_mat, SVector{N}(L...), H) end # ────────────────────────────────────────────── @@ -190,7 +190,7 @@ scattering(G::SLH) = G.scattering """ jump_operator(G::SLH) -Return the jump-operator vector `J` of an SLH object. +Return the jump-operator vector `L` 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}, J::SVector{N}) where {N} +@generated function _slh_matvec(S::SMatrix{N,N}, L::SVector{N}) where {N} if N == 1 - return :(SVector{1}(_mul(S[1, 1], J[1]))) + return :(SVector{1}(_mul(S[1, 1], L[1]))) end exprs = [] for i = 1:N first_name = Symbol("tmp_$(i)_1") - terms = [:($first_name = _mul(S[$i, 1], J[1]))] + terms = [:($first_name = _mul(S[$i, 1], L[1]))] acc = first_name for j = 2:N tname = Symbol("tmp_$(i)_$(j)") - push!(terms, :($tname = _add($acc, _mul(S[$i, $j], J[$j])))) + push!(terms, :($tname = _add($acc, _mul(S[$i, $j], L[$j])))) acc = tname end push!(exprs, Expr(:block, terms..., acc)) @@ -245,13 +245,13 @@ end return :(SVector($(exprs...))) end -@generated function _slh_dot(J1::SVector{N}, J2::SVector{N}) where {N} +@generated function _slh_dot(L1::SVector{N}, L2::SVector{N}) where {N} if N == 1 - return :(_mul(J1[1], J2[1])) + return :(_mul(L1[1], L2[1])) end - expr = :(_mul(J1[1], J2[1])) + expr = :(_mul(L1[1], L2[1])) for i = 2:N - expr = :(_add($expr, _mul(J1[$i], J2[$i]))) + expr = :(_add($expr, _mul(L1[$i], L2[$i]))) end return expr end @@ -265,27 +265,27 @@ end Cascade two SLH triples: -``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))`` +``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))`` Unicode `\\triangleright`. See also [`cascade`](@ref). """ function ▷(G1::SLH{N}, G2::SLH{N}) where {N} - S1, J1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) - S2, J2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) + S1, L1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) + S2, L2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) S_t = _post.(S2 * S1) - S2J1 = _slh_matvec(S2, J1) - J_t = SVector{N}(ntuple(i -> _post(_add(J2[i], S2J1[i])), Val(N))) + S2L1 = _slh_matvec(S2, L1) + L_t = SVector{N}(ntuple(i -> _post(_add(L2[i], S2L1[i])), Val(N))) - J2_adj = SVector{N}(ntuple(i -> _adj(J2[i]), Val(N))) - cross1 = _slh_dot(J2_adj, S2J1) + L2_adj = SVector{N}(ntuple(i -> _adj(L2[i]), Val(N))) + cross1 = _slh_dot(L2_adj, S2L1) 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, J_t, H_t, op_hint) + return _build_slh(S_t, L_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} J_1 \\\\ J_2 \\end{pmatrix},\\; H_1 + H_2\\right)`` +\\begin{pmatrix} L_1 \\\\ L_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, J1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) - S2, J2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) + S1, L1, H1 = scattering(G1), jump_operator(G1), hamiltonian(G1) + S2, L2, H2 = scattering(G2), jump_operator(G2), hamiltonian(G2) S_t = SMatrix{$N,$N}($(s_exprs...)) - J_t = vcat(J1, J2) + L_t = vcat(L1, L2) H_t = _add(H1, H2) op_hint = _op_type(G1) op_hint === nothing && (op_hint = _op_type(G2)) - return _build_slh(S_t, J_t, H_t, op_hint) + return _build_slh(S_t, L_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(J::SVector{N}, idx::Int, ::Val{M}) where {N,M} - SVector{M}(ntuple(i -> J[i >= idx ? i + 1 : i], Val(M))) +function _drop_index(L::SVector{N}, idx::Int, ::Val{M}) where {N,M} + SVector{M}(ntuple(i -> L[i >= idx ? i + 1 : i], Val(M))) end function _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) - J = jump_operator(G) + L = 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) - 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)) + 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)) S_col_y_full = SVector{N}(ntuple(i -> S[i, y], Val(N))) - 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)) + L_adj = SVector{N}(ntuple(i -> _adj(L[i]), Val(N))) + term = _mul(_slh_dot(L_adj, S_col_y_full), _mul(loop_gain, L_x)) H_red = _post(_add(H, _mul(1 / (2im), _add(term, _mul(-1, _adj(term)))))) - return _build_slh(S_red, J_red, H_red, _op_type(G)) + return _build_slh(S_red, L_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, J_QO)`, where `J_QO` is a vector holding one translated operator +Returns the tuple `(H_QO, L_QO)`, where `L_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...) - J_QO = [SQA.to_numeric(J_, b; kwargs...) for J_ in jump_operator(G)] - return H_QO, J_QO + L_QO = [SQA.to_numeric(L_, b; kwargs...) for L_ in jump_operator(G)] + return H_QO, L_QO end diff --git a/src/correlations.jl b/src/correlations.jl index ec226c91..8a73be55 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -1,23 +1,23 @@ """ - correlation_matrix(T, ρt, f, Js; kwargs...) - correlation_matrix(T, ρt, H, J, Js; kwargs...) + correlation_matrix(T, ρt, f, Ls; kwargs...) + correlation_matrix(T, ρt, H, J, Ls; kwargs...) Compute the two-time correlation matrix -``g^{(1)}(t_1, t_2) = \\langle J_s^\\dagger(t_1) J_s(t_2) \\rangle`` +``g^{(1)}(t_1, t_2) = \\langle L_s^\\dagger(t_1) L_s(t_2) \\rangle`` on the time grid `T`. Writes directly into output matrix. Supply the dynamics either as a `master_dynamic`-style function `f(t, ρ)`, or as operators passed straight to the solver: a time-dependent `H` (e.g. the `TimeDependentSum` from [`to_numeric`](@ref)) with jump operators `J`, or a constant `H` with constant `J`. The operator form is much faster for time-dependent problems (the integrator is built once). -`Js` is either a constant operator or a function `Js(t)` returning the operator at `t`. +`Ls` is either a constant operator or a function `Ls(t)` returning the operator at `t`. The returned matrix is a `Hermitian` wrapper. To extract the dominant temporal modes, diagonalize it with `eigen(g1_m)`. When only the leading modes are needed, te cheaper eigenvalue-range method can be used, e.g. `eigen(g1_m, (n-4):n)` for the five dominant modes, where `n = size(g1_m, 1)`. """ -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 +function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) end end @@ -27,18 +27,18 @@ function correlation_matrix( ρt::Vector, H::QuantumOpticsBase.AbstractTimeDependentOperator, J::AbstractVector, - Js; + Ls; 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 + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, copy(H), [copy(j) for j in J]; kwargs...) end end -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 +function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Ls; kwargs...) + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) + _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master(T_slice, ρ0, H, J; kwargs...) end end @@ -53,8 +53,8 @@ function _sample_operator_and_adjoint( ) throw( ArgumentError( - "`Js` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " * - "function `Js(t)` returning the concrete operator at time `t`.", + "`Ls` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " * + "function `Ls(t)` returning the concrete operator at time `t`.", ), ) end @@ -63,22 +63,22 @@ function _sample_operator_and_adjoint(T::Vector, op) return fill(op, nt), fill(dagger(op), nt) end -function _correlation_loop(solve_fn, T, ρt, Js_vec, Js_dag_vec) +function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_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 = Js_vec[it] * ρt[it] + ρ0_it = Ls_vec[it] * ρt[it] τ_, ρ_bar_τ = solve_fn(@view(T[it:end]), ρ0_it) @inbounds for i in eachindex(ρ_bar_τ) - val = expect(Js_dag_vec[it+i-1], ρ_bar_τ[i]) + val = expect(Ls_dag_vec[it+i-1], ρ_bar_τ[i]) g1_m[it, it+i-1] = val g1_m[it+i-1, it] = conj(val) end end - g1_m[l_T, l_T] = expect(Js_dag_vec[l_T], Js_vec[l_T] * ρt[l_T]) + g1_m[l_T, l_T] = expect(Ls_dag_vec[l_T], Ls_vec[l_T] * ρt[l_T]) return Hermitian(g1_m) end diff --git a/test/test_SLH.jl b/test/test_SLH.jl index f8218047..1ac400e7 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_J = jump_operator(G_c) + G_c_L = jump_operator(G_c) G_c_H = hamiltonian(G_c) @test G_c_S == G_c.scattering - @test G_c_J == G_c.jump_operator + @test G_c_L == 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_J isa SVector{1} + @test G_c_L 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) - J_s = sparse(sqrt(1.0) * a_op) + L_s = sparse(sqrt(1.0) * a_op) gu_f(t) = exp(-t^2) * sparse(a_op) gv_f(t) = exp(-(t - 2)^2) * sparse(a_op) @testset "static SLH concreteness" begin - G = SLH(1, J_s, H_s) - @test eltype(jump_operator(G)) === typeof(J_s) + G = SLH(1, L_s, H_s) + @test eltype(jump_operator(G)) === typeof(L_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, J_s, H_s) ▷ SLH(1, gu_f, H_s) + G_cas = SLH(1, L_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, J_s, H_s) ⊞ SLH(1, gu_f, H_s) - JT = eltype(jump_operator(G_cat)) - @test JT <: FunctionWrapper - @test JT !== Any + 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 @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, J_s, H_s) ⊞ SLH(1, J_s, H_s) - @test eltype(jump_operator(G_cat)) === typeof(J_s) + G_cat = SLH(1, L_s, H_s) ⊞ SLH(1, L_s, H_s) + @test eltype(jump_operator(G_cat)) === typeof(L_s) @test !(eltype(jump_operator(G_cat)) <: FunctionWrapper) end @testset "FunctionWrapper call is inferred" begin G_td = SLH(1, gu_f, H_s) - j = jump_operator(G_td)[1] - @inferred j(0.5) + l = jump_operator(G_td)[1] + @inferred l(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, J_s, H_s) + G_s = SLH(1, L_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 09246e90..d36d5c6a 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) - Jop = jump_operator(G_t) - J_R = Jop[1] - J_L = Jop[2] + L = jump_operator(G_t) + L_R = L[1] + L_L = L[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) - 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) + 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) σ_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 = [J_R_QO(t), J_L_QO(t), J_add...] + J = [L_R_QO(t), L_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(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)] + I_R_1 = [real(expect(L_R_QO(ti)' * L_R_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] + I_L_1 = [real(expect(L_L_QO(ti)' * L_L_QO(ti), ρt1[i])) for (i, ti) in enumerate(t1)] # -------- Example 05-2 style (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) - J_qo = jump_operator(G_t_qo) - J_R_qo = J_qo[1] - J_L_qo = J_qo[2] + L_qo = jump_operator(G_t_qo) + L_R_qo = L_qo[1] + L_L_qo = L_qo[2] # FunctionWrapper is callable but not <: Function _callable(x) = x isa Union{Function,FunctionWrapper} Hf = _callable(H_qo) ? H_qo : (t -> H_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) + L_R_f = _callable(L_R_qo) ? L_R_qo : (t -> L_R_qo) + L_L_f = _callable(L_L_qo) ? L_L_qo : (t -> L_L_qo) J_add_qo = [√(γ_add[i]) * σ_qds(i, 1, 2) for i = 1:N] function input_output_2(t, ρ) Ht = Hf(t) - J = [J_R_f(t), J_L_f(t), J_add_qo...] + J = [L_R_f(t), L_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(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)] + 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)] @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 d9b58647..f992d09b 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] - Js = J[1] + Ls = 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, Js) + g1_dyn = correlation_matrix(T, ρt_dyn, f_const, Ls) _, ρt_static = timeevolution.master(T, ψ0, H, J) - g1_static = correlation_matrix(T, ρt_static, H, J, Js) + g1_static = correlation_matrix(T, ρt_static, H, J, Ls) @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 Js rejected" begin - td_Js = TimeDependentSum([t -> 1.0 + 0im], J) - @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Js) + @testset "time-dependent Ls rejected" begin + td_Ls = TimeDependentSum([t -> 1.0 + 0im], J) + @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Ls) end @testset "time-dependent operator path matches static" begin H_td = TimeDependentSum([t -> 0.0 + 0im], [one(ba)]) J_td = [TimeDependentSum([t -> sqrt(γ) + 0im], [σm])] _, ρt_td = timeevolution.master_dynamic(T, ψ0, H_td, J_td) - g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Js) + g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Ls) @test maximum(abs.(g1_td .- g1_static)) < 1e-8 end end diff --git a/test/test_example_cavity_scattering.jl b/test/test_example_cavity_scattering.jl index 138b5918..d74c61ce 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) - Jop = jump_operator(G_cas)[1] + L = 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) - J_QO = to_numeric(Jop, b; parameter = dict_p, time_parameter = dict_p_t) + L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) function input_output_1(t, ρ) Ht = H_QO(t) - J = [J_QO(t)] + J = [L_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 - Js(t) = gu_t(t)*au_qo + √(γ_)*c_qo - g1_m = correlation_matrix(T, ρt, input_output_1, Js) + Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo + g1_m = correlation_matrix(T, ρt, input_output_1, Ls) 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) - J_QO_2 = to_numeric(Jop, 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) function input_output_2(t, ρ) H = H_QO_2(t) - J = [J_QO_2(t)] + J = [L_QO_2(t)] return H, J, dagger.(J) end diff --git a/test/test_feedback.jl b/test/test_feedback.jl index 5dac0e72..024d73dd 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") - # j1_r = rnumber("j1_r") - # j2_r = rnumber("j2_r") + # l1_r = rnumber("l1_r") + # l2_r = rnumber("l2_r") # # # s11_i = rnumber("s11_i") # s12_i = rnumber("s12_i") # s21_i = rnumber("s21_i") # s22_i = rnumber("s22_i") - # j1_i = rnumber("j1_i") - # j2_i = rnumber("j2_i") + # l1_i = rnumber("l1_i") + # l2_i = rnumber("l2_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 - # j1 = j1_r + 1im*j1_i - # j2 = j2_r + 1im*j2_i + # l1 = l1_r + 1im*l1_i + # l2 = l2_r + 1im*l2_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 - j1 = 0.18 + 1im*0.21 - j2 = 0.15 + 1im*0.25 + l1 = 0.18 + 1im*0.21 + l2 = 0.15 + 1im*0.25 h0 = 1.0 - G = SLH([s11 s12; s21 s22], [j1, j2], h0) + G = SLH([s11 s12; s21 s22], [l1, l2], h0) G_red = feedback(G, 1, 1) loop_gain = (1 - s11)^(-1) expected_S = simplify(s22 + s21 * loop_gain * s12) - expected_J = simplify(j2 + s21 * loop_gain * j1) - expected_term = simplify((j1' * s11 + j2' * s21) * loop_gain * j1) + expected_L = simplify(l2 + s21 * loop_gain * l1) + expected_term = simplify((l1' * s11 + l2' * s21) * loop_gain * l1) expected_H = simplify(h0 + (expected_term - expected_term') / (2im)) @test 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_J) < 1e-10 + @test abs(jump_operator(G_red)[1] - expected_L) < 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) - J_s = sparse(sqrt(1.0) * a_op) + L_s = sparse(sqrt(1.0) * a_op) gu_f(t) = exp(-t^2) * sparse(a_op) gv_f(t) = exp(-(t - 2)^2) * sparse(a_op) G_cat = SLH(1, gu_f, H_s) ⊞ SLH(1, gv_f, H_s) G_fb = feedback(G_cat, 1, 1) - JT = eltype(jump_operator(G_fb)) - @test JT <: FunctionWrapper - @test JT !== Any + LT = eltype(jump_operator(G_fb)) + @test LT <: FunctionWrapper + @test LT !== 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 ba76b13c..9db62b54 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) - Jop = jump_operator(G_cas)[1] + L = 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)) - J_int_sym = simplify(substitute(Jop, int_dict)) + L_int_sym = simplify(substitute(L, 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) - J_QO = to_numeric(J_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) function input_output_I(t, ρ) Ht = H_int_QO(t) - J = [J_QO(t)] + J = [L_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) - J_QO_S = to_numeric(Jop, 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) function input_output_S(t, ρ) Ht = H_QO(t) - J = [J_QO_S(t)] + J = [L_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) - J_QO_small = - to_numeric(J_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) function input_output_I_small(t, ρ) Ht = H_int_QO_small(t) - J = [J_QO_small(t)] + J = [L_QO_small(t)] return Ht, J, dagger.(J) end diff --git a/test/test_translate.jl b/test/test_translate.jl index fe33fed7..5fe1168c 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, J_QO = to_numeric( + H_QO, L_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(J_QO) == 1 - @test J_QO[1] isa TimeDependentSum + @test length(L_QO) == 1 + @test L_QO[1] isa TimeDependentSum @test sum(abs.((mat(H_QO, 0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8 - @test sum(abs.((mat(J_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < + @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 end