Compare timeevolution of coherent states created with analytic method and displacement method

coherentstates.jl

using QuantumOptics
using PyPlot

# Parameters
const ω = 1.89 # Frequency of driving laser
const ωc = 2.13 # Cavity frequency
const η = 0.76 # Pump strength
const κ = 0.34 # Decay rate
const δc = ωc - ω # Detuning
const T = [0, 10]
const α0 = 1.3 - 1.1im # Initial coherent state
const αT = (α0 + η/δc)*exp(-1im*δc*T[end]) - η/δc # Final coherent state
println("αT: ", αT)

# "Exact" solution at t=0 and t=T[end]
b_check = FockBasis(500)
psi_t0_check1 = coherentstate(b_check, α0)
psi_t0_check2 = displace(b_check, α0)*fockstate(b_check, 0)
println("Error of 'exact' solution at t=t0: ", norm(psi_t0_check1 - psi_t0_check2))
@assert 1e-14 > norm(psi_t0_check1 - psi_t0_check2)

psi_t1_check1 = coherentstate(b_check, αT)*exp(1im*1.7149805399866265)
psi_t1_check2 = displace(b_check, αT)*fockstate(b_check, 0)*exp(1im*1.7149805399866265)
println("Error of 'exact' solution at t=T[end]: ", norm(psi_t1_check1 - psi_t1_check2))
@assert 1e-14 > norm(psi_t1_check1 - psi_t1_check2)

psi_t0_check = (psi_t0_check1 + psi_t0_check2)/2
psi_t1_check = (psi_t1_check1 + psi_t1_check2)/2

# Difference between "exact" solution and solution obatined with a certain cutoff
function error(psi_large::Ket, psi_small::Ket)
    psi_extended = Ket(b_check)
    psi_extended.data[1:length(psi_small.data)] = psi_small.data
    norm(psi_large - psi_extended)
end

error_analytic_t0 = Float64[]
error_analytic_t1 = Float64[]
error_displace_t0 = Float64[]
error_displace_t1 = Float64[]

function simulate(Ncutoff)
    b = FockBasis(Ncutoff)

    a = destroy(b)
    at = create(b)
    n = number(b)
    Hint = δc*n + η*(a + at)

    # Initial state
    psi0_analytic = coherentstate(b, α0)
    psi0_displace = displace(b, α0)*fockstate(b, 0)

    # No decay
    tout, psit_analytic = timeevolution.schroedinger(T, psi0_analytic, Hint, abstol=1e-10, reltol=1e-10)
    tout, psit_displace = timeevolution.schroedinger(T, psi0_displace, Hint, abstol=1e-10, reltol=1e-10)

    push!(error_analytic_t0, error(psi_t0_check, psit_analytic[1]))
    push!(error_analytic_t1, error(psi_t1_check, psit_analytic[end]))
    push!(error_displace_t0, error(psi_t0_check, psit_displace[1]))
    push!(error_displace_t1, error(psi_t1_check, psit_displace[end]))
end

Ncutoffs = [1:1:200;]
for Ncutoff in Ncutoffs
    simulate(Ncutoff)
end

println(error_analytic_t0)
println(error_analytic_t1)
println(error_displace_t0)
println(error_displace_t1)


# Visualization
subplot(2, 1, 1)
semilogy(Ncutoffs, error_analytic_t0, label="analytic")
semilogy(Ncutoffs, error_displace_t0, label="displace")
xlabel("cutoff")
ylabel("Error t=0")
legend()
grid()

subplot(2, 1, 2)
semilogy(Ncutoffs, error_analytic_t1, label="analytic")
semilogy(Ncutoffs, error_displace_t1, label="displace")
xlabel("cutoff")
ylabel("Error t=10")
legend()
grid()

tight_layout()

show()