Correlation Functions

Contents

• Two-Time Correlation Functions
• Steady State Correlation Functions
• Emission Spectrum
• Non-Steady State Correlation Function

Two-Time Correlation Functions

With the QuTiP time-evolution functions (for example mesolve and mcsolve), a state vector or density matrix can be evolved from an initial state at :math:t_0 to an arbitrary time $t$, $\rho(t)=V(t, t_0)\left\{\rho(t_0)\right\}$, where $V(t, t_0)$ 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 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 $\rho(t)=V(t, 0)\left\{\rho(0)\right\}$ using one of the QuTiP evolution solvers with $\rho(0)$ as initial state, and then again use the same solver to calculate $V(t+\tau, t)\left\{B\rho(t)\right\}$ using $B\rho(t)$ as the initial state. Note that if the intial 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 $\tau$.

QuTiP provides a family of functions that assists in the process of calculating two-time correlation functions. The available functions and their usage is show 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.

QuTiP Function	Correlation Function Type
`correlation` or `correlation_2op_2t`	$\left$ or $\left$.
`correlation_ss` or `correlation_2op_1t`	$\left$ or $\left$.
`correlation_3op_1t`	$\left$.
`correlation_3op_2t`	$\left$.
`correlation_4op_1t` (Depreciated)	$\left$
`correlation_4op_2t` (Depreciated)	$\left$

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

Steady State Correlation Function

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

Emission Spectrum

Given a correlation function $\left<A(\tau)B(0)\right>$ 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(\omega)$ using either spectrum, which first calculates the correlation function using the essolve solver and then performs the Fourier transform semi-analytically, or we can use the function 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.

Non-Steady State Correlation Function

More generally, we can also calculate correlation functions of the kind $\left<A(t_1+t_2)B(t_1)\right>$, i.e., the correlation function of a system that is not in its steadystate. In QuTiP, we can evoluate 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 ($t_1$ and $t_2$).

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

Ex: First-Order Optical Coherence Function

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

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

Example: Second-Order Optical Coherence Function

The second-order optical coherence function, with time-delay $\tau$, 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^{(2)}(\tau) = 1$, for a thermal state $g^{(2)}(\tau=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^{(2)}(\tau = 0) = n(n - 1)/n^2 < 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 could use correlation_4op_1t, which computes a correlation function of the form $\left<A(0)B(\tau)C(\tau)D(0)\right>$ (four operators, one delay-time vector). However, the middle pair of operators are evaluated at the same time $\tau$, and thus can be simplified to a single operator $E(\tau)=B(\tau)C(\tau)$, and we can instead call the correlation_3op_1t function to compute $\left<A(0)E(\tau)D(0)\right>$. This simplification is done automatically inside the depreciated correlation_4op_1t function that calls correlation_3op_1t internally.

The following code calculates and plots $g^{(2)}(\tau)$ as a function of $\tau$ for coherent, thermal and fock states.

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