View source for Derivation of the optical Bloch equations

The master equation is an equation of motion for a density
matrix describing an open quantum system, much like the Schrodinger
equation describes the evolution of a closed quantum system.  This
section provides a derivation of the master equation for a
spontaneously emitting atom, driven by a classical field, which is
known as the Optical Bloch Equation.

<categorytree mode=pages style="float:right; clear:right; margin-left:1ex; border:1px solid gray; padding:0.7ex; background-color:white;" hideprefix=auto>8.422</categorytree>

=== Classical model of atom and field ===
A good starting point, to appreciate the problem of open quantum
systems, is the classical model for a two-level atom coupled to a
black-body electromagnetic field, the Einstein rate equations
:<math>\begin{array}{rcl}  
	\frac{dN_g}{dt} &=& A_{eg} N_e + u(\omega_{eg})   \left[    { B_{eg}N_e -
	B_{ge} N_g } \right] 
\\
	\frac{dN_e}{dt} &=& - A_{eg} N_e + u(\omega_{eg})   \left[    { B_{ge}N_g -
	B_{eg} N_e } \right] 
\,,
\end{array}</math>
where <math>A_{eg}</math> is the spontaneous emission rate, <math>B_{ge}</math> and <math>B_{eg}</math>
are stimulated emission rates, and <math>u(\omega_{eg})</math> is the field
energy density at atomic frequency <math>\omega_{eg}</math>, for levels denoted
by <math>|e{\rangle}</math> and <math>|g{\rangle}</math>.

Is there a straightforward quantum analogue of this?  We might be
tempted to simply add a damping term to the Schr&ouml;dinger equation,
like
:<math> 
	i\hbar \partial_t |\psi \rangle  = H |\psi \rangle  - i\Gamma |\psi{\rangle}
\,,
</math>
but this is not physically allowed by quantum mechanics!
How, then, can we construct a fully quantum-mechanical description of
open system dynamics?  The key concept is that we must properly
account for noise:

::[[Image:chapter4-obe-part-1-obe-fig2.png|thumb|400px|none|]]
From classical thermodynamics, we know that any time there is energy
transfer from a system to the environment, there is entropy exchanged
back from the environment to the system.  This is a simple
illustration of the very basic fluctuation--dissipation
principle: there can be no relaxation without con-commitment noise!
To study quantum open system, we must model the appropriate quantum
noise contribution which goes along with relaxation.

=== Density matrices and closed system dynamics ===
The main tool we shall use to model open quantum systems is the
density matrix representation for quantum states, so it is helpful to
begin with a [[Decoherence Processes, Mixed States, Density Matrices: Review|review of density matrices]] and how they evolve under
Hamiltonian dynamics.

==== Review of density matrices: definition, properties, and unravlings ====

Recall that a density matrix <math>\rho</math> for a pure state <math>|\psi{\rangle}</math> is the
matrix <math>|\psi{\rangle}{\langle}\psi|</math>.  Thus, for example
:<math>\begin{array}{rcl}  
	|0 \rangle  &\rightarrow&  \left[	\begin{array}{cc}{1}&{0}\\{0}&{0}\end{array}\right] 
\\ &&~
\\	|1 \rangle  &\rightarrow&  \left[	\begin{array}{cc}{0}&{0}\\{0}&{1}\end{array}\right] 
\\ &&~
\\	\frac{|0{\rangle}+|1{\rangle}}{\sqrt{2}} &\rightarrow& \frac{1}{2}  \left[	\begin{array}{cc}{1}&{1}\\{1}&{1}\end{array}\right] 
\,.
\end{array}</math>
Density matrices may also represent statistical mixtures of pure
states; this state
:<math> 
	\frac{1}{2} \left[	\begin{array}{cc}{1}&{0}\\{0}&{1}\end{array}\right] 
</math>
can be interpreted as a 50/50 mixture of <math>|0{\rangle}</math> and <math>|1{\rangle}</math>.  However,
one must be careful, because there are infinitely many ways to ''unravel'' a density matrix into statistical mixtures of pure states.

For example,
:<math> 
	\rho = \frac{1}{4} \left[	\begin{array}{cc}{3}&{0}\\{0}&{1}\end{array}\right] 
</math>
is
:<math> 
	\rho = \frac{3}{4} |0 \rangle  \langle 0| + \frac{1}{4} |1 \rangle  \langle 1|
\,,
</math>
but it is also
:<math> 
	\rho = \frac{1}{2} |a \rangle  \langle a| + \frac{1}{2} |b \rangle  \langle b|
\,,
</math>
where
:<math>\begin{array}{rcl}  
	|a \rangle  &=& \frac{\sqrt{3}|0{\rangle}+|1{\rangle}}{2}
\\
	|b \rangle  &=& \frac{\sqrt{3}|0{\rangle}-|1{\rangle}}{2}
\,.
\end{array}</math>
In general, a density matrix <math>\rho</math> may always be written as a
statistical mixture of pure states,
:<math> 
	\rho = \sum_k p_k |\psi_k{\rangle}
\,,
</math>
where <math>p_k</math> are probabilities, such that <math>\sum_k p_k = 1</math>.

A matrix <math>\rho</math> is a valid density matrix if and only if the
eigenvalues of <math>\rho</math> are non-negative, and sum to one, such that
<math>{\rm Tr}(\rho)=1</math>.  <math>\rho</math> represents a pure state if and only if
<math>{\rm Tr}(\rho^2)=1</math>.

==== Density matrices: closed system evolution ====

How does a density matrix evolve in a closed system?  From the
Schrodinger equation
:<math> 
	i\hbar\partial_t |\psi \rangle  = H|\psi{\rangle}
\,,
</math>
it follows that a pure state 
density matrix <math>\rho=|\psi{\rangle}{\langle}\psi|</math>
evolves as
:<math>\begin{array}{rcl}  
	\dot{\rho} &=& |\dot{\psi}{\rangle}{\langle}\psi| + |\psi{\rangle}{\langle}\dot{\psi}|
\\
	&=& -\frac{i}{\hbar}[H,\rho]
\,.
\end{array}</math>
For example, if <math>\rho</math> is a two-level atom evolving under the
classical field Jaynes-Cummings Hamiltonian
:<math> 
	H = \frac{\hbar\omega_0}{2}(|e \rangle  \langle e|-|g \rangle  \langle g|) + \Omega (|g \rangle  \langle e|+|e \rangle  \langle g|)
\,,
</math>
which we may express using Pauli matrices as
:<math> 
	H = \frac{\hbar\omega_0}{2} Z + \Omega X
\,,
</math>
then for
:<math> 
	\rho =  \left[	\begin{array}{cc}{\rho_{ee}}&{\rho_{eg}}\\{\rho_{ge}}&{\rho_{gg}}\end{array}\right] 
\,,
</math>
the equation of motion for <math>\rho</math> is
:<math> 
	\dot{\rho} =  \left[	\begin{array}{cc}{i\Omega(\rho_{eg}-\rho_{ge})}&{i\omega_0\rho_{ge}-i\Omega(\rho_{ee}-\rho_{gg})}\\{-i\omega_0\rho_{eg}+i\Omega(\rho_{ee}-\rho_{gg})}&{i\Omega(\rho_{eg}-\rho_{ge})}\end{array}\right] 
\,.
</math>
We can recognize this as a rotation of the Bloch sphere about the axis
defined by
:<math> 
	\vec{n} = \frac{\hbar\omega_0}{2} \hat{z} + \Omega \hat{x}
\,,
</math>
by using the Bloch sphere representation for a density matrix,
:<math> 
	\rho = \frac{I+\vec{r}\cdot{\sigma}}{2}
\,.
</math>

=== Density matrices and open system dynamics: approach ===
We can build a mathematical model for quantum open system dynamics,
based on four basic ideas:
<table border=1>
<tr>
<td>
Density matrix evolution
</td>
<td>
Instead of pure states (eg <math>|\psi{\rangle}</math>), we describe the system
state using a density matrix <math>\rho</math>.  The equation of motion
for <math>\rho</math> is
:<math> 
	\dot{\rho} = { \mathcal  L}[ \rho ]
	\,,
</math>
where <math>{ \mathcal  L}</math> is known as the Liouvillian operator (or a
"superoperator").  For example, for Hamiltonian evolution,
we have:
:<math> 
		i\hbar \partial_t |\psi \rangle  = H |\psi{\rangle}
</math>
	which, for <math>\rho = |\psi{\rangle}{\langle}\psi|</math> gives
:<math>\begin{array}{rcl}  
		\partial_t \rho &=& |\dot\psi{\rangle}{\langle}\psi| + |\psi{\rangle}{\langle}\dot\psi|
	\\	&=& -\frac{i}{\hbar} H |\psi{\rangle}{\langle}\psi| +  |\psi{\rangle}{\langle}\psi|
			\frac{i}{\hbar} H
	\\	&=& -\frac{i}{\hbar} [H,\rho]
	\,.
\end{array}</math>
This is unitary evolution, but in general, the differential
equation for <math>\rho</math> can describe non-unitary evolution.
Such differential equations are known as "master equations."
They are nontrivial to construct, because they must restrict
<math>\rho</math> to be a legitimate density matrix at all times.
</td></tr>
<tr>
<td>
[[Entangled_Photons#Partial_Traces|Partial trace]]
</td>
<td>
We are interested in the state of the system alone, and want
to disregard the state of the environment.  If <math>\rho_{\rm
total}</math> is the state of the whole system + environment, then
the state of the system alone is
:<math>\begin{array}{rcl}  
		\rho_{sys} &=& {\rm Tr}_{env}   \left[    { \rho_{\rm total} } \right] 
	\\	&=& \sum_{\rm env} {\langle}{\rm env}| \rho_{\rm total} |{\rm env}{\rangle}
\,.
\end{array}</math>
</td></tr>
<tr>
<td>
Assumptions about the environment
</td><td>
The environment is also known as a "bath" or a "reservoir"
(cf API).  We model it as being an ensemble of oscillators, of
a variety of frequencies, which are weakly coupled to the
system.  It has several important properties:
* Large and unchanging -- Born approximation
* Short correlation time <math>\tau_c</math> -- Markov	approximation
</td></tr>
<tr>
<td>
Two (very different) timescales
</td><td>
There are two important timescales in this model:
* <math>T_{relax}</math>: A slow evolution of the system
* <math>\tau_c</math>: The fast fluctuations of the environment

We will build equations of motion which have a timescale
<math>\Delta t</math>, chosen such that <math>\tau_c \ll \Delta t \ll
T_{relax}</math>.
</td></tr>
</table>

Our goal is to construct a model dynamical equation of motion for
<math>\rho_{sys}</math> of the form
:<math> 
	\frac{\Delta \rho_{sys}}{\Delta t} = M \rho_{sys}
\,,
</math>
where <math>M</math> is time independent.  This is known as a "coarse grained"
evolution equation.  It is desirable to obtain <math>M</math> for a variety of
scenarios, including interactions where the system + environment are
atom + light, light + light, and atom + motion, for example.  Below,
we construct a master equation for the atom + vacuum using two
different approaches.

=== Beamsplitter model of the master equation ===
The physical intuition behind the master equation we desire to
construct can be captured with a simple example, which builds on the
beamsplitter we studied in previous lectures.  Consider a single
photon state <math>|\psi \rangle  =  \alpha  |0 \rangle  +  \beta  |1{\rangle}</math> (for simplicity, let the
coefficients be real-valued) entering a beamsplitter of angle
<math>\theta</math>:
::[[Image:chapter4-obe-part-1-obe-fig3.png|thumb|400px|none|]]

A vacuum state <math>|0{\rangle}</math> is input to the other port, whose output we
discard.  Let us consider the photon as being our system, and the
other (initially vacuum) mode as being our environment.  What is the
quantum state of the undiscarded output?  Naively, we might argue that
a single photon is discarded into the environment with probability
<math>p_1 =  \beta ^2 \sin^2 \theta</math>, so that we might expect the output to be
<math>|0{\rangle}</math> with probability <math>p_1</math>, and <math>|\psi{\rangle}</math> with probability <math>1-p_1</math>.

However, that (semi-)classical argument is incorrect.
The output state <math>|\phi{\rangle}</math> of the system + environment is
:<math>\begin{array}{rcl}  
	|\phi \rangle  &=& e^{i\theta ( a^\dagger  b +  b^\dagger  a)}
			  \left[    { |\psi \rangle  \otimes |0 \rangle  } \right] 
\\	&=&  \alpha  |00 \rangle  + \beta   \left[    {  \cos\theta  |10 \rangle  +  \sin\theta  |01 \rangle  } \right] 
\\	&=&   \left[    {  \alpha  |0 \rangle  + \beta  \cos\theta  |1 \rangle  } \right]  \otimes |0{\rangle}
		+   \left[    { \beta  \sin\theta  |0 \rangle  } \right]  \otimes |1{\rangle}
\,.
\end{array}</math>
Thus, the correct result is that the output is <math>|\psi_0 \rangle  =
( \alpha |0 \rangle  +  \beta  \cos\theta  |1 \rangle )/\sqrt{ \alpha ^2 +  \beta ^2 \cos^2\theta}</math>, with
probability <math>1-p_1</math>, and <math>|\psi_1 \rangle  = |0{\rangle}</math> with probability <math>p_1</math>.

These states can conveniently be written as density matrices.  The
input state is
:<math> 
	\rho_{in} = |\psi{\rangle}{\langle}\psi| =  \left[	\begin{array}{cc}{ \alpha ^2}&{ \alpha  \beta }\\{ \alpha  \beta }&{ \beta ^2}\end{array}\right] 
\,,
 
</math>
and the output state is
:<math>\begin{array}{rcl}  
	\rho_{out} = p_1 |\psi_1{\rangle}{\langle}\psi_1| + p_0 |\psi_0{\rangle}{\langle}\psi_0|
 =  \left[	\begin{array}{cc}{ \alpha ^2 +  \beta ^2 \sin^2\theta}&{ \alpha  \beta  \cos\theta }\\{ \alpha  \beta  \cos\theta }&{ \beta ^2\cos^2\theta}\end{array}\right] 
\,.
\end{array}</math>

Note that <math>\rho_{out}</math> cannot be written as <math>|\chi{\rangle}{\langle}\chi|</math> for any
pure state <math>|\chi{\rangle}</math>, because it is not pure (it is a statistical mixture).
The change in <math>\rho</math> is
:<math> 
	\Delta \rho = \rho_{out} - \rho_{in}
	=  \left[	\begin{array}{cc}{-\beta^2 (\cos 2\theta-1)/2}&{ \alpha  \beta ( \cos\theta -1)}\\{ \alpha  \beta ( \cos\theta -1)}&{ \beta ^2(\cos 2\theta-1)/2}\end{array}\right] 
\,.
</math>

Now imagine that we send the single photon state through many
beamsplitters in a sequence, each with some small tap angle <math>\theta =
\sqrt{\Gamma \Delta t/2}</math>:
::[[Image:chapter4-obe-part-1-obe-fig4.png|thumb|400px|none|]]

We make two assumptions: the environment modes always begin in the
vacuum <math>|0{\rangle}</math> (this is the Born approximation), and the environment is
completely different and uncorrelated between scattering events (this
is the Markov approximation).

This allows us to write a coarse-grained differential equation for the
photon state
:<math> 
	\frac{\Delta \rho}{\Delta t} \approx
	 \left[	\begin{array}{cc}{\dot{\rho}_{00}}&{\dot{\rho}_{01}}\\{\dot{\rho}_{10}}&{\dot{\rho}_{11}}\end{array}\right] 
	= - \Gamma  \left[	\begin{array}{cc}{- \beta ^2}&{ \alpha  \beta  /2}\\{ \alpha  \beta /2}&{ \beta ^2}\end{array}\right] 
\,.
 
</math>
Expressed as differential equations for each of the independent matrix
elements, we get (using Eq.(\ref{eq:rhoin}) for <math>\rho(t=0)</math>)
:<math>\begin{array}{rcl}  
	 \frac{d}{dt}  \rho_{00} &=& +\Gamma \rho_{11}
\\	 \frac{d}{dt}  \rho_{11} &=& -\Gamma \rho_{11}
\end{array}</math>
for the diagonal elements.  These describe the evolutions of the
probabilities of finding the photon in the <math>|0{\rangle}</math> and <math>|1{\rangle}</math> states,
and are analogous to the Einstein rate equations.  And for the
off-diagonal elements we get
:<math>\begin{array}{rcl}  
	 \frac{d}{dt}  \rho_{01} &=& -\frac{\Gamma}{2} \rho_{01}
\\	 \frac{d}{dt}  \rho_{10} &=& -\frac{\Gamma}{2} \rho_{10}
\end{array}</math>
which show the decay of the quantum coherence of the state.

The form of these differential equations, which are master equations,
is very general, and almost exactly the same result is obtained for a
two-level atom interacting with the vacuum.  In that situation, the
solution differs essentially only in that the coherences evolve as
:<math>\begin{array}{rcl}  
	 \frac{d}{dt}  \rho_{01} &=& - \left(   { i\Delta + \frac{\Gamma}{2}} \right)  \rho_{01}
\,,
\end{array}</math>
where <math>\Delta</math> is a frequency shift of the system known as the "Lamb
shift," which is due to virtual excitations to higher atomic levels.

In the atomic master equation, <math>\Gamma</math> is the spontaneous emission
rate, given by Fermi's golden rule
:<math> 
	\Gamma = \frac{2\pi}{\hbar} \sum_{\kappa \epsilon}
		| \langle g;\kappa\epsilon| V | e; 0{\rangle}|^2 \delta(\hbar\omega
		- \hbar\omega_{ge})
\,,
</math>
as derived [[Lineshape|elsewhere]].

=== Full derivation -- walk-through ===
We now turn to a full derivation of the general master equation.

==== Master equation in integral form for arbitrary interaction ====

Following the notation used in API, Chapter 4, let <math>A</math> denote the
system, and <math>R</math> the environment (known as the reservoir in API).  The
full Hamiltonian is
:<math> 
	H = H_A + H_R + V
\,,
</math>
where <math>V</math> is the system-reservoir interaction potential.  In the
interaction picture defined by <math>H_A</math> and <math>H_R</math>, the equation of motion
for the full system + reservoir density matrix is
{{Eq
|math=
:<math> 
	i\hbar \dot{\rho} = [V,\rho]
\,.
</math>
|num=3.2
}}
Integrating this once gives
:<math> 
	i\hbar \rho(t) = \int_0^t [ V(t'),\rho(t') ]\, dt'
\,.
</math>
Substituting this back into Eq.(3.2) gives
:<math> 
	\dot{\rho} = -\frac{1}{\hbar^2} \int_0^t [V(t),[V(t'),\rho(t')]]\,dt'
\,.
</math>
If we assume that the system and reservoir are initially uncorrelated,
and make the approximation that the reservoir stays unchanged (the
Born approximation), then
:<math> 
	\rho(t) \approx \rho_A(t) \otimes \rho_R(0)
	=  \rho_A(t) \otimes \rho_R
</math>
This gives us our starting point for a general master equation:
:<math> 
	\dot{\rho}_{A}
	= -\frac{1}{\hbar^2} {\rm Tr}_R
	 \int_0^t [V(t),[V(t'),\rho_A(t') \otimes \rho_R ]]\,dt'
\,.
 
</math>
An example is helpful in seeing how this equation works.  Generally,
we will take system + reservoir interactions of the form <math>V = -A R</math>,
where <math>A</math> acts only on the system, and <math>R</math> acts only on the
reservoir.  Specifically, let the system be a two-level atom, and the
environment be a single electromagnetic mode initially in the vacuum
state <math>|0{\rangle}</math>.  The atom interacts with the usual dipole interaction,
:<math> 
	V = -\frac{\Gamma}{2} (  b^\dagger  |g \rangle  \langle e| + b |e \rangle  \langle g| )
\,,
</math>
which is conveniently written using Pauli raising and lowering
operators
:<math> 
	\sigma_{\pm} = \frac{\sigma_x \pm i\sigma_y}{2}
</math>
where
:<math> 
	\sigma_- = |g \rangle  \langle e| ~~~~~~{\rm and}~~~~~~ \sigma_+ = |e \rangle  \langle g|
\,.
</math>
Insert this now into Eq.(\ref{eq:rhome}), but disregard the integral
over time (this lets us see what the essential dynamics are, at the
expense of not obtaining the correct specific rates).  The
relevant commutators are
:<math> 
	  \left[    {V, \rho_A\otimes |0 \rangle  \langle 0| } \right]  = ( \sigma_- \rho_A) \otimes |1 \rangle  \langle 0|
	 			- (\rho_A \sigma_+ ) \otimes |0 \rangle  \langle 1|
</math>
and
:<math> 
	  \left[    { V,   \left[    { V,\rho_A\otimes |0 \rangle  \langle 0| } \right]  } \right]  =
		( \sigma_+  \sigma_- \rho_A)\otimes |0 \rangle  \langle 0|
		- 2 ( \sigma_- \rho_A \sigma_+ )\otimes |1 \rangle  \langle 1|
		+ (\rho_A \sigma_+  \sigma_- ) |0 \rangle  \langle 0|
		+ {\rm other}
\,,
</math>
where the "other" terms are not diagonal in the electromagnetic mode
states, and thus disappear in the partial trace over the reservoir.
We find, finally:
:<math> 
	\dot{\rho}_{A} = -\frac{\Gamma}{2}   \left[    {  \sigma_+  \sigma_-  \rho_A - 2 \sigma_- \rho_A \sigma_+ 
		+ \rho_A \sigma_+  \sigma_- } \right] 
\,.
</math>
This is the master equation for a two-level atom dipole coupled to a
single electromagnetic mode initially in the vacuum state.  It is
written in a form known as the "Lindblad" form, which is very
common.  In atomic physics, you will often see master equations like
this.  The Lindblad form has the special property that it ensures
<math>\rho</math> is a legitimate density matrix at all times; not only does
<math>{\rm Tr}(\rho)=1</math> always, but also, its eigenvalues remain
non-negative.  And more importantly, the map from <math>\rho(t)</math> to
<math>\rho(t')</math> is completely positive, meaning that if the map
operates on just part of a larger system, the state of the larger
system remains described by a valid, positive density matrix.
Using the definitions for <math>\sigma_{\pm}</math>, if we express <math>\rho_A</math> as
:<math> 
	\rho_A =  \left[	\begin{array}{cc}{a}&{b}\\{c}&{d}\end{array}\right] 
\,,
</math>
then we find
:<math> 
	\dot{\rho}_A = -\Gamma  \left[	\begin{array}{cc}{ 2 a}&{b}\\{c}&{-2a}\end{array}\right] 
\,,
</math>
which is identical to the master equation we constructed for the
beamsplitter example, Eq.(\ref{eq:bsme}), up to a relabeling of <math>|g{\rangle}</math>
and <math>|e{\rangle}</math>.
As shown by this example, the physical picture behind the master
equation is not so complicated, even though the mathematics (used in
all its glory) can be overwhelming.
The equation of motion for <math>\dot{\rho}</math> we have obtained is very close
to the classical Einstein equations we began with, as we can see by
writing out equations of motion for the individual components of
<math>\rho</math>.  Explicitly, and including Hamiltonian evolution under the
classical field Jaynes-Cummings interaction, we find
:<math>\begin{array}{rcl}  
	\dot{\rho_{ee}} &=& i\Omega(\rho_{eg}-\rho_{ge}) - \Gamma \rho_{ee}
\\
	\dot{\rho_{ge}} &=& i\omega_0\rho_{ge} - i\Omega(\rho_{ee}-\rho_{gg})
	- \frac{\Gamma}{2} \rho_{ge}
\,,
\end{array}</math>
where the other two components are given by <math>\rho_{gg}=1-\rho{ee}</math> and
<math>\rho_{eg} = \rho_{ge}^*</math>.  These differential equations are known as
the Optical Bloch Equations, and we will base a great deal of
our study of atoms and light forces on this quantum description of
open system dynamics of a spontaneously emitting atom driven by a
classical field.

=== References ===

[[Category:Optical Bloch Equations]]