View source for Open Quantum Systems/The Quantum Optical Master Equation

{{lecture}}

==== Weak-coupling approximation ====

In the following, we will consider a system <math>S</math> that is weakly coupled to a bath <math>B</math>. The Hamiltonian of the combined system is given by
:<math>H = H_S + H_B + H_I,</math>
where <math>H_S</math> denotes the part of Hamiltonian only acting on <math>S</math>, <math>H_B</math> only acts on <math>B</math>, and <math>H_I</math> accounts for the interaction between the two. We first perform a unitary transformation <math> \rho = U(t) \rho' U(t)^\dagger</math> into the interaction picture, i.e. (assuming <math>\hbar = 1</math>),
:<math> U(t) = \exp[-i(H_S+H_B)t].</math>
Inserting the transformed density matrix into the Liouville von Neumann equation, we obtain
:<math> \frac{d}{dt} \left(U(t) \rho' U(t)^\dagger\right) = -i[H, U(t)\rho' U(t)^\dagger],</math>
which can be brought into the convenient form
:<math> \frac{d}{dt} \rho = -i[H_I(t),\rho].</math>
The interaction Hamiltonian is now time-dependent according to
:<math> H_I(t) = U^\dagger(t)H_I U(t) = \exp[i (H_S+H_B)t] H_I \exp[-i(H_S+H_B)t].</math>
We can formally integrate the Liouville--Von Neumann equation and obtain
:<math>\rho(t) = \rho(0) - i\int\limits_0^t ds[H_I(s),\rho(s)].</math>
This expression can be inserted back into the Liouville von Neumann equation and after taking the trace over the bath, we arrive at <ref name="Breuer2002">{{Cite book
| publisher = Oxford University Press
| isbn = 9780198520634
| last1 = Breuer
| first1 = Heinz-Peter
| last2 = Petruccione
| first2  = F 
| title = The theory of open quantum systems
| location = Oxford; New York
| date = 2002
}}</ref>
:<math>\frac{d}{dt}\rho_S(t) = -\int\limits_0^t ds \text{Tr}_B\{[H_I(t),[H_I(s),\rho(s)]]\},</math>
where we have assumed that the initial state is such that the interaction does not generate any (first-order) dynamics in the bath, i.e.,
:<math>\text{Tr}_B\{[H_I(t),\rho(0)]\} = 0.</math>
So far, we have made little progress, as the right-hand side of the equation of motion still contains the density operator of the full system. To obtain a closed equation of motion for <math>\rho_S</math> only, we assume that the interaction is weak such that the influence on the bath is small. This is also known as the Born approximation. Then, we may treat the bath as approximately constant and write for the total density operator
:<math>\rho(t) = \rho_S(t) \otimes \rho_B.</math>
Then, we obtain a closed integro-differential equation for the density operator <math>\rho_S</math>,
:<math>\frac{d}{dt}\rho_S(t) = -\int\limits_0^t ds \text{Tr}_B\{[H_I(t),[H_I(s),\rho_S(s)\otimes \rho_B]]\}.</math>
Such an integro-differential equation is very difficult to handle as the dynamics at time <math>t</math> depends on the state of the system in all previous times. The equation of motion can be brought into a time-local form by replacing <math>\rho_S(s)</math> by <math>\rho_S(t)</math>. As we will see later on, this is not really an approximation, but already implicit in the weak-coupling assumption. This step brings us to an equation known as the Redfield master equation <ref name="Breuer2002" />,
:<math>\frac{d}{dt}\rho_S(t) = -\int\limits_0^t ds \text{Tr}_B\{[H_I(t),[H_I(s),\rho_S(t)\otimes \rho_B]]\}.</math>
This equation is still a non-Markovian master equation and does not guarantee to conserve positivity of the density matrix due to approximations we have made.

==== Markov approximation ====

To obtain a Markovian master equation, we first substitute <math>s</math> by <math>t-s</math> in the integrand, which does not change the bounds of the integration. Then, we can understand the parameter <math>s</math> as indicating how far we go backwards in time to account for memory effects, which can be characteristic timescale <math>\tau_B</math>, over which correlations in the bath decay. Under the Markov approximation, these memory effects are short-lived and therefore the integrand decays very quickly for <math> s \gg \tau_B</math>. Then, we may replace the upper bound of the integration by infinity, and obtain a Markovian master equation,
:<math>\frac{d}{dt}\rho_S(t) = -\int\limits_0^\infty ds \text{Tr}_B\{[H_I(t),[H_I(t-s),\rho_S(t)\otimes \rho_B]]\}.</math>
The typical timescale of the dynamics generated by this master equation is characterized by <math>\tau_R</math>, which is the relaxation time of the system due to the interaction with the bath. For the Markov approximation to be valid, this relaxation time has to be long compared to the bath correlation time, i.e, <math>\tau_R \gg \tau_B</math>. For quantum optical systems, <math>\tau_B</math> is the inverse of the optical frequency, i.e., several inverse THz, while the lifetime of an optical excitation is in the inverse MHz range. Therefore, the Markov approximation is well justified in quantum optical systems.

The two approximation we have made so far are often grouped together as the Born-Markov approximation. However, they still do not guarantee that the resulting master equation generates a quantum dynamical semigroup and hence cannot be cast into a Lindblad form. For this, another approximation is necessary.

==== Secular approximation ====

The secular approximation involves discarding fast oscillating terms in the Markovian master equation. It is therefore similar to the rotating-wave approximation (RWA) used in NMR and quantum optics. However, applying the RWA directly on the level on the interaction Hamiltonian can cause problems such as an incorrect renormalization of the system Hamiltonian <ref name="Fleming2010">{{Cite journal
| doi = 10.1088/1751-8113/43/40/405304
| issn = 1751-8121
| volume = 43
| issue = 40
| pages = 405304
| last = Fleming
| first = Chris
| coauthors = N. I. Cummings, Charis Anastopoulos, B. L. Hu
| title = The rotating-wave approximation: consistency and applicability from an open quantum system analysis
| journal = Journal of Physics A: Mathematical and Theoretical
| accessdate = 2014-05-22
| date = 2010-10-08
| url = http://iopscience.iop.org/1751-8121/43/40/405304
}}</ref>. The secular approximation (which is sometimes also called RWA), however, is carried out on the level of the quantum master equation.

To be explicit, let us write the interaction Hamiltonian in the form
:<math>H_I = \sum_\alpha A_\alpha \otimes B_\alpha,</math>
where the Hermitian operators <math>A_\alpha</math> and <math>B_\alpha</math> only act on system and bath, respectively. Assume, we have already diagonalized the system Hamiltonian <math>H_S</math>, so we know its eigenvalues <math>\varepsilon</math> and its projectors onto eigenstates <math>\Pi(\varepsilon)</math>. Then, we can project the operators <math>A_\alpha</math> on subspaces with a fixed energy difference <math>\omega</math> <ref name="Breuer2002" />,
:<math>A_\alpha(\omega) = \sum\limits_{\varepsilon'-\varepsilon=\omega}\Pi(\varepsilon)A_\alpha \Pi(\varepsilon').</math>
As the eigenvectors of <math>H_S</math> form a complete set, we can recover <math>A_\alpha</math> by summing over all frequencies,
:<math>\sum\limits_\omega A_\alpha(\omega) = \sum\limits_\omega A_\alpha(\omega)^\dagger = A_\alpha.</math>
Then, we can write the interaction Hamiltonian as
:<math>H_I = \sum\limits_{\alpha,\omega} A_\alpha(\omega) \otimes B_\alpha = \sum\limits_{\alpha,\omega} A_\alpha(\omega)^\dagger \otimes B_\alpha^\dagger.</math>
Using the relation
:<math>\exp(iH_s t)A_\alpha(\omega)\exp(-iH_s t) = \exp(-i\omega t)A_\alpha(\omega)</math>
and its Hermitian conjugate, we can write the interaction Hamiltonian in the interaction picture as
:<math>H_I(t) = \sum\limits_{\alpha,\omega} \exp(-i\omega t) A_\alpha(\omega) \otimes B_\alpha(t) = \sum\limits_{\alpha,\omega} \exp(i\omega t) A_\alpha(\omega)^\dagger \otimes B_\alpha(t)^\dagger,</math>
where we have introduced the interaction picture operators of the bath,
:<math>B_\alpha(t) = \exp(iH_B t) B_\alpha \exp(-i H_B t).</math>
Inserting this expression into the Markovian master equation leads us to <ref name="Breuer2002" />
:<math>\begin{align}\frac{d}{dt}\rho&=\int\limits_0^\infty ds \text{Tr}_B\left\{H_I(t-s)\rho_S(t)\rho_B H_I(t) -H_I(t)H_I(t-s)\rho_S(t)\rho_B\right\} + \text{H.c.}\\&=\sum\limits_{\omega,\omega'}\sum\limits_{\alpha\beta}e^{i(\omega'-\omega)t}\Gamma_{\alpha\beta}(\omega)\left[A_\beta(\omega)\rho_S(t)A_\alpha(\omega')^\dagger-A_\alpha(\omega')^\dagger A_\beta(\omega)\rho_S(t)\right] + \text{H.c.},\end{align}</math>
where we have used the one-sided Fourier transform of the bath correlation functions,
:<math>\Gamma_{\alpha\beta} = \int\limits_0^\infty ds \exp(i\omega s) \text{Tr}_B\left\{B_\alpha(t)^\dagger B_\beta(t-s)\rho_B\right\}.</math>
In the case where the state of the bath <math>\rho_B</math> is an eigenstate of the bath Hamiltonian <math>H_B</math>, the bath correlations do not depend on time.

The secular approximation is then referred to as the omission of all terms with <math>\omega' \ne \omega</math>, as these terms oscillate fast and average out. Again, in quantum optical systems, this comes from the fact that optical transition frequencies are much larger than the decay rates of excited states, i.e., <math>\tau_R \gg \tau_S</math>. Hence, we obtain
:<math>\frac{d}{dt}\rho = \sum\limits_{\omega}\sum\limits_{\alpha\beta}\Gamma_{\alpha\beta}(\omega)\left[A_\beta(\omega)\rho_S(t)A_\alpha(\omega)^\dagger-A_\alpha(\omega)^\dagger A_\beta(\omega)\rho_S(t)\right] + \text{H.c.}</math>

This master equation can now be cast into a Lindblad form. For this, we split real and imaginary parts of the coefficients <math>\Gamma_{\alpha\beta}</math> according to
:<math>\Gamma_{\alpha\beta}(\omega) = \frac{1}{2}\gamma_{\alpha\beta}(\omega) + iS_{\alpha\beta}(\omega),</math>
where the real part can be written as
:<math>\gamma_{\alpha\beta}(\omega) = \Gamma_{\alpha\beta}(\omega) + \Gamma_{\alpha\beta}(\omega)^* = \int\limits_{-\infty}^\infty ds \exp(i\omega s) \langle B_\alpha(s)^\dagger B_\beta(0)\rangle</math>
and forms a positive matrix <ref name="Breuer2002" />. Then, diagonalizing the coefficient matrix yields the Lindblad form for the dynamics
:<math>\frac{d}{dt}\rho = -i\left[H_{LS},\rho_S(t)\right] +\sum\limits_{\omega,k} \gamma_k(\omega)\left(A_k(\omega)\rho_S A_k(\omega)^\dagger - \frac{1}{2} \left\{A_k(\omega)^\dagger A_k(\omega),\rho_S\right\}\right),</math>
where the Lamb shift Hamiltonian <math>H_{LS}</math> commutes with the system Hamiltonian and hence results in a renormalization of the energy levels. This master equation is still expressed in the interaction picture, it can be transformed back to the Schrödinger picture by adding the system Hamiltonian <math>H_S</math> to the coherent part of the dynamics.

==== Interaction with the radiation field ====

In the case of an atom interacting with the quantized radiation field, the Hamiltonian for the latter is given by a sum of harmonic oscillators,
:<math>H_B = \sum\limits_{\vec{k}\lambda} \omega_\vec{k} a_{\vec{k}\lambda}^\dagger a_{\vec{k}\lambda},</math>
where <math>\lambda</math> is the polarization index <ref name="Walls1994">{{Cite book
| publisher = Springer-Verlag
| isbn = 9780387588315
| last1 = Walls
| first1 = D. F.
| last2 = Milburn
| first2 = Gerard J. 
| title = Quantum Optics
| date = 1994
}}</ref>. If the wavelength of the radiation field is much larger than the spatial extent of the atomic wavefunction, the dominant interaction term is given by the electric dipole term,
:<math>H_I = -\vec{d}\vec{E},</math>
where <math>\vec{d}</math> is the dipole operator of the atom and <math>\vec{E}</math> is the quantized electric field,
:<math>\vec{E} = i\sum\limits_{\vec{k}\lambda} \sqrt{\frac{h\omega_k}{V}} \vec{e}_\lambda\left(a_{\vec{k}\lambda}-a_{\vec{k}\lambda}^\dagger\right).</math>

We now assume the reference state of the bath, <math>\rho_B</math> is the vacuum without any photons. Then, we can use the following relations for the reservoir correlations <ref name="Breuer2002" />
:<math>\begin{align} \langle a_{\vec{k}\lambda} a_{\vec{k}'\lambda'}\rangle &= \langle a_{\vec{k}\lambda}^\dagger a_{\vec{k}'\lambda'}^\dagger \rangle = \langle a_{\vec{k}\lambda}^\dagger a_{\vec{k}'\lambda'}\rangle = 0\\
\langle a_{\vec{k}\lambda} a_{\vec{k}'\lambda'}^\dagger\rangle &= \delta_{\vec{k}\vec{k'}}\delta_{\lambda\lambda'}\\
\sum\limits_\lambda e_{\vec{k}\lambda}^i e_{\vec{k}\lambda}^j &= \delta_{ij} - \frac{k_ik_j}{k^2}.\end{align}</math>
Using this, we find for the spectral correlation tensor
:<math>\Gamma_{ij} = \sum\limits_\vec{k} \frac{h\omega_k}{V} \left(\delta_{ij}-\frac{k_ik_j}{k^2}\right)\int\limits_0^\infty ds \exp[-i(\omega_k-\omega)s].</math>
Going to the continuum limit, we can use
:<math>\begin{align}\frac{1}{V}\sum\limits_\vec{k} \to \int \frac{dk}{(2\pi)^3} &= \frac{1}{(2\pi)^3c^3}\int\limits_0^\infty d\omega_k \omega_k^2 \int d\Omega\\
\int d\Omega\left(\delta_{ij} - \frac{k_ik_j}{k^2}\right) &= \frac{8\pi}{3}\delta_{ij}.\end{align}</math>
Fortunately, this means that the spectral correlation tensor is already in a diagonal form and we can directly obtain a Lindblad form for the master equation. For the <math>s</math>-integration, we make use of the relation
:<math>\int\limits_0^\infty ds \exp[-i(\omega_k-\omega)s] = \pi\delta(\omega_k-\omega) + i\text{PV}\frac{1}{\omega_k-\omega},</math>
where <math>\text{PV}</math> denotes the Cauchy principal value. We can thus split the spectral correlation tensor into its real and imaginary part to obtain the decay rate and the Lamb shift, respectively. Note. however, that the expression for the Lamb shift is divergent as we have neglected relativistic effects that become relevant for large values of <math>\omega_k</math>. Since it only gives a slight renormalization of the energy levels of the atom, we neglect it in the following. For simplicity, let us consider the case where we are interested in the transition between only two atomic levels differing by the frequency <math>\omega</math>. Then, we finally arrive at the master equation in Lindblad form,
:<math>\frac{d}{dt}\rho = \gamma\left(\sigma_-\rho\sigma_+ - \frac{1}{2}\left\{\sigma_+\sigma_-,\rho\right\}\right),</math>
where we have used the spontaneous emission rate
:<math>\gamma = \frac{4\omega^3 d^2}{3\hbar c^3}</math>
and the spin flip operators
:<math>\sigma_\pm = \frac{1}{2}(\sigma_x\pm i\sigma_y).</math>

The master equation can be solved by expanding it into Pauli matrices. The vector <math>\langle \vec{\sigma} \rangle</math> is known as the Bloch vector and the equations of motion for its component decouple,
:<math>\begin{align}\frac{d}{dt}\langle \sigma_x\rangle &= -\frac{\gamma}{2}\langle \sigma_x\rangle\\
\frac{d}{dt}\langle \sigma_y\rangle &= -\frac{\gamma}{2}\langle \sigma_y\rangle\\
\frac{d}{dt}\langle \sigma_z\rangle &= -\gamma(\langle\sigma_z\rangle + 1).\end{align}</math>
From this, we see that the off-diagonal elements of the density matrix decay exponentially with the rate <math>\gamma/2</math>, a phenomenon which is called decoherence. The <math>z</math>-component of the Bloch vector decays exponentially with the rate <math>\gamma</math> to a steady state of <math>\langle \sigma_z\rangle = -1</math>, i.e., the atom ends up in the state with lower energy.

==== Resonance fluorescence ====

Let us assume now that in addition to the radiation field being in the vacuum state, there is a single mode that is being driven by an external laser. As the laser only affects a single mode, we can ignore its effect on the dissipative dynamics, i.e., the only consequence results from the Hamiltonian of the form
:<math>H_L = dE \cos(\omega t),</math>
where we have assumed that the laser is on resonance with the atom. Denoting the two levels of the atom by <math>|g\rangle</math>, and <math>|e\rangle</math>, respectively, we can go into the rotating frame of the laser driving by making the unitary transformation
:<math>|e\rangle = \exp(i\omega t)|e\rangle.</math>
Inserting back into the master equation, we see that the rotation exactly compensates the energy difference between <math>|e\rangle</math> and <math>|g\rangle</math>, and the laser term becomes
:<math> H_L = \frac{dE}{2}[1+\exp(2i\omega t)]</math>
In the rotating wave approximation we neglect the fast oscillating term. Then, we can write the laser Hamiltonian as
:<math>H_L = \frac{\Omega}{2}(\sigma_+ + \sigma_-),</math>
where we have introduced the Rabi frequency <math>\Omega = dE</math>. The total quantum master equation then reads
:<math>\frac{d}{dt}\rho = -i\left[H_L,\rho\right] + \gamma\left(\sigma_-\rho\sigma^+ - \frac{1}{2}\left\{\sigma_+\sigma_-,\rho\right\}\right).</math>
The interesting aspect about this master equation is that it exhibits a competition between the coherent dynamics generated by the laser and the dissipative dynamics arising from the decay into the vacuum of the radiation field. As before, we can write the master equation as an equation of motion for the Bloch vector <ref name="Breuer2002" />,
:<math>\frac{d}{dt}\langle \vec{\sigma}\rangle = G\langle \vec{\sigma}\rangle + \vec{b},</math>
using the matrix <math>G</math> given by
:<math>G = \left(\begin{array}{ccc}-\gamma/2 & 0 & 0\\0 &-\gamma/2 & -\Omega\\0 & \Omega & -\gamma\end{array}\right)</math>
and the vector <math>\vec{b}</math>,
:<math>\vec{b} = \left(\begin{array}{c} 0\\ 0\\ -\gamma\end{array}\right).</math>
This equation of motion is called the optical Bloch equation. Its stationary state <math>\langle \vec{\sigma}\rangle_s</math> can be found from the condition <math>d/dt \langle \vec{\sigma} \rangle = 0</math> and is given by
:<math>\begin{align}\langle \sigma_z\rangle_s &= -\frac{\gamma^2}{\gamma^2 + 2\Omega^2}\\
\langle \sigma_+\rangle_s = \langle \sigma_-\rangle_s^* &= -\frac{\Omega\gamma}{\gamma^2+2\Omega^2}.\end{align}</math>
Note that the population of the excited state,
:<math>p_e = \frac{1}{2}(1+\langle \sigma_z\rangle_s) = \frac{\Omega^2}{\gamma^2+2\Omega^2}</math>
is always less than <math>1/2</math> even in the limit of strong driving, i.e., <math>\Omega \gg \gamma</math>. Thus, it is not possible to create population inversion in a two level system in the stationary state by coherent driving. The population will merely saturate at <math>p_e = 1/2</math>.

It is also interesting to look at the relaxation dynamics towards the stationary state. For this, we introduce another vector expressing the difference between the Bloch vector and the stationary solution <ref name="Breuer2002" />, i.e.,
:<math> \delta\vec{\sigma} = \langle \vec{\sigma}\rangle - \langle \vec{\sigma} \rangle_s.</math>
The dynamics of this vector is described by a homogeneous differential equation,
:<math>\frac{d}{dt}\delta\vec{\sigma} = G \delta\vec{\sigma}.</math>
We find the eigenvalues of <math>G</math> to be
:<math>\begin{align}\lambda_1 &= -\frac{\gamma}{2}\\
\lambda_2 &= -\frac{3}{4}\gamma + i\mu\\
\lambda_3 &= -\frac{3}{4}\gamma - i\mu.\end{align}</math>
where <math>\mu</math> is given by
:<math>\mu = \sqrt{\Omega^2 -\left(\frac{\gamma}{4}\right)^2}.</math>
Since all eigenvalues have negative real parts, all coefficients of the <math>\delta\vec{\sigma}</math> vector will eventually decay and the stationary state <math>\langle \vec{\sigma}\rangle_s</math> is reached. If the atom is initially in the state <math>|g\rangle</math>, the resulting dynamics is given by
:<math>\begin{align}p_e &= \frac{\Omega^2}{\gamma^2+2\Omega^2}\left[1-\exp\left(-\frac{3}{4}\gamma t\right)\left(\cos \mu t +  \frac{3}{4}\frac{\gamma}{\mu} \sin \mu t\right)\right]\\
\langle\sigma_\pm\rangle &= \mp \frac{\Omega\gamma}{\gamma^2+2\Omega^2}\left[1-\exp\left(-\frac{3}{4}\gamma t\right)\left(\cos \mu t + \left\{\frac{\gamma}{4\mu}-\frac{\Omega^2}{\gamma\mu}\right\} \sin \mu t\right)\right].\end{align}</math>

==== References ====

{{reflist}}