From 0fdf1daa9cd223f259a1c1599e1d9d4bbd348724 Mon Sep 17 00:00:00 2001 From: ChristophHotter Date: Mon, 24 Aug 2026 11:53:35 -0400 Subject: [PATCH 1/2] Hermitian in g1 --- src/QuantumInputOutput.jl | 2 +- src/correlations.jl | 6 +++++- test/test_correlations.jl | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl index 58412d9d..dd1a6bfd 100644 --- a/src/QuantumInputOutput.jl +++ b/src/QuantumInputOutput.jl @@ -8,7 +8,7 @@ using Symbolics: Symbolics using SpecialFunctions: erf using DataInterpolations: LinearInterpolation, ExtrapolationType using NumericalIntegration: cumul_integrate -using LinearAlgebra: LinearAlgebra, I, mul! +using LinearAlgebra: LinearAlgebra, I, mul!, Hermitian using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, Tsit5, solve using StaticArrays: StaticArrays, SMatrix, SVector using FunctionWrappers: FunctionWrappers, FunctionWrapper diff --git a/src/correlations.jl b/src/correlations.jl index 80f1a05c..ec226c91 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -11,6 +11,9 @@ passed straight to the solver: a time-dependent `H` (e.g. the `TimeDependentSum` [`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`. + +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) @@ -76,5 +79,6 @@ function _correlation_loop(solve_fn, T, ρt, Js_vec, Js_dag_vec) g1_m[it+i-1, it] = conj(val) end end - return g1_m + g1_m[l_T, l_T] = expect(Js_dag_vec[l_T], Js_vec[l_T] * ρt[l_T]) + return Hermitian(g1_m) end diff --git a/test/test_correlations.jl b/test/test_correlations.jl index 81c30934..d9b58647 100644 --- a/test/test_correlations.jl +++ b/test/test_correlations.jl @@ -30,7 +30,7 @@ using Test @test maximum(abs.(g1_dyn .- g1_static)) < 1e-8 - F = eigen(Hermitian((g1_static + g1_static') / 2)) + F = eigen(g1_static) n_modes = real.(F.values) * ΔT expected_n = 1 - exp(-γ * T_end) From 3fe856332cc4d22a5af9852639d6e0c978881b38 Mon Sep 17 00:00:00 2001 From: ChristophHotter Date: Mon, 24 Aug 2026 11:58:59 -0400 Subject: [PATCH 2/2] eigen(g1_m, (n-4):n) in example 1 --- .../01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 3141fc82..85fddf91 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 @@ -129,7 +129,8 @@ p # The eigenvalues and corresponding eigenvectors are sorted in ascending order, which means the last eigenvalue corresponds to the highest populated temporal mode. -F = eigen(g1_m) +n = size(g1_m, 1) +F = eigen(g1_m, (n-4):n) n_avg = round.(real.(F.values)*ΔT; digits = 3) modes = F.vectors v_mode = (modes[:, end]) / sqrt(ΔT)