GitHub

Theoretical background

In this section, we will describe the fundamental theoretical concepts used within QuantumCumulants.jl.

The Quantum Langevin equation

In an open quantum system, the equation of motion of an operator $\mathcal{O}$ is given by the Quantum Langevin equation. Let a system be described by the Hamiltonian $H$, and subject to a number of decay channels with the rates $\gamma_n$ and corresponding damping operators $c_n$. The equation for $\mathcal{O}$ is then given by

\[\dot{\mathcal{O}} = \frac{i}{\hbar}[H,\mathcal{O}] + \sum_n \frac{\gamma_n}{2}\left(2c_n^\dagger \mathcal{O}c_n - c_n^\dagger c_n \mathcal{O} - \mathcal{O}c_n^\dagger c_n\right) + \text{noise}.\]

We did not specify the noise term. For the unconditional dynamics (the average over all measurement outcomes) the white-noise term does not contribute to averages, so we neglect it in the following. This is no longer the case when a decay channel is continuously monitored: the measurement backaction reintroduces a noise term, and QuantumCumulants can derive those stochastic equations too (see Noise & measurement backaction). The above equation is an operator equation, i.e. solving it directly has the same numerical complexity as solving a stochastic master equation. However, averaging over the above we obtain a c-number equation, which, in principle, is easy to solve. This is the basic idea in QuantumCumulants.jl: derive equations of motions of operators, and then convert them to easily solvable c-number differential equations. However, as we will see, there is another crucial step required, namely the cumulant expansion.

A brief example

To illustrate this, let us consider a simple example, the Jaynes-Cummings model. This model describes a two-level atom that is coupled to an optical resonator. The Hamiltonian reads

\[H_\mathrm{JC} = \hbar\Delta a^\dagger a + \hbar g\left(a^\dagger \sigma^{ge} + a\sigma^{eg}\right),\]

where $\Delta = \omega_\mathrm{c} - \omega_\mathrm{a}$ is the detuning between the cavity resonance frequency $\omega_\mathrm{c}$ and the atomic transition frequency $\omega_\mathrm{a}$. The frequency $g$ describes the strength of the dipole-coupling between the atom and the cavity. The operators $a$ and $a^\dagger$ are the photonic annihilation and creation operators of the cavity mode. The atom has a ground state $|g\rangle$ and an excited state $|e\rangle$, and its dynamics are described by the projection operators $\sigma^{ij} = |i\rangle\langle j|$, where $i,j \in\{g,e\}$. For simplicity, we will here assume that the system is closed, i.e. it is not subject to decay, such that the operator dynamics are described by the Heisenberg equation (i.e. only the first part of the Quantum Langevin equation).

Now, say we want to compute the field dynamics, i.e. we want to derive the equation for $\dot{a}$. Using the fundamental relations

\[[a,a^\dagger] = 1, \quad \sigma^{ij}\sigma^{kl} = \delta_{jk}\sigma^{il},\]

we derive

\[\begin{align*} \dot{a} &= -i\Delta a - ig \sigma^{ge}, \\ \dot{\sigma}^{ge} &= 2ig a \sigma^{ee} - ig a, \\ \dot{\sigma^{ee}} &= ig\left(a^\dagger\sigma^{ge} - a\sigma^{eg}\right). \end{align*}\]

Since $\dot{a}$ couples to $\sigma^{ge}$, and $\dot{\sigma}^{ge}$ to $\sigma^{ee}$, we needed to derive a total of three equations to arrive at a complete set. In order to make them easy to handle, we average over the above system of equations to obtain c-number equations. We find

\[\begin{align*} \langle\dot{a}\rangle &= -i\Delta \langle a\rangle - ig \langle\sigma^{ge}\rangle, \\ \langle\dot{\sigma}^{ge}\rangle &= 2ig \langle a\sigma^{ee}\rangle - ig \langle a \rangle, \\ \langle \dot{\sigma^{ee}}\rangle &= -2g\text{Im}\left\{\langle a^\dagger\sigma^{ge}\rangle\right\}. \end{align*}\]

The above system can, however, not be solved since we encounter terms such as $\langle a\sigma^{ee}\rangle$, meaning the set of equations is incomplete, since in general $\langle a\sigma^{ee}\rangle \neq \langle a\rangle\langle\sigma^{ee}\rangle$. A naive approach would be to derive the equations for all missing average values. Unfortunately, these equations will couple to averages of ever longer operator products. A complete set of equations can therefore not be derived, since it would consist of infinitely many equations.

The cumulant hierarchy

The coupling we just hit is generic. Every equation of motion sends an operator product through the commutator with $H$, and for an interacting Hamiltonian (the dipole coupling $g$ here) the commutator is a longer product, so an average of order $n$ couples to averages of order $n+1$. Those couple to order $n+2$, and so on. The equations form an infinite tower, the cumulant hierarchy:

order   averages at this level         what couples them
─────   ────────────────────────────   ──────────────────────────────────
  1     ⟨a⟩   ⟨σ^ge⟩   ⟨σ^ee⟩           the dipole coupling g
        │                             raises the order by one
        ▼
  2     ⟨a σ^ee⟩   ⟨a† σ^ge⟩   …
        │
        ▼
  3     ⟨a† a σ^ge⟩   …
        ⋮   infinite: the coupling never stops climbing

For an open system the two parts of the Quantum Langevin equation play opposite roles in this tower. The interaction terms in $H$ are what climb it (the vertical edges above), while the dissipator $\sum_n \gamma_n(\ldots)$ and the free, quadratic part of $H$ couple an average only to others of the same or lower order. Dissipation damps the dynamics but does not grow the hierarchy. A purely quadratic Hamiltonian with linear damping therefore never climbs and closes exactly at order $2$; the interesting case is the non-quadratic one, where the tower is genuinely infinite and has to be cut. That cut is the cumulant expansion.

Cumulant expansion

To obtain a closed set of c-number equations, we truncate the in principle infinite set of equations at a certain order. By order, we essentially mean the length of an operator product, e.g. $\langle a \rangle$ is of order $1$, $\langle a^\dagger a \rangle$ and $\langle a^\dagger \sigma^{ge}\rangle$ are of the order $2$. The order of a system determines its size and the accuracy of the underlying approximation. It is therefore an essential concept in QuantumCumulants.jl.

The way in which we truncate a system of equations is called the generalized cumulant expansion (see also R. Kubo, Generalized Cumulant Expansion Method). The joint cumulant, which we denote by $\langle\cdot\rangle_c$ of a product of operators $X_1 X_2 ... X_n$ of order $n$ is given by

\[\langle X_1 X_2 ... X_n \rangle_c := \sum_{p \in P(\mathcal{I})} \left(|p| - 1\right)! (-1)^{|p|-1} \prod_{B \in p} \langle \prod_{i\in B} X_i\rangle.\]

In the above, $\mathcal{I}=\{1,2,...,n\}$, $P(\mathcal{I})$ is the set of all partitions of $\mathcal{I}$, $|p|$ denotes the length of the partition $p$, and $B$ runs over the blocks of each partition. For example, in the case of $n=3$ we find

\[\begin{align*} \langle X_1X_2X_3 \rangle_c &= \langle X_1X_2X_3\rangle -\langle X_1X_2\rangle\langle X_3\rangle - \langle X_1X_3\rangle\langle X_2\rangle \\ &- \langle X_1\rangle\langle X_2X_3\rangle + 2\langle X_1\rangle\langle X_2\rangle\langle X_3\rangle. \end{align*}\]

Note that the joint cumulant of order $n$ is proportional to averages of order $\leq n$. Furthermore, the average of order $n$ occurs precisely once on the right-hand-side.

The joint cumulant can be thought of as a general measure for the correlation between operators. According to Theorem I from R. Kubo, Generalized Cumulant Expansion Method, the joint cumulant is zero iff any of the operators is statistically independent of the others. The key assumption we are making is to essentially invert this statement: instead of computing the joint cumulant of a given order to see if it is zero, we assume that it is. Since the average value of the same order occurs only once in the definition of the joint cumulant, we may invert the relation to express the average in terms of lower-order terms; i.e., if we assume $\langle X_1 X_2 ... X_n\rangle_c = 0$, we can write

\[\langle X_1X_2...X_n \rangle = \sum_{p \in P(\mathcal{I})\backslash \mathcal{I}} \left(|p| - 1\right)! (-1)^{|p|} \prod_{B \in p} \langle \prod_{i\in B} X_i\rangle,\]

where now $P(\mathcal{I})\backslash\mathcal{I}$ is the set of all partitions of $\mathcal{I}$ that does not contain $\mathcal{I}$ itself. In the example of $n=3$, we have

\[\langle X_1X_2X_3\rangle = \langle X_1X_2\rangle\langle X_3\rangle + \langle X_1X_3\rangle\langle X_2\rangle + \langle X_1\rangle\langle X_2X_3\rangle - 2\langle X_1\rangle\langle X_2\rangle\langle X_3\rangle.\]

In other words, by neglecting the cumulant of order $n$ we can express all averages of order $n$ in terms of averages of order $n-1$ and below. By applying this expansion recursively, we can reduce the order of any term to one as low as we choose.

Returning to the example of the Jaynes-Cummings Hamiltonian, we could use the cumulant expansion to express all the terms of order $2$ in first order only. This is also called the mean-field approach which neglects all quantum correlations of a system. For the Jaynes-Cummings model, we would then have

\[\begin{align*} \langle\dot{a}\rangle &= -i\Delta \langle a\rangle - ig \langle\sigma^{ge}\rangle, \\ \langle\dot{\sigma}^{ge}\rangle &= 2ig \langle a\rangle\langle\sigma^{ee}\rangle - ig \langle a \rangle, \\ \langle \dot{\sigma^{ee}}\rangle &= -2g\text{Im}\left\{\langle a^\dagger\rangle\langle\sigma^{ge}\rangle\right\}. \end{align*}\]

Now this system of equations forms a closed set, and can readily be implemented and solved numerically. Of course, solving the Jaynes-Cummings model in mean field is not very interesting. We could make it more interesting by doing a second-order cumulant expansion. However, this would involve deriving the equations for all the second-order averages, which is quite tedious. So instead, let's just be lazy and use QuantumCumulants.jl to do it,

using QuantumCumulants

# Symbolic parameters
@variables Δ::Real g::Real

# Hilbert space
hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = hf⊗ha

# Operators
@qnumbers a::Destroy(h)
σ(i, j) = Transition(h, :σ, i, j)

# Hamiltonian
H = Δ*a'*a + g*(a'*σ(:g,:e) + a*σ(:e,:g))

# List of first and second-order operators
ops = find_operators(h,2)

# Derive equations
meanfield(ops,H;order=2) # expand to second order

\[\begin{align*} \frac{d}{dt} \langle a\rangle =& -1.0 i \langle a\rangle \Delta -1.0 i \langle {\sigma}^{{ge}}\rangle g \\ \frac{d}{dt} \langle {\sigma}^{{ge}}\rangle =& -1.0 i \langle a\rangle g + 2.0 i g \langle a {\sigma}^{{ee}}\rangle \\ \frac{d}{dt} \langle {\sigma}^{{ee}}\rangle =& 1.0 i g \langle a^\dagger {\sigma}^{{ge}}\rangle -1.0 i g \langle a {\sigma}^{{eg}}\rangle \\ \frac{d}{dt} \langle a {\sigma}^{{ge}}\rangle =& 2.0 i \left( \langle {\sigma}^{{ee}}\rangle \langle a a\rangle + 2 \langle a\rangle \langle a {\sigma}^{{ee}}\rangle -2 \langle {\sigma}^{{ee}}\rangle \langle a\rangle ^{2} \right) g -1.0 i g \langle a a\rangle -1.0 i \Delta \langle a {\sigma}^{{ge}}\rangle \\ \frac{d}{dt} \langle a {\sigma}^{{ee}}\rangle =& -1.0 i \left( \langle {\sigma}^{{eg}}\rangle \langle a a\rangle + 2 \langle a\rangle \langle a {\sigma}^{{eg}}\rangle -2 \langle {\sigma}^{{eg}}\rangle \langle a\rangle ^{2} \right) g + 1.0 i \left( \langle a^\dagger\rangle \langle a {\sigma}^{{ge}}\rangle + \langle a\rangle \langle a^\dagger {\sigma}^{{ge}}\rangle + \langle {\sigma}^{{ge}}\rangle \langle a^\dagger a\rangle -2 \langle a^\dagger\rangle \langle a\rangle \langle {\sigma}^{{ge}}\rangle \right) g -1.0 i \Delta \langle a {\sigma}^{{ee}}\rangle \\ \frac{d}{dt} \langle a^\dagger a\rangle =& -1.0 i g \langle a^\dagger {\sigma}^{{ge}}\rangle + 1.0 i g \langle a {\sigma}^{{eg}}\rangle \\ \frac{d}{dt} \langle a {\sigma}^{{eg}}\rangle =& -1.0 i \langle {\sigma}^{{ee}}\rangle g -2.0 i \left( \langle a^\dagger\rangle \langle a {\sigma}^{{ee}}\rangle + \langle a\rangle \langle a^\dagger {\sigma}^{{ee}}\rangle + \langle {\sigma}^{{ee}}\rangle \langle a^\dagger a\rangle -2 \langle a^\dagger\rangle \langle a\rangle \langle {\sigma}^{{ee}}\rangle \right) g + 1.0 i g \langle a^\dagger a\rangle -1.0 i \Delta \langle a {\sigma}^{{eg}}\rangle \\ \frac{d}{dt} \langle a a\rangle =& -2.0 i g \langle a {\sigma}^{{ge}}\rangle -2.0 i \Delta \langle a a\rangle \end{align*}\]

Note that QuantumCumulants.jl automatizes the derivation of equations and the cumulant expansion. The final step of numerical implementation is handled by converting to the ModelingToolkitBase.jl framework via System.

References