With the QuTiP time-evolution functions (for example :func:`qutip.mesolve` and :func:`qutip.mcsolve`), a state vector or density matrix can be evolved from an initial state at 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 [see, e.g., Gardineer and Zoller, Quantum Noise, Springer, 2004] 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 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.
.. tabularcolumns:: | p{4cm} | L |
| QuTiP function | Correlation function |
|---|---|
| :func:`qutip.correlation.correlation` or :func:`qutip.correlation.correlation_2op_2t` | \left<A(t+\tau)B(t)\right> or \left<A(t)B(t+\tau)\right>. |
| :func:`qutip.correlation.correlation_ss` or :func:`qutip.correlation.correlation_2op_1t` | \left<A(\tau)B(0)\right> or \left<A(0)B(\tau)\right>. |
| :func:`qutip.correlation.correlation_4op_1t` | \left<A(0)B(\tau)C(\tau)D(0)\right>. |
| :func:`qutip.correlation.correlation_4op_2t` | \left<A(t)B(t+\tau)C(t+\tau)D(t)\right>. |
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 :func:`qutip.correlation.correlation_2op_1t` function. These correlation function sovlers return a vector or matrix (in general complex) with the correlations as a function of the delays times.
The following code demonstrates how to calculate the \left<x(t)x(0)\right> correlation for a leaky cavity with three different relaxation rates.
.. plot:: guide/scripts/correlation_ex1.py :width: 4.0in :include-source:
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 :func:`qutip.correlation.spectrum_ss`, which first calculates the correlation function using the :func:`qutip.essolve.essolve` solver and then performs the Fourier transform semi-analytically, or we can use the function :func:`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.
.. plot:: guide/scripts/spectrum_ex1.py :width: 4.0in :include-source:
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 :func:`qutip.correlation.correlation`. 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).
.. plot:: guide/scripts/correlation_ex2.py :width: 4.0in :include-source:
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 :func:`qutip.correlation.correlation` function, if we pass the density matrix at time t_1 as second argument, and None as third argument. The :func:`qutip.correlation.correlation` function then returns a vector with the correlation values corresponding to the times in taulist (the fourth argument).
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.
.. plot:: guide/scripts/correlation_ex3.py :width: 4.0in :include-source:
For convenience, the steps for calculating the first-order coherence function have been collected in the function :func:`qutip.correlation.coherence_function_g1`.
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 can use :func:`qutip.correlation.correlation_4op_1t`, which computes a correlation function on the form \left<A(0)B(\tau)C(\tau)D(0)\right> (four operators, one delay-time vector).
The following code calculates and plots g^{(2)}(\tau) as a function of \tau for a coherent, thermal and fock state.
.. plot:: guide/scripts/correlation_ex4.py :width: 4.0in :include-source:
For convenience, the steps for calculating the second-order coherence function have been collected in the function :func:`qutip.correlation.coherence_function_g2`.