introduce correlation_3op_1t

Author: ytdHuang

docs/src/resources/api.md

@@ -230,6 +230,7 @@ lr_mesolve
 
 ```@docs
 correlation_3op_2t
+correlation_3op_1t
 correlation_2op_2t
 correlation_2op_1t
 spectrum_correlation_fft

src/correlations.jl

@@ -1,4 +1,10 @@
-export correlation_3op_2t, correlation_2op_2t, correlation_2op_1t
+export correlation_3op_2t, correlation_3op_1t, correlation_2op_2t, correlation_2op_1t
+
+function _check_correlation_time_list(tlist::AbstractVector)
+    any(t -> t == 0, tlist) || throw(ArgumentError("The time list for calculating correlation function must contain the element `0`"))
+    all(>=(0), tlist) || throw(ArgumentError("All the elements in the time list for calculating correlation function must be positive."))
+    return nothing
+end
 
 @doc raw"""
     correlation_3op_2t(H::AbstractQuantumObject,
@@ -34,6 +40,10 @@ function correlation_3op_2t(
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
+    # check tlist and τlist
+    _check_correlation_time_list(tlist)
+    _check_correlation_time_list(τlist)
+
     L = liouvillian(H, c_ops)
     if ψ0 isa Nothing
         ψ0 = steadystate(L)
@@ -47,7 +57,44 @@ function correlation_3op_2t(
 
     corr = map((t, ρ) -> mesolve(L, C * ρ * A, τlist .+ t, e_ops = [B]; kwargs...).expect[1, :], tlist, ρt)
 
-    return corr
+    return reduce(hcat, corr)
+end
+
+@doc raw"""
+    correlation_3op_1t(H::AbstractQuantumObject,
+        ψ0::Union{Nothing,QuantumObject},
+        τlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple},
+        A::QuantumObject,
+        B::QuantumObject,
+        C::QuantumObject;
+        kwargs...)
+
+Returns the one-time correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(0) \hat{B}(\tau) \hat{C}(0) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+
+If the initial state `ψ0` is given as `nothing`, then the [`steadystate`] will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
+"""
+function correlation_3op_1t(
+    H::AbstractQuantumObject{DataType,HOpType},
+    ψ0::Union{Nothing,QuantumObject{<:AbstractArray{T1},StateOpType}},
+    τlist::AbstractVector,
+    c_ops::Union{Nothing,AbstractVector,Tuple},
+    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
+    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
+    C::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject};
+    kwargs...,
+) where {
+    DataType,
+    T1,
+    T2,
+    T3,
+    T4,
+    HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
+}
+    corr = correlation_3op_2t(H, ψ0, [0], τlist, c_ops, A, B, C; kwargs...)
+
+    return corr[:, 1]
 end
 
 @doc raw"""
@@ -92,7 +139,7 @@ function correlation_2op_2t(
         corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, C, A, B; kwargs...)
     end
 
-    return reduce(hcat, corr)
+    return corr
 end
 
 @doc raw"""

test/core-test/correlations_and_spectrum.jl

@@ -1,11 +1,13 @@
 @testset "Correlations and Spectrum" begin
+    Id = qeye(10)
     a = destroy(10)
     H = a' * a
     c_ops = [sqrt(0.1 * (0.01 + 1)) * a, sqrt(0.1 * (0.01)) * a']
 
     t_l = range(0, 333 * π, length = 1000)
-    corr = correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
-    ω_l1, spec1 = spectrum_correlation_fft(t_l, corr)
+    corr1 = correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
+    corr2 = correlation_3op_1t(H, nothing, t_l, c_ops, Id, a', a; progress_bar = Val(false))
+    ω_l1, spec1 = spectrum_correlation_fft(t_l, corr1)
 
     ω_l2 = range(0, 3, length = 1000)
     spec2 = spectrum(H, ω_l2, c_ops, a', a)
@@ -21,14 +23,25 @@
     idxs2 = test_func2 .> 0.05
     @test sum(abs2.(spec1[idxs1] .- test_func1[idxs1])) / sum(abs2.(test_func1[idxs1])) < 0.01
     @test sum(abs2.(spec2[idxs2] .- test_func2[idxs2])) / sum(abs2.(test_func2[idxs2])) < 0.01
+    @test all(corr1 .≈ corr2)
     @test all(spec2 .≈ spec3)
 
     @testset "Type Inference spectrum" begin
         @inferred correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
-        @inferred spectrum_correlation_fft(t_l, corr)
+        @inferred spectrum_correlation_fft(t_l, corr1)
         @inferred spectrum(H, ω_l2, c_ops, a', a)
         @inferred spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
     end
 
+    t_wrong1 = [1, 2, 3]
+    t_wrong2 = [-1, 0, 1]
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_l, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_l, t_wrong2, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_l, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_l, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong1, t_wrong2, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_wrong1, c_ops, Id, a', a)
+    @test_throws ArgumentError correlation_3op_2t(H, nothing, t_wrong2, t_wrong2, c_ops, Id, a', a)
     @test_throws ErrorException FFTCorrelation()
 end