guide-correlation.rst

Two-time correlation functions

With the QuTiP time-evolution functions (for example qutip.mesolve and qutip.mcsolve), a state vector or density matrix can be evolved from an initial state at t₀ to an arbitrary time t, ρ(t) = V(t, t₀){ρ(t₀)}, where V(t, t₀) 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 $\left<A(t+\tau)B(t)\right>$, we can use the quantum regression theorem (see, e.g., [Gar03]) to write

$$\left<A(t+\tau)B(t)\right> = {\rm Tr}\left[A V(t+\tau, t)\left\{B\rho(t)\right\}\right] = {\rm Tr}\left[A V(t+\tau, t)\left\{BV(t, 0)\left\{\rho(0)\right\}\right\}\right]$$

We therefore first calculate ρ(t) = V(t, 0){ρ(0)} using one of the QuTiP evolution solvers with ρ(0) as initial state, and then again use the same solver to calculate V(t + τ, t){Bρ(t)} using Bρ(t) as initial state.

Note that if the initial state is the steady state, then $\rho(t)=V(t, 0)\left\{\rho_{\rm ss}\right\}=\rho_{\rm ss}$ and

$$\left<A(t+\tau)B(t)\right> = {\rm Tr}\left[A V(t+\tau, t)\left\{B\rho_{\rm ss}\right\}\right] = {\rm Tr}\left[A V(\tau, 0)\left\{B\rho_{\rm ss}\right\}\right] = \left<A(\tau)B(0)\right>,$$

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

QuTiP 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. Each of these functions can use one of the following evolution solvers: Master-equation, Exponential series and the Monte-Carlo. The choice of solver is defined by the optional argument solver.

table-striped

QuTiP function	Correlation function
qutip.correlation.correlation_2op_2t	$\left<A(t+\tau)B(t)\right>$ or $\left<A(t)B(t+\tau)\right>$.
qutip.correlation.correlation_2op_1t	$\left<A(\tau)B(0)\right>$ or $\left<A(0)B(\tau)\right>$.
qutip.correlation.correlation_3op_1t	$\left<A(0)B(\tau)C(0)\right>$.
qutip.correlation.correlation_3op_2t	$\left<A(t)B(t+\tau)C(t)\right>$.

The most common use-case is to calculate correlation functions of the kind $\left&lt;A(\tau)B(0)\right&gt;$, in which case we use the correlation function solvers that start from the steady state, e.g., the qutip.correlation.correlation_2op_1t function. These correlation function solvers return a vector or matrix (in general complex) with the correlations as a function of the delays times.

Steadystate correlation function

The following code demonstrates how to calculate the $\left&lt;x(t)x(0)\right&gt;$ correlation for a leaky cavity with three different relaxation rates.

times = np.linspace(0,10.0,200) a = destroy(10) x = a.dag() + a H = a.dag() * a

corr1 = correlation_2op_1t(H, None, times, [np.sqrt(0.5) * a], x, x) corr2 = correlation_2op_1t(H, None, times, [np.sqrt(1.0) * a], x, x) corr3 = correlation_2op_1t(H, None, times, [np.sqrt(2.0) * a], x, x)

plt.figure() plt.plot(times, np.real(corr1), times, np.real(corr2), times, np.real(corr3)) plt.legend(['0.5','1.0','2.0']) plt.xlabel(r'Time $t$') plt.ylabel(r'Correlation $left&lt;x(t)x(0)right&gt;$') plt.show()

Emission spectrum

Given a correlation function $\left&lt;A(\tau)B(0)\right&gt;$ we can define the corresponding power spectrum as

$$S(\omega) = \int_{-\infty}^{\infty} \left<A(\tau)B(0)\right> e^{-i\omega\tau} d\tau.$$

In QuTiP, we can calculate S(ω) using either qutip.correlation.spectrum_ss, which first calculates the correlation function using one of the time-dependent solvers and then performs the Fourier transform semi-analytically, or we can use the function qutip.correlation.spectrum_correlation_fft to numerically calculate the Fourier transform of a given correlation data using FFT.

The following example demonstrates how these two functions can be used to obtain the emission power spectrum.

guide/scripts/spectrum_ex1.py

Non-steadystate correlation function

More generally, we can also calculate correlation functions of the kind $\left&lt;A(t_1+t_2)B(t_1)\right&gt;$, i.e., the correlation function of a system that is not in its steady state. In QuTiP, we can evaluate such correlation functions using the function qutip.correlation.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 (t₁ and t₂).

guide/scripts/correlation_ex2.py

However, in some cases we might be interested in the correlation functions on the form $\left&lt;A(t_1+t_2)B(t_1)\right&gt;$, but only as a function of time coordinate t₂. In this case we can also use the qutip.correlation.correlation_2op_2t function, if we pass the density matrix at time t₁ as second argument, and None as third argument. The qutip.correlation.correlation_2op_2t function then returns a vector with the correlation values corresponding to the times in taulist (the fourth argument).

Example: first-order optical coherence function

This example demonstrates how to calculate a correlation function on the form $\left&lt;A(\tau)B(0)\right&gt;$ 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 $g^{(1)}(\tau) = \frac{\left&lt;a^\dagger(\tau)a(0)\right&gt;}{\sqrt{\left&lt;a^\dagger(\tau)a(\tau)\right&gt;\left&lt;a^\dagger(0)a(0)\right&gt;}}$. For a coherent state |g⁽¹⁾(τ)| = 1, and for a completely incoherent (thermal) state g⁽¹⁾(τ) = 0. The following code calculates and plots g⁽¹⁾(τ) as a function of τ.

guide/scripts/correlation_ex3.py

For convenience, the steps for calculating the first-order coherence function have been collected in the function qutip.correlation.coherence_function_g1.

Example: second-order optical coherence function

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

$$\displaystyle g^{(2)}(\tau) = \frac{\langle a^\dagger(0)a^\dagger(\tau)a(\tau)a(0)\rangle}{\langle a^\dagger(0)a(0)\rangle^2}$$

For a coherent state g⁽²⁾(τ) = 1, for a thermal state g⁽²⁾(τ = 0) = 2 and it decreases as a function of time (bunched photons, they tend to appear together), and for a Fock state with n photons g⁽²⁾(τ = 0) = n(n − 1)/n² < 1 and it increases with time (anti-bunched photons, more likely to arrive separated in time).

To calculate this type of correlation function with QuTiP, we can use qutip.correlation.correlation_3op_1t, which computes a correlation function on the form $\left&lt;A(0)B(\tau)C(0)\right&gt;$ (three operators, 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 a^†(τ)a(τ) = (a^†a)(τ).

The following code calculates and plots g⁽²⁾(τ) as a function of τ for a coherent, thermal and Fock state.

guide/scripts/correlation_ex4.py

For convenience, the steps for calculating the second-order coherence function have been collected in the function qutip.correlation.coherence_function_g2.