Add function correlation_dynamic

src/timecorrelations.jl

@@ -1,6 +1,6 @@
 module timecorrelations
 
-export correlation, spectrum, correlation2spectrum
+export correlation, spectrum, correlation2spectrum, correlation_dynamic
 
 using QuantumOpticsBase
 using ..timeevolution, ..steadystate
@@ -62,6 +62,39 @@ function correlation(rho0::Operator, H::AbstractOperator, J,
 end
 
 
+"""
+    timecorrelations.correlation_dynamic(tspan, rho0, f, op1, op2; <keyword arguments>)
+
+
+Calculate two time correlation values ``⟨A(t)B(0)⟩`` for time-dependent Liouvillian
+
+The calculation is done by multiplying the initial density operator
+with ``B`` performing a time evolution according to a master equation
+and then calculating the expectation value ``\\mathrm{Tr} \\{A ρ\\}``
+
+
+# Arguments
+* `tspan`: Points of time at which the correlation should be calculated.
+* `rho0`: Initial density operator.
+* `f`: Function `f(t, rho) -> (H, J, Jdagger)` or `f(t, rho) -> (H, J, Jdagger, rates)`
+* `op1`: Operator at time `t`.
+* `op2`: Operator at time `t=0`.
+* `rates=ones(N)`: Vector or matrix specifying the coefficients (decay rates)
+        for the jump operators.
+* `kwargs...`: Further arguments are passed on to the ode solver.
+"""
+function correlation_dynamic(tspan, rho0::Operator, f, op1, op2;
+                             rates=nothing, kwargs...)
+
+    function fout(t, rho)
+        expect(op1, rho)
+    end
+
+    t,u = timeevolution.master_dynamic(tspan, op2*rho0, f, rates=rates,fout=fout, kwargs...)
+    u
+end
+
+
 """
     timecorrelations.spectrum([omega_samplepoints,] H, J, op; <keyword arguments>)
 

test/test_timecorrelations.jl

@@ -30,6 +30,10 @@ Ja2 = embed(basis, 1, sqrt(0.5*γ)*sp)
 Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
 J = [Ja, Ja2, Jc]
 
+# time-dependent Hamiltonian in rotating frame of cavity mode,
+# which should not change the time correlation for spin operators.
+f_HJ(t,rho) = [ Ha + exp(im*ωc*t) * sm ⊗ create(fockbasis) + exp(-im*ωc*t) * sp ⊗ destroy(fockbasis), J, dagger.(J) ]
+
 Ψ₀ = basisstate(spinbasis, 2) ⊗ fockstate(fockbasis, 5)
 ρ₀ = dm(Ψ₀)
 
@@ -41,12 +45,16 @@ exp_values = timecorrelations.correlation(tspan, ρ₀, H, J, dagger(op), op)
 ρ₀ = dm(Ψ₀)
 
 tout, exp_values2 = timecorrelations.correlation(ρ₀, H, J, dagger(op), op; tol=1e-5)
+
+exp_values3 = timecorrelations.correlation_dynamic(tspan, ρ₀, f_HJ, dagger(op), op)
 
 @test length(exp_values) == length(tspan)
 @test length(exp_values2) == length(tout)
+@test length(exp_values3) == length(tspan)
 @test norm(exp_values[1]-exp_values2[1]) < 1e-15
 @test norm(exp_values[end]-exp_values2[end]) < 1e-4
-
+@test all(norm.(exp_values .- exp_values3) .< 1e-4)
+
 n = length(tspan)
 omega_sample = mod(n, 2) == 0 ? [-n/2:n/2-1;] : [-(n-1)/2:(n-1)/2;]
 omega_sample .*= 2pi/tspan[end]