View source for Solutions of the optical Bloch equations

The optical Bloch equations
:<math> 
	\dot{\rho} = -\frac{i}{\hbar}[H,\rho]
		-\frac{\Gamma}{2}   \left[    {  \sigma_+  \sigma_-  \rho - 2 \sigma_- \rho \sigma_+ 
		+ \rho \sigma_+  \sigma_- } \right] 
\,,
</math>
provide a time-dependent quantum
description of a spontaneously emitting atom driven by a classical
electromagnetic field.  Considerable insight into the physical
processes involved can be gained by studying these equations in the
transient excitation limit, as well as the steady-state limit, as we
see in this section.  We begin by considering the coherent part of the
evolution, then extend this to re-visit the Bloch sphere picture of
the optical Bloch equations, which provides useful visualizations of
transient responses and steady state solutions.  Finally, we return to [[Atoms_and_cavities#Vacuum_Rabi_Oscillations|vacuum Rabi oscillations]] and investigate how cavity loss leads to damping of the oscillations, as an illustration of master equations which are more general than 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>

=== Eigenstates of the Jaynes-Cummings Hamiltonian ===

The Hamiltonian involved in the optical Bloch equations,
:<math> 
	H = \frac{\hbar\omega_0}{2}(|e \rangle  \langle e|-|g \rangle  \langle g|)
			+ \frac{\hbar\Omega}{2} (|g \rangle  \langle e|+|e \rangle  \langle g|)
\,,
</math>
describes the evolution of an atom in a classical field.
We begin here by reviewing the coherent
evolution under <math>H</math>.

<math>H</math> arises from the Jaynes-Cummings interaction we have previously
considered in the context of [[Atoms and cavities|cavity QED]], describing a single two-level
atom interacting with a single mode of the electromagnetic field:
:<math> 
	H_{JC} = \frac{\hbar\omega_0}{2}  \left[    { |e \rangle  \langle e|-|g \rangle  \langle g| } \right] 
	+ \hbar\omega  a^\dagger  a + \frac{\hbar\Omega}{2}
			  \left[    { |g \rangle  \langle e|+|e \rangle  \langle g| } \right] 
			  \left[    { a +  a^\dagger  } \right] 
\,.
</math>
The last term in this expression is <math>H_I</math>, the dipole interaction
between atom and field.  By defining <math> \sigma_+  = |e \rangle  \langle g|</math> and
<math> \sigma_- =|g \rangle  \langle e|</math>, we may write this interaction as
:<math> 
	H_I = \frac{\hbar\Omega}{2}
			  \left[    {  \sigma_+ + \sigma_-  } \right] 
			  \left[    { a +  a^\dagger  } \right] 
\,.
</math>
In the frame of reference of the atom and field, recall that
:<math>\begin{array}{rcl}  
	 \sigma_+  a^\dagger  &\rightarrow& e^{i(\omega+\omega_0)t}    \sigma_+  a^\dagger 
\\	 \sigma_-  a    &\rightarrow& e^{-i(\omega+\omega_0)t}   \sigma_-  a
\\	 \sigma_+  a    &\rightarrow& e^{i(\omega_0-\omega)t}    \sigma_+  a
\\	 \sigma_-  a^\dagger  &\rightarrow& e^{-i(\omega_0-\omega)t} \sigma_-  a^\dagger 
\,.
\end{array}</math>
When near resonance, <math>\omega\approx\omega_0</math>, and because the
<math> \sigma_+  a^\dagger </math> and <math> \sigma_-  a</math> terms oscillate at nearly twice the frequency
of <math>\omega</math>, those terms can be dropped.  Doing so is known as the
rotating wave approximation, and it gives us a simplified
interaction Hamiltonian
:<math> 
	H_I = \frac{\hbar\Omega}{2}
			  \left[    {  a^\dagger   \sigma_- + a  \sigma_+  } \right] 
\,.
</math>
Under this approximation, it is useful to note that this interaction
merely exchanges one quantum of excitation from atom to field, and
back, so that the total number of excitations <math>N= a^\dagger  a + |e \rangle  \langle e|</math>
is a constant of the motion.  We may thus write the total Hamiltonian,
in the rotating frame with the rotating wave approximation, as
:<math> 
	H = \hbar\omega N - \delta \sigma_z +
			\frac{\hbar\Omega}{2}( a^\dagger   \sigma_- +a \sigma_+ )
\,,
</math>
where we have defined <math>\sigma_z = |e \rangle  \langle e|-|g \rangle  \langle g|</math>, and <math>\delta =
\omega-\omega_0</math>.  Below, we may use <math>g = \hbar\Omega/2</math> to simplify
writing.

What are the eigenstates of this Hamiltonian?  It describes a
two-level system coupled to a simple harmonic oscillator; when
uncoupled, if <math>\delta=0</math>, then the eigenstates are simply those of
<math>|e{\rangle}</math>, <math>|g{\rangle}</math> and <math>|n{\rangle}</math>, as shown here:
::[[Image:Transient_and_steady_state_solutions_of_the_optical_Bloch_equations-obe-Hjc.png|thumb|400px|none|]]

When coupled, degenerate energy levels split (this is sometimes called ''dynamimc Stark splitting''), with harmonic oscillator
levels <math>|n{\rangle}</math> and <math>|n+1{\rangle}</math> splitting into two energy levels separated
by <math>2g\sqrt{n+1}</math>.  Since the coupling only pairs levels separated by
one quantum of excitation, it is straightforward to show that the
eigenstates of the Jaynes-Cummings Hamiltonian fall into well defined
pairs of states, which we may label as <math>|\pm,n{\rangle}</math>; these are
:<math> 
	|\pm,n \rangle  = \frac{1}{\sqrt{2}}  \left[    { |e,n \rangle  \pm |g,n+1 \rangle  } \right] 
\,,
</math>
and they have energies
:<math> 
	E_{\pm,n} = \hbar\omega(n+1) \pm g\sqrt{n+1}
\,.
</math>
When <math>\delta\neq 0</math>, similar physics result, but with slightly more
complicated expressions describing the eigenstates, as we shall see
when we later return to the "dressed states" picture.

==== Strongly driven atom: Mollow triplet ====

An atom strongly coupled to a single mode electromagnetic field, or an
atom driven strongly by a single mode field, will thus have an
emission spectrum described by the coupled energy level diagram:
::[[Image:Transient_and_steady_state_solutions_of_the_optical_Bloch_equations-obe-mollow-energy.png|thumb|306px|none|]]
where, to good approximation, the energy level differences are
<math>\omega</math> and <math>\omega\pm\Omega</math>.  These three lines which appear in the
spectrum are known as the Mollow triplet:
::[[Image:Transient_and_steady_state_solutions_of_the_optical_Bloch_equations-obe-mollow.png|thumb|400px|none|]]

The Mollow triplet is experimentally observed in a wide variety of
systems.  However, while our energy eigenstate analysis has predicted
the number and frequencies of the emission lines, it fails to explain
a key characteristic: the widths are not the same.  If the central
peak at <math>\omega_0</math> has width HWHM <math>\Gamma/2</math>, the two sidebands each have a
HWHM of <math>3\Gamma/4</math>.  To explain this, we need the optical Bloch
equations.

=== Optical Bloch equation evolution on the Bloch sphere ===

Let's examine the case of a two-level atom with an excited state <math>|e\rangle</math> with linewdith <math>\Gamma</math> that
can decay to the ground state <math>|g\rangle</math>.  We'll use the following conventions
:<math> 
	|0{\rangle}=|e{\rangle} = \left[ \begin{array}{cc}{1}\\{0}\end{array}\right], |1{\rangle}=|g{\rangle} = \left[ \begin{array}{cc}{0}\\{1}\end{array}\right],
 \rho =  \left[ \begin{array}{cc}{\rho_{ee}}&{\rho_{eg}}\\{\rho_{ge}}&{\rho_{gg}}\end{array} \right] 
\\    X = \left[ \begin{array}{cc}{0}&{1}\\{1}&{0}\end{array} \right] , Y = \left[ \begin{array}{cc}{0}&{-i}\\{i}&{0}\end{array} \right] , Z = \left[ \begin{array}{cc}{1}&{0}\\{0}&{-1}\end{array} \right], 
\sigma_+ = \left[ \begin{array}{cc}{0}&{1}\\{0}&{0}\end{array} \right], \sigma_- = \left[ \begin{array}{cc}{0}&{0}\\{1}&{0}\end{array} \right]
\,,
</math>
Recall that the optical Bloch equation is 
:<math> 
	\dot{\rho} = -\frac{i}{\hbar}[H,\rho]
		-\frac{\Gamma}{2}   \left[    {  \sigma_+  \sigma_-  \rho - 2 \sigma_- \rho \sigma_+ 
		+ \rho \sigma_+  \sigma_- } \right] 
\,,
</math>
and the rotating frame Hamiltonian (as derived in [[Resonances|AMO I Resonances]]) is
:<math> 
	H = \frac{\hbar}{2}\left(-\delta Z + \Omega X \right) = \frac{\hbar}{2} \left[ \begin{array}{cc}{-\delta}&{\Omega}\\{\Omega}&{\delta}\end{array}\right] 
\,.
</math>
Substituting these matrices into the optical Bloch equation and evaluating the right had side gives us the equations of motion
:<math>\begin{array}{rcl}  
	\dot{\rho}_{ee} &=&  i\frac{\Omega}{2}(\rho_{eg}-\rho_{ge})  - \Gamma\rho_{ee}
\\	\dot{\rho}_{eg} &=& i\delta\rho_{eg} + i\frac{\Omega}{2}(\rho_{ee}-\rho_{gg}) - \frac{\Gamma}{2} \rho_{eg}
\\	\dot{\rho}_{ge} &=& -i\delta\rho_{ge} - i\frac{\Omega}{2}(\rho_{ee}-\rho_{gg}) - \frac{\Gamma}{2} \rho_{ge}
\\	\dot{\rho}_{gg} &=& -i\frac{\Omega}{2}(\rho_{eg}-\rho_{ge})  + \Gamma\rho_{ee}
\,.
\end{array}</math>
As a quick sanity check, we see that <math>\dot{\rho}_{ee}=-\dot{\rho}_{gg}</math> and <math>\dot{\rho}_{ge}=\dot{\rho}^*_{eg}</math>
as you would expect from taking the time derivative of the normalization condition <math>\rho_{ee}=1-\rho_{gg}</math> and the hermicity
requirement <math>\rho_{ge}=\rho^*_{eg}</math>.  Furthermore <math>i(\rho_{eg}-\rho_{ge})=i(\rho_{eg}-\rho^*_{eg})</math> is purely
real, which means that the time derivatives of <math>\rho_{ee}</math> and <math>\rho_{gg}</math> are purely real, which in turn means
that <math>\rho_{ee}</math> and <math>\rho_{gg}</math> stay purely real as the density matrix evolves.

=== Visualization: Bloch vector ===

This seemingly complicated evolution of the density matrix is easier to understand when visualized on the Bloch sphere.
We have previously seen that an arbitrary qubit state <math>|\psi \rangle  =
\cos\frac{\theta}{2}|0 \rangle  + e^{i\phi}\sin\frac{\theta}{2}|1{\rangle}</math> can be represented as being a
point on a unit sphere, located at <math>(\theta,\phi)</math> in polar
coordinates.  Similarly, a density matrix <math>\rho</math> may be depicted as
being a point inside or on the unit sphere, using
:<math> 
	\rho = \frac{I+\vec{r}\cdot{\sigma}}{2}
\,,
</math>
where <math>\vec{r}</math> is the Bloch vector representation of <math>\rho</math>.
Writing out both sides of that expression allows us to see how the Bloch vector components are related to the elements of the density matrix
:<math>\begin{array}{rcl}  
	\rho_{ee} &=& \frac{1}{2}\left(1+r_z\right)
\\	\rho_{eg} &=& \frac{1}{2}\left(r_x-ir_y\right)
\\	\rho_{ge} &=& \frac{1}{2}\left(r_x+ir_y\right)
\\	\rho_{gg} &=& \frac{1}{2}\left(1-r_z\right)
\,.
\end{array}</math>
And these relations can easily be reversed to find
:<math>\begin{array}{rcl}  
	r_x &=& \langle X \rangle &=& \rho_{ge}+\rho_{eg}
\\	r_y &=& \langle Y \rangle &=& (\rho_{ge}-\rho_{eg})/i
\\	r_z &=& \langle Z \rangle &=& \rho_{ee}-\rho_{gg}
\,.
\end{array}</math>
We can then use these relations to rewrite the above differential equations for the elements of the density matrix into differential equations for
the components of the Bloch vector.  This gives
:<math>\begin{array}{rcl}  
	\dot{r}_x &=& \delta\, r_y - \frac{\Gamma}{2} r_x
\\	\dot{r}_y &=& -\delta\, r_x - \Omega\, r_z - \frac{\Gamma}{2} r_y
\\	\dot{r}_z &=& \Omega\, r_y - \Gamma (r_z +1)
\,.
\end{array}</math>
Note how these equations of motion provide a simple set of flows on
the Bloch sphere: the <math>\delta</math> terms correspond to a rotation in the
<math>\hat{x}-\hat{y}</math> plane, <math>\Omega</math> corresponds to a rotation in the
<math>\hat{y}-\hat{z}</math> plane, and <math>\Gamma</math> drives a relaxation process
which shrinks <math>\hat{x}</math> and <math>\hat{y}</math> components of the Bloch vector,
while moving the <math>\hat{z}</math> component toward <math>r_z=-1</math>.

=== Physics: in-phase and quadrature components ===

Physically, what is the meaning of <math>r_x</math>, <math>r_y</math>, and <math>r_z</math>?  <math>r_z</math> is
manifestly the population difference between the excited and ground
states.  The other two components may be interpreted by recognizing
that the average dipole moment of the atom is
:<math>\begin{array}{rcl}  
	 \langle d \rangle  &=& d_{ge} {\rm Tr}(\rho(|e \rangle  \langle g|+|g \rangle  \langle e|))
\\
	&=& d_{ge} (\rho_{ge} e^{i\omega t} + \rho_{eg} e^{-i\omega t})
\\
	&=& d_{ge} (r_x \cos\omega t - r_y \sin\omega t)
\,.
\end{array}</math>
Thus, <math>r_x</math> and <math>r_y</math> correspond to the phase components of the atomic
dipole moment which are in-phase and in quadrature with the incident
electromagnetic field.

=== Animation of full solutions ===

Here are plots of the full solutions of the optical Bloch equations on the Bloch sphere, for various cases of <math>\Omega</math>, <math>\delta</math>, and <math>\Gamma</math> (sample matlab files to solve the differential equations: [http://cua.mit.edu/wikipost/20090304-131009/obefun1.m.txt obefun1.m]
[http://cua.mit.edu/wikipost/20090304-131009/plotobe2.m.txt plotobe2.m]):
<blockquote>
<jwplayer width="560" height="440" repeat="true" displayheight="420" 
image="http://feynman.mit.edu/8.422/plotobe6.png"
autostart="false">http://feynman.mit.edu/8.422/plotobe6.flv</jwplayer>
</blockquote>

=== Solutions of the optical Bloch equations ===

The optical Bloch equations, describing the evolution of a single atom coupled to the vacuum, 
:<math> 
	\dot{\rho} = -\frac{i}{\hbar}[H,\rho]
		-\frac{\Gamma}{2}   \left[    {  \sigma_+  \sigma_-  \rho - 2 \sigma_- \rho \sigma_+ 
		+ \rho \sigma_+  \sigma_- } \right] 
\,,
</math>
are a set of three coupled differential equations,
:<math>\begin{array}{rcl}  
	\dot{r}_x &=& \delta\, r_y - \frac{\Gamma}{2} r_x
\\	\dot{r}_y &=& -\delta\, r_x - \Omega\, r_z - \frac{\Gamma}{2} r_y
\\	\dot{r}_z &=& \Omega\, r_y - \Gamma (r_z +1)
\,.
\end{array}</math>
written in the rotating frame of the atom's Hamiltonian.  These equations can be solved analytically, but a great deal of intuition can be obtained from just studying the solutions in several limits.  Here, we consider the steady state, transient, and weak excitation limits.

==== Transient response of the optical Bloch equations ====

The optical Bloch equations allow us to study the internal state of
the atom as it changes due to the external driving field, and due to
spontaneous emission.

Starting from the time-independent form of the equations,
:<math>\begin{array}{rcl}  
	\dot{r}_x &=& \delta\, r_y - \frac{\Gamma}{2} r_x
\\	\dot{r}_y &=& -\delta\, r_x - \Omega\, r_z - \frac{\Gamma}{2} r_y
\\	\dot{r}_z &=& \Omega\, r_y - \Gamma (r_z +1)
\,.
\end{array}</math>
we may note that when <math>\Omega\rightarrow 0</math> and at resonance, <math>\delta=0</math>,
the Bloch vector exhibits pure damping behavior, towards <math>r_z=1</math>, and
<math>r_x=r_y=0</math>.  Note that a convenient way to write these differential equations is in matrix form,
:<math> 
\dot{\vec{r}} = \left[ \begin{array}{ccc}
        {-\Gamma/2}& {\delta} & {0} \\
        {-\delta} & -\Gamma/2 & -\Omega \\
        0 & \Omega & -\Gamma
        \end{array} \right] \vec{r}
        + \left[ \begin{array}{c}
        0 \\ 0 \\ -\Gamma \end{array} \right]
</math>

When <math>\Omega\gg \Gamma</math>, Rabi oscillations occur, represnted by rapid
rotations of the Bloch vector about <math>\hat{x}</math>.  Since the relaxation
along <math>\hat{y}</math> occurs at rate <math>\Gamma/2</math>, and the relaxation about
<math>\hat{z}</math> occurs at rate <math>\Gamma</math>, we might expect that the average
relaxation rate of the rotating components under such a strong driving
field would be <math>(\Gamma+\Gamma/2)/2 = 3\Gamma/4</math>.  The remaining
component <math>r_x</math> does not rotate, because it sits along <math>\hat{x}</math>, the
axis of rotation.  Thus, it relaxes with rate <math>\Gamma/2</math>.  [http://cua.mit.edu/wikipost/20090303-132628/transient.pdf Computation of the eigenvalues] of the equations of motion verify this qualitative
picture, and show that for <math>\Omega\gg \Gamma</math>, and <math>\delta=0</math>, the
eigenvalues of motion are <math>\pm i\Omega + 3\Gamma/4</math> and <math>-\Gamma/2</math>. 

These
correspond to a main peak at <math>\omega_0</math> with width <math>\Gamma/2</math>, and two
sidebands at <math>\omega_0\pm \Omega</math>, with widths <math>3\Gamma/4</math>, thus
explaining the widths of the observed Mollow triplet lines.

==== Steady-state solution of the optical Bloch equations ====

The steady state solution of the optical Bloch equations are found by
setting all the time derivatives to zero, giving a set of three
simultaneous equations,
Starting from the time-independent form of the equations,
:<math>\begin{array}{rcl}  
	0 &=& \delta\, r_y - \frac{\Gamma}{2} r_x
\\	0 &=& -\delta\, r_x - \Omega\, r_z - \frac{\Gamma}{2} r_y
\\	0 &=& \Omega\, r_y - \Gamma (r_z +1)
\,.
\end{array}</math>
The solutions are:
:<math>\begin{align}
	r_x &= (\Omega\delta) \frac{1}{\delta^2 + \Omega^2/2 + (\Gamma/2)^2}
\\
	r_y &= \left(\frac{\Omega\Gamma}{2}\right) \frac{1}{\delta^2 + \Omega^2/2 + (\Gamma/2)^2}
\\
	r_z &= -\left(\delta^2+\frac{\Gamma^2}{4}\right) \frac{1}{\delta^2 + \Omega^2/2 + (\Gamma/2)^2}
\end{align}</math>

Physically, these are Lorentzians; the <math>r_y</math> solution (the component
in quadrature with the dipole) corresponds to an absorption curve with
half-width
:<math> 
	\sqrt{\frac{\Gamma^2}{4}+\frac{\Omega^2}{2}}
\,,
</math>
and the <math>r_x</math> solution (the component in-phase with the dipole)
corresponds to a dispersion curve.  And under a strong driving field,
as <math>\Omega\rightarrow\infty</math>, <math>r_z\rightarrow 0</math>, indicating that the
populations in the excited and ground states are equalizing.  The
steady-state population in the excited state is
:<math> 
	\rho_{ee} = \frac{1+r_z}{2} = \frac{\Omega^2/4}{\delta^2+\Omega^2/2+\Gamma^2/4}
\,,
</math>
an important result that will later be used in studying light forces.

<!--
Here are plots of these solutions, as a function of detuning, for a range of values of <math>\Gamma</math>:
<blockquote>
<jwplayer width="877" height="700" repeat="true" displayheight="680" 
image="http://feynman.mit.edu/8.422/steadyobe3.png"
autostart="false">http://feynman.mit.edu/8.422/steadyobe3.flv</jwplayer>
</blockquote>
-->

These solutions can be re-expressed in a simplified manner by defining
the saturation parameter
:<math> 
	s = \frac{\Omega^2/2}{\delta^2+\Gamma^2/4}
\,,
</math>
in terms of which we find
:<math>\begin{align}
	r_x &= \frac{2\delta}{\Omega} \frac{s}{1+s}
\\
	r_y &= \frac{\Gamma}{\Omega} \frac{s}{1+s}
\\
	r_z &= -\frac{1}{1+s}
\,.
\end{align}</math>
As <math>s\rightarrow \infty</math>, the atomic transitions become [[Saturation of atomic transitions|saturated]], and the linewidth of the transition broadens from its
natural value <math>\Gamma</math>, becoming <math>\Gamma' = \Gamma\sqrt{1+s}</math> on
resonance, at <math>\delta=0</math>.

Here are plots of the steady state solution for varying <math>\Omega/\Gamma</math>:
<blockquote>
<jwplayer width="560" height="440" repeat="true" displayheight="420" 
image="http://feynman.mit.edu/8.422/steady4.png"
autostart="false">http://feynman.mit.edu/8.422/steady4.flv</jwplayer>
</blockquote>

==== Weak excitation limit solution of the optical Bloch equations ====

Once again, the time-independent form of the optical Bloch equations is
:<math>
\begin{align}
        \dot{\rho}_{ee} &= i\frac{\Omega}{2}(\rho_{eg}-\rho_{ge})
                                - \Gamma\rho_{ee}
\\
        \dot{\rho}_{ge} &= -i\delta\rho_{ge} - 
                i\frac{\Omega}{2}(\rho_{ee}-\rho_{gg})
                                - \frac{\Gamma}{2} \rho_{ge}
\,,
\end{align}
</math>
where the remaining two components of the density matrix are given by
<math>\rho_{gg} = 1-\rho_{ee}</math>, and <math>\rho_{eg} = \rho_{ge}^*</math>.  It is
insightful to study these equations in the limit of weak excitation,
and for short evolution times.

It turns out that the solution of these equations to lowest order
in <math>|\Omega|</math>, and in the limit <math>|\Omega|\ll\Gamma</math>, with the initial
conditions <math>\rho_{ee}=0</math> and <math>\rho_{ge}=0</math>, gives
:<math>
        \rho_{ee} =
        \frac{\frac{1}{4}|\Omega|^2}{\delta^2+\left(\frac{\Gamma}{2}\right)^2} 
        \left[{ 1+e^{-\Gamma t}-2\cos(\delta t)e^{-\Gamma t/2}}\right]
\,.
</math>

Moreover the solution of these equations to lowest order in
<math>|\Omega|</math> in the limit <math>|\Omega|t\ll 1</math>, with the initial conditions <math>\rho_{ee}=0</math> and <math>\rho_{ge}=0</math>, gives
:<math>
        \rho_{ee} = \frac{1}{4}|\Omega|^2 t^2
\,,
</math>
irrespective of the values of <math>(\omega_0-\omega_L)</math> and <math>\Gamma</math>.

=== Damped vacuum Rabi oscillations ===

{{:Damped vacuum Rabi oscillations}}

== References ==

[[Category:Optical Bloch Equations]]