Homework 3¶
Due: Monday, 2/9/26, 23:59:59 \begin{equation*} \newcommand{\mat}[1]{\mathbb{#1}} \newcommand{\vb}[1]{\mathbf{#1}} \newcommand{\dd}[1]{\,\mathrm{d}#1} \end{equation*}
Solution¶
First, a word from our sponsor: the exponential series. For real $x$, \begin{equation*} e^{x} = 1 + x + \frac{x^{2}}{2!} + \frac{x^{3}}{3!} + \cdots \qquad\text{and}\qquad e^{-x} = 1 - x + \frac{x^{2}}{2!} - \frac{x^{3}}{3!}+ \cdots \end{equation*} By the rules of exponentials, $e^{x}e^{-x} = e^{x-x} = e^{0} = 1$. We can also see this result arising from the product of the two exponential series: \begin{equation*} e^{x} e^{-x} = 1 + (x - x) + \bigg(\frac{x^{2}}{2!} - x^{2} + \frac{x^{2}}{2!} \bigg) + \cdots = 1 \end{equation*} If we replace $x$ by $\mat{A}$, the same algebra yields the same result (although we have to represent 1 as the identity matrix).
Turning now to the problem at hand, By definition, a unitary matrix satisfies $\mat{U}^\dagger = \mat{U}^{-1}$. In our case, $\mat{U} = e^{\mat{A}}$, so \begin{equation} \mat{U} = \mat{1} + \mat{A} + \frac{\mat{A}^2}{2!} + \frac{\mat{A}^3}{3!} + \cdots \end{equation} To take the adjoint (conjugate transpose) of this equation, it sure would be nice if we could just replace every $\mat{A}$ with $\mat{A}^{\dagger}$, wouldn't it! Can we justify it?
Let's start with $\mat{A}^{2}$ and see about its adjoint. \begin{equation*} (\mat{A}^{2})_{ik} = \sum_{j} A_{ij}A_{jk} = A_{ij}A_{jk} \end{equation*} where I am using Einstein summation notation in the final expression: any index in a product that appears more than once is automatically summed over. What is the conjugate transpose of this matrix? I just have to exchange $i$ and $k$ and take the complex conjugate: \begin{equation*} (A_{ik})^{\dagger} = A^{*}_{ki} = (A^{*}_{ij} A^{*}_{jk})^{\rm T} = A^{*}_{kj}A^{*}_{ji} = (\mat{A}^{\dagger} \mat{A}^{\dagger})_{ki} \end{equation*} In other words, \begin{equation*} (\mat{A}^2)^{\dagger} = (\mat{A}^{\dagger})^{2} \end{equation*} We can extend this result by induction to show that $(\mat{A}^{n})^{\dagger} = (\mat{A}^{\dagger})^{n}$. Now, armed with this insight, we can take the adjoint of \myref{eq:Udef}: \begin{equation} \mat{U}^\dagger = \mat{1} + \mat{A}^\dagger + \frac{(\mat{A}^\dagger)^2}{2!} + \cdots = e^{\mat{A}^\dagger} \end{equation}
For $\mat{U}^\dagger = \mat{U}^{-1}$, the product of the two exponential series must be the identity. If $\mat{A}^{\dagger} = -\mat{A}$, we're all set. Can we show that this is the only solution? Assume that $\mat{A}^{\dagger} = -\mat{A} + \mat{B}$, for some nonzero matrix $\mat{B}$. Then the algebra above for the product of exponential series yields \begin{equation*} \mat{U} \mat{U}^{\dagger} = e^{\mat{B}} \end{equation*} which must be the identity. The only way that is possible is if $\mat{B}^{n} = 0$ for all positive powers $n$. I think that means that $\mat{B} = 0$.
The operator that rotates the state of a spin-1/2 particle around the $x$ axis through angle $\phi$ is \begin{equation*} \hat{R}_{x}(\phi) = e^{-i \hat{S}_{x} \phi/\hbar} \end{equation*} where $\hat{S}_{x}$ is the angular momentum operator along the $x$ direction. We say that angular momentum is the generator of rotation.
The basis of the eigenstates of spin along the $z$ axis has two orthonormal states: $\ket{\uparrow}$ and $\ket{\downarrow}$. In this basis, the matrix representing $\hat{S}_{x}$ has the form \begin{equation*} \newcommand{\mat}[1]{\mathbb{#1}} \hat{S}_{x} \xrightarrow{\hat{S}_{z}\text{ basis}} \frac{\hbar}{2} \mat{X} \qquad\text{where}\qquad \mat{X} \equiv \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \end{equation*}
a. Describe in words what the matrix $\mat{X}$ does. What happens when you apply it twice?
b. Using the Taylor series for the exponential function, write out enough terms of the matrix representation of $\hat{R}_{x}(\phi)$ that you see the pattern and can sum it to yield a simple $2\times2$ matrix.
c. Show that multiplying the matrix representation of $\hat{R}_{x}(\phi)$ by its adjoint (conjugate transpose) yields the identity. This shows that $\hat{R}_{x}(\phi)$ is a unitary operator.
Solution¶
a. The matrix $\mat{X}$ exchanges the two components of a state vector. The state representing spin-up along $z$, $\begin{pmatrix} 1 \\ 0\end{pmatrix}$, becomes the state representing spin-down along $z$, $\begin{pmatrix}0 \\ 1\end{pmatrix}$, and vice versa. So, if you apply it twice, you get back to where you started; $\mat{X}^{2} = \mat{I}$.
b. Using the Taylor series for $e^{x} = 1 + x + x^{2}/2! + x^{3}/3! + \cdots $, we get \begin{align} \hat{R}_{x}(\phi) \xrightarrow{\hat{S}_{z}\text{ basis}} \mat{R}_{x}(\phi) &= \mat{I} + \frac{-i\phi}{2} \mat{X} + + \bigg( \frac{-i\phi}{2} \bigg)^{2} \frac{\mat{X}^{2}}{2!} + \bigg( \frac{-i\phi}{2} \bigg)^{3} \frac{\mat{X}^{3}}{3!} + \cdots \\ &= \mat{I} \bigg[1 - \frac{(\phi/2)^{2}}{2!} + \frac{(\phi/2)^{4}}{4!} - \cdots \bigg] + \mat{X} \bigg[ (-i \phi/2) + \frac{(-i\phi/2)^{3}}{3!} + \frac{(-i\phi/2)^{5}}{5!} + \cdots \bigg] \\ \mat{R}_{x}(\phi) &= \cos(\phi/2) \begin{pmatrix}1 & 0 \\ 0 & 1 \end{pmatrix} - i \sin(\phi/2) \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} \cos\frac{\phi}{2} & -i \sin\frac{\phi}{2} \\ -i \sin\frac{\phi}{2} & \cos\frac{\phi}{2} \end{pmatrix} \end{align}
c. The adjoint matrix is the conjugate transpose, so \begin{equation*} \mat{R}^{\dagger}_{x}(\phi) = \begin{pmatrix} \cos\frac{\phi}{2} & i \sin\frac{\phi}{2} \\ i \sin\frac{\phi}{2} & \cos\frac{\phi}{2} \end{pmatrix} \end{equation*} Multiplying the adjoint times the rotation matrix gives \begin{equation*} \mat{R}^{\dagger}_{x}(\phi) \mat{R}_{x}(\phi) = \begin{pmatrix} \cos\frac{\phi}{2} & i \sin\frac{\phi}{2} \\ i \sin\frac{\phi}{2} & \cos\frac{\phi}{2} \end{pmatrix} \begin{pmatrix} \cos\frac{\phi}{2} & -i \sin\frac{\phi}{2} \\ -i \sin\frac{\phi}{2} & \cos\frac{\phi}{2} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} = \mat{I} \end{equation*}
Therefore, $\mat{R}_{x}^{\dagger}(\phi) = \mat{R}_{x}^{-1}(\phi)$, which shows that $\hat{R}_{x}(\phi)$ is a unitary operator.
Solution¶
$\newcommand\dd[1]{\,\mathrm{d}#1}$ Starting with the constant $a_0$, we have $$ 1 = \int_0^1 a_0^2 \dd{x} = a_0^2 x \big|_0^1 = a_0^2 \qquad\text{so}\qquad\boxed{a_0 = 1}$$
The second function, $b_0 + b_1 t$, must be orthogonal to $a_0$ and it must be normalized: \begin{align*} 0 &= \int_0^1 a_0 (b_0 + b_1 t) \dd{t} = \left( b_0 t + b_1 \frac{t^2}{2} \right)_0^1 = b_0 + \frac12 b_1 = 0 \\ 1 &= \int_0^1 (b_0 + b_1 t)^2 \dd{t} = \left( b_0^2 t + 2 b_0 b_1 \frac{t^2}{2} + b_1^2 \frac{t^3}{3} \right)_0^1 = b_0^2 + b_0 b_1 + \frac{b_1^2}{3} = 1 \end{align*}
The first equation gives $b_0 = -b_1 / 2$, which we substitute into the second equation to get \begin{equation*} \frac{b_1^2}{4} - \frac{b_1^2}{2} + \frac{b_1^2}{3} = \frac{1}{12} b_1^2 = 1 \qquad\text{so}\qquad \boxed{ b_1 = -2\sqrt{3} } \qquad\text{and}\qquad \boxed{ b_0 = \sqrt{3} } \end{equation*} where I have chosen to make the constant term positive. The third function, $c_0 + c_1 t + c_2 t^2$, must be orthogonal to the first two polynomials and also normalized. Let’s take this step by step. $$ \int_0^1 a_0 (c_0 + c_1 t + c_2 t^2)\dd{t} = c_0 + \frac{c_1}{2} + \frac{c_2}{3} = 0 $$
$$ \int_0^1 \sqrt{3} (1 - 2t)(c_0 + c_1 t + c_2 t^2)\dd{t} = \sqrt{3} \left( 0 c_0 - \frac{c_1}{6} - \frac{c_2}{6} \right) = 0 \qquad\text{so}\qquad c_2 = -c_1 $$
Substituting this result into the first equation gives $c_0 + \frac{c_1}{6} = 0$, so $c_1 = -6 c_0$. Now, to normalize, we must have $$ 1 = \int_0^1 c_0^2(1 - 6t + 6t^2)^2\dd{t} = \frac{c_0^2}{5} \qquad\text{and}\qquad c_0 = \sqrt{5}$$ So, the third polynomial is $\sqrt{5}(1 - 6t + 6t^2)$. In summary, the three orthonormal polynomials are $$ 1, \sqrt{3}(1 - 2 t), \sqrt{5}(1 -6t + 6t^2) $$
import matplotlib.pyplot as plt
import numpy as np
%matplotlib widget
plt.close('all')
fig, ax = plt.subplots(figsize=(4,3))
xv = np.linspace(0, 1, 51)
ax.plot(xv, xv*0+1)
ax.plot(xv, (1 - 2*xv)*np.sqrt(3))
ax.plot(xv, (1 - 6*xv +6 * xv**2) * np.sqrt(5))
ax.grid();
For which values of $\alpha$ are the following matrices invertible? Find the inverses whenever possible. (You may do this analytically or by using code, or both!) \begin{equation*} \newcommand{\mat}[1]{\mathbb{#1}} \mat{A} = \begin{pmatrix} 1 & \alpha & 0 \\ \alpha & 1 & \alpha \\ 0 & \alpha & 1 \end{pmatrix} \qquad\qquad \mat{B} = \begin{pmatrix} 0 & 1 & \alpha \\ 1 & \alpha & 0 \\ \alpha & 0 & 1 \end{pmatrix} \end{equation*}
Solution¶
For a matrix to be invertible, its row vectors must be linearly independent, which means that its determinant must not vanish. The determinants are
\begin{equation} \det(\mat{A}) = 1 (1 - \alpha^2) - \alpha (\alpha) = 1 - 2 \alpha^2 \qquad\text{and}\qquad \det(\mat{B}) = -1 (1) + \alpha (-\alpha^2) = -(1 + \alpha^3) \end{equation} Hence, for $\mat{A}$ we must have $\alpha_{\mat{A}} \ne \pm 1/\sqrt{2}$ and $\alpha_{\mat{B}} \ne e^{i \pi n/3}$, for $n \in \{1, 2, 3\}$.
We can use the Gauss-Jordan procedure to determine the inverse of $\mat{A}$. Start with \begin{equation*} \left(\begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ \alpha & 1 & \alpha & 0 & 1 & 0 \\ 0 & \alpha & 1 & 0 & 0 & 1 \end{array}\right) \xrightarrow{R_2 \leftarrow R_2 - \alpha R_1} \left(\begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 - \alpha^2 & \alpha & -\alpha & 1 & 0 \\ 0 & \alpha & 1 &0 & 0 & 1 \end{array}\right) \end{equation*} Now divide the middle row ($R_2$) by $1-\alpha^2$ to produce a 1 in the pivot position: \begin{equation} \left(\begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 & \frac{\alpha}{1 - \alpha^2} & -\frac{\alpha}{1-\alpha^2} & \frac{1}{1-\alpha^2} & 0 \\ 0 & \alpha & 1 &0 & 0 & 1 \end{array}\right) \xrightarrow{R_3 \leftarrow R_3 - \alpha R_2} \left(\begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 & \frac{\alpha}{1 - \alpha^2} & -\frac{\alpha}{1-\alpha^2} & \frac{1}{1-\alpha^2} & 0 \\ 0 & 0 & 1 - \frac{\alpha^2}{1-\alpha^2} & \frac{\alpha^2}{1-\alpha^2} & -\frac{\alpha}{1-\alpha^2} & 1 \end{array}\right) \end{equation}
Noting that $1 - \frac{\alpha^2}{1-\alpha^2} = \frac{1 - 2 \alpha^2}{1 - \alpha^2}$, we multiply $R_3$ by $\frac{1-\alpha^2}{1 - 2\alpha^2}$ to get \begin{equation*} \left( \begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 & \frac{\alpha}{1 - \alpha^2} & -\frac{\alpha}{1-\alpha^2} & \frac{1}{1-\alpha^2} & 0 \\ 0 & 0 & 1 & \frac{\alpha^2}{1-2\alpha^2} & -\frac{\alpha}{1-2\alpha^2} & \frac{1-\alpha^2}{1-2\alpha^2} \end{array}\right) \xrightarrow{R_2 \leftarrow R_2 - \frac{\alpha}{1-\alpha^2} R_3} \left( \begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & -\frac{\alpha}{1-\alpha^2} \big(1 + \frac{\alpha^2}{1-2\alpha^2} \big) & \frac{1}{1-\alpha^2} \big(1 + \frac{\alpha^2}{1-2\alpha^2}\big) & -\frac{\alpha}{1-2\alpha^2} \\ 0 & 0 & 1 & \frac{\alpha^2}{1-2\alpha^2} & -\frac{\alpha}{1-2\alpha^2} & \frac{1-\alpha^2}{1-2\alpha^2} \end{array}\right) \end{equation*} which simplifies to $\newcommand{\DEN}{1 - 2 \alpha^2}$ \begin{equation*} \left( \begin{array}{ccc|ccc} 1 & \alpha & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & -\frac{\alpha}{\DEN} & \frac{1}{\DEN} & -\frac{\alpha}{\DEN} \\ 0 & 0 & 1 & \frac{\alpha^2}{\DEN} & -\frac{\alpha}{\DEN} & \frac{1-\alpha^2}{\DEN} \end{array} \right) \xrightarrow{R_1 \leftarrow R_1 - \alpha R_2} \left( \begin{array}{ccc|ccc} 1 & 0 & 0 & 1 + \frac{\alpha^2}{\DEN} & -\frac{\alpha}{\DEN} & \frac{\alpha^2}{\DEN} \\ 0 & 1 & 0 & -\frac{\alpha}{\DEN} & \frac{1}{\DEN} & -\frac{\alpha}{\DEN} \\ 0 & 0 & 1 & \frac{\alpha^2}{\DEN} & -\frac{\alpha}{\DEN} & \frac{1-\alpha^2}{\DEN} \end{array} \right) \end{equation*} So, \begin{equation*} \mat{A}^{-1} = \frac{1}{\DEN} \begin{pmatrix} 1 - \alpha^2 & - \alpha & \alpha^2 \\ - \alpha & 1 & -\alpha \\ \alpha^2 & - \alpha & 1 - \alpha^2 \end{pmatrix} \end{equation*} Unsurprisingly, the prefactor diverges for the values of $\alpha$ that cause the determinant to vanish.
The procedure for $\mat{B}$ is similar; I will quote the result: \begin{equation*} \mat{B}^{-1} = \frac{1}{1+\alpha^3} \begin{pmatrix} -\alpha & 1 & \alpha^2 \\ 1 & \alpha^2 & -\alpha \\ \alpha^2 & -\alpha & 1 \end{pmatrix} \end{equation*}
By midsemester in Physics 52, you will show that the energy levels of a simple harmonic oscillator of angular frequency $\omega$ are $$ \epsilon_n = \hbar \omega \left(n + \frac{1}{2} \right) \qquad\text{where}\qquad n = 0, 1, 2, ...$$ That is, the energy of the oscillator is quantized in units of $\hbar \omega$, with a minimum energy (zero-point offset) of $\hbar\omega/2$.
When such an oscillator is in contact with a reservoir held at absolute temperature $T$, its energy is not fixed. Rather, the probability of observing it to have any particular energy depends on the ratio of that energy to thermal energy, defined as $k_{\mathrm{B}} T,$ where $k_{\mathrm{B}}$ is Boltzmann's constant ($1.38\times 10^{-23}\;\mathrm{J/K}$). Specifically, the probability that the oscillator will have energy $\epsilon_n$ is given by $$ P(n) = \frac{e^{-\beta \epsilon_n}}{Z} $$ where $\beta \equiv \frac{1}{k_{\mathrm{B}} T}$ and $$ Z = \sum_{n=0}^{\infty} e^{-\beta \epsilon_n} $$ is called the partition function. (It’s just the sum of all the Boltzmann factors; dividing by $Z$ produces a normalized probability.) Therefore, the average energy in the oscillator is $$ \bar{E} = \sum_{n=0}^{\infty} P(n) \epsilon_n = \sum_{n=0}^{\infty} \frac{\epsilon_n e^{-\beta \epsilon_n}}{Z} $$
Write a function to compute the (approximate) thermal average energy in an oscillator of angular frequency $\omega$ at temperature $T$ by summing $N$ terms in the series. Make a plot of $\bar{E}(x)$, where $x = k_{\mathbf{B}} T / \hbar \omega$ for $0.025 \le x \le 2$. (Note that the “right” way to handle dimensions such as temperature and energy is to define a variable $x$ that is formed by a dimensionless ratio, as we do here.]
NB: Should you wish to compare your result to the exact expression, the following results are derived early in Physics 117 (Statistical Mechanics and Thermodynamics), $$ \bar{E} = - \frac{\partial(\ln Z)}{\partial\beta} \qquad\text{and}\qquad Z = \frac{1}{ 2 \sinh(\beta \hbar \omega/2) } $$ Here’s the template:
def energy(N:int, x:float)->float:
"""
Using the first N terms of the sum that gives the
energy of an oscillator in contact with a heat reservoir
at absolute temperature T, estimate its average energy.
Inputs:
N: the number of terms to include in the sum
b: the dimensionless ratio of thermal energy to the energy,
quantum, k_B T / hbar ω.
Output:
average energy / hbar omega
"""
ε = np.arange(0.5, 0.5 + N, 1.0) # energy in units of hbar
p = np.exp(- ε / x)
Z = np.sum(p)
return np.sum(p * ε) / Z
Solution¶
T = np.arange(0.025,2.01,0.025)
eavg = [energy(20, x) for x in T]
fig, ax = plt.subplots()
ax.plot(T, eavg, label="Numerical")
ax.set_xlabel(r"$k_\mathrm{B}T / \hbar\omega$")
ax.set_ylabel(r"$\bar{E} / \hbar\omega$");
The analytic expression is not so hard to work out: \begin{align*} -\ln Z &= \ln 2 + \ln[\sinh(\beta \hbar \omega / 2)] \\ \frac{\partial (-\ln Z)}{\partial \beta} &= \frac{1}{\sinh(\beta \hbar \omega / 2)} \cosh(\beta \hbar \omega / 2) \frac{\hbar \omega}{2} \\ \bar{E} &= \frac{\hbar \omega}{2} \coth\left(\frac{\hbar\omega}{2 k_{\mathrm{B}} T} \right) \end{align*} Let $x = \frac{2 k_{\mathrm{B}}T}{\hbar\omega}$. Then, $$ \bar{E} = \frac{\hbar \omega}{2} \frac{\cosh(1/x)}{\sinh(1/x)} $$ As $x \to 0$, both the $\cosh$ and the $\sinh$ tend to $\frac12 e^{-x}$, so their ratio tends to one and the average energy settles down to the ground-state energy, $\hbar \omega/2$. In the opposite limit, where $x \gg 1$, $\cosh(1/x) \to 1$ and $\sinh(1/x) \to 1/x$, so $\bar{E} \to \frac{\hbar\omega}{2} \times \frac{2 k_{\mathrm{B}} T}{\hbar \omega} = k_\mathrm{B} T$.
T = np.arange(0.1,2.01,0.1)
Eanalytic = 0.5 / np.tanh(0.5 / T)
ax.plot(T, Eanalytic, 'ro', label="Analytic")
ax.plot(T, T, 'k--', label="Asymptotic")
ax.legend();