19. The Density Matrix#

Name:

Total Points: –/100 pts

### UNCOMMENT AND RUN THIS CELL IF USING GOOGLE COLAB
# !pip install ipympl -q
# from google.colab import output
# output.enable_custom_widget_manager()

19.1. Review#

In the last notebook, we arrived at the equation of motion for the stationary states that make up a wavefunction. To summarize, if the wavefunction is

\[ \begin{align} \Psi(r,t) = \sum_nb_n|n\rangle e^{-j\omega_n t} \end{align} \]

and the Hamiltonian is

\[ \begin{align} \hat{H} = \hat{H_0} + V \end{align} \]

then the equation of motion for the \(b\) coefficients is

\[ \begin{align} \dot{b}_n = -\frac{j}{\hbar}\sum_l V_{nl}e^{j\omega_{nl}t} b_l \end{align} \]

where \(V_{ln} = \langle n|\hat{V}|l \rangle\) is the perturbation matrix element.

We noted that for nonlinear optics, we ar most interested in wavefunctions of electron clouds around atoms. In this case the perturbing Hamiltonian \(V\) accounts for the slight polarization of the atom when an external field is applied:

\[ \begin{align} \hat{V} = e\hat{\mathbf{r}}\cdot{\mathbf{E}} &= -\hat{\mathbf{\mu}}\cdot\mathbf{E} \end{align} \]

Here’s the explicit derivation:

\[\begin{split} \begin{align} \frac{\partial |\Psi \rangle}{\partial t} &= \frac{-j}{\hbar}\hat{H} | \Psi \rangle \\ \frac{\partial}{\partial t}\sum_nb_n|n\rangle e^{-j\omega_n t}&= \frac{-j}{\hbar} \left(\hat{H_0} + V \right ) \sum_nb_n|n\rangle e^{-j\omega_n t} \\ \sum_n \left ( \dot{b}_n -j\omega_n b_n \right)|n\rangle e^{-j\omega_n t}&= \sum_n -j\omega_n b_n |n\rangle e^{-j\omega_n t} - \frac{j}{\hbar} V \sum_nb_n|n\rangle e^{-j\omega_n t} \\ \sum_n \dot{b}_n|n\rangle e^{-j\omega_n t}&= - \frac{j}{\hbar} V\ \sum_nb_n|n\rangle e^{-j\omega_n t} \\ \end{align} \end{split}\]

Projecting both sides onto a state \(\langle l |e^{j\omega_l t}\) we obtain

\[\begin{split} \begin{align} \langle l |e^{-j\omega_l t} \sum_n \dot{b}_n|n\rangle e^{-j\omega_n t}&= - \frac{j}{\hbar} \langle l |e^{j\omega_l t} V \sum_nb_n|n\rangle e^{-j\omega_n t} \\ \sum_n \dot{b}_n \langle l |n\rangle e^{-j(\omega_n-\omega_l) t}&= - \frac{j}{\hbar} \sum_nb_n \langle l|V|n\rangle e^{-j(\omega_n-\omega_l) t} \\ \dot{b}_n &= -\frac{j}{\hbar}\sum_l V_{nl}e^{j\omega_{nl}t} b_l \end{align} \end{split}\]

where in the last step we have used the fact that the stationary states are orthogonal:

\[ \begin{align} \langle l | n \rangle = \int_\mathbf{r} \psi_l^*(r) \psi_n(r)d\mathbf{r} = \delta_{nl} \end{align} \]

Let’s explore what this rate equation is telling us. We can think of it this way: When the perturbation part of the Hamiltonion \(V\) operates on the wavefunction \(\psi\), it modifies each of the stationary states that make up the “bare” wavefunctions. Each of the bare states are now in a superposition of the other bare states, with \(V_{nl}\) telling us how much of the state \(|n\rangle\) is mapped onto the state \(|l \rangle\) by the perturbation \(V\). But the new state \(|l\rangle\) has a different time evolution \(e^{-j\omega_l t}\), so as the state evolves, the result of a small perturbation can actually result in a large re-mapping in the original state into the other stationary states over time. Does this sound like what we’ve been studying in nonlinear optics? I hope so!

19.2. The Density Matrix#

When treating nonlinear optics phenomena with the machinery of quantum mechanics, there are significant advantages to employing what is called the “density operator”. We already know that state vector \(|\psi \rangle\) can be written as a superpostion of the stationary states: \(| \psi \rangle = \sum_n b_n |n \rangle\). Consider now an operator \(\hat{\rho}\), which we’ll call the density matrix:

\[\begin{split} \begin{align} \hat{\rho} &= |\psi \rangle \langle \psi | \\ & = \left(\sum_m b_m |m\rangle \right) \left( \sum_n b_n^*\langle n | \right ) \end{align} \end{split}\]

To get a better feeling for this, let’s write it our in matrix form. Our state vectors can be written:

\[\begin{split} \begin{align} |\psi\rangle = \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots\\ b_N \\ \end{bmatrix} \\ \end{align} \end{split}\]
\[ \begin{align} \langle \psi | = \begin{bmatrix} b_0^* & b_1^* & b_2^* & \cdots & b_N^* \end{bmatrix} \end{align} \]

Which makes our density matrix

\[\begin{split} \begin{align} \hat{\rho} = |\psi \rangle \langle \psi | &= \begin{bmatrix} b_0^* & b_1^* & b_2^* & \cdots & b_N^* \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots\\ b_N \\ \end{bmatrix} \\ &= \begin{bmatrix} b_0b_0^* & b_0b_1^* & b_0b_2^* & \cdots & b_0b_N^* \\ b_1b_0^* & b_1b_1^* & b_1b_2^* & \cdots & b_1b_N^* \\ b_2b_0^* & b_2b_1^* & b_2b_2^* & \cdots & b_2b_N^* \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ b_Nb_0^* & b_Nb_1^* & b_Nb_2^* & \cdots & b_Nb_N^* \\ \end{bmatrix} \end{align} \end{split}\]

The density matrix contains the same exact information as the state vector itself–no more and no less. In fact, starting with the Schrodinger equation, we can write down the time evolution of the density matrix:

\[ \begin{align} \frac{\partial \hat{\rho}}{\partial t} = \frac{-j}{\hbar}[\hat{H}, \hat{\rho}] \end{align} \]

where \([\hat{H}, \hat{\rho}] = \hat{H} \hat{\rho} - \hat{\rho} \hat{H} \) is the commutator of \(\hat{H}\) and \(\hat{\rho}\). This equation contains the identical information as the Schrodinger equation.

Each element \(\rho_{mn}\) of the density matrix may be written

\[ \begin{align} \rho_{nm} = \langle n |\hat{\rho} | m \rangle = b_n b_m^* \end{align} \]

Take a look at the diagonal elements of \(\hat{\rho}\). They are of the form \(b_nb_n^* = |b_n|^2\). This is just the probability of finding the wavevector in state \(|n\rangle\). Since all the probabilities must add up to 1, the sum of the diagonal elements of \(\hat{\rho}\), also called the trace of \(\hat{\rho}\) must be 1.

19.3. Two Level Atom#

To explore how to use the density matrix, let’s consider the simplest possible system: an atom with just two states. In this case we have

\[\begin{split} \begin{align} \hat{\rho} = \begin{bmatrix} |b_0|^2 & b_0b_1^* \\ b_1b_0^* & |b_1|^2 \\ \end{bmatrix} \end{align} \end{split}\]

Again, the diagonal elements represent the probability of find the atom in that state. It turns out the off-diagonal elements represent the polarization of the atom, but we’ll save that for a minute. First, let’s practice putting the atom in the ground state but setting \(b_0 = 1\) and \(b_1 = 0\). Our density matrix becomes:

\[\begin{split} \begin{align} \hat{\rho}_{ground} = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ \end{bmatrix} \end{align} \end{split}\]

Ok, in order find the time evolution of the state vector probabilities we need to find the commutator \([\hat{H}, \hat{\rho}]\):

\[\begin{split} \begin{align} [\hat{H}_0 + V, \hat{\rho}] &= \hat{H}_0\hat{\rho} - \hat{\rho}\hat{H}_0 + V\hat{\rho} - \hat{\rho}V \\ &= \hat{H}_0|\psi \rangle \langle \psi | - |\psi \rangle \langle \psi |\hat{H}_0 + V|\psi \rangle \langle \psi - |\psi \rangle \langle \psi| V \\ \end{align} \end{split}\]

We can write this in matrix form:

\[\begin{split} \begin{align} [\hat{H}, \hat{\rho}] &= \begin{bmatrix} E_0|b_0|^2 & E_0b_0b_1^* \\ E_1b_1b_0^* & E_1|b_1|^2 \\ \end{bmatrix} - \begin{bmatrix} E_0|b_0|^2 & E_1b_0b_1^* \\ E_0b_1b_0^* & E_1|b_1|^2 \\ \end{bmatrix} + \begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \\ \end{bmatrix} \begin{bmatrix} |b_0|^2 & b_0b_1^* \\ b_1b_0^* & |b_1|^2 \\ \end{bmatrix} - \begin{bmatrix} |b_0|^2 & b_0b_1^* \\ b_1b_0^* & |b_1|^2 \\ \end{bmatrix} \begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \\ \end{bmatrix}\\ &= \begin{bmatrix} 0 & (E_0 - E_1)b_0b_1^* \\ (E_1-E_0)b_1b_0^* & 0 \\ \end{bmatrix} + \begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \\ \end{bmatrix} \begin{bmatrix} |b_0|^2 & b_0b_1^* \\ b_1b_0^* & |b_1|^2 \\ \end{bmatrix} - \begin{bmatrix} |b_0|^2 & b_0b_1^* \\ b_1b_0^* & |b_1|^2 \\ \end{bmatrix} \begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \\ \end{bmatrix} \end{align} \end{split}\]

Using the definition of matrix multiplication, we can find the time derivative of any element of the density matrix to be:

\[ \begin{align} \dot{\rho}_{mn} = -j\omega_{mn}\rho_{mn} - \frac{j}{\hbar}\sum_q\left(V_{mq}\rho_{qn} - \rho_{mq}V_{qn} \right) \end{align} \]

Since the rule for matrix multiplication doesn’t depend on the size the of the matrices, this is actually a general result.

19.4. Example#

Let’s consider again the Hydrogen atom, but let’s assume for a moment that it only has two levels \(\psi_{100}\) and \(\psi_{200}\) and that there is no perturbation to these levels (\(V=0\)). Now let’s put it in the superposition state

\[\begin{split} \begin{align} \Psi(r,t) &= \frac{\psi_{100}e^{\frac{-j}{\hbar}E_1 t} + \psi_{200}e^{\frac{-j}{\hbar}E_2 t}}{\sqrt{2}} = \frac{1}{\sqrt{2}}\begin{bmatrix} e^{-j\omega_1 t} \\ e^{-j\omega_2 t} \\ \end{bmatrix} \end{align} \end{split}\]

and write down the density matrix:

\[\begin{split} \begin{align} \hat{\rho} &= |\Psi \rangle \langle \Psi | \\ & = \begin{bmatrix} \frac{1}{2}e^{-j\omega_1 t} & \frac{1}{2}e^{-j\omega_{12} t} \\ \frac{1}{2}e^{j\omega_{12} t} & \frac{1}{2}e^{-j\omega_{2} t} \\ \end{bmatrix} \end{align} \end{split}\]

We take the time derivative to find

\[\begin{split} \begin{align} \frac{\partial \hat{\rho}}{\partial t} = \begin{bmatrix} \frac{-j\omega_1}{2}e^{-j\omega_1 t} & \frac{-j\omega_{12}}{2}e^{-j\omega_{12} t} \\ \frac{j\omega_{12}}{2}e^{j\omega_{12} t} & \frac{-j\omega_2}{2}e^{-j\omega_{2} t} \\ \end{bmatrix} \end{align} \end{split}\]

and we can verify that it follows the time evolution equation we derived.

19.5. Exercises#

  1. (30 pts) For the superposition state considered in this example, plot the time evolution of the probability of finding the atom in the ground state and the excited state. Use \(\omega_1 = 13.6\) and \(\omega_2 = 13.6/4\).

  2. (20 pts)Now treat all four elements of the density matrix as complex numbers \(x(t) + jp(t)\) and plot them in the complex plane with \(x\) and \(p\) as the coordinate axes.

  3. (10 pts) Check that \(r^2 = x^2 + y^2\) corresponds to the probabilities you found in part (1).

  4. (10 pts) Some people call the off diagonal elements the “coherence” between two states. Does your plot in the complex plane give any insight into what the off-diagonal elements may be describing?

  5. (30 pts) Add a very small perturbation $\(\begin{align} V = \begin{bmatrix} 0 & V_{12} \\ V_{21} & 0 \end{bmatrix} \end{align}\)\( where \)V_{12} = V_{21}\( is a variable that you can adjust. Since \)E_n = \hbar\omega_n\(, try making \)V_{12}\( a fraction of this number, and setting \)\hbar = 1$. Plot the new time evolution of the density matrix elements and the state populations.

It is often much more convenient to define a new density matrix \(\sigma\) such that every element is modified according to

\[ \begin{align} \sigma_{mn} = \rho_{mn}e^{j\omega_{mn}t} \end{align} \]

This takes away the natural time dependence of the bare states and justs looks at the effect of the oscillations induced by the perturbations. Plot the elements of \(\sigma\) and the state populations. Animate your plot so that you can adjust \(V_{12}\).