G1 hermitian

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)

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

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

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)