Two-time Correlation Functions

Introduction

With the QuantumToolbox.jl time-evolution function mesolve, a state vector (Ket) or density matrix (Operator) can be evolved from an initial state at to an arbitrary time , namely

where   is the propagator defined by the equation of motion. The resulting density matrix can then be used to evaluate the expectation values of arbitrary combinations of same-time operators.

To calculate two-time correlation functions on the form , we can use the quantum regression theorem [see, e.g., DocumenterCitations.CitationSiteNode("Gardiner-Zoller2004-cite-1")

] to write

We therefore first calculate   using mesolve with as initial state, and then again use mesolve to calculate   using as initial state.

Note that if the initial state is the steady state, then    and

which is independent of , so that we only have one time coordinate .

QuantumToolbox.jl provides a family of functions that assists in the process of calculating two-time correlation functions. The available functions and their usage is shown in the table below.

Function call	Correlation function
correlation_2op_2t	or
correlation_2op_1t	or
correlation_3op_1t
correlation_3op_2t

The most common used case is to calculate the two time correlation function , which can be done by correlation_2op_1t.

Steadystate correlation function

The following code demonstrates how to calculate the correlation for a leaky cavity with three different relaxation rates .

tlist = LinRange(0, 10, 200)
a = destroy(10)
x = a' + a
H = a' * a

# if the initial state is specified as `nothing`, the steady state will be calculated and used as the initial state.
corr1 = correlation_2op_1t(H, nothing, tlist, [sqrt(0.5) * a], x, x)
corr2 = correlation_2op_1t(H, nothing, tlist, [sqrt(1.0) * a], x, x)
corr3 = correlation_2op_1t(H, nothing, tlist, [sqrt(2.0) * a], x, x)

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = L"Time $t$", ylabel = L"\langle \hat{x}(t) \hat{x}(0) \rangle")
lines!(ax, tlist, real(corr1), label = L"\gamma = 0.5", linestyle = :solid)
lines!(ax, tlist, real(corr2), label = L"\gamma = 1.0", linestyle = :dash)
lines!(ax, tlist, real(corr3), label = L"\gamma = 2.0", linestyle = :dashdot)

axislegend(ax, position = :rt)

fig

Emission spectrum

Given a correlation function , we can define the corresponding power spectrum as

In QuantumToolbox.jl, we can calculate using either spectrum, which provides several solvers to perform the Fourier transform semi-analytically, or we can use the function spectrum_correlation_fft to numerically calculate the fast Fourier transform (FFT) of a given correlation data.

The following example demonstrates how these methods can be used to obtain the emission (  and  ) power spectrum.

N = 4             # number of cavity fock states
ωc = 1.0 * 2 * π  # cavity frequency
ωa = 1.0 * 2 * π  # atom frequency
g  = 0.1 * 2 * π  # coupling strength
κ  = 0.75         # cavity dissipation rate
γ  = 0.25         # atom dissipation rate

# Jaynes-Cummings Hamiltonian
a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))
H = ωc * a' * a + ωa * sm' * sm + g * (a' * sm + a * sm')

# collapse operators
n_th = 0.25
c_ops = [
    sqrt(κ * (1 + n_th)) * a,
    sqrt(κ *      n_th)  * a',
    sqrt(γ)              * sm,
];

# calculate the correlation function using mesolve, and then FFT to obtain the spectrum.
# Here we need to make sure to evaluate the correlation function for a sufficient long time and
# sufficiently high sampling rate so that FFT captures all the features in the resulting spectrum.
tlist = LinRange(0, 100, 5000)
corr = correlation_2op_1t(H, nothing, tlist, c_ops, a', a; progress_bar = Val(false))
ωlist1, spec1 = spectrum_correlation_fft(tlist, corr)

# calculate the power spectrum using spectrum
# using Exponential Series (default) method
ωlist2 = LinRange(0.25, 1.75, 200) * 2 * π
spec2 = spectrum(H, ωlist2, c_ops, a', a; solver = ExponentialSeries())

# calculate the power spectrum using spectrum
# using Pseudo-Inverse method
spec3 = spectrum(H, ωlist2, c_ops, a', a; solver = PseudoInverse())

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], title = "Vacuum Rabi splitting", xlabel = "Frequency", ylabel = "Emission power spectrum")
lines!(ax, ωlist1 / (2 * π), spec1, label = "mesolve + FFT", linestyle = :solid)
lines!(ax, ωlist2 / (2 * π), spec2, label = "Exponential Series", linestyle = :dash)
lines!(ax, ωlist2 / (2 * π), spec3, label = "Pseudo-Inverse", linestyle = :dashdot)

xlims!(ax, ωlist2[1] / (2 * π), ωlist2[end] / (2 * π))
axislegend(ax, position = :rt)

fig

Non-steadystate correlation function

More generally, we can also calculate correlation functions of the kind , i.e., the correlation function of a system that is not in its steady state. In QuantumToolbox.jl, we can evaluate such correlation functions using the function correlation_2op_2t. The default behavior of this function is to return a matrix with the correlations as a function of the two time coordinates ( and ).

t1_list = LinRange(0, 10.0, 200)
t2_list = LinRange(0, 10.0, 200)

N = 10
a = destroy(N)
x = a' + a
H = a' * a

c_ops = [sqrt(0.25) * a]

α = 2.5
ρ0 = coherent_dm(N, α)

corr = correlation_2op_2t(H, ρ0, t1_list, t2_list, c_ops, x, x; progress_bar = Val(false))

# plot by CairoMakie.jl
fig = Figure(size = (500, 400))

ax = Axis(fig[1, 1], title = L"\langle \hat{x}(t_1 + t_2) \hat{x}(t_1) \rangle", xlabel = L"Time $t_1$", ylabel = L"Time $t_2$")

heatmap!(ax, t1_list, t2_list, real(corr))

fig

Example: first-order optical coherence function

This example demonstrates how to calculate a correlation function on the form for a non-steady initial state. Consider an oscillator that is interacting with a thermal environment. If the oscillator initially is in a coherent state, it will gradually decay to a thermal (incoherent) state. The amount of coherence can be quantified using the first-order optical coherence function

For a coherent state  , and for a completely incoherent (thermal) state  . The following code calculates and plots as a function of :

τlist = LinRange(0, 10, 200)

# Hamiltonian
N = 15
a = destroy(N)
H = 2 * π * a' * a

# collapse operator
G1 = 0.75
n_th = 2.00  # bath temperature in terms of excitation number
c_ops = [
    sqrt(G1 * (1 + n_th)) * a,
    sqrt(G1 *      n_th)  * a'
]

# start with a coherent state of α = 2.0
ρ0 = coherent_dm(N, 2.0)

# first calculate the occupation number as a function of time
n = mesolve(H, ρ0, τlist, c_ops, e_ops = [a' * a], progress_bar = Val(false)).expect[1,:]
n0 = n[1] # occupation number at τ = 0

# calculate the correlation function G1 and normalize with n to obtain g1
g1 = correlation_2op_1t(H, ρ0, τlist, c_ops, a', a, progress_bar = Val(false))
g1 = g1 ./ sqrt.(n .* n0)

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], title = "Decay of a coherent state to an incoherent (thermal) state", xlabel = L"Time $\tau$")
lines!(ax, τlist, real(g1), label = L"g^{(1)}(\tau)", linestyle = :solid)
lines!(ax, τlist, real(n), label = L"n(\tau)", linestyle = :dash)

axislegend(ax, position = :rt)

fig

Example: second-order optical coherence function

The second-order optical coherence function, with time-delay , is defined as

For a coherent state  , for a thermal state    and it decreases as a function of time (bunched photons, they tend to appear together), and for a Fock state with -photons      and it increases with time (anti-bunched photons, more likely to arrive separated in time).

To calculate this type of correlation function with QuantumToolbox.jl, we can use correlation_3op_1t, which computes a correlation function on the form (three operators and one delay-time vector). We first have to combine the central two operators into one single one as they are evaluated at the same time, e.g. here we do   .

The following code calculates and plots as a function of for a coherent, thermal and Fock state:

τlist = LinRange(0, 25, 200)

# Hamiltonian
N = 25
a = destroy(N)
H = 2 * π * a' * a

κ = 0.25
n_th = 2.0  # bath temperature in terms of excitation number
c_ops = [
    sqrt(κ * (1 + n_th)) * a,
    sqrt(κ *      n_th)  * a'
]

cases = [
    Dict("state" => coherent_dm(N, sqrt(2)), "label" => "coherent state", "lstyle" => :solid),
    Dict("state" => thermal_dm(N, 2), "label" => "thermal state", "lstyle" => :dash),
    Dict("state" => fock_dm(N, 2), "label" => "Fock state", "lstyle" => :dashdot),
]

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = L"Time $\tau$", ylabel = L"g^{(2)}(\tau)")

for case in cases
    ρ0 = case["state"]

    # calculate the occupation number at τ = 0
    n0 = expect(a' * a, ρ0)

    # calculate the correlation function g2
    g2 = correlation_3op_1t(H, ρ0, τlist, c_ops, a', a' * a, a, progress_bar = Val(false))
    g2 = g2 ./ n0^2

    lines!(ax, τlist, real(g2), label = case["label"], linestyle = case["lstyle"])
end

axislegend(ax, position = :rt)

fig